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:
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.