# NAG Library Routine Document

## 1Purpose

d02pqf is a setup routine which must be called prior to the first call of either of the integration routines d02pef, d02pff and d02pgf.

## 2Specification

Fortran Interface
 Subroutine d02pqf ( n, tend, tol,
 Integer, Intent (In) :: n, method Integer, Intent (Inout) :: ifail Integer, Intent (Out) :: iwsav(130) Real (Kind=nag_wp), Intent (In) :: tstart, tend, yinit(n), tol, thresh(n), hstart Real (Kind=nag_wp), Intent (Out) :: rwsav(32*n+350)
#include nagmk26.h
 void d02pqf_ (const Integer *n, const double *tstart, const double *tend, const double yinit[], const double *tol, const double thresh[], const Integer *method, const double *hstart, Integer iwsav[], double rwsav[], Integer *ifail)

## 3Description

d02pqf and its associated routines (d02pef, d02pff, d02pgf, d02phf, d02pjf, d02prf, d02psf, d02ptf and d02puf) solve the initial value problem for a first-order system of ordinary differential equations. The routines, based on Runge–Kutta methods and derived from RKSUITE (see Brankin et al. (1991)), integrate
 $y′=ft,y given yt0=y0$
where $y$ is the vector of $n$ solution components and $t$ is the independent variable.
The integration proceeds by steps from the initial point ${t}_{0}$ towards the final point ${t}_{f}$. An approximate solution $y$ is computed at each step. For each component ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$, the error made in the step, i.e., the local error, is estimated. The step size is chosen automatically so that the integration will proceed efficiently while keeping this local error estimate smaller than a tolerance that you specify by means of arguments tol and thresh.
d02pef can be used to solve the ‘usual task’, namely integrating the system of differential equations to obtain answers at points you specify. d02pff is used for more ‘complicated’ tasks where $f\left(t,y\right)$ can readily be coded within a routine argument and high-order interpolation is not required. d02pgf is used for the most ‘complicated’ tasks where $f\left(t,y\right)$ is best evaluated outside the integrator or where high-order interpolation is required.
You should consider carefully how you want the local error to be controlled. Essentially the code uses relative local error control, with tol being the desired relative accuracy. For reliable computation, the code must work with approximate solutions that have some correct digits, so there is an upper bound on the value used for tol. It is impossible to compute a numerical solution that is more accurate than the correctly rounded value of the true solution, so you are not allowed to specify tol too small for the precision you are using. The magnitude of the local error in ${y}_{i}$ on any step will not be greater than ${\mathbf{tol}}×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mu }_{i},{\mathbf{thresh}}\left(i\right)\right)$ where ${\mu }_{i}$ is an average magnitude of ${y}_{i}$ over the step. If ${\mathbf{thresh}}\left(i\right)$ is smaller than the current value of ${\mu }_{i}$, this is a relative error test and tol indicates how many significant digits you want in ${y}_{i}$. If ${\mathbf{thresh}}\left(i\right)$ is larger than the current value of ${\mu }_{i}$, this is an absolute error test with tolerance ${\mathbf{tol}}×{\mathbf{thresh}}\left(i\right)$. Relative error control is the recommended mode of operation, but pure relative error control, ${\mathbf{thresh}}\left(i\right)=0.0$, is not permitted. See Section 9 for further information about error control.
d02pef, d02pff and d02pgf control local error rather than the true (global) error, the difference between the numerical and true solution. Control of the local error controls the true error indirectly. Roughly speaking, the code produces a solution that satisfies the differential equation with a discrepancy bounded in magnitude by the error tolerance. What this implies about how close the numerical solution is to the true solution depends on the stability of the problem. Most practical problems are at least moderately stable, and the true error is then comparable to the error tolerance. To judge the accuracy of the numerical solution, you could reduce tol substantially, e.g., use $0.1×{\mathbf{tol}}$, and solve the problem again. This will usually result in a rather more accurate solution, and the true error of the first integration can be estimated by comparison. Alternatively, a global error assessment can be computed automatically using the argument method. Because indirect control of the true error by controlling the local error is generally satisfactory and because both ways of assessing true errors cost twice, or more, the cost of the integration itself, such assessments are used mostly for spot checks, selecting appropriate tolerances for local error control, and exploratory computations.
d02pef, d02pff and d02pgf each implement three Runge–Kutta formula pairs, and you must select one for the integration. The best choice for method depends on the problem. The order of accuracy is $3$, $5$ and $8$ respectively. As a rule, the smaller tol is, the larger you should take the order of the method. If the components thresh are small enough that you are effectively specifying relative error control, experience suggests
 tol efficient method ${10}^{-2}-{10}^{-4}$ order $2$ and $3$ pair ${10}^{-3}-{10}^{-6}$ order $4$ and $5$ pair ${10}^{-5}-$ order $7$ and $8$ pair
