NAG Toolbox: nag_ode_ivp_2nd_rkn_setup (d02lx)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{h}$ – double scalar

If $\mathbf{start}=\mathit{true}$, h may specify an initial step size to be attempted in nag_ode_ivp_2nd_rkn (d02la).

If $\mathbf{start}=\mathit{false}$, h may specify a step size to override the choice of next step attempted made internally to nag_ode_ivp_2nd_rkn (d02la).

The sign of h is not important, as the absolute value of h is chosen and the appropriate sign is selected by nag_ode_ivp_2nd_rkn (d02la).

If this option is not required then you must set $\mathbf{h}=0.0$.

2:     $\mathrm{tol}$ – double scalar

Must be set to a relative tolerance for controlling the error in the integration by nag_ode_ivp_2nd_rkn (d02la). nag_ode_ivp_2nd_rkn (d02la) has been designed so that, for most problems, a reduction in tol leads to an approximately proportional reduction in the error in the solution. However the actual relation between tol and the accuracy of the solution cannot be guaranteed. You are strongly recommended to repeat the integration with a smaller value of tol and compare the results. See the description of thres and thresp for further details of how tol is used.

Constraint: $10\times \epsilon \le \mathbf{tol}\le 1.0$ ($\epsilon$ is the machine precision, see nag_machine_precision (x02aj)).

3:     $\mathrm{thres}\left(\mathbf{neq}\right)$ – double array

4:     $\mathrm{thresp}\left(\mathbf{neq}\right)$ – double array

thres and thresp may be set to thresholds for use in the error control of nag_ode_ivp_2nd_rkn (d02la). At each step in the numerical integration estimates of the local errors $\mathrm{E1}\left(\mathit{i}\right)$ and $\mathrm{E2}\left(\mathit{i}\right)$ in the solution, $y_{\mathit{i}}$, and its derivative, $y_{\mathit{i}}^{\prime }$, respectively are computed, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$. For the step to be accepted conditions of the following type must be satisfied:

$$\begin{array}{l}\underset{1\le i\le \mathbf{neq}}{\max }\left(\frac{\mathrm{E1}\left(i\right)}{\max \left(\mathbf{thres}\left(i\right),\left|y_i\right|\right)}\right)\le \mathbf{tol}\text{,}\\ \\ \\ \underset{1\le i\le \mathbf{neq}}{\max }\left(\frac{\mathrm{E2}\left(i\right)}{\max \left(\mathbf{thresp}\left(i\right),\left|y_i^{\prime }\right|\right)}\right)\le \mathbf{tol}\text{.}\end{array}$$

If one or both of these is not satisfied then the step size is reduced and the solution is recomputed.

If $\mathbf{thres}\left(1\right)\le 0.0$ on entry, then a value of $50.0\times \epsilon$ is used for $\mathbf{thres}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$, where $\epsilon$ is machine precision. Similarly for thresp.

Constraints:

• $\mathbf{thres}\left(1\right)\le 0.0$ or $\mathbf{thres}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$;
• $\mathbf{thresp}\left(1\right)\le 0.0$ or $\mathbf{thresp}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,\mathbf{neq}$.

5:     $\mathrm{maxstp}$ – int64int32nag_int scalar

A bound on the number of steps attempted in any one call of nag_ode_ivp_2nd_rkn (d02la).

If $\mathbf{maxstp}\le 0$ on entry, a value of $1000$ is used.

6:     $\mathrm{start}$ – logical scalar

Specifies whether or not the call of nag_ode_ivp_2nd_rkn (d02la) is for a new problem. $\mathbf{start}=\mathit{true}$ indicates that a new problem is to be solved. $\mathbf{start}=\mathit{false}$ indicates the call of nag_ode_ivp_2nd_rkn_setup (d02lx) is prior to a continuation call of nag_ode_ivp_2nd_rkn (d02la).

7:     $\mathrm{onestp}$ – logical scalar

The mode of operation for nag_ode_ivp_2nd_rkn (d02la).

$\mathbf{onestp}=\mathit{true}$

