Best strategies to solve an ODE with piecewise functions using OrdinaryDiffEq

[ChrisRackauckas]

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

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 or

Handling Instability When Solving ODE Problems

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?