Integer, Intent (In)	::	neq, lrwork
Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), Intent (In)	::	tend
Real (Kind=nag_wp), Intent (Inout)	::	t, y(neq), yp(neq), ydp(neq), rwork(lrwork)
External	::	fcn

C Header Interface

#include <nag.h>

void	d02laf_ ( void (NAG_CALL fcn)(const Integer neq, const double t, const double y[], double f[]), const Integer neq, double t, const double tend, double y[], double yp[], double ydp[], double rwork[], const Integer lrwork, Integer ifail)

The routine may be called by the names d02laf or nagf_ode_ivp_2nd_rkn.

3 Description

Given the initial values

x, y_{1}, y_{2}, \dots, y_{neq}, y_{1}^{'}, y_{2}^{'}, \dots, y_{neq}^{'}

d02laf integrates a non-stiff system of second-order differential equations of the type

y_{i}^{''} = f_{i} (x, y_{1}, y_{2}, \dots, y_{neq}), i = 1, 2, \dots, neq,

from

x = t

x = tend

using a Runge–Kutta–Nystrom formula pair. The system is defined by fcn, which evaluates

f_{i}

in terms of

x

and

y_{1}, y_{2}, \dots, y_{neq}

, where

y_{1}, y_{2}, \dots, y_{neq}

are supplied at

x

There are two Runge–Kutta–Nystrom formula pairs implemented in this routine. The lower order method is intended if you have moderate accuracy requirements and may be used in conjunction with the interpolation routine d02lzf to produce solutions and derivatives at user-specified points. The higher order method is intended if you have high accuracy requirements.

In one-step mode the routine returns approximations to the solution, derivative and

f_{i}

at each integration point. In interval mode these values are returned at the end of the integration range. You select the order of the method, the mode of operation, the error control and various optional inputs by a prior call to d02lxf.

For a description of the Runge–Kutta–Nystrom formula pairs see Dormand et al. (1986a) and Dormand et al. (1986b) and for a description of their practical implementation see Brankin et al. (1989).

4 References

Brankin R W, Dormand J R, Gladwell I, Prince P J and Seward W L (1989) Algorithm 670: A Runge–Kutta–Nystrom Code ACM Trans. Math. Software 15 31–40

Dormand J R, El–Mikkawy M E A and Prince P J (1986a) Families of Runge–Kutta–Nystrom formulae Mathematical Report TPMR 86-1 Teesside Polytechnic

Dormand J R, El–Mikkawy M E A and Prince P J (1986b) High order embedded Runge–Kutta–Nystrom formulae Mathematical Report TPMR 86-2 Teesside Polytechnic

5 Arguments

1: $fcn$ – Subroutine, supplied by the user. External Procedure

fcn must evaluate the functions

f_{i}

(that is the second derivatives

y_{i}^{''}

) for given values of its arguments

x

y_{1}, y_{2}, \dots, y_{neq}

The specification of fcn is:

Fortran Interface

Subroutine fcn (

neq, t, y, f)

Integer, Intent (In)	::	neq
Real (Kind=nag_wp), Intent (In)	::	t, y(neq)
Real (Kind=nag_wp), Intent (Out)	::	f(neq)

C Header Interface

void	fcn (const Integer neq, const double t, const double y[], double f[])

1: $neq$ – Integer Input: On entry: the number of differential equations.
2: $t$ – Real (Kind=nag_wp) Input: On entry: $x$ , the value of the argument.
3: $y (neq)$ – Real (Kind=nag_wp) array Input: On entry: $y_{i}$ , for $i = 1, 2, \dots, neq$ , the value of the argument.
4: $f (neq)$ – Real (Kind=nag_wp) array Output: On exit: the value of $f_{i}$ , for $i = 1, 2, \dots, neq$ .

fcn must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02laf is called. Arguments denoted as Input must not be changed by this procedure.

Note: fcn should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02laf. If your code inadvertently does return any NaNs or infinities, d02laf is likely to produce unexpected results.

2: $neq$ – Integer Input

On entry: the number of second-order ordinary differential equations to be solved by d02laf. It must contain the same value as the argument neq used in a prior call to d02lxf.

Constraint:

neq \geq 1

3: $t$ – Real (Kind=nag_wp) Input/Output

On entry: the initial value of the independent variable

x

Constraint:

t \neq tend

On exit: the value of the independent variable, which is usually tend, unless an error has occurred or the code is operating in one-step mode. If the integration is to be continued, possibly with a new value for tend, t must not be changed.

4: $tend$ – Real (Kind=nag_wp) Input

On entry: the end point of the range of integration. If

tend < t

on initial entry, integration will proceed in the negative direction. tend may be reset, in the direction of integration, before any continuation call.

5: $y (neq)$ – Real (Kind=nag_wp) array Input/Output

On entry: the initial values of the solution

y_{1}, y_{2}, \dots, y_{neq}

On exit: the computed values of the solution at the exit value of t. If the integration is to be continued, possibly with a new value for tend, these values must not be changed.

6: $yp (neq)$ – Real (Kind=nag_wp) array Input/Output

On entry: the initial values of the derivatives

y_{1}^{'}, y_{2}^{'}, \dots, y_{neq}^{'}

On exit: the computed values of the derivatives at the exit value of t. If the integration is to be continued, possibly with a new value for tend, these values must not be changed.

7: $ydp (neq)$ – Real (Kind=nag_wp) array Input/Output

On entry: must be unchanged from a previous call to d02laf.

On exit: the computed values of the second derivative at the exit value of t, unless illegal input is detected, in which case the elements of ydp may not have been initialized. If the integration is to be continued, possibly with a new value for tend, these values must not be changed.

8: $rwork (lrwork)$ – Real (Kind=nag_wp) array Communication Array

This must be the same argument rwork as supplied to d02lxf. It is used to pass information from d02lxf to d02laf, and from d02laf to both d02lyf and d02lzf. Therefore, the contents of this array must not be changed before the call to d02laf or calling either of the routines d02lyf and d02lzf.

9: $lrwork$ – Integer Input

On entry: the dimension of the array rwork as declared in the (sub)program from which d02laf is called.

This must be the same argument lrwork as supplied to d02lxf.

10: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

−1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

−1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

−1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

−1

is recommended since useful values can be provided in some output arguments even when

ifail \neq 0

on exit. When the value $- 1$ or $1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

−1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

Note: in some cases d02laf may return useful information.

$ifail = 1$: On entry, $lrwork = ⟨ value ⟩$ , but in the previous call to the setup routine $lrwork = ⟨ value ⟩$ .

On entry, $neq = ⟨ value ⟩$ , but in the previous call to the setup routine $neq = ⟨ value ⟩$ .

On entry, $t = ⟨ value ⟩$ and $tend = ⟨ value ⟩$ .
Constraint: $t \neq tend$ .

tend has been reset such that the direction of integration is reversed.

The previous call to the setup routine resulted in the error exit $⟨ value ⟩$ .

The setup routine d02lxf has not been called.

$ifail = 2$: The maximum number of steps, $⟨ value ⟩$ , has been attempted.

$ifail = 3$: To satisfy the accuracy requirements the step size, $⟨ value ⟩$ , at $t = ⟨ value ⟩$ , is too small for the machine precision.

$ifail = 4$: Two successive errors detected at the current value of $t = ⟨ value ⟩$ .

$ifail = 5$: Inefficiency detected in integrating exactly to values of tend.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The accuracy of integration is determined by the arguments tol, thres and thresp in a prior call to d02lxf. Note that only the local error at each step is controlled by these arguments. The error estimates obtained are not strict bounds but are usually reliable over one step. Over a number of steps the overall error may accumulate in various ways, depending on the system. The code is designed so that a reduction in tol should lead to an approximately proportional reduction in the error. You are strongly recommended to call d02laf with more than one value for tol and to compare the results obtained to estimate their accuracy.

The accuracy obtained depends on the type of error test used. If the solution oscillates around zero a relative error test should be avoided, whereas if the solution is exponentially increasing an absolute error test should not be used. For a description of the error test see the specifications of the arguments tol, thres and thresp in routine document d02lxf.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.

d02laf is not thread safe and should not be called from a multithreaded user program. Please see Section 1 in FL Interface Multithreading for more information on thread safety.

d02laf is not threaded in any implementation.

9 Further Comments

If d02laf fails with

ifail = 3

then the value of tol may be so small that a solution cannot be obtained, in which case the routine should be called again with a larger value for tol. If the accuracy requested is really needed then you should consider whether there is a more fundamental difficulty. For example:

(a)in the region of a singularity the solution components will usually be of a large magnitude. d02laf could be used in one-step mode to monitor the size of the solution with the aim of trapping the solution before the singularity. In any case numerical integration cannot be continued through a singularity, and analytical treatment may be necessary;
(b)if the solution contains fast oscillatory components, the routine will require a very small step size to preserve stability. This will usually be exhibited by excessive computing time and sometimes an error exit with $ifail = 3$ . The Runge–Kutta–Nystrom methods are not efficient in such cases and you should consider reposing your problem as a system of first-order ordinary differential equations and then using a routine from Sub-chapter D02M–N with the Blend formulae (see d02mvf).

d02laf can be used for producing results at short intervals (for example, for tabulation), in two ways. By far the less efficient is to call d02laf successively over short intervals,

t + (i - 1) \times h

t + i \times h

, although this is the only way if the higher order method has been selected and precisely not what it is intended for. A more efficient way, only for use when the lower order method has been selected, is to use d02laf in one-step mode. The output values of arguments y, yp, ydp, t and rwork are set correctly for a call to d02lzf to compute the solution and derivative at the required points.

10 Example

This example solves the following system (the two body problem)

\begin{matrix} y_{1}^{''} & = & - y_{1} / {(y_{1}^{2} + y_{2}^{2})}^{3 / 2} \\ y_{2}^{''} & = & - y_{2} / {(y_{1}^{2} + y_{2}^{2})}^{3 / 2} \end{matrix}

over the range

[0, 20]

with initial conditions

y_{1} = 1.0 - ε

y_{2} = 0.0

y_{1}^{'} = 0.0

and

y_{2}^{'} = \sqrt{(\frac{1 + ε}{1 - ε})}

where

ε

, the eccentricity, is

0.5

. The system is solved using the lower order method with relative local error tolerances

1.0E−4

and

1.0E−5

and default threshold tolerances. d02laf is used in one-step mode (

onestp = .TRUE.

) and d02lzf provides solution values at intervals of

2.0

d02la: FL CL CPP AD PY MB

NAG FL Interfaced02laf (ivp_​2nd_​rkn)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
d02laf (ivp_2nd_rkn)