nag_ode_ivp_2nd_rkn (d02la) will operate in one-step mode, that is it will return after each successful step.

$\mathbf{onestp}=\mathit{false}$

nag_ode_ivp_2nd_rkn (d02la) will operate in interval mode, that is it will return at the end of the integration interval.

8:     $\mathrm{high}$ – logical scalar

If $\mathbf{high}=\mathit{true}$, a high-order method will be used, whereas if $\mathbf{high}=\mathit{false}$, a low-order method will be used. (See the specification of nag_ode_ivp_2nd_rkn (d02la) for further details.)

9:     $\mathrm{rwork}\left(\mathit{lrwork}\right)$ – double array

lrwork, the dimension of the array, must satisfy the constraint

• if $\mathbf{high}=\mathit{true}$, $\mathit{lrwork}\ge 16+20\times \mathbf{neq}$;
• if $\mathbf{high}=\mathit{false}$, $\mathit{lrwork}\ge 16+11\times \mathbf{neq}$.

This must be the same argument rwork supplied to nag_ode_ivp_2nd_rkn (d02la). It is used to pass information to nag_ode_ivp_2nd_rkn (d02la) and therefore the contents of this array must not be changed before calling nag_ode_ivp_2nd_rkn (d02la).

Optional Input Parameters

1:     $\mathrm{neq}$ – int64int32nag_int scalar

Default: the dimension of the arrays thres, thresp. (An error is raised if these dimensions are not equal.)

The number of second-order ordinary differential equations to be solved by nag_ode_ivp_2nd_rkn (d02la).

Constraint: $\mathbf{neq}\ge 1$.

Output Parameters

1:     $\mathrm{start}$ – logical scalar

$\mathbf{start}=\mathit{false}$.

2:     $\mathrm{rwork}\left(\mathit{lrwork}\right)$ – double array

3:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

   $\mathbf{ifail}=1$

$\mathbf{thres}\left(1\right)>0.0$ and for some $i$ $\mathbf{thres}\left(i\right)\le 0.0$, $1\le i\le \mathbf{neq}$, and/or, $\mathbf{thresp}\left(1\right)>0.0$ and for some $i$ $\mathbf{thresp}\left(i\right)\le 0.0$, $1\le i\le \mathbf{neq}$.

   $\mathbf{ifail}=2$

lrwork is too small.

   $\mathbf{ifail}=3$

tol does not satisfy $10\times \epsilon \le \mathbf{tol}\le 1.0$ ($\epsilon$ is the machine precision, see nag_machine_precision (x02aj))

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: d02lx_example

function d02lx_example


fprintf('d02lx example results\n\n');

neq     = 2;

% Parameters values for integrator setup
h       = 0;
reltol  = 1e-4;
thresh  = zeros(neq,1);
threshp = zeros(neq,1);
maxstp  = int64(0); 
start   = true;
onestp  = true;
high    = false;
rwork   = zeros(16+20*neq,1);

% Call setup
[start, rwork, ifail] = ...
    d02lx(h, reltol, thresh, threshp, maxstp, start, onestp, high, rwork);

% Initializations for integration
t     = 0;
tend  = 20;
y     = [0.5; 0];
yp    = [0; sqrt(3)];
ydp   = zeros(neq, 1);

% take one time-step towards tend
[t, y, yp, ydp, rwork, ifail] = d02la(@fcn, t, tend, y, yp, ydp, rwork);

fprintf('Solution at first time step, t = %6.3f:\n', t);
fprintf('    y   = %6.3f, %6.3f\n',y);
fprintf('    y''  = %6.3f, %6.3f\n',yp);
fprintf('    y'''' = %6.3f, %6.3f\n',ydp);



function ydp = fcn(neq, t, y)
% Evaluate second derivatives.
r = sqrt(y(1)^2+y(2)^2)^3;
f = zeros(2,1);
ydp(1) = -y(1)/r;
ydp(2) = -y(2)/r;

d02lx example results

Solution at first time step, t =  0.062:
    y   =  0.492,  0.108
    y'  = -0.246,  1.705
    y'' = -3.848, -0.841