NAG Library Routine Document

D02RAF

1  Purpose

2  Specification

3  Description

4  References

5  Parameters

1:     N – INTEGERInput

On entry: $\mathit{n}$, the number of differential equations.

Constraint: $\mathbf{N}>0$.

2:     MNP – INTEGERInput

On entry: MNP must be set to the maximum permitted number of points in the finite difference mesh. If LWORK or LIWORK are too small then internally MNP will be replaced by the maximum permitted by these values. (A warning message will be output if on entry IFAIL is set to obtain monitoring information.)

Constraint: $\mathbf{MNP}\ge 32$.

3:     NP – INTEGERInput/Output

On entry: must be set to the number of points to be used in the initial mesh.

Constraint: $4\le \mathbf{NP}\le \mathbf{MNP}$.

On exit: the number of points in the final mesh.

4:     NUMBEG – INTEGERInput

On entry: the number of left-hand boundary conditions (that is the number involving $y\left(a\right)$ only).

Constraint: $0\le \mathbf{NUMBEG}<\mathbf{N}$.

5:     NUMMIX – INTEGERInput

On entry: the number of coupled boundary conditions (that is the number involving both $y\left(a\right)$ and $y\left(b\right)$).

Constraint: $0\le \mathbf{NUMMIX}\le \mathbf{N}-\mathbf{NUMBEG}$.

6:     TOL – REAL (KIND=nag_wp)Input

On entry: a positive absolute error tolerance. If

$$a=x_1<x_2<\cdots <x_{\mathbf{NP}}=b$$

is the final mesh, $z_j\left(x_i\right)$ is the $j$th component of the approximate solution at $x_i$, and $y_j\left(x\right)$ is the $j$th component of the true solution of (1) and (2), then, except in extreme circumstances, it is expected that

$$\left|z_j\left(x_i\right)-y_j\left(x_i\right)\right|\le \mathbf{TOL}\text{, \hspace{1em}}i=1,2,\dots ,\mathbf{NP}\text{ and }j=1,2,\dots ,n\text{.}$$

(5)

Constraint: $\mathbf{TOL}>0.0$.

7:     INIT – INTEGERInput

On entry: indicates whether you wish to supply an initial mesh and approximate solution ($\mathbf{INIT}=1$) or whether default values are to be used, ($\mathbf{INIT}=0$).

Constraint: $\mathbf{INIT}=0$ or $1$.

8:     X(MNP) – REAL (KIND=nag_wp) arrayInput/Output

On entry: you must set $\mathbf{X}\left(1\right)=a$ and $\mathbf{X}\left(\mathbf{NP}\right)=b$. If $\mathbf{INIT}=0$ on entry a default equispaced mesh will be used, otherwise you must specify a mesh by setting $\mathbf{X}\left(\mathit{i}\right)=x_{\mathit{i}}$, for $\mathit{i}=2,3,\dots ,\mathbf{NP}-1$.

Constraints:

• if $\mathbf{INIT}=0$, $\mathbf{X}\left(1\right)<\mathbf{X}\left(\mathbf{NP}\right)$;
• if $\mathbf{INIT}=1$, $\mathbf{X}\left(1\right)<\mathbf{X}\left(2\right)<\cdots <\mathbf{X}\left(\mathbf{NP}\right)$.

On exit: $\mathbf{X}\left(1\right),\mathbf{X}\left(2\right),\dots ,\mathbf{X}\left(\mathbf{NP}\right)$ define the final mesh (with the returned value of NP) and $\mathbf{X}\left(1\right)=a$ and $\mathbf{X}\left(\mathbf{NP}\right)=b$.

9:     Y(LDY,MNP) – REAL (KIND=nag_wp) arrayInput/Output

On entry: if $\mathbf{INIT}=0$, then Y need not be set.

If $\mathbf{INIT}=1$, then the array Y must contain an initial approximation to the solution such that $\mathbf{Y}\left(j,i\right)$ contains an approximation to

$$y_j\left(x_i\right)\text{, \hspace{1em}}i=1,2,\dots ,\mathbf{NP}\text{ and }j=1,2,\dots ,n\text{.}$$

