I am solving an ode (~40 variables) for many different initial conditions. Mostly TRBDF2 gives fastest results. But sometimes it shows that the system is unstable. In such cases I switched to Rodas5P or RadauIIA5 which gave results but was much slower.
What criteria could I define to catch a solution going unstable and switch solver to Radau and then switch back once the region of extreme stiffness has passed? I tried integrator.dt<1e-10 in the choice function of composite algorithm as done in the documented example but that doesnāt work.
On a side note, do I mess up the numerical errors if I do such switch?
What do you mean by ābut that doesnāt work.ā?
Sorry, I meant that it doesnāt switch. It still gives āunstableā. While when called with just radau (instead of the composite) it works fine.
interestingā¦ How long is your tspan? Any chance you can post an example? Iām pretty sure that for recent versions of OrdinaryDiffEq it can only be unstable if dt<eps(t)
The time span is something like 1e3. Posting an example would be difficult without posting the whole code. There are a lot of things happening to calculate the du vector, so itād be hard to come up with a MWE.
Interestingly, FBDF and QNDF always failed for me irrespective of the initial conditions. I am running with reltol=1e-6. (And same abstol)
I am able to get everything work when I increase the tolerance.
To rephrase my question, what I want to implement is something like:
#=
Default solver = Tsit5()
If stiff, solver = TRBDF2(autodiff=false)
If going unstable, solver = RadauIIA5(autodiff=false)
=#
myalg = CompositeAlgorithm((AutoTsit5(TRBDF2(autodiff=false)),
RadauIIA5(autodiff=false)),choice_function)
What should the choice function be to catch when solution goes unstable and return to previous algo when not?
yeah, the thing thatās surprising is that theoretically, the thing that youāre doing should work. Iām somewhat unsure how you are getting instability without dt going to 0. Have you checked to see whether the integrator switches and then goes unstable or whether it never actually switches? Another thing worth considering is whether you have discontinuities that could be solved by callbacks.
How do I check if the solver switched or not?
It gets unstable because ādt was forced below floating point epsilonā.
I define my choice_function like this:
choice_function(integrator) = (Int(integrator.dt < 1e-9) + 1)
Now, calling the solution with solver:
Autotsit5(TRBDF2) ā unstable
Autotsit5(RadauIIA5) ā stable, but slow
RadauIIA5 ā stable, but slow
(Autotsit5(TRBDF2),RadauIIA5),choice_function ā unstable
(TRBDF2,RadauIIA5),choice_function ā unstable
sol.alg_choice will tell you which alg was used at each timestep. It might be worth trying AutoTsit5(Rosenbrock23) since it seem like higher order solvers arenāt working well here.
On checking sol.alg_choice, it never switched (from trbdf2 to radau).
Iāll check with rosenbrock23.
Edit: Rosenbrock23 seems to be working. But it is also slower than TRBDF2. So again leaving me with the same issue.
Also, what time does it fail? If it never switches that might just mean you need to switch at a bigger dt (like 1e-6?)
Another thing thatās worth checking is solving it with TRBDF2, finding the last step before dt<1e-9 and seeing if another solver can solve from there. It is possible that your equations are such that by the time TRBDF runs into trouble, the system is inherently headed towards instability.
I am feeding in the point where trbdf2 failed for all further calculations.
I donāt see dt<1e-9 in any output points so maybe integrator.dt only updates after a successfull step?
I am afraid that if I switch at a bigger dt, then wouldnāt it always go to radau whenever it does stiff? Iāll try doing it though.
TRBDF already can handle pretty massive stiffness. You also should be able to switch back if the dt raises sufficiently high.
Iāve run it now with a switch at dt<1e-5. It doesnāt say unstable yet, so hopefully itāll run to completion. Iāll then check alg_choice to see when the switches happened.
Thanks
Just an update. Switching at dt<1e-5 does the job. The algorithms switches back to TRBDF2 later and reduces the computation time.
The problem was there was no successful step where dt<1e-9 and hence integrator.dt never crossed that threshold (successfully) before dropping below floating point epsilon and giving the unstable error.
It might be useful to be able to access the integratorās attempted dt to handle sharp discontinuities that cannot be handled via callbacks.
1 Like