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.

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.

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.