# Switching stiff solvers when one is unstable

**URL:** <https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609>\
**Category:** Modelling & Simulations\
**Tags:** differentialequation\
**Created:** [September 19, 2024, 2:56pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609 "2024-09-19T14:56:41Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 19, 2024, 2:56pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/1 "2024-09-19T14:56:42Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 19, 2024, 3:44pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/2 "2024-09-19T15:44:18Z")

</div>

What do you mean by “but that doesn’t work.”?

---

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 19, 2024, 4:01pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/3 "2024-09-19T16:01:50Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 19, 2024, 4:04pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/4 "2024-09-19T16:04:41Z")

</div>

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)`

---

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 19, 2024, 4:32pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/5 "2024-09-19T16:32:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 19, 2024, 6:31pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/6 "2024-09-19T18:31:16Z")

</div>

To rephrase my question, what I want to implement is something like:

```julia
#=
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?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 19, 2024, 6:40pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/7 "2024-09-19T18:40:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 19, 2024, 6:52pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/8 "2024-09-19T18:52:12Z")

</div>

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:

```julia
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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 19, 2024, 6:55pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/9 "2024-09-19T18:55:55Z")

</div>

`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.

---

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 19, 2024, 7:05pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/10 "2024-09-19T19:05:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 19, 2024, 7:13pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/11 "2024-09-19T19:13:53Z")

</div>

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?)

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 19, 2024, 7:21pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/12 "2024-09-19T19:21:27Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 19, 2024, 7:28pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/13 "2024-09-19T19:28:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 19, 2024, 7:40pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/15 "2024-09-19T19:40:04Z")

</div>

TRBDF already can handle pretty massive stiffness. You also should be able to switch back if the dt raises sufficiently high.

---

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 19, 2024, 7:44pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/16 "2024-09-19T19:44:00Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![pulk\_jain](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pulk_jain/32/204900_2.png) [@pulk\_jain](https://discourse.julialang.org/u/pulk_jain)\
**Post date:** [September 25, 2024, 6:03pm UTC](https://discourse.julialang.org/t/switching-stiff-solvers-when-one-is-unstable/119609/17 "2024-09-25T18:03:59Z")

</div>

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.