# Best strategies to solve an ODE with piecewise functions using OrdinaryDiffEq
**URL:** https://discourse.julialang.org/t/best-strategies-to-solve-an-ode-with-piecewise-functions-using-ordinarydiffeq/87750
**Category:** Numerics
**Tags:** question, package, diffeq
**Created:** [September 24, 2022, 4:23pm UTC](https://discourse.julialang.org/t/best-strategies-to-solve-an-ode-with-piecewise-functions-using-ordinarydiffeq/87750 "2022-09-24T16:23:25Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@RayleighLord](https://discourse.julialang.org/u/RayleighLord)
#### Post date: [September 24, 2022, 4:23pm UTC](https://discourse.julialang.org/t/best-strategies-to-solve-an-ode-with-piecewise-functions-using-ordinarydiffeq/87750/1 "2022-09-24T16:23:26Z")
I want to solve an ODE problem involving some piecewise functions and I would like to ask for advice on what could be the best strategy to solve them.
Physically, the ODEs I want to solve represent a contact element sticking and sliding according to Coulomb’s friction law and with variable normal load. The code for the forces is the following
```julia
using OrdinaryDiffEq, Plots
normal(xn) = xn > 0.0 ? xn : 0.0
function tangential!(w, xt, N)
w_init = w[1]
if (N > 0.0) && (abs(xt - w_init) < N)
T = xt - w_init
return T
elseif (N > 0.0) && (abs(xt - w_init) ≥ N)
sg = sign(xt - w_init)
@. w = xt - sg * N
T = sg * N
return T
else
@. w = xt
T = 0.0
return T
end
end
```
where `w` is a mutating 1-element vector that takes into account the “memory” of the friction state. Then, the ODE is defined and solved as
```julia
function f!(du, u, p, t)
w = p
N = normal(u[2])
T = tangential!(w, u[1], N)
du[1] = u[3]
du[2] = u[4]
du[3] = 0.75sin(t) - T
du[4] = -N + 0.5
end
begin
u0 = [0.0, 2.0, 0.0, 0.0]
p = [0.0]
tspan = (0, 40π)
prob = ODEProblem(f!, u0, tspan, p)
sol = solve(prob, Midpoint())
plot(sol)
end
```
where the Midpoint method has been used. The idea is to find the best strategy to solve this simple problem to, in the end, use multiple of these contacts in a more complex setup. I have thought of the following strategies using the OrdinaryDiffEq package:
1. Low order explicit methods
2. High order explicit methods using callbacks
3. Implicit methods
The use of low order explicit methods is easy to implement since the code above just runs perfectly.
High order explicit methods cannot be used directly, since our function does not have continuous derivatives. Therefore, the only way I can see the use of a high order scheme is to add callbacks to the problem. Using callbacks I could stop the integration when we move from one section to the piecewise function to another and, since each section is infinitely differentiable, the use of high order schemes could give good results.
Finally, I found that I could only make implicit methods work by implementing them myself specifically suited for this problem. If I try to execute an implicit scheme with my system I get errors when it tries to automatically differentiate, which I suspect is expected due to the mutations that are going on in my function. Therefore, I do not know if using an implicit method is viable in my problem with the interface of DifferentialEquations and I would like if someone could give some advice on this issue.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 24, 2022, 5:47pm UTC](https://discourse.julialang.org/t/best-strategies-to-solve-an-ode-with-piecewise-functions-using-ordinarydiffeq/87750/2 "2022-09-24T17:47:08Z")
You probably want to look into the methods for dynamical ODE problems as this is a partitioned (second order) ODE.
> [@RayleighLord](#):
>
> Finally, I found that I could only make implicit methods work by implementing them myself specifically suited for this problem. If I try to execute an implicit scheme with my system I get errors when it tries to automatically differentiate, which I suspect is expected due to the mutations that are going on in my function. Therefore, I do not know if using an implicit method is viable in my problem with the interface of DifferentialEquations and I would like if someone could give some advice on this issue.
Set autodiff=false, like `Rodas5P(autodiff=false)`.
[https://diffeq.sciml.ai/stable/basics/faq/#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods](https://diffeq.sciml.ai/stable/basics/faq/#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods)
Though with your equation I don’t see why one would use an implicit method. Nothing sticks out to me as stiff there.
> [@RayleighLord](#):
>
> High order explicit methods cannot be used directly, since our function does not have continuous derivatives. Therefore, the only way I can see the use of a high order scheme is to add callbacks to the problem. Using callbacks I could stop the integration when we move from one section to the piecewise function to another and, since each section is infinitely differentiable, the use of high order schemes could give good results.
See [https://diffeq.sciml.ai/stable/basics/faq/#My-Problem-Has-Discontinuities-and-is-Unstable-/-Slow](https://diffeq.sciml.ai/stable/basics/faq/#My-Problem-Has-Discontinuities-and-is-Unstable-/-Slow) or
> [@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…
For your case, you callbacks and a decreased `interp_points` probably makes sense. Did you try it?
---
### Author:  [@RayleighLord](https://discourse.julialang.org/u/RayleighLord)
#### Post date: [September 25, 2022, 10:24am UTC](https://discourse.julialang.org/t/best-strategies-to-solve-an-ode-with-piecewise-functions-using-ordinarydiffeq/87750/3 "2022-09-25T10:24:55Z")
> [@ChrisRackauckas](#):
>
> Set autodiff=false, like `Rodas5P(autodiff=false)`.
I tried doing this, and it works fine, though it gives very different results. I do not know if, when solving the system implicitly, the mutations I have in my system can induce some error by getting some unexpected mutation.
You are right that this equation is not stiff by itself. This system will actually be part of a much larger system trying to model the oscillations of a structure with friction at its contact interfaces. By our experience with previous algorithms in Python, solving the whole system implicitly gave as the fastest integration times, while using IMEX methods had a more restrictive `dt` so that is why I am interested in the possibility of fully solving the system implicitly.
> [@ChrisRackauckas](#):
>
> You probably want to look into the methods for dynamical ODE problems as this is a partitioned (second order) ODE.
Thank you for pointing this out. I tried to specify the problem as a dynamical ODE and got a x3 faster solving times using the same numerical scheme. Interestingly, defining it as a SecondOrderODEProblem gives just slightly slower times (0.1 ms seconds more).
> [@ChrisRackauckas](#):
>
> For your case, you callbacks and a decreased `interp_points` probably makes sense. Did you try it?
I haven’t tried it, since I think it could get a little messy once you have several elements that behave in this piecewise manner. Inspired by your link above, I have tried to find a smooth version of the normal and tangential forces, but it does not seem to be easy since one piecewise function depends on the other. Now I have
```julia
normal_smooth(xn) = 1 / 65 * log(1 + exp(65xn))
```
for the normal force and
```julia
function tangential_smooth!(w, xt, N)
w_init = w[1]
if (N > 0.0) && (abs(xt - w_init) < N)
T = xt - w_init
return T
elseif (N > 0.0) && (abs(xt - w_init) ≥ N)
sg = sign(xt - w_init)
@. w = xt - sg * N
T = sg * N
return T
end
end
```
for the tangential force, where the last condition of the if-else block was removed but I do not find an easy way to approximate the rest of `tangential_smooth!` to a better behaved function. Using this I get very close results to my original setup but it probably would be optimal to have this other function as a differentiable function.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 28, 2022, 11:45am UTC](https://discourse.julialang.org/t/best-strategies-to-solve-an-ode-with-piecewise-functions-using-ordinarydiffeq/87750/4 "2022-09-28T11:45:54Z")
> [@RayleighLord](#):
>
> You are right that this equation is not stiff by itself. This system will actually be part of a much larger system trying to model the oscillations of a structure with friction at its contact interfaces. By our experience with previous algorithms in Python, solving the whole system implicitly gave as the fastest integration times, while using IMEX methods had a more restrictive `dt` so that is why I am interested in the possibility of fully solving the system implicitly.
Sounds like it’s worth trying RKC methods like ROCK2/ROCK4?
> [@RayleighLord](#):
>
> By our experience with previous algorithms in Python
Note that Python’s slowness tends to cause a bias towards “more vectorized” methods, giving a bias towards implicit methods simply because they use more linear algebra which has comparatively less overhead than BLAS 1 operations. I wouldn’t extrapolate too far.
---
### Author:  [@RayleighLord](https://discourse.julialang.org/u/RayleighLord)
#### Post date: [September 28, 2022, 2:59pm UTC](https://discourse.julialang.org/t/best-strategies-to-solve-an-ode-with-piecewise-functions-using-ordinarydiffeq/87750/5 "2022-09-28T14:59:26Z")
> [@ChrisRackauckas](#):
>
> Sounds like it’s worth trying RKC methods like ROCK2/ROCK4?
I actually did not know about the existence of this methods, so thank you for pointing it out. I will try them out to see what I can get.