The overlap in the ranges of tolerances appropriate for a given method merely reflects the dependence of efficiency on the problem being solved. Making tol smaller will normally make the integration more expensive. However, in the range of tolerances appropriate to a method, the increase in cost is modest. There are situations for which one method, or even this kind of code, is a poor choice. You should not specify a very small value for ${\mathbf{thresh}}\left(i\right)$, when the $i$th solution component might vanish. In particular, you should not do this when ${y}_{i}=0.0$. If you do, the code will have to work hard with any value for method to compute significant digits, but the lowest order method is a particularly poor choice in this situation. All three methods are inefficient when the problem is ‘stiff’. If it is only mildly stiff, you can solve it with acceptable efficiency with the order $2$ and $3$ pair, but if it is moderately or very stiff, a code designed specifically for such problems will be much more efficient. The higher the order the more smoothness is required of the solution in order for the method to be efficient.
When assessment of the true (global) error is requested, this error assessment is updated at each step. Its value can be obtained at any time by a call to d02puf. The code monitors the computation of the global error assessment and reports any doubts it has about the reliability of the results. The assessment scheme requires some smoothness of $f\left(t,y\right)$, and it can be deceived if $f$ is insufficiently smooth. At very crude tolerances the numerical solution can become so inaccurate that it is impossible to continue assessing the accuracy reliably. At very stringent tolerances the effects of finite precision arithmetic can make it impossible to assess the accuracy reliably. The cost of this is roughly twice the cost of the integration itself with the 5th and 8th order methods, and three times with the 3rd order method.
The first step of the integration is critical because it sets the scale of the problem. The integrator will find a starting step size automatically if you set the argument hstart to $0.0$. Automatic selection of the first step is so effective that you should normally use it. Nevertheless, you might want to specify a trial value for the first step to be certain that the code recognizes the scale on which phenomena occur near the initial point. Also, automatic computation of the first step size involves some cost, so supplying a good value for this step size will result in a less expensive start. If you are confident that you have a good value, provide it via the argument hstart.

## 4References

Brankin R W, Gladwell I and Shampine L F (1991) RKSUITE: A suite of Runge–Kutta codes for the initial value problems for ODEs SoftReport 91-S1 Southern Methodist University

## 5Arguments

