d02rac solves a two-point boundary value problem with general boundary conditions for a system of ordinary differential equations, using a deferred correction technique and Newton iteration.
The function may be called by the names: d02rac or nag_ode_bvp_fd_nonlin_gen.
3Description
d02rac solves a two-point boundary value problem for a system of $n$ ordinary differential equations in the interval $[a,b]$ with $b>a$. The system is written in the form
The functions ${g}_{i}$ are evaluated by g. The solution is computed using a finite difference technique with deferred correction allied to a Newton iteration to solve the finite difference equations. The technique used is described fully in Pereyra (1979).
You must supply an absolute error tolerance and may also supply an initial mesh for the finite difference equations and an initial approximate solution (alternatively a default mesh and approximation are used). The approximate solution is corrected using Newton iteration and deferred correction. Then, additional points are added to the mesh and the solution is recomputed with the aim of making the error everywhere less than your tolerance and of approximately equidistributing the error on the final mesh. The solution is returned on this final mesh.
If the solution is required at a few specific points then these should be included in the initial mesh. If, on the other hand, the solution is required at several specific points then you should use the interpolation functions provided in Chapter E01 if these points do not themselves form a convenient mesh.
These may be supplied through jacobf for $\left(\frac{\partial {f}_{i}}{\partial {y}_{j}}\right)$ and jacobg for the others. Alternatively the Jacobians may be calculated by numerical differentiation using the algorithm described in Curtis et al. (1974).
For problems of the type (1) and (2) for which it is difficult to determine an initial approximation from which the Newton iteration will converge, a continuation facility is provided. You must set up a family of problems
where $f={[{f}_{1},{f}_{2},\dots ,{f}_{n}]}^{\mathrm{T}}$ etc., and where $\epsilon $ is a continuation parameter. The choice $\epsilon =0$ must give a problem (3) which is easy to solve and $\epsilon =1$ must define the problem whose solution is actually required. The function solves a sequence of problems with $\epsilon $ values
The number $p$ and the values ${\epsilon}_{i}$ are chosen by the function so that each problem can be solved using the solution of its predecessor as a starting approximation. Jacobians $\frac{\partial f}{\partial \epsilon}$ and $\frac{\partial g}{\partial \epsilon}$ are required and they may be supplied by you via jaceps and jacgep respectively or may be computed by numerical differentiation.
4References
Curtis A R, Powell M J D and Reid J K (1974) On the estimation of sparse Jacobian matrices J. Inst. Maths. Applics.13 117–119
Pereyra V (1979) PASVA3: An adaptive finite-difference Fortran program for first order nonlinear, ordinary boundary problems Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science (eds B Childs, M Scott, J W Daniel, E Denman and P Nelson) 76 Springer–Verlag
5Arguments
1: $\mathbf{neq}$ – IntegerInput
On entry:
$\mathit{n}$, the number of differential equations.
Constraint:
${\mathbf{neq}}>0$.
2: $\mathbf{deleps}$ – double *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{machineprecision}}$ 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
the NAG defined null function pointer NULLFN.
3: $\mathbf{fcn}$ – function, supplied by the userExternal Function
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).
On exit: the values of the derivatives
${f}_{\mathit{i}}$ evaluated at $x$ given $\epsilon $, for $\mathit{i}=1,2,\dots ,\mathit{n}$.
6: $\mathbf{comm}$ – Nag_User *
Pointer to structure of type Nag_User.
p – Pointer
The type Pointer will be void *. Before calling d02rac these pointers may be allocated memory and initialized with various quantities for use by comm when called from d02rac.
Note:fcn should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02rac. If your code inadvertently does return any NaNs or infinities, d02rac is likely to produce unexpected results.
4: $\mathbf{numbeg}$ – IntegerInput
On entry: the number of left-hand boundary conditions (that is the number involving $y\left(a\right)$ only).
On exit: the values
${g}_{\mathit{i}}(y\left(a\right),y\left(b\right),\epsilon )$, 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{neq}}-{\mathbf{numbeg}}-{\mathbf{nummix}}$).
6: $\mathbf{comm}$ – Nag_User *
Pointer to structure of type Nag_User.
p – Pointer
The type Pointer will be void *. Before calling d02rac these pointers may be allocated memory and initialized with various quantities for use by comm when called from d02rac.
Note:g should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02rac. If your code inadvertently does return any NaNs or infinities, d02rac is likely to produce unexpected results.
7: $\mathbf{init}$ – Nag_MeshSetInput
On entry: indicates whether you wish to supply an initial mesh and approximate solution (${\mathbf{init}}=\mathrm{Nag\_UserInitMesh}$) or whether default values are to be used, (${\mathbf{init}}=\mathrm{Nag\_DefInitMesh}$).
Constraint:
${\mathbf{init}}=\mathrm{Nag\_DefInitMesh}$ or $\mathrm{Nag\_UserInitMesh}$.
8: $\mathbf{mnp}$ – IntegerInput
On entry:
mnp must be set to the maximum permitted number of points in the finite difference mesh.
Constraint:
${\mathbf{mnp}}\ge 32$.
9: $\mathbf{np}$ – Integer *Input/Output
On entry: must be set to the number of points to be used in the initial mesh.
On entry: you must set ${\mathbf{x}}\left[0\right]=a$ and ${\mathbf{x}}\left[{\mathbf{np}}-1\right]=b$. If ${\mathbf{init}}=\mathrm{Nag\_DefInitMesh}$ on entry a default equispaced mesh will be used, otherwise you must specify a mesh by setting
${\mathbf{x}}\left[\mathit{i}-1\right]={x}_{\mathit{i}}$, for $\mathit{i}=2,3,\dots ,{\mathbf{np}}-1$.
Constraints:
if ${\mathbf{init}}=\mathrm{Nag\_DefInitMesh}$, ${\mathbf{x}}\left[0\right]<{\mathbf{x}}\left[{\mathbf{np}}-1\right]$;
if ${\mathbf{init}}=\mathrm{Nag\_UserInitMesh}$, ${\mathbf{x}}\left[0\right]<{\mathbf{x}}\left[1\right]<\cdots <{\mathbf{x}}\left[{\mathbf{np}}-1\right]$.
On exit: ${\mathbf{x}}\left[0\right],{\mathbf{x}}\left[1\right],\dots ,{\mathbf{x}}\left[{\mathbf{np}}-1\right]$ define the final mesh (with the returned value of np) and ${\mathbf{x}}\left[0\right]=a$ and ${\mathbf{x}}\left[{\mathbf{np}}-1\right]=b$.
On entry: if ${\mathbf{init}}=\mathrm{Nag\_DefInitMesh}$, then y need not be set.
If ${\mathbf{init}}=\mathrm{Nag\_UserInitMesh}$, then the array y must contain an initial approximation to the solution such that ${\mathbf{y}}\left[\left(j-1\right)\times {\mathbf{mnp}}+i-1\right]$ contains an approximation to
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.
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
On exit:
${\mathbf{abt}}\left[\mathit{i}-1\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.
14: $\mathbf{jacobf}$ – function, supplied by the userExternal Function
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.
On exit: ${\mathbf{f}}\left[\left(\mathit{j}-1\right)\times {\mathbf{neq}}+\mathit{i}-1\right]$ must be set to the value of $\frac{\partial {f}_{\mathit{i}}}{\partial {y}_{\mathit{j}}}$, evaluated at the point $(x,y)$, for $\mathit{i}=1,2,\dots ,n$ and $\mathit{j}=1,2,\dots ,n$.
6: $\mathbf{comm}$ – Nag_User *
Pointer to structure of type Nag_User.
p – Pointer
The type Pointer will be void *. Before calling d02rac these pointers may be allocated memory and initialized with various quantities for use by comm when called from d02rac.
Note:jacobf should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02rac. If your code inadvertently does return any NaNs or infinities, d02rac is likely to produce unexpected results.
Note that if jacobf is supplied then jacobg must also be supplied. Note that if jacobf is supplied and continuation is requested then jaceps and jacgep must also be supplied.
15: $\mathbf{jacobg}$ – function, supplied by the userExternal Function
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.
On exit: ${\mathbf{aj}}\left[\left(\mathit{i}-1\right)\times {\mathbf{neq}}+\mathit{j}-1\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$.
On exit: ${\mathbf{bj}}\left[\left(\mathit{i}-1\right)\times {\mathbf{neq}}+\mathit{j}-1\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$.
7: $\mathbf{comm}$ – Nag_User *
Pointer to structure of type Nag_User.
p – Pointer
The type Pointer will be void *. Before calling d02rac these pointers may be allocated memory and initialized with various quantities for use by comm when called from d02rac.
Note:jacobg should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02rac. If your code inadvertently does return any NaNs or infinities, d02rac is likely to produce unexpected results.
Note that if jacobg is supplied then jacobf must also be supplied.
16: $\mathbf{jaceps}$ – function, supplied by the userExternal Function
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.
On exit: ${\mathbf{f}}\left[\mathit{i}-1\right]$ must contain the value $\frac{\partial {f}_{\mathit{i}}}{\partial \epsilon}$ at the point $(x,y)$, for $\mathit{i}=1,2,\dots ,n$.
6: $\mathbf{comm}$ – Nag_User *
Pointer to structure of type Nag_User.
p – Pointer
The type Pointer will be void *. Before calling d02rac these pointers may be allocated memory and initialized with various quantities for use by comm when called from d02rac.
Note:jaceps should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02rac. If your code inadvertently does return any NaNs or infinities, d02rac is likely to produce unexpected results.
Note that if jaceps is defined then jacgep must also be defined.
17: $\mathbf{jacgep}$ – function, supplied by the userExternal Function
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.
On exit: ${\mathbf{bcep}}\left[\mathit{i}-1\right]$ must contain the value of $\frac{\partial {g}_{\mathit{i}}}{\partial \epsilon}$, for $\mathit{i}=1,2,\dots ,n$.
6: $\mathbf{comm}$ – Nag_User *
Pointer to structure of type Nag_User.
p – Pointer
The type Pointer will be void *. Before calling d02rac these pointers may be allocated memory and initialized with various quantities for use by comm when called from d02rac.
Note:jacgep should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02rac. If your code inadvertently does return any NaNs or infinities, d02rac is likely to produce unexpected results.
Note that if jacgep is defined then jaceps must also be defined.
18: $\mathbf{comm}$ – Nag_User *
Pointer to structure of type Nag_User.
p – Pointer
The type Pointer will be void *. Before calling d02rac these pointers may be allocated memory and initialized with various quantities for use by comm when called from d02rac.
19: $\mathbf{fail}$ – NagError *Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).
6Error Indicators and Warnings
NE_2_INT_ARG_ZERO
On entry, ${\mathbf{numbeg}}=0$ and ${\mathbf{nummix}}=0$. These arguments must not both be zero.
NE_2_REAL_ARG_LE
On entry, ${\mathbf{x}}\left[0\right]=\u27e8\mathit{\text{value}}\u27e9$ and ${\mathbf{x}}\left[{\mathbf{np}}-1\right]=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{x}}\left[0\right]<{\mathbf{x}}\left[{\mathbf{np}}-1\right]$.
Convergence failure. There are a number of possible causes.
(a)Faulty coding of the Jacobian calculation functions.
(b)If Jacobians have not been supplied then inaccurate Jacobians have been calculated internally (not very likely).
(c)A poor choice of initial mesh or initial starting conditions either by the user or by default. Try using the continuation facility.
The Newton iteration has failed to converge.
This could be due to there being too few points in the initial mesh or to
the initial approximate solution being too inaccurate.
If this latter reason is suspected or you cannot make changes to prevent
this error, you should use the function with a continuation facility instead.
NE_CONV_CONT_DELEPS
deleps is required to be less than machine precision for continuation to proceed. It is likely that either the problem has no solution for some value near the current value of $\epsilon $ or that the problem is so difficult that even with continuation it is unlikely to be solved using this function. Using more mesh points may help.
The continuation step is required to be less than machine precision for
continuation to proceed. It is likely that either the problem has no
solution for some value of the continuation parameter near the current
value or that the problem is so difficult that even with continuation
it is unlikely to be solved using this function. In the latter case
using more mesh points initially may help.
NE_CONV_CONT_DEP
There is no dependence on the continuation parameter when continuation
is being used. This can be due to faulty coding of derivatives with
respect to the continuation parameter or to a zero initial choice of
approximate solution.
There is no dependence on $\epsilon $ when continuation is being used. This may 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}}=\mathrm{Nag\_DefInitMesh}$).
NE_CONV_JACOBG
The Jacobian calculated by jacobg (or the equivalent matrix calculated by numerical differentiation) is singular. This may be 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}}=\mathrm{Nag\_DefInitMesh}$).
The Jacobian for the boundary conditions is singular.
This may occur due to faulty coding of the Jacobian or, in some
circumstances, to a zero initial choice of approximate solution.
NE_CONV_MESH
A finer mesh is required for the accuracy requested; that is,
${\mathbf{mnp}}=\u27e8\mathit{\text{value}}\u27e9$ is not large enough.
NE_CONV_ROUNDOFF
Newton iteration has reached round-off level.
If desired accuracy has not been reached, then tol is too
small for this problem and this machine precision.
Solution cannot be improved due to roundoff error. Too much accuracy might have been requested.
NE_INT_ARG_LT
On entry, ${\mathbf{mnp}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{mnp}}\ge 32$.
On entry, ${\mathbf{neq}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{neq}}\ge 1$.
On entry, ${\mathbf{np}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{np}}\ge 4$.
On entry, ${\mathbf{numbeg}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{numbeg}}\ge 0$.
On entry, ${\mathbf{nummix}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{nummix}}\ge 0$.
NE_INT_RANGE_CONS
On entry, ${\mathbf{np}}=\u27e8\mathit{\text{value}}\u27e9$ and ${\mathbf{mnp}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{np}}\le {\mathbf{mnp}}$.
On entry, ${\mathbf{numbeg}}=\u27e8\mathit{\text{value}}\u27e9$ and ${\mathbf{neq}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{numbeg}}<{\mathbf{neq}}$.
On entry, ${\mathbf{numbeg}}=\u27e8\mathit{\text{value}}\u27e9$, ${\mathbf{nummix}}=\u27e8\mathit{\text{value}}\u27e9$
and ${\mathbf{neq}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{numbeg}}+{\mathbf{nummix}}\le {\mathbf{neq}}$.
NE_INTERNAL_ERROR
A continuation error occurred, but continuation is not being used.
Please contact NAG.
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
A serious error occurred in a call to the internal integrator.
The error code internally was $\u27e8\mathit{\text{value}}\u27e9$.
Please contact NAG.
NE_INVALID_FUN_JAC
Only one of jacobf or jacobg has been set to non-null possibly implying user-defined Jacobians. Both must be non-null.
NE_INVALID_FUN_JAC_CONT
deleps has been set to $\u27e8\mathit{\text{value}}\u27e9$ implying continuation and both jacobf and jacobg have been set to non-null implying user-defined Jacobians. Hence the functions jaceps and jacgep must also be non-null.
NE_INVALID_FUN_JAC_NO_CONT
deleps has been set to $\u27e8\mathit{\text{value}}\u27e9$ implying no continuation and both jacobf and jacobg have been set to non-null implying user-defined Jacobians. Hence the functions jaceps and jacgep must be NULL.
NE_INVALID_FUN_NO_JAC_CONT
deleps has been set to $\u27e8\mathit{\text{value}}\u27e9$ implying continuation and both jacobf and jacobg have been set to NULL implying no user-defined Jacobians. Hence the functions jaceps and jacgep must also be NULL.
NE_NOT_STRICTLY_INCREASING
On entry the mesh points are not in strictly ascending order.
For $i=\u27e8\mathit{\text{value}}\u27e9$, mesh point $i=\u27e8\mathit{\text{value}}\u27e9$,
but mesh point $i+1=\u27e8\mathit{\text{value}}\u27e9$.
NE_REAL_ARG_LE
On entry, ${\mathbf{tol}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{tol}}>0.0$.
7Accuracy
The solution returned by the function will be accurate to your tolerance as defined by the relation (5) except in extreme circumstances. The final error estimate over the whole mesh for each component is given in the array abt. If too many points are specified in the initial mesh, the solution may be more accurate than requested and the error may not be approximately equidistributed.
8Parallelism and Performance
d02rac is not threaded in any implementation.
9Further Comments
There are too many factors present to quantify the timing. The time taken by d02rac is negligible only on very simple problems.
In the case where you wish to solve a sequence of similar problems, the use of the final mesh and solution from one case as the initial mesh is strongly recommended for the next.
to an accuracy specified by ${\mathbf{tol}}=\text{1.0e\u22124}$. The continuation facility is used with the continuation parameter $\epsilon $ introduced as in the differential equation above and with ${\mathbf{deleps}}=0.1$ initially. (The continuation facility is not needed for this problem and is used here for illustration.)