# 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:** 10
**Page:** 1
### Author:  [@francisco](https://discourse.julialang.org/u/francisco)
#### Post date: [July 5, 2018, 4:56pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/1 "2018-07-05T16:56:41Z")
Hello,
I am new at Julia, I would like to solve this system
\frac{dx}{dt} = k1y \
\frac{dy}{dt} = k2y+I
where k1 and k2 are constant parameters. however, I=0 when y,0 or Ky otherwise, where k is a constant value.
I followed the tutorial about ODE in The question is how it is possible to solve this piecewise differential equation in DifferentialEquations.jl however I did not find anything about piecewise functions
Thank in advance
Francisco
I
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 5, 2018, 5:24pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/2 "2018-07-05T17:24:46Z")
This is discussed here and I haven’t added it to the FAQ yet:
> [@Handling Instability When Solving ODE Problems](https://discourse.julialang.org/t/handling-instability-when-solving-ode-problems/9019/5):
>
> Thanks for the example. I will share as little as possible but just a snippet in order to explain what’s going on to others. This is a good teaching opportunity. I was able to find out that the problem is that, in the statement k9\*X3/(k10\*k11/(k12\*(1.-min(1.,k13+k14\*X5))) + X3) when X3 and X5 are dual numbers, the value of k11/(k12\*(1.-min(1.,k13+k14\*X5))) is NaN for the derivative components when the condition makes it the constant 1.0 (since the derivative of the minimum of in the denominat…
---
### Author:  [@francisco](https://discourse.julialang.org/u/francisco)
#### Post date: [July 5, 2018, 8:30pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/3 "2018-07-05T20:30:20Z")
Dear Chris,
Thank you for your answer. I know that a discretecallback would be the solution of my problem. however I can realize how to program it because in all example, as I far as understood,the discontinuity is given by some time point and in my case the parameter change is due to the actual value of the function itself. Could you help me a bit please?
Thank you very much
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 5, 2018, 9:02pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/4 "2018-07-05T21:02:22Z")
Just use a `ContinuousCallback` that triggers when the zeroing happens. Or use the smoothed version of the operation. I’ll see if I get a chance to do a detailed response sooner rather than later. I see you posted here as well:
>
---
### Author:  [@francisco](https://discourse.julialang.org/u/francisco)
#### Post date: [July 5, 2018, 9:10pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/5 "2018-07-05T21:10:52Z")
Thank you very much, all your help would be very appreciated! 😄 I am not such a good programer, I am going to read about callback functions
---
### Author:  [@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")
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?)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 6, 2018, 12:25am UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/7 "2018-07-06T00:25:57Z")
> [@dawbarton](#):
>
> 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.
That shouldn’t be a problem. It won’t interpolate over the derivative discontinuity. It has the discontinuity as one of its interval endpoints.
> [@dawbarton](#):
>
> 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
DAEs and ODEs solve the same way. If your system is being solved by a solver that already handles DAEs (i.e. Rosenbrock with mass matrices, IDA, etc.) I don’t see it causing a problem.
---
### Author:  [@dawbarton](https://discourse.julialang.org/u/dawbarton)
#### Post date: [July 6, 2018, 9:14am UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/8 "2018-07-06T09:14:00Z")
> [@ChrisRackauckas](#):
>
> That shouldn’t be a problem. It won’t interpolate over the derivative discontinuity. It has the discontinuity as one of its interval endpoints.
Unfortunately it seems to be a problem. Try my code above and then
```julia
t = linspace(0, 100, 10001)
solt = sol(t)
plot(solt[1,:], solt[2,:])
```
and you’ll see an extra kink in the solution around x=1, x'=0.7 as it switches between the different vector fields. I did this initially without restricting the time step and the results are even more pronounced (and confusing). I don’t know if it is because just because I’m abusing the event detection routines somehow.
> [@ChrisRackauckas](#):
>
> DAEs and ODEs solve the same way. If your system is being solved by a solver that already handles DAEs (i.e. Rosenbrock with mass matrices, IDA, etc.) I don’t see it causing a problem.
That’s a good point; I’ll have to try that some time.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 6, 2018, 12:55pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/9 "2018-07-06T12:55:43Z")
> [@dawbarton](#):
>
> and you’ll see an extra kink in the solution around x=1x=1, x′=0.7x’=0.7 as it switches between the different vector fields. I did this initially without restricting the time step and the results are even more pronounced (and confusing). I don’t know if it is because just because I’m abusing the event detection routines somehow.
Oh that’s a “bug” due to the Verners doing a lazy interpolant grow. I need to document why this is happening and offer a workaround flag.
```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, DP8(), callback=cb, dtmax=0.1)
t = linspace(0, 100, 10001)
solt = sol(t)
plot(solt[1,:], solt[2,:])
```
Let me edit the docs and put a fix into the coming v0.7 version. This is a very specific case that shows up only in the case where:
1. You use a method with a lazy interpolant (currently and probably forever, this is only the 4 Verner methods and `BS5`)
2. You change a parameter using a callback
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 6, 2018, 1:52pm UTC](https://discourse.julialang.org/t/piecewise-differential-equations/12187/10 "2018-07-06T13:52:37Z")
That laziness caveat is now mentioned in the docs ( [https://github.com/JuliaDiffEq/DiffEqDocs.jl/commit/eda9b90942ca84dbef1f7bcc8bbfcee73b537738](https://github.com/JuliaDiffEq/DiffEqDocs.jl/commit/eda9b90942ca84dbef1f7bcc8bbfcee73b537738) ) and the non-lazy mode has a PR which will merge when tests pass ( [https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/418](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/418) ). I still want to keep the lazy default since it has about half as many steps if you aren’t interpolating every step, so you really only want the extra cost in very specific cases (and if you’re using a continuous callback it’s not an extra cost anyways since you have to be interpolating).
Ugh I don’t like the details this adds though. This side note is invading this thread now