1:     $\mathbf{n}$ – IntegerInput
On entry: $n$, the number of ordinary differential equations in the system to be solved by the integration routine.
Constraint: ${\mathbf{n}}\ge 1$.
2:     $\mathbf{tstart}$ – Real (Kind=nag_wp)Input
On entry: the initial value of the independent variable, ${t}_{0}$.
3:     $\mathbf{tend}$ – Real (Kind=nag_wp)Input
On entry: the final value of the independent variable, ${t}_{f}$, at which the solution is required. tstart and tend together determine the direction of integration.
Constraint: ${\mathbf{tend}}$ must be distinguishable from tstart for the method and the precision of the machine being used.
4:     $\mathbf{yinit}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: ${y}_{0}$, the initial values of the solution, ${y}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
5:     $\mathbf{tol}$ – Real (Kind=nag_wp)Input
On entry: a relative error tolerance. The actual tolerance used is ; that is, the minimum tolerance is set at $10$ times machine precision and the maximum tolerance is set at $0.01$.
6:     $\mathbf{thresh}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput
On entry: a vector of thresholds. For the $i$th component, the actual threshold used is $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(\sqrt{\mathbf{safe range}},{\mathbf{thresh}}\left(i\right)\right)$, where $\mathbf{safe range}$ is the value returned by x02amf.
7:     $\mathbf{method}$ – IntegerInput
On entry: the Runge–Kutta method to be used.
${\mathbf{method}}=1$ or $-1$
A $2\left(3\right)$ pair is used.
${\mathbf{method}}=2$ or $-2$
A $4\left(5\right)$ pair is used.
${\mathbf{method}}=3$ or $-3$
A $7\left(8\right)$ pair is used.
Constraint: ${\mathbf{method}}=1$, $-1$, $2$, $-2$, $3$ or $-3$. Note: if ${\mathbf{method}}>0$ then global error assessment is to be computed with the main integration; if ${\mathbf{method}}<0$ then global error assessment will not be computed.
8:     $\mathbf{hstart}$ – Real (Kind=nag_wp)Input
On entry: a value for the size of the first step in the integration to be attempted. The absolute value of hstart is used with the direction being determined by tstart and tend. The actual first step taken by the integrator may be different to hstart if the underlying algorithm determines that hstart is unsuitable. If ${\mathbf{hstart}}=0.0$ then the size of the first step is computed automatically.
Suggested value: ${\mathbf{hstart}}=0.0$.
9:     $\mathbf{iwsav}\left(130\right)$ – Integer arrayCommunication Array
10:   $\mathbf{rwsav}\left(32×{\mathbf{n}}+350\right)$ – Real (Kind=nag_wp) arrayCommunication Array
On exit: the contents of the communication arrays must not be changed prior to calling one of the integration routines.
11:   $\mathbf{ifail}$ – IntegerInput/Output
On entry: ifail must be set to $0$, $-1\text{​ or ​}1$. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{​ or ​}1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is $0$. When the value $-\mathbf{1}\text{​ or ​}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

## 6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
On entry, ${\mathbf{method}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{method}}=-3$, $-2$, $-1$, $1$, $2$ or $3$.
On entry, ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{n}}\ge 1$.
On entry, too much workspace required.
Workspace provided was $〈\mathit{\text{value}}〉$, workspace required is $〈\mathit{\text{value}}〉$.
On entry, ${\mathbf{tstart}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{tstart}}\ne {\mathbf{tend}}$.
On entry, tstart is too close to tend.
$\left|{\mathbf{tstart}}-{\mathbf{tend}}\right|=〈\mathit{\text{value}}〉$, but this quantity should be at least $〈\mathit{\text{value}}〉$.
${\mathbf{ifail}}=-99$
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

Not applicable.

## 8Parallelism and Performance

d02pqf is not threaded in any implementation.

It is often the case that a solution component (the $i$th, say) is of no interest when it is smaller in magnitude than a certain threshold. You can inform the code of this by setting the $i$th component of thresh to this threshold. In this way you avoid the cost of computing significant digits in the $i$th component of $y$ when it is smaller than the threshold of interest. This matter is important when a component of $y$ vanishes, and in particular, when the initial value is zero. An appropriate threshold depends on the general size of $y$ in the course of the integration. Physical reasoning may help you select suitable threshold values. If you do not know what to expect of $y$, you can find out by a preliminary integration using d02pef with nominal values of thresh. As d02pef integrates by steps in time, it stores, for each component, the largest magnitude of solution computed so far; these values are output in the array ymax. This can help determine more appropriate values for thresh for an accurate integration. For example, the values in thresh could be set to  times the final value of ymax.