RKN1210 12th/10th order Runge-Kutta-Nyström integrator
RKN1210() is a 12th/10th order variable-step numerical integrator for second-order ordinary differential equations of the form
y'' = f(t, y) (1)
with initial conditions
y(t0) = y0
y'(t0) = yp0 (2)
This second-order differential equation is integrated with a Runge-Kutta-Nyström method using 17 function evaluations per step. RKN12(10) is a very high-order method, to be used in problems with *extremely* stringent error tolerances.
The RKN-class of integrators is especially suited for problems of type (1). Compared to a classic Runge-Kutta integration scheme, the same accuracy can be obtained with fewer function evaluations. Also, it has been shown in various studies that this particular integration method is overall more efficient than (symplectic) multi-step or extrapolation methods that give the same accuracy.
RKN1210's behavior is very similar MATLAB's ODE-integrator suite; you can set options via ODESET, and input and output values are also practically the same.
Both output functions and event functions are fully supported.
The construction of RKN12(10) is described in
High-Order Embedded Runge-Kutta-Nyström Formulae
J. R. DORMAND, M. E. A. EL-MIKKAWY, AND P. J. PRINCE
IMA Journal of Numerical Analysis (1987) 7, 423-430
Coefficients obtained from
These are also available in any format on request to these authors.
If you find this work useful, please consider a small donation:
@Adam: you're right. I've made a few changes that correct the problem with the demo.
Runge/Kutta/Nyström integrators can't deal with equations involving the first derivative. They are a specialization of the more general Runge/Kutta integrators, exploiting an efficiency gain that is only possible if there is no first derivative to consider. For your case, use ODE113.
Also, is it possible to use this integrator on a system of ODEs that are both second order and first order? For example, in the motion of a satellite that that both influenced by gravity (second order ODE) and expending mass (first order mass flow rate).
Looks like there's a problem with the demo in R2018b. It never terminates. It runs the event detection demo, gets just to the point where it would normally stop (including plotting), and the hangs. Any ideas, Rody?
@Sammy Olivaw I uploaded today, so should be tomorrow or so.
@Rody Oldenhuis Thank you! This integrator is really good. Any idea when the update will be up?
@Sammy Olivaw: In principle, yes. There was a bug relating to this issue, but this has been fixed with the pending update.
Does it support integrating backward in time? (ex: tSpan = linspace(0, -10);)
@Tarek: thank you for the heads up, it has been corrected and a fixed version will be up as soon as they've approved it!
There seems to be a bug when using tspan as a vector to display results at certain time steps. An error message is generated for the index of tout as
Undefined function or variable "index".
Error in ==> rkn1210>finalize at 768
tout = tout (1:index,:);
The codes works fine with specifying tspan as [t0 tf]
- implemented NormControl
See GitHub for way more details (255 char limit here)
fixed divide-by-zero error for isolated cases where f turns all-zero
- Misc. bugfixes (minor)
- Found & removed some stray debugging code
- Fixed 2 bugs (thanks everyone!)
Improvements regarding performance and memory usage.
bugfix for missing "index" variable on numel(tspan)>2 as noted by Tarek
- Removed bug inherited from RKN86; the error estimates there were all wrong. Now RKN1210 integrates at least 50x faster!
- Fixed 2 bugs: 1) exitflag = 3 was never assigned or described in the argument [output] 2) intermediate results were erased by event functions and only the event's zero was returned
- Interactive demonstration included
Inspired by: rkn86