On exit: the approximate solution $z_j\left(x_i\right)$ satisfying (5) on the final mesh, that is

$$\mathbf{Y}\left(j,i\right)=z_j\left(x_i\right)\text{, \hspace{1em}}i=1,2,\dots ,\mathbf{NP}\text{ and }j=1,2,\dots ,n\text{,}$$

where NP is the number of points in the final mesh. If an error has occurred then Y contains the latest approximation to the solution. The remaining columns of Y are not used.

10:   LDY – INTEGERInput

On entry: the first dimension of the array Y as declared in the (sub)program from which D02RAF is called.

Constraint: $\mathbf{LDY}\ge \mathbf{N}$.

11:   ABT(N) – REAL (KIND=nag_wp) arrayOutput

On exit: $\mathbf{ABT}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,n$, holds the largest estimated error (in magnitude) of the $i$th component of the solution over all mesh points.

12:   FCN – SUBROUTINE, supplied by the user.External Procedure

FCN must evaluate the functions $f_i$ (i.e., the derivatives $y_i^{\prime }$) at a general point $x$ for a given value of $\epsilon$, the continuation parameter (see Section 3).

The specification of FCN is:

SUBROUTINE FCN ( X, EPS, Y, F, N) INTEGER  N REAL (KIND=nag_wp)  X, EPS, Y(N), F(N)

1:     X – REAL (KIND=nag_wp)Input

On entry: $x$, the value of the independent variable.

2:     EPS – REAL (KIND=nag_wp)Input

On entry: $\epsilon$, the value of the continuation parameter. This is $1$ if continuation is not being used.

3:     Y(N) – REAL (KIND=nag_wp) arrayInput

On entry: $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, the values of the dependent variables at $x$.

4:     F(N) – REAL (KIND=nag_wp) arrayOutput

On exit: the values of the derivatives $f_{\mathit{i}}$ evaluated at $x$ given $\epsilon$, for $\mathit{i}=1,2,\dots ,\mathit{n}$.

5:     N – INTEGERInput

On entry: $n$, the number of equations.

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

13:   G – SUBROUTINE, supplied by the user.External Procedure

G must evaluate the boundary conditions in equation (3) and place them in the array BC.

The specification of G is:

SUBROUTINE G ( EPS, YA, YB, BC, N) INTEGER  N REAL (KIND=nag_wp)  EPS, YA(N), YB(N), BC(N)

1:     EPS – REAL (KIND=nag_wp)Input

On entry: $\epsilon$, the value of the continuation parameter. This is $1$ if continuation is not being used.

2:     YA(N) – REAL (KIND=nag_wp) arrayInput

On entry: the value $y_{\mathit{i}}\left(a\right)$, for $\mathit{i}=1,2,\dots ,n$.

3:     YB(N) – REAL (KIND=nag_wp) arrayInput

On entry: the value $y_{\mathit{i}}\left(b\right)$, for $\mathit{i}=1,2,\dots ,n$.

4:     BC(N) – REAL (KIND=nag_wp) arrayOutput

On exit: the values $g_{\mathit{i}}\left(y\left(a\right),y\left(b\right),\epsilon \right)$, for $\mathit{i}=1,2,\dots ,n$. These must be ordered as follows:

(i)	first, the conditions involving only $y\left(a\right)$ (see NUMBEG);
(ii)	next, the NUMMIX coupled conditions involving both $y\left(a\right)$ and $y\left(b\right)$ (see NUMMIX); and,
(iii)	finally, the conditions involving only $y\left(b\right)$ ($\mathbf{N}-\mathbf{NUMBEG}-\mathbf{NUMMIX}$).

5:     N – INTEGERInput

On entry: $\mathit{n}$, the number of equations.

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

14:   IJAC – INTEGERInput

On entry: indicates whether or not you are supplying Jacobian evaluation routines.

$\mathbf{IJAC}\ne 0$

You must supply JACOBF and JACOBG and also, when continuation is used, JACEPS and JACGEP.

$\mathbf{IJAC}=0$

