# Help with hyper-stiff ODE system

**URL:** <https://discourse.julialang.org/t/help-with-hyper-stiff-ode-system/22681>\
**Category:** Modelling & Simulations\
**Tags:** diffeq\
**Created:** [April 3, 2019, 8:43am UTC](https://discourse.julialang.org/t/help-with-hyper-stiff-ode-system/22681 "2019-04-03T08:43:41Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jcook](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcook/32/18211_2.png) [@jcook](https://discourse.julialang.org/u/jcook)\
**Post date:** [April 3, 2019, 8:43am UTC](https://discourse.julialang.org/t/help-with-hyper-stiff-ode-system/22681/1 "2019-04-03T08:43:41Z")

</div>

Dear All,

I’ve come against a pretty tricky problem that I can’t find an accurate and performant solution to in Julia. There is a closed source FORTRAN code that can solve this problem quickly and accurately. So at least it’s possible to get a good solution.

I’m hoping that the community can shed some light on it.

The ODE is defined [here](https://github.com/ukaea/decay-rate-ode-solver), along with initial conditions and expected results.

At the expense of cross-posting, the corresponding issue in DifferentialEquations.jl is [this](https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/444).

Cheers,  
James

---

<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:** [April 3, 2019, 9:28am UTC](https://discourse.julialang.org/t/help-with-hyper-stiff-ode-system/22681/2 "2019-04-03T09:28:54Z")

</div>

> <https://github.com/ukaea/decay-rate-ode-solver/blob/master/gist.jl#L79-L83>

That’s a little odd because it happens to not pick out a single one of the recommended algorithms for highly stiff equations… ImplicitMidpoint for example is not a good choice because it is symplectic and thus not L-stable. The autoswitch algorithm are relying on autoswitching which can have issues with asymtopically large stiffness. Instead, the candidate “standard algorithms” would be:

- Rosenbrock23
- Rodas5
- KenCarp4

Now if you want this exactness

[https://github.com/ukaea/decay-rate-ode-solver/blob/master/gist.jl#L97](https://github.com/ukaea/decay-rate-ode-solver/blob/master/gist.jl#L97)

You’ll need to set a tolerance on the ODE solver.

[https://github.com/ukaea/decay-rate-ode-solver/blob/master/gist.jl#L60](https://github.com/ukaea/decay-rate-ode-solver/blob/master/gist.jl#L60)

reltol=1e-3 and abstol=1e-6 isn’t going to cut it.

But I think the issue is that using a nonlinear ODE solver for a linear problem isn’t the best method to use anyways. If they are getting ~0 abstol difference (how do you know the exact solution? What this computed independently using a high precision matrix exponential?), they are likely specializing on the linear form. If that’s the case, make `A = DiffEqArrayOperator(mat)` and then use one of the exponential integrators.

[http://docs.juliadiffeq.org/latest/solvers/ode\_solve.html#Exponential-Methods-for-Linear-and-Affine-Problems-1](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Exponential-Methods-for-Linear-and-Affine-Problems-1)

That will be a Krylov method for the matrix exponential, and then its extension to nonlinear problems will be things like

[http://docs.juliadiffeq.org/latest/solvers/ode\_solve.html#Exponential-Propagation-Iterative-Runge-Kutta-Methods-(EPIRK)-1](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Exponential-Propagation-Iterative-Runge-Kutta-Methods-(EPIRK)-1)

Those would be the methods that give exact solutions on linear problems, but I wouldn’t expect that property on any of the methods you chose.

---

<div class="post-metadata">

**Author:** ![jcook](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jcook/32/18211_2.png) [@jcook](https://discourse.julialang.org/u/jcook)\
**Post date:** [April 3, 2019, 10:22am UTC](https://discourse.julialang.org/t/help-with-hyper-stiff-ode-system/22681/3 "2019-04-03T10:22:19Z")

</div>

Yes a lot of the commented out integrators aren’t listed [here](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Stiff-Problems-1), but `AutoVern9(Rodas5())` is listed [here](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Unknown-Stiffness-Problems-1). I was just experimenting large time steps etc.

I’m afraid I don’t know how the expected results were generated, hence this problem. I am told that results with tiny numbers like `10^-150` are not to be ignored.

Where is `DiffEqArrayOperator` ? I’ve found in `DiffEqArray` in `DifferentialEquations` but not `DiffEqArrayOperator`.

I’ll give the Exponential methods a go.