# Help with hyper-stiff ODE system

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, along with initial conditions and expected results.

At the expense of cross-posting, the corresponding issue in DifferentialEquations.jl is this.

Cheers,
James

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

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

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.

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

Yes a lot of the commented out integrators aren’t listed here, but `AutoVern9(Rodas5())` is listed here. 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`.