Numerical differentiation is used to calculate the Jacobian and the routines D02GAW, D02GAX, D02GAY and D02GAZ respectively may be used as the dummy parameters.

15:   JACOBF – SUBROUTINE, supplied by the NAG Library or the user.External Procedure

JACOBF evaluates the Jacobian $\left(\frac{\partial f_{\mathit{i}}}{\partial y_{\mathit{j}}}\right)$, for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,n$, given $x$ and $y_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,n$.

If all Jacobians are to be approximated internally by numerical differentiation then it must be replaced by the NAG defined null function pointer NULLFN.

The specification of JACOBF is:

SUBROUTINE JACOBF ( X, EPS, Y, F, N) INTEGER  N REAL (KIND=nag_wp)  X, EPS, Y(N), F(N,N)

1:     X – REAL (KIND=nag_wp)Input

On entry: $x$, the value of the independent variable.

2:     EPS – REAL (KIND=nag_wp)Input

On entry: $\epsilon$, the value of the continuation parameter. This is $1$ if continuation is not being used.

3:     Y(N) – REAL (KIND=nag_wp) arrayInput

On entry: $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, the values of the dependent variables at $x$.

4:     F(N,N) – REAL (KIND=nag_wp) arrayOutput

On exit: $\mathbf{F}\left(\mathit{j},\mathit{i}\right)$ must be set to the value of $\frac{\partial f_{\mathit{i}}}{\partial y_{\mathit{j}}}$, evaluated at the point $\left(x,y\right)$, for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,n$.

5:     N – INTEGERInput

On entry: $\mathit{n}$, the number of equations.

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

16:   JACOBG – SUBROUTINE, supplied by the NAG Library or the user.External Procedure

JACOBG evaluates the Jacobians $\left(\frac{\partial g_i}{\partial y_j\left(a\right)}\right)$ and $\left(\frac{\partial g_i}{\partial y_j\left(b\right)}\right)$. The ordering of the rows of AJ and BJ must correspond to the ordering of the boundary conditions described in the specification of G.

If all Jacobians are to be approximated internally by numerical differentiation then it must be replaced by the NAG defined null function pointer NULLFN.

The specification of JACOBG is:

SUBROUTINE JACOBG ( EPS, YA, YB, AJ, BJ, N) INTEGER  N REAL (KIND=nag_wp)  EPS, YA(N), YB(N), AJ(N,N), BJ(N,N)

1:     EPS – REAL (KIND=nag_wp)Input

On entry: $\epsilon$, the value of the continuation parameter. This is $1$ if continuation is not being used.

2:     YA(N) – REAL (KIND=nag_wp) arrayInput

On entry: the value $y_{\mathit{i}}\left(a\right)$, for $\mathit{i}=1,2,\dots ,n$.

3:     YB(N) – REAL (KIND=nag_wp) arrayInput

On entry: the value $y_{\mathit{i}}\left(b\right)$, for $\mathit{i}=1,2,\dots ,n$.

4:     AJ(N,N) – REAL (KIND=nag_wp) arrayOutput

On exit: $\mathbf{AJ}\left(\mathit{i},\mathit{j}\right)$ must be set to the value $\frac{\partial g_{\mathit{i}}}{\partial y_{\mathit{j}}\left(a\right)}$, for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,n$.

5:     BJ(N,N) – REAL (KIND=nag_wp) arrayOutput

On exit: $\mathbf{BJ}\left(\mathit{i},\mathit{j}\right)$ must be set to the value $\frac{\partial g_{\mathit{i}}}{\partial y_{\mathit{j}}\left(b\right)}$, for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,n$.

6:     N – INTEGERInput

On entry: $\mathit{n}$, the number of equations.

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

17:   DELEPS – REAL (KIND=nag_wp)Input/Output

On entry: must be given a value which specifies whether continuation is required. If $\mathbf{DELEPS}\le 0.0$ or $\mathbf{DELEPS}\ge 1.0$ then it is assumed that continuation is not required. If $0.0<\mathbf{DELEPS}<1.0$ then it is assumed that continuation is required unless $\mathbf{DELEPS}<\sqrt{\mathit{machine precision}}$ when an error exit is taken. DELEPS is used as the increment ${\epsilon }_2-{\epsilon }_1$ (see (4)) and the choice $\mathbf{DELEPS}=0.1$ is recommended.

