Repeated occurrence of Runge-Kutta-Nyström issues concerning their accuracy

[HenryLangner]

So far there where a few issues with the Runge-Kutta-Nyström solvers concerning their accuracy:

• https://github.com/SciML/DifferentialEquations.jl/issues/726
• https://github.com/SciML/DifferentialEquations.jl/issues/751
• https://github.com/SciML/OrdinaryDiffEq.jl/issues/1030
• https://github.com/SciML/OrdinaryDiffEq.jl/issues/1372

All Nyström solvers used, in these issues, are designed for 2nd order ODEs which are not dependent on the first derivative. However most of the issues feature this type of 2nd order ODE. The lack of accuracy may just boil down to the solvers not being designed for theses types of problems. Any ideas on how to prevent these issues in the future?

[ranocha]

Thanks for listing these issues! Could you please check whether RKN methods allowing a second derivative depending on the first derivative work for the problems described there?

@chrisrackauckas Would it make sense to create a new ODE problem type for this class of problems to avoid these issues in the future? RKN methods allowing a dependency on the first derivative can continue to use the old setup but the other ones would require a new problem type without possibility to access the first derivative du when computing the second derivative ddu? Or do you have another suggestion how to handle this?

[HenryLangner]

For the testing I used the Nystrom4 solver which allows a dependence on the first derivative. The code is attached.

The following table lists the difference between the numeric solution given by Tsit5, DPRKN6 and Nystrom4 for each issue in the last time step. Since Nystrom4 can only be used with fixed time steps I chose to use fixed time steps on the other solvers too, so adaptive is set false. Results for dt = 1e-04:

Issue	difference between Tsit5 - Nystrom4	difference between Tsit5 - DPRKN6
726	dx = 3.552713678800501e-15	dx = 0.00014262397589970277
726	dy = -1.4210854715202004e-14	dy = 4.2214684228270016e-5
726	x = 3.501554601825774e-11	x = 0.035219797821355314
726	y = -8.003553375601768e-11	y = 0.021683669002868555
751	u = 2.673914423212409e-10	u = 5.44898921289132
751	du = 1.2369127944111824e-10	du = 2.4798490935190785
1030	dx = 7.086646831623969e-5	dx = -0.49883487250072506
1030	dy = 1.388565170623579e-7	dy = -0.0018123043314393072
1030	x = 2.1609006668477093e-5	x = -1.5611792619641611
1030	y = 4.349838630445598e-7	y = -0.07033738943846785
1372	dx = 3.552713678800501e-15	dx = -1.3117096002694987e-5
1372	dy = -1.4210854715202004e-14	dy = -7.0716174730312e-5
1372	x = 3.501554601825774e-11	x = 0.00011958319965943431
1372	y = -8.003553375601768e-11	y = 0.009709477552178214

As seen above the solutions from Tsit5 and Nystrom4 are really close together, so Nystrom4 works fine.

In https://github.com/SciML/OrdinaryDiffEq.jl/issues/1372 and in https://github.com/SciML/DifferentialEquations.jl/issues/751 DPRKN6 with adaptive step size gives even worse solutions. From which we can conclude that fixed time steps (dt = 1e-04) favor accuracy in our case (at least for tolerances under 1e-12) That being said it is odd that ¸DPRKN6gives different results than Tsit5 and Nystrom4 do.

Digging deeper into it I find that plenty of the implemented Runge-Kutta-Nyström solvers, which are constructed for 2nd order ODEs without dependence on the first derivative, are inaccurate and show a similar error. Taking https://github.com/SciML/DifferentialEquations.jl/issues/751 for example. The analytical solution in the second argument yields 134.50945587649028 at t = 3.0. Solving it numerically (again with fixed timesteps dt=1e-04) with different RKN solvers we get

Solver	numerical solution at t=3.0	error at t=3.0
Nystrom4	134.50945587597616	5.1412e-10
ERKN4	134.4840611591785	0.025394717311769455
ERKN5	134.48406115917854	0.025394717311741033
ERKN7	134.48406115917854	0.025394717311741033
Nystrom4Velocity Independent	134.4840611591781	0.02539471731216736
Nystrom5Velocity Independent	134.48406115917854	0.025394717311741033
DPRKN4	134.4840611591785	0.025394717311769455
DPRKN5	134.48406162812904	0.025394248361237715
DPRKN6	134.4819611976302	0.027494678860080057
DPRKN6FM	134.48406115917854	0.025394717311741033
DPRKN8	134.4840611591785	0.025394717311769455
DPRKN12	134.48406115917854	0.025394717311741033

The error is the difference between the numerical and analytical solution at t=3.0 Note that Nystrom4 is the only solver which allows the 2nd order ODE to depend on the first derivative. All the other solvers listed are designed for first derivative independent 2nd order ODEs.

In comparison to Nystrom4 all other solvers are similarly inaccurate which indicates that it is not a DPRKN6 (or as seen in https://github.com/SciML/OrdinaryDiffEq.jl/issues/1030 a DPRKN12) issue rather that these solvers in general have trouble solving 2nd order ODEs which depend on the first derivative.

The code used, it is just slightly changed from the Issues listed: test_726.txt test_751.txt test_1030.txt test_1372.txt

[ranocha]

Thanks a lot for these tests! This really tells us that it may be too easy for users to accidentally use the wrong method. Let's see what @ChrisRackauckas (and other main developers) think about this.

@HenryLangner Could you please comment briefly in the related issues, explaining the problem briefly and pointing to this issues? This should make it easier for new users running into the same problems to identify the cause.

[ChrisRackauckas]

Yeah this problem is tricky. The reason why we've kept these a bit open before is because they are still convergent, just at a lower order. But that's true of all sorts of methods. With highly stiff equations for example, you have the phenomena of stage order, where for example most ESDIRK methods converge at order 2 in such cases even though they are nominally orders 2-5. Similarly in the SDE case, where if you use a method for additive noise in a non-additive noise case you get order 0.5 convergence instead of the order 1.5.

We've tried to make sure that these are documented as much as possible with the solvers (https://docs.sciml.ai/DiffEqDocs/stable/solvers/dynamical_solve/#Runge-Kutta-Nystr%C3%B6m-Integrators), though I can see how "Allows acceleration to depend on velocity." slipped through telling people that others don't allow it for higher order convergence. The methods still achieve first order convergence of course, and that means they are bad.

So what to do?

@ChrisRackauckas Would it make sense to create a new ODE problem type for this class of problems to avoid these issues in the future? RKN methods allowing a dependency on the first derivative can continue to use the old setup but the other ones would require a new problem type without possibility to access the first derivative du when computing the second derivative ddu? Or do you have another suggestion how to handle this?

We've tried to avoid that because it can be very difficult to support so many problem types hanging around. One easy thing that we could do is do a single call to the acceleration function with velocity set as NaN values, and then error with these methods if the derivative is NaN. That of course can have some weird caveats though. I was hoping to someday use compiler plugin tools and do these kinds of proofs via abstract interpretation, but you can then only do that to Julia codes so it's limiting. I have though about doing this kind of stuff at problem construction time and tagging traits for things like this, and is_autonomous (which would allow Rosenbrock methods to simplify), etc.

But maybe doing the new problem type is easiest. What would it be named? VelocityFreeSecondOrderODEProblem?

[ranocha]

We've tried to make sure that these are documented as much as possible with the solvers (https://docs.sciml.ai/DiffEqDocs/stable/solvers/dynamical_solve/#Runge-Kutta-Nystr%C3%B6m-Integrators), though I can see how "Allows acceleration to depend on velocity." slipped through telling people that others don't allow it for higher order convergence. The methods still achieve first order convergence of course, and that means they are bad.

I guess a first step could be to add explicit warnings that the second derivative must not depend on the first derivative to the docs, following the docstring updates @HenryLangner did

[ranocha]

Adding a type parameter to DynamicalODEProblem and SecondOrderODEProblem may be okay - either assuming that the second derivative depends on the first derivative by default or trying to check it as you described (for floating point types). This parameter could also be used to improve the efficiency of some methods since less computations need to be performed.