# Can OrdinaryDiffEq.jl solve Fillippov system?

**URL:** <https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818>\
**Category:** Specific Domains\
**Tags:** question\
**Created:** [August 15, 2023, 8:28am UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818 "2023-08-15T08:28:04Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![xiaoming](https://avatars.discourse-cdn.com/v4/letter/x/a6a055/32.png) [@xiaoming](https://discourse.julialang.org/u/xiaoming)\
**Post date:** [August 15, 2023, 8:28am UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/1 "2023-08-15T08:28:04Z")

</div>

Fillippov systems are piecewise smooth ODEs and the solutions can slide along the switch surfaces. Here is an excellent illustration of such systems, taken from this page [[Events and Discontinuities in Differential Equations—Wolfram Language Documentation](https://reference.wolfram.com/language/tutorial/NDSolveWhenEvents.html.en?source=footer)]

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/d/ad3a20a61de0dd2bb80474eadb1570cafb6872aa.png)  
Here is a typical Fillippov system:

\begin{aligned} \dot{x}=&y, \\ \dot{y}=&-cy-kx+F(y-v\_{s})+\epsilon \cos(\omega t),\\ \end{aligned}

where F is given by

F(x)=-A\ \text{Sign}(x)(\frac{\alpha}{1+\gamma|x|}+\lambda+\eta x^{2}).

So I use OrdinaryDiffEq.jl to solve this ODE. Here is my code:

```julia
using StaticArrays, OrdinaryDiffEq, Plots

function F(x)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    -A*sign(x)*(α/(1+γ*abs(x))+λ+η*x^2)
end

function vec(x,p,t)
    SVector(x[2],-p[1]*x[2]-p[2]*x[1]+F(x[2]-p[3])+p[4]*cos(p[5]*t)) # p=[c,k,v_s,ϵ,ω]
end

prob=ODEProblem(vec, SVector(1.0,3.0),(0,20),SVector(0.2,1,1,0.2,1.067))

sol=solve(prob, Vern6())

plot(sol, idxs = (1, 2))

```

Here is the solution:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/c/cc5cdb97125a5e67aaa427d05e6beaf78b08b0ff.png)  
As you can see from the picture, in the switch surface y=1, the solution is oscillating slightly and the integrator spent lots of time on the switch surface. It’s clearly not what I want.

I have tried using the event handler, but it just gives the wrong solution:

```julia
using StaticArrays, OrdinaryDiffEq, Plots

function F(x,s)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    -A*s*(α/(1+γ*abs(x))+λ+η*x^2)
end

function vec(x,p,t)
    SVector(x[2],-0.2*x[2]-x[1]+F(x[2]-1, p[1])+0.2*cos(1.067*t)) 
end

function h(u, t, integrator)
    u[2]-integrator.p[1]
end

function affect!(integrator)
    integrator.p=-integrator.p
end

cb=ContinuousCallback(h, affect!)

prob=ODEProblem(vec, SVector(1.0,3.0),(0,20),[1])

sol=solve(prob, Vern6(), callback=cb, dtmax=0.1)

plot(sol, idxs = (1, 2))

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/8/98d98fe1673bd530009a08f765962b6c8cf915ee.png)  
So, does any solution available?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 15, 2023, 1:08pm UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/2 "2023-08-15T13:08:33Z")

</div>

Those aren’t the same equation. In the first case you have `-A*sign(x)* ....` In the second case, you have `-A *s*...`, so presumably `s = sign(x)` where `s = p[1]`. So does `p[1] = sign(x)`? In the second case you have the condition `u[2]-integrator.p[1]`, so when `u[2] = integrator.p[1]` you switch the parameters. When you start, `p[1] = 1`, so when `u[2] = 1` you swap `p[1] = -1`. Now next time around, the sign change would occur when `u[2] = -1`. In your original model, the point where the sign changed did not move after the first crossing, so these are not the same model.

> [@xiaoming](#):
>
> the solution is oscillating slightly and the integrator spent lots of time on the switch surface. It’s clearly not what I want.

If I do:

```julia
function F(x)
    @show sign(x)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    -A*sign(x)*(α/(1+γ*abs(x))+λ+η*x^2)
end

```

you see:

```julia
sign(x) = 1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0
sign(x) = 1.0
sign(x) = -1.0

```

Those oscillations are real. They are your system switching over and over and over. It seems that there are just “infinitely many” events going on in the small interval.

---

<div class="post-metadata">

**Author:** ![xiaoming](https://avatars.discourse-cdn.com/v4/letter/x/a6a055/32.png) [@xiaoming](https://discourse.julialang.org/u/xiaoming)\
**Post date:** [August 15, 2023, 1:45pm UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/3 "2023-08-15T13:45:20Z")

</div>

> [@ChrisRackauckas](#):
>
> Those aren’t the same equation. In the first case you have `-A*sign(x)* ....` In the second case, you have `-A *s*...`, so presumably `s = sign(x)` where `s = p[1]`. So does `p[1] = sign(x)`? In the second case you have the condition `u[2]-integrator.p[1]`, so when `u[2] = integrator.p[1]` you switch the parameters. When you start, `p[1] = 1`, so when `u[2] = 1` you swap `p[1] = -1`. Now next time around, the sign change would occur when `u[2] = -1`. In your original model, the point where the sign changed did not move after the first crossing, so these are not the same model.

Thank you for pointing out why my second code was wrong! Yes, the solution was oscillating around u[2]=1, since the vector fields in the two sides of u[2]=1 are squeezing towards u[2]=1. Just like this picture:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/2/822695345527db12eab4c1018281d60eb553d0e7.png)  
In non-smooth nonlinear dynamics, this phenomenon is called the sliding mode. My question is whether the existing solvers can handle this sliding mode more accurately. For example, detect when the solution enters the sliding mode or when it exits, and the solution can just slide along the switch surface, without oscillating.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [August 15, 2023, 1:54pm UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/4 "2023-08-15T13:54:07Z")

</div>

This seems more like something for a bifurcation analysis tool to me. Have you considered something like BifurcationKit?

[https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/)

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [August 15, 2023, 2:20pm UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/5 "2023-08-15T14:20:44Z")

</div>

For this type of thing, you have to simulate differently at the switching surface to eliminate infinite switch. This series of lectures explain this [https://www.youtube.com/watch?v=MNbc69pcio0&list=PLp8ijpvp8iCvFDYdcXqqYU5Ibl\_aOqwjr&index=37](https://www.youtube.com/watch?v=MNbc69pcio0&list=PLp8ijpvp8iCvFDYdcXqqYU5Ibl_aOqwjr&index=37)

---

<div class="post-metadata">

**Author:** ![xiaoming](https://avatars.discourse-cdn.com/v4/letter/x/a6a055/32.png) [@xiaoming](https://discourse.julialang.org/u/xiaoming)\
**Post date:** [August 15, 2023, 3:10pm UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/6 "2023-08-15T15:10:23Z")

</div>

OK. I will check that.

---

<div class="post-metadata">

**Author:** ![xiaoming](https://avatars.discourse-cdn.com/v4/letter/x/a6a055/32.png) [@xiaoming](https://discourse.julialang.org/u/xiaoming)\
**Post date:** [August 15, 2023, 3:12pm UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/7 "2023-08-15T15:12:52Z")

</div>

I am afraid not. The parameters of the system do not change during the simulation. I am just looking for a suitable solver for this type of system.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [August 15, 2023, 7:00pm UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/8 "2023-08-15T19:00:10Z")

</div>

I only have a small amount of experience with these kinds of problems.

I would try to use `x[1]` as a parameter (instead of a dynamical variable) in the bifurcation analysis. It should then yield this ‘sliding mode’ as a fix point manifold.

---

<div class="post-metadata">

**Author:** ![xiaoming](https://avatars.discourse-cdn.com/v4/letter/x/a6a055/32.png) [@xiaoming](https://discourse.julialang.org/u/xiaoming)\
**Post date:** [August 17, 2023, 1:32am UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/9 "2023-08-17T01:32:14Z")

</div>

Hi, Dear Chris @ChrisRackauckas. I am designing a function to switch the sliding mode manually. In the document of DifferentialEquation.jl, the function ContinuousCallback has a keyword `repeat_nudge` to avoid repeating find zeros. So, if I have an initial time `t0`, can I set a fixed time, say `1//100`? Only after `t0+1//100`, the event handler starts to work.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 17, 2023, 2:54am UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/10 "2023-08-17T02:54:37Z")

</div>

> [@xiaoming](#):
>
> So, if I have an initial time `t0`, can I set a fixed time, say `1//100`? Only after `t0+1//100`, the event handler starts to work.

You can change the parameter. It’s just a keyword argument. But this has so many nearby events that I would suggest @josuagrw 's approach as the issue is going to be down to floating point error

---

<div class="post-metadata">

**Author:** ![xiaoming](https://avatars.discourse-cdn.com/v4/letter/x/a6a055/32.png) [@xiaoming](https://discourse.julialang.org/u/xiaoming)\
**Post date:** [August 19, 2023, 12:09pm UTC](https://discourse.julialang.org/t/can-ordinarydiffeq-jl-solve-fillippov-system/102818/11 "2023-08-19T12:09:52Z")

</div>

Hello everyone, thank you for all your advice. Now I have a solution for this type of problem. I will just give the code, and refer to [this paper](https://www.sciencedirect.com/science/article/abs/pii/S0378475410003186) for basic knowledge about the Fillippov systems.

```julia
using StaticArrays, OrdinaryDiffEq, Plots, LinearAlgebra

struct ODESol
    solution # ODE solution solved by OrdinaryDiffEq.jl
    state::Symbol # only two symbols are used: :regularmode and :slidemode
end

function fillippov_solve(f1, f2, h, ∇h, v0, p0, t0, tend)
    # it is assumed that f1 is at h(x)<0 and f2 is at h(x)>0
    # f1, f2, h, ∇h are all type of f(u, p, t), the output is a vector
    # initial condition: v0 
    # time interval: t0,tend
    # parameter p0
    # 
    # first define the vector field on the switch surface
    function sv(x, p, t)
        α = dot(∇h(x, p, t), f1(x, p, t)) / dot(∇h(x, p, t), f1(x, p, t) - f2(x, p, t))
        (1 - α) * f1(x, p, t) + α * f2(x, p, t)
    end

    # the hyper surface 
    function newh(x, t, integrator)
        h(x, integrator.p, t)
    end
    # stop the integrator
    function stop(integrator)
        if (integrator.t-integrator.sol.t[1])< 1e-2 # to avoid the integrator stopping at the beginning when it exits the sliding surface
            nothing
        else
            terminate!(integrator)
        end
    end

    # exit the sliding mode when slide_exit=0
    function slide_exit(x, t, integrator)
        dot(∇h(x, integrator.p, t), f2(x, integrator.p, t)) * dot(∇h(x, integrator.p, t), f1(x, integrator.p, t))
    end

    # main data consist of different solutions at h>0, h<0 and h=0
    sols = Vector{ODESol}(undef, 1)

    if abs(h(v0, p0, t0)) < 1e-10 # if this happen, we assume that the initial condition is on h=0
        prob = ODEProblem(sv, v0, (t0, tend), p0)
        cb = ContinuousCallback(slide_exit, stop)
        firstsol = ODESol(solve(prob, Vern6(), callback=cb), :slidemode)
        sols[1] = firstsol
    elseif h(v0, p0, t0) > 0
        prob = ODEProblem(f2, v0, (t0, tend), p0)
        cb = ContinuousCallback(newh, stop)
        firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
        sols[1] = firstsol
    else
        prob = ODEProblem(f1, v0, (t0, tend), p0)
        cb = ContinuousCallback(newh, stop)
        firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
        sols[1] = firstsol
    end
    while sols[end].solution.t[end] != tend # if the last of sols's time ≠ tend
        v1 = sols[end].solution.u[end]
        t1 = sols[end].solution.t[end]
        if sols[end].state == :slidemode # if the last solution is in slidemode, we have to determine which vector field it enter
            if abs(dot(∇h(v1, p0, t1), f2(v1, p0, t1))) < 1e-8 && dot(∇h(v1, p0, t1), f1(v1, p0, t1)) > 1e-5
                prob = ODEProblem(f2, v1, (t1, tend), p0)
                cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
                firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
                append!(sols, [firstsol])
            elseif abs(dot(∇h(v1, p0, t1), f1(v1, p0, t1))) < 1e-8 && dot(∇h(v1, p0, t1), f2(v1, p0, t1)) < -1e-5
                prob = ODEProblem(f1, v1, (t1, tend), p0)
                cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
                firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
                append!(sols, [firstsol])
            else
                break
                print("The solution is not defined on state $v1 and time $t1")
            end
        else
            if dot(∇h(v1, p0, t1), f1(v1, p0, t1)) * dot(∇h(v1, p0, t1), f2(v1, p0, t1)) > 0 # crossing h=0
                if dot(∇h(v1, p0, t1), f1(v1, p0, t1)) > 0
                    prob = ODEProblem(f2, v1, (t1, tend), p0)
                    cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
                    firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
                    append!(sols, [firstsol])
                else
                    prob = ODEProblem(f1, v1, (t1, tend), p0)
                    cb = ContinuousCallback(newh, stop, repeat_nudge=0.2)
                    firstsol = ODESol(solve(prob, Vern6(), callback=cb), :regularmode)
                    append!(sols, [firstsol])
                end
            else # entering h=0
                prob = ODEProblem(sv, v1, (t1, tend), p0)
                cb = ContinuousCallback(slide_exit, stop)
                firstsol = ODESol(solve(prob, Vern6(), callback=cb), :slidemode)
                append!(sols, [firstsol])
            end
        end
    end
    sols
end

```

The main function is `fillippov_solve`, and its output is a vector consist of `ODESol`.  
Now we can test this function using the system I gave before. First, define two vector fields:

```julia
function F1(x)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    A * (α / (1 + γ * abs(x)) + λ + η * x^2)
end

function F2(x)
    A, γ, λ, η, α = 10.0, 3.0, 0.1, 0.01, 0.3
    -A * (α / (1 + γ * abs(x)) + λ + η * x^2)
end

function vec1(x, p, t)
    SVector(x[2], -p[1] * x[2] - p[2] * x[1] + F1(x[2] - p[3]) + p[4] * cos(p[5] * t)) # p=[c,k,vs,ϵ,ω]
end

function vec2(x, p, t)
    SVector(x[2], -p[1] * x[2] - p[2] * x[1] + F2(x[2] - p[3]) + p[4] * cos(p[5] * t)) # p=[c,k,vs,ϵ,ω]
end

```

And then define the hypersurface and its grad:

```julia
function H(x, p, t)
    x[2] - p[3]
end

function ∇H(x, p, t)
    SVector(0.0, 1.0)
end

```

Note that we have to make sure that when H\>0 the vector field is `vec2` and when H\<0 the vector field is `vec1`.  
Now we can use `fillippov_solve` to solve this ODE. Just run

```julia
sol = fillippov_solve(vec1, vec2, H, ∇H, SVector(1.0, 3.0), SVector(0.2, 1, 1, 0.2, 1.067), 0.0, 20.0)

```

You will get a vector.

```julia
sol[2].state

```

gives

```julia
:slidemode

```

```julia
sol[2].solution

```

gives a ODE solution

```julia
retcode: Terminated
Interpolation: specialized 
6th order lazy interpolation
t: 3-element Vector{Float64}:
 0.5146478247426896
 2.117126414261962
 2.117126414261962
u: 3-element Vector{StaticArraysCore.SVector{2, Float64}}:
 [2.07049529563982, 1.0000000000000895]
 [3.6729738851570852, 1.000000000000007]
 [3.6729738851570852, 1.000000000000007]

```

Since we define a new vector field `sv` on the switch surface when the vector fields are squeezing towards the switch surface, the solution is now stable.  
We can define a simple function to plot all the solutions:

```julia
function myplot(sols)
    plotlyjs()
    figure = plot()
    for i in sols
        plot!(i.solution, idxs=(1,2))
    end
    figure
end

```

Now we run

```julia
myplot(sol)

```

it gives this picture:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/f/4fac2af124818f8e44d2898ecf8c789d6ad9fe70.png)  
Since I also consider the case of crossing, the codes should work on the usual piecewise smooth case.