Help with hyper-stiff ODE system

Specific Domains Modelling & Simulations
1

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

2

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

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

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

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

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.

2 Likes
3

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.

I’ll give the Exponential methods a go.