12th Order Runge-Kutta-Nystrom Integrator (RKN1210)

This is RKN1210 12th/10th order Runge-Kutta-Nystrom integrator. RKN1210() is a 12th/10th order numerical integrator for ordinary differential equations of the special form (1)

y” = f(t, y)                  (1)

With initial conditions (2)

y (t0) = y0             (2)
y'(t0) = yp0

This second-order differential equation is integrated with a Runge-Kutta-Nystrom method, with 17 function evaluations per step. The RKN-class of integrators is especially suited for this purpose, since compared to a classic Runge-Kutta integration scheme the same accuracy can be obtained with less function evaluations.

This RKN12(10) method is a very high-order method, to be used in problems  with *extremely* stringent error tolerances. As the name implies, the error should be less than O(h^13). In verious studies, it has been shown  that this particular integration technique is overall more efficient for ODE’s of the form (1) than multi-step or extrapolation methods that give the same accuracy.

RKN1210 USAGE:
[t, y, yp] = RKN1210(funfcn, tspan, y0, yp0)
[t, y, yp] = RKN1210(funfcn, tspan, y0, yp0, options)
[t, y, yp, exitflag, output] = RKN1210(…)
[t, y, yp, TE, YE, YPE, IE, exitflag, output] = RKN1210(…)
sol = RKN1210(funfcn, tspan, …)
RKN1210()’s behavior is very similar MATLAB’s ODE-integrator suite.

INPUT ARGUMENTS:

• funfcn – definition of the second-derivative function f(t, y) (See (1)). It should accept a scalar value [t] and a column vector [y] which has the same number of elements as the initial values [y0] and [dy0] provided.
• tspan – time interval over which to integrate. It can be a two-element vector, in which case it will be interpreted as an interval. In case [tspan] has more than two elements, the integration is carried out to all times in [tspan]. Only the values for those times are then returned in [y] and [yp].
•  y0, yp0 – initial values as in (2). Both should contain the same number of elements.
• options – options structure, created with ODESET(). Used options are MaxStep, InitialStep, AbsTol, Stats, Event, OutputFcn, OutputSel, and Refine. See the help for ODESET() for more information.

How to use Event and/or Output functions is described in the documentation on ODESET(). There is one difference: RKN1210() now also passes the first derivative [yp] at each step as an argument:

•  status = outputFcn(t, y, yp, flag)
  [value, isterminal, direction] = event(t, y, yp)

where [t] is scalar, and [y] and [yp] are column vectors, as with f(t,y).

OUTPUT ARGUMENTS:

• t, y, yp – The approximate solutions for [y] and [y’] at times [t]. All are concatenated row-wise, that is
  t = N-by-1
  y = N-by-numel(y0)
  y’ = N-by-numel(y0)
  with N the number of sucessful steps taken during the integration.
• exitflag – A scalar value, indicating the termination conditions of the integration:
  -2: a non-finite function value was encountered during the integration (INF of NaN); the integration was stopped.
  -1: the step size [h] fell below the minimum acceptable value at some time(s) [t]; results may be inaccurate.
  0: nothing was done; initial state.
  +1: sucessful integration, normal exit.
  +2: integration was stopped by one of the output functions.
  +3: One or more events were detected, and their corresponding [isterminal] condition also evaluated to [true].
  TE,YE, – These arguments are only returned when one or more event YPE,IE functions are used. [TE] contains the times at which events were detected. [YE] and [YPE] lists the corresponding values of the solution [y] and the first derivative [yp] at these times. [IE] contains indices to the event-functions with  which these events were detected. Use a smaller value for AbsTol (in [options]) to increase the accuracy of these roots when required.
• output – structure containing additional information about the integration. It has the fields:
  output.h step size (sucesful steps only) at each time [tn]
  output.rejected amount of rejected steps
  output.accepted amount of accepted steps
  output.delta estimate of the largest possible error at each time [tn]
  output.message Short message describing the termination conditions

Note that these fields contain the information of ALL steps taken, even for cases where [tspan] contains more than 2 elements.

• sol – a structure that can be passed to deval() in order to evaluate the solution at any point in the interval [tspan]. The structure [sol] always includes these fields:
  sol.x Steps chosen by the solver.
  sol.y Each column sol.y(:,ii) contains the solution of the function at sol.x(ii).
  sol.yp Each column sol.yp(:,ii) contains the solution of the first derivative at sol.x(ii).
  sol.solver Solver name (‘rkn1210’).

If you specify the ‘Events’ option and events are detected, [sol] also includes these fields:

• sol.xe Points at which events, if any, occurred.
  sol.xe(end) contains the exact point of a terminal event, if any.
  sol.ye Solutions for the function corresponding to events in sol.xe.
  sol.ype Solutions for the derivative corresponding to events in sol.xe.
  sol.ie Indices into the vector returned by the function specified in the Events option. The values indicate  which event the solver detected.

The construction of RKN12(10) is described in:

High-Order Embedded Runge-Kutta-Nystrom 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:
http://www.tampa.phys.ucl.ac.uk/rmat/test/rknint.f
These are also available in any format on request to these authors.

Example on Using The Code:
Equations of motion for a circular orbit in 2D
f1 = @(t, y) -y/sqrt(y’*y)^3;
f2 = @(t, y) [y(3:4); -y(1:2)/sqrt(y(1:2).’*y(1:2))^3];

% RKN1210 with moderate accuracy setting
disp(‘RKN1210, with AbsTol = RelTol = 1e-6:’)
options = odeset(‘RelTol’, 1e-6, ‘AbsTol’, 1e-6);
tic
[t1, y1] = rkn1210(f1, [0 1000], [1; 0], [0; 1], options);
toc
disp(‘Maximum absolute error:’)
max( (y1(:,1)-cos(t1)).^2 + (y1(:,2)-sin(t1)).^2 )

This is how much ODE45 will have to be tuned to achieve similar accuracy
fprintf(‘\n\n’)
disp(‘Compare with ODE45, which needs AbsTol = RelTol = 1e-8:’)
options = odeset(‘RelTol’, 1e-8, ‘AbsTol’, 1e-8);
tic
[t2, y2] = ode45(f2, [0, 1000], [1; 0; 0; 1], options);
toc
disp(‘Maximum absolute error:’)
max( (y2(:,1)-cos(t2)).^2 + (y2(:,2)-sin(t2)).^2