On exit: an overestimate of the increment ${\epsilon }_p-{\epsilon }_{p-1}$ (in fact the value of the increment which would have been tried if the restriction ${\epsilon }_p=1$ had not been imposed). If continuation was not requested then $\mathbf{DELEPS}=0.0$.

If continuation is not requested then JACEPS and JACGEP may each be replaced by dummy actual parameters in the call to D02RAF. (D02GAW and D02GAX respectively may be used as the dummy parameters.)

18:   JACEPS – SUBROUTINE, supplied by the NAG Library or the user.External Procedure

JACEPS evaluates the derivative $\frac{\partial f_i}{\partial \epsilon }$ given $x$ and $y$ if continuation is being used.

If all Jacobians (derivatives) are to be approximated internally by numerical differentiation, or continuation is not being used, then it must be replaced by the NAG defined null function pointer NULLFN.

The specification of JACEPS is:

SUBROUTINE JACEPS ( X, EPS, Y, F, N) INTEGER  N REAL (KIND=nag_wp)  X, EPS, Y(N), F(N)

1:     X – REAL (KIND=nag_wp)Input

On entry: $x$, the value of the independent variable.

2:     EPS – REAL (KIND=nag_wp)Input

On entry: $\epsilon$, the value of the continuation parameter.

3:     Y(N) – REAL (KIND=nag_wp) arrayInput

On entry: the solution values $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$, at the point $x$.

4:     F(N) – REAL (KIND=nag_wp) arrayOutput

On exit: $\mathbf{F}\left(\mathit{i}\right)$ must contain the value $\frac{\partial f_{\mathit{i}}}{\partial \epsilon }$ at the point $\left(x,y\right)$, for $\mathit{i}=1,2,\dots ,n$.

5:     N – INTEGERInput

On entry: $\mathit{n}$, the number of equations.

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

19:   JACGEP – SUBROUTINE, supplied by the NAG Library or the user.External Procedure

JACGEP evaluates the derivatives $\frac{\partial g_i}{\partial \epsilon }$ if continuation is being used.

If all Jacobians (derivatives) are to be approximated internally by numerical differentiation, or continuation is not being used, then it must be replaced by the NAG defined null function pointer NULLFN.

The specification of JACGEP is:

SUBROUTINE JACGEP ( EPS, YA, YB, BCEP, N) INTEGER  N REAL (KIND=nag_wp)  EPS, YA(N), YB(N), BCEP(N)

1:     EPS – REAL (KIND=nag_wp)Input

On entry: $\epsilon$, the value of the continuation parameter.

2:     YA(N) – REAL (KIND=nag_wp) arrayInput

On entry: the value of $y_{\mathit{i}}\left(a\right)$, for $\mathit{i}=1,2,\dots ,n$.

3:     YB(N) – REAL (KIND=nag_wp) arrayInput

On entry: the value of $y_{\mathit{i}}\left(b\right)$, for $\mathit{i}=1,2,\dots ,n$.

4:     BCEP(N) – REAL (KIND=nag_wp) arrayOutput

On exit: $\mathbf{BCEP}\left(\mathit{i}\right)$ must contain the value of $\frac{\partial g_{\mathit{i}}}{\partial \epsilon }$, for $\mathit{i}=1,2,\dots ,n$.

5:     N – INTEGERInput

On entry: $\mathit{n}$, the number of equations.

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

20:   WORK(LWORK) – REAL (KIND=nag_wp) arrayWorkspace

21:   LWORK – INTEGERInput

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

Constraint: $\mathbf{LWORK}\ge \mathbf{MNP}\times \left(3{\mathbf{N}}^2+6\mathbf{N}+2\right)+4{\mathbf{N}}^2+3\mathbf{N}$.

22:   IWORK(LIWORK) – INTEGER arrayWorkspace

23:   LIWORK – INTEGERInput

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

Constraints:

• if $\mathbf{IJAC}\ne 0$, $\mathbf{LIWORK}\ge \mathbf{MNP}\times \left(2\times \mathbf{N}+1\right)+\mathbf{N}$;
• if $\mathbf{IJAC}=0$, $\mathbf{LIWORK}\ge \mathbf{MNP}\times \left(2\times \mathbf{N}+1\right)+{\mathbf{N}}^2+4\times \mathbf{N}+2$.

24:   IFAIL – INTEGERInput/Output

For this routine, the normal use of IFAIL is extended to control the printing of error and warning messages as well as specifying hard or soft failure (see Section 3.3 in the Essential Introduction).

On entry: IFAIL must be set to a value with the decimal expansion $\mathit{cba}$, where each of the decimal digits $c$, $b$ and $a$ must have a value of $0$ or $1$.

$a=0$	specifies hard failure, otherwise soft failure;
$b=0$	suppresses error messages, otherwise error messages will be printed (see Section 6);
$c=0$	suppresses warning messages, otherwise warning messages will be printed (see Section 6).

The recommended value for inexperienced users is $110$ (i.e., hard failure with all messages printed).

On exit: $\mathbf{IFAIL}=\mathbf{0}$ unless the routine detects an error or a warning has been flagged (see Section 6).

6  Error Indicators and Warnings

$\mathbf{IFAIL}=1$

One or more of the parameters N, MNP, NP, NUMBEG, NUMMIX, TOL, DELEPS, LWORK or LIWORK is incorrectly set, or $\mathbf{X}\left(1\right)\ge \mathbf{X}\left(\mathbf{NP}\right)$ or the mesh points $\mathbf{X}\left(i\right)$ are not in strictly ascending order.

$\mathbf{IFAIL}=2$

A finer mesh is required for the accuracy requested; that is MNP is not large enough. This error exit normally occurs when the problem being solved is difficult (for example, there is a boundary layer) and high accuracy is requested. A poor initial choice of mesh points will make this error exit more likely.

$\mathbf{IFAIL}=3$

The Newton iteration has failed to converge. There are several possible causes for this error:

(i)	faulty coding in one of the Jacobian calculation routines;
(ii)	if $\mathbf{IJAC}=0$ then inaccurate Jacobians may have been calculated numerically (this is a very unlikely cause); or,
(iii)	a poor initial mesh or initial approximate solution has been selected either by you or by default or there are not enough points in the initial mesh. Possibly, you should try the continuation facility.

$\mathbf{IFAIL}=4$

The Newton iteration has reached round-off error level. It could be however that the answer returned is satisfactory. The error is likely to occur if too high an accuracy is requested.

$\mathbf{IFAIL}=5$

The Jacobian calculated by JACOBG (or the equivalent matrix calculated by numerical differentiation) is singular. This may occur due to faulty coding of JACOBG or, in some circumstances, to a zero initial choice of approximate solution (such as is chosen when $\mathbf{INIT}=0$).

$\mathbf{IFAIL}=6$

There is no dependence on $\epsilon$ when continuation is being used. This can be due to faulty coding of JACEPS or JACGEP or, in some circumstances, to a zero initial choice of approximate solution (such as is chosen when $\mathbf{INIT}=0$).

$\mathbf{IFAIL}=7$

DELEPS is required to be less than machine precision for continuation to proceed. It is likely that either the problem (3) has no solution for some value near the current value of $\epsilon$ (see the advisory print out from D02RAF) or that the problem is so difficult that even with continuation it is unlikely to be solved using this routine. If the latter cause is suspected then using more mesh points initially may help.

$\mathbf{IFAIL}=8$

$\mathbf{IFAIL}=9$

A serious error has occurred in an internal call. Check all array subscripts and subroutine parameter lists in calls to D02RAF. Seek expert help.

7  Accuracy

8  Further Comments

9  Example

9.1  Program Text

Program Text (d02rafe.f90)

9.2  Program Data

Program Data (d02rafe.d)

9.3  Program Results

Program Results (d02rafe.r)