NAG CL Interface
d02rac (bvp_​fd_​nonlin_​gen)

1 Purpose

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.

2 Specification

#include <nag.h>
void  d02rac (Integer neq, double *deleps,
void (*fcn)(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm),
Integer numbeg, Integer nummix,
void (*g)(Integer neq, double eps, const double ya[], const double yb[], double bc[], Nag_User *comm),
Nag_MeshSet init, Integer mnp, Integer *np, double x[], double y[], double tol, double abt[],
void (*jacobf)(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm),
void (*jacobg)(Integer neq, double eps, const double ya[], const double yb[], double aj[], double bj[], Nag_User *comm),
void (*jaceps)(Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm),
void (*jacgep)(Integer neq, double eps, const double ya[], const double yb[], double bcep[], Nag_User *comm),
Nag_User *comm, NagError *fail)
The function may be called by the names: d02rac or nag_ode_bvp_fd_nonlin_gen.

3 Description

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
yi = f i x, y 1 , y 2 ,, y n ,   i=1,2,,n (1)
and the derivatives fi are evaluated by fcn. With the differential equations (1) must be given a system of n (nonlinear) boundary conditions
gi ya,yb = 0 ,   i=1,2,,n ,  
where
y x = y 1 x , y 2 x ,, y n x T . (2)
The functions gi 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.
The Newton iteration requires Jacobian matrices
fi yj , gi yja  and  gi yjb .  
These may be supplied through jacobf for fi yj 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
y = f x,y,ε ,   g ya,yb,ε = 0 , (3)
where f = f 1 , f 2 ,, f n T etc., and where ε is a continuation parameter. The choice ε=0 must give a problem (3) which is easy to solve and ε=1 must define the problem whose solution is actually required. The function solves a sequence of problems with ε values
0 = ε1 < ε2 < < εp = 1 . (4)
The number p and the values εi are chosen by the function so that each problem can be solved using the solution of its predecessor as a starting approximation. Jacobians f ε and g ε are required and they may be supplied by you via jaceps and jacgep respectively or may be computed by numerical differentiation.

4 References

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

5 Arguments

1: neq Integer Input
On entry: n, the number of differential equations.
Constraint: neq>0.
2: deleps double * Input/Output
On entry: must be given a value which specifies whether continuation is required. If deleps0.0 or deleps1.0 then it is assumed that continuation is not required. If 0.0<deleps<1.0 then it is assumed that continuation is required unless deleps<machine precision when an error exit is taken. deleps is used as the increment ε2-ε1 (see (4)) and the choice deleps=0.1 is recommended.
On exit: an overestimate of the increment εp-εp-1 (in fact the value of the increment which would have been tried if the restriction εp=1 had not been imposed). If continuation was not requested then 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: fcn function, supplied by the user External Function
fcn must evaluate the functions fi (i.e., the derivatives yi) at a general point x for a given value of ε, the continuation parameter (see Section 3).
The specification of fcn is:
void  fcn (Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm)
1: neq Integer Input
On entry: n, the number of equations.
2: x double Input
On entry: x, the value of the independent variable.
3: eps double Input
On entry: ε, the value of the continuation parameter. This is 1 if continuation is not being used.
4: y[neq] const double Input
On entry: yi, for i=1,2,,n, the values of the dependent variables at x.
5: f[neq] double Output
On exit: the values of the derivatives fi evaluated at x given ε, for i=1,2,,n.
6: comm Nag_User *
Pointer to structure of type Nag_User.
pPointer 
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: numbeg Integer Input
On entry: the number of left-hand boundary conditions (that is the number involving ya only).
Constraint: 0numbeg<neq.
5: nummix Integer Input
On entry: the number of coupled boundary conditions (that is the number involving both ya and yb).
Constraint: 0nummixneq-numbeg.
6: g function, supplied by the user External Function
g must evaluate the boundary conditions in equation (3) and place them in the array bc.
The specification of g is:
void  g (Integer neq, double eps, const double ya[], const double yb[], double bc[], Nag_User *comm)
1: neq Integer Input
On entry: n, the number of equations.
2: eps double Input
On entry: ε, the value of the continuation parameter. This is 1 if continuation is not being used.
3: ya[neq] const double Input
On entry: the value yia, for i=1,2,,n.
4: yb[neq] const double Input
On entry: the value yib, for i=1,2,,n.
5: bc[neq] double Output
On exit: the values gi ya,yb,ε , for i=1,2,,n. These must be ordered as follows:
  1. (i)first, the conditions involving only ya (see numbeg);
  2. (ii)next, the nummix coupled conditions involving both ya and yb (see nummix); and,
  3. (iii)finally, the conditions involving only yb (neq-numbeg-nummix).
6: comm Nag_User *
Pointer to structure of type Nag_User.
pPointer 
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: init Nag_MeshSet Input
On entry: indicates whether you wish to supply an initial mesh and approximate solution (init=Nag_UserInitMesh) or whether default values are to be used, (init=Nag_DefInitMesh).
Constraint: init=Nag_DefInitMesh or Nag_UserInitMesh.
8: mnp Integer Input
On entry: mnp must be set to the maximum permitted number of points in the finite difference mesh.
Constraint: mnp32.
9: np Integer * Input/Output
On entry: must be set to the number of points to be used in the initial mesh.
Constraint: 4npmnp.
On exit: the number of points in the final mesh.
10: x[mnp] double Input/Output
On entry: you must set x[0]=a and x[np-1]=b. If init=Nag_DefInitMesh on entry a default equispaced mesh will be used, otherwise you must specify a mesh by setting x[i-1]=xi, for i=2,3,,np-1.
Constraints:
  • if init=Nag_DefInitMesh, x[0]<x[np-1];
  • if init=Nag_UserInitMesh, x[0]<x[1]<<x[np-1].
On exit: x[0],x[1],,x[np-1] define the final mesh (with the returned value of np) and x[0]=a and x[np-1]=b.
11: y[neq×mnp] double Input/Output
On entry: if init=Nag_DefInitMesh, then y need not be set.
If init=Nag_UserInitMesh, then the array y must contain an initial approximation to the solution such that y[j-1×mnp+i-1] contains an approximation to
yjxi,  i=1,2,,np​ and ​j=1,2,,n.  
On exit: the approximate solution zjxi satisfying (5) on the final mesh, that is
y[j-1×mnp+i-1]=zjxi,  i=1,2,,np​ and ​j=1,2,,n,  
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.
12: tol double Input
On entry: a positive absolute error tolerance. If
a=x1<x2<<xnp=b  
is the final mesh, zjxi is the jth component of the approximate solution at xi, and yjx is the jth component of the true solution of (1) and (2), then, except in extreme circumstances, it is expected that
zjxi-yjxitol,  i=1,2,,np​ and ​j=1,2,,n. (5)
Constraint: tol>0.0.
13: abt[neq] double Output
On exit: abt[i-1], for i=1,2,,n, holds the largest estimated error (in magnitude) of the ith component of the solution over all mesh points.
14: jacobf function, supplied by the user External Function
jacobf evaluates the Jacobian fi yj , for i=1,2,,n and j=1,2,,n, given x and yj, for j=1,2,,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:
void  jacobf (Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm)
1: neq Integer Input
On entry: n, the number of equations.
2: x double Input
On entry: x, the value of the independent variable.
3: eps double Input
On entry: ε, the value of the continuation parameter. This is 1 if continuation is not being used.
4: y[neq] const double Input
On entry: yi, for i=1,2,,n, the values of the dependent variables at x.
5: f[neq×neq] double Output
On exit: f[j-1×neq+i-1] must be set to the value of fi yj , evaluated at the point x,y, for i=1,2,,n and j=1,2,,n.
6: comm Nag_User *
Pointer to structure of type Nag_User.
pPointer 
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: jacobg function, supplied by the user External Function
jacobg evaluates the Jacobians gi yja and gi yjb . 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:
void  jacobg (Integer neq, double eps, const double ya[], const double yb[], double aj[], double bj[], Nag_User *comm)
1: neq Integer Input
On entry: n, the number of equations.
2: eps double Input
On entry: ε, the value of the continuation parameter. This is 1 if continuation is not being used.
3: ya[neq] const double Input
On entry: the value yia, for i=1,2,,n.
4: yb[neq] const double Input
On entry: the value yib, for i=1,2,,n.
5: aj[neq×neq] double Output
On exit: aj[i-1×neq+j-1] must be set to the value gi yja , for i=1,2,,n and j=1,2,,n.
6: bj[neq×neq] double Output
On exit: bj[i-1×neq+j-1] must be set to the value gi yjb , for i=1,2,,n and j=1,2,,n.
7: comm Nag_User *
Pointer to structure of type Nag_User.
pPointer 
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: jaceps function, supplied by the user External Function
jaceps evaluates the derivative fi ε 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:
void  jaceps (Integer neq, double x, double eps, const double y[], double f[], Nag_User *comm)
1: neq Integer Input
On entry: n, the number of equations.
2: x double Input
On entry: x, the value of the independent variable.
3: eps double Input
On entry: ε, the value of the continuation parameter.
4: y[neq] const double Input
On entry: the solution values yi, for i=1,2,,n, at the point x.
5: f[neq] double Output
On exit: f[i-1] must contain the value fi ε at the point x,y, for i=1,2,,n.
6: comm Nag_User *
Pointer to structure of type Nag_User.
pPointer 
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: jacgep function, supplied by the user External Function
jacgep evaluates the derivatives gi ε 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:
void  jacgep (Integer neq, double eps, const double ya[], const double yb[], double bcep[], Nag_User *comm)
1: neq Integer Input
On entry: n, the number of equations.
2: eps double Input
On entry: ε, the value of the continuation parameter.
3: ya[neq] const double Input
On entry: the value of yia, for i=1,2,,n.
4: yb[neq] const double Input
On entry: the value of yib, for i=1,2,,n.
5: bcep[neq] double Output
On exit: bcep[i-1] must contain the value of gi ε , for i=1,2,,n.
6: comm Nag_User *
Pointer to structure of type Nag_User.
pPointer 
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: comm Nag_User *
Pointer to structure of type Nag_User.
pPointer 
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: fail NagError * Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

6 Error Indicators and Warnings

NE_2_INT_ARG_ZERO
On entry, numbeg=0 and nummix=0 . These arguments must not both be zero.
NE_2_REAL_ARG_LE
On entry, x[0]=value and x[np-1]=value.
Constraint: x[0]<x[np-1].
NE_ALLOC_FAIL
Dynamic memory allocation failed.
NE_BAD_PARAM
On entry, argument init had an illegal value.
NE_CONV_CONT
Convergence failure. There are a number of possible causes.
  1. (a)Faulty coding of the Jacobian calculation functions.
  2. (b)If Jacobians have not been supplied then inaccurate Jacobians have been calculated internally (not very likely).
  3. (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 ε 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 ε 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 init=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 init=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, mnp=value 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, mnp=value.
Constraint: mnp32.
On entry, neq=value.
Constraint: neq1.
On entry, np=value.
Constraint: np4.
On entry, numbeg=value.
Constraint: numbeg0.
On entry, nummix=value.
Constraint: nummix0.
NE_INT_RANGE_CONS
On entry, np=value and mnp=value.
Constraint: npmnp.
On entry, numbeg=value and neq=value.
Constraint: numbeg<neq.
On entry, numbeg=value, nummix=value
and neq=value.
Constraint: numbeg+nummixneq.
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 value.
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 value 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 value 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 value 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=value, mesh point i=value, but mesh point i+1=value.
NE_REAL_ARG_LE
On entry, tol=value.
Constraint: tol>0.0.

7 Accuracy

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.

8 Parallelism and Performance

d02rac is not threaded in any implementation.

9 Further 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.

10 Example

This example solves the differential equation
y=-yy-2ε1-y2  
with ε=1 and boundary conditions
y0=y0=0,  y10=1  
to an accuracy specified by tol=1.0e−4. The continuation facility is used with the continuation parameter ε introduced as in the differential equation above and with deleps=0.1 initially. (The continuation facility is not needed for this problem and is used here for illustration.)

10.1 Program Text

Program Text (d02race.c)

10.2 Program Data

None.

10.3 Program Results

Program Results (d02race.r)
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 0 2 4 6 8 10 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 y y' and y'' x Example Program Solution of Third-order BVP y y' y'' gnuplot_plot_1 gnuplot_plot_2 gnuplot_plot_3