# Piecewise differential equations

**URL:** <https://discourse.julialang.org/t/piecewise-differential-equations/12187>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [July 5, 2018, 4:56pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187 "2018-07-05T16:56:41Z")\
**Posts on this page:** 1\
**Showing post:** 6

<div class="post-metadata">

**Author:** ![dawbarton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dawbarton/32/215461_2.png) [@dawbarton](https://discourse.julialang.org/u/dawbarton)\
**Post date:** [July 5, 2018, 11:26pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/6 "2018-07-05T23:26:21Z")

</div>

Here is a (mildly) interesting example x''+x'+x=\pm p\_1 where the sign of p\_1 changes when a switching manifold is encountered at x=p\_2. To make things more interesting, consider hysteresis in the switching manifold such that p\_2\mapsto -p\_2 whenever the switching manifold is crossed.

The code is relatively straightforward; the StaticArrays/SVector/MVector can be ignored, they are only for speed.

```julia
using OrdinaryDiffEq
using StaticArrays

f(x, p, t) = SVector(x[2], -x[2]-x[1]+p[1]) # x'' + x' + x = ±p₁
h(u, t, integrator) = u[1]-integrator.p[2] # switching surface x = ±p₂;
g(integrator) = (integrator.p .= -integrator.p) # impact map (p₁, p₂) = -(p₁, p₂)

prob = ODEProblem(f, # RHS
                  SVector(0.0, 1.0), # initial value
                  (0.0, 100.0), # time interval
                  MVector(1.0, 1.0)) # parameters
cb = ContinuousCallback(h, g)
sol = solve(prob, Vern6(), callback=cb, dtmax=0.1)

```

Then plot `sol[2,:]` against `sol[1,:]` to see the phase plane - a nice non-smooth limit cycle in this case.

Note that if you try to use interpolation of the resulting solution (i.e., `sol(t)`) you need to be very careful around the points that have a discontinuous derivative as the interpolant goes a little awry. That’s why I’ve used `dtmax=0.1` to get a smoother solution output in this case. (I’m probably not using the most appropriate integrator either but it’s the one that I was using in a previous piece of code that I copied-and-pasted 🙂)

If you are interested in the dynamics of piecewise differential equations there are plenty of gotchas to watch out for, particularly if your system is a Filippov system; effects like sliding require a bit more care than shown in the code above. (E.g., I’m not sure if DifferentialEquations.jl can handle a system that changes from an ODE to a DAE part way through solving, such as happens with the sliding vector fields that come from Filippov systems; ever tried it @ChrisRackauckas?)

---

_[View the full topic](https://discourse.julialang.org/t/piecewise-differential-equations/12187)._