# Simple approachable example highlighting the usage of a stiff ODE solver in DiffEq

**URL:** <https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, teaching\
**Created:** [August 12, 2023, 8:28am UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741 "2023-08-12T08:28:51Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [August 12, 2023, 8:28am UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/1 "2023-08-12T08:28:51Z")

</div>

Hi all, I will be presenting a Julia tutorial on an international PHD school, and OrdinaryDiffEq.jl will be one of the packages I will be highlighting. I want to demonstrate some of the various solvers of the ecosystem via approachable examples. I will showcase two examples: 1) using Tsit5, then going from Tsit5 to Vern9 for better performance when high accuracy is in demand. 2) An example where a stiff solver leads to a better (more stable / more performant) solution. I couldn’t find a simple stiff example in the online docs. Ideally I would like a low dimensional ODE where it makes a big impact to use a stiff solver vs. a non-stiff. Does anyone have an example in mind?

---

<div class="post-metadata">

**Author:** ![vettert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vettert/32/30599_2.png) [@vettert](https://discourse.julialang.org/u/vettert)\
**Post date:** [August 12, 2023, 10:28am UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/2 "2023-08-12T10:28:54Z")

</div>

> **[GitHub - SciML/DiffEqProblemLibrary.jl: A library of premade problems for...](https://github.com/SciML/DiffEqProblemLibrary.jl)**
>
> A library of premade problems for examples and testing differential equation solvers and other SciML scientific machine learning tools - GitHub - SciML/DiffEqProblemLibrary.jl: A library of premade...

Should contain suitable examples, which are even already preprogrammed 🙂

Specifically, the Robertson or the van der Pol examples may be what you are looking for.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 12, 2023, 12:05pm UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/3 "2023-08-12T12:05:38Z")

</div>

I tend to like user ROBER (Robertson).

[https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/ROBER/](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/ROBER/)

It is very clearly not possible to solve with explicit methods, just try it. And you can easily see the time scale separation from the coefficients being orders of magnitude apart.

[https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/VanDerPol/](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/VanDerPol/)

Van Der Pol is indeed a great example too because it has tunable stiffness. As you change `μ` you change the time scales and you can see the behavior get more and more smooth. The stuff benchmark we have is from Hairer’s book and uses `1e6`, while MATLAB has an example of nonstiff ODEs

> **[Solve Nonstiff ODEs
- MATLAB & Simulink](https://www.mathworks.com/help/matlab/math/solve-nonstiff-odes.html)**
>
> This page contains two examples of solving nonstiff ordinary differential equations using ode45.

which is Van Der Pol at `μ=1`. This highlights that the difference can be not just structural but also dependent on parameter values.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [August 12, 2023, 12:21pm UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/4 "2023-08-12T12:21:28Z")

</div>

I would say slow fast systems with relaxation oscillations are likely to show case this.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [August 12, 2023, 5:00pm UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/5 "2023-08-12T17:00:19Z")

</div>

> [@ChrisRackauckas](#):
>
> Van Der Pol is indeed a great example too because it has tunable stiffness. As you change `μ` you change the time scales and you can see the behavior get more and more smooth

The Van der Pol example is very appealing to me because the audience I am teaching to will likely be familiar with it, and they are also interested in timescale separation. However, I am having difficulty constructing an example where the stiff solver is “more appropriate”. I’ve cooked up:

```julia
function vanderpol_rule(u, μ, t)
    x, y = u
    dx = y
    dy = μ*((1-x^2)*y - x)
    return SVector(dx, dy)
end

fig = Figure()
axs = [Axis(fig[i, :]) for i in 1:2]
algs = (Tsit5(), Rodas5P()) # second is stiff solver

for (i, alg) in enumerate(algs)
    for μ in (1, 1e3, 1e6)
        prob = ODEProblem(vanderpol_rule, SVector(1.0, 1.0), (0.0, 10.0), μ)
        sol = solve(prob; alg, abstol = 1e-12, reltol = 1e-12, saveat = 0.01, maxiters = typemax(Int))
        t = 0:0.01:10
        lines!(axs[i], t, sol[1, :])
    end
end
fig

```

I was hoping that in this MWE the solutions using Tsit5 would fail or diverge but they don’t, provided that I use large enough `maxiters`. Is this really all there is to it? that the stiff solvers are just much more efficient in solving these problems and hence take less steps overall?

EdIT: output figure is:

![image](https://global.discourse-cdn.com/julialang/original/3X/9/2/92c17c025b570ce80425d7baace3aa9e983daae5.png)

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [August 12, 2023, 5:03pm UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/6 "2023-08-12T17:03:18Z")

</div>

Actually it gets even weirded: I started printing the number of rejected steps, `sol.stats.nreject` and saw that the stiff solver rejected much more steps:

Alg: Rodas5P, rejected steps: 1551

vs

Alg: Tsit5, rejected steps: 0

for same large `μ` and same tolerances.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [August 12, 2023, 5:05pm UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/7 "2023-08-12T17:05:16Z")

</div>

ah okay, but the Rodas solver is **much more performant** :

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/6/e6eb6977d0a0f0902bb7155dc809a4eab899f786.png)

thanks a lot @ChrisRackauckas , I’ll use the Van der Pol as an example!!!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 12, 2023, 5:14pm UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/8 "2023-08-12T17:14:38Z")

</div>

Yes, and you can see this for example in the Filament PDE example.

[https://docs.sciml.ai/SciMLBenchmarksOutput/html/MOLPDE/Filament.html](https://docs.sciml.ai/SciMLBenchmarksOutput/html/MOLPDE/Filament.html)

 ![Screenshot 2023-08-12 131227](https://global.discourse-cdn.com/julialang/original/3X/8/b/8ba74e601cadc5d6c9225b5fdb754698163ef5f8.png)

An explicit method can take a billion steps to try and solve it, it just takes 50x as long.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 12, 2023, 5:15pm UTC](https://discourse.julialang.org/t/simple-approachable-example-highlighting-the-usage-of-a-stiff-ode-solver-in-diffeq/102741/9 "2023-08-12T17:15:13Z")

</div>

Though note that some cases are about stability. ROBER is pretty hard (if possible) to make stable with an explicit method.