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
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),
|
|
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
ordinary differential equations in the interval
with
. The system is written in the form
and the derivatives
are evaluated by
fcn. With the differential equations
(1) must be given a system of
(nonlinear) boundary conditions
where
The functions
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
These may be supplied through
jacobf for
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
etc., and where
is a continuation parameter. The choice
must give a problem
(3) which is easy to solve and
must define the problem whose solution is actually required. The function solves a sequence of problems with
values
The number
and the values
are chosen by the function so that each problem can be solved using the solution of its predecessor as a starting approximation. Jacobians
and
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:
– Integer
Input
-
On entry:
, the number of differential equations.
Constraint:
.
-
2:
– double *
Input/Output
-
On entry: must be given a value which specifies whether continuation is required. If
or
then it is assumed that continuation is not required. If
then it is assumed that continuation is required unless
when an error exit is taken.
deleps is used as the increment
(see
(4)) and the choice
is recommended.
On exit: an overestimate of the increment
(in fact the value of the increment which would have been tried if the restriction
had not been imposed). If continuation was not requested then
.
If continuation is not requested then
jaceps and
jacgep may each be replaced by
the NAG defined null function pointer NULLFN.
-
3:
– function, supplied by the user
External Function
-
fcn must evaluate the functions
(i.e., the derivatives
) at a general point
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:
– Integer
Input
-
On entry: , the number of equations.
-
2:
– double
Input
-
On entry: , the value of the independent variable.
-
3:
– double
Input
-
On entry: , the value of the continuation parameter. This is if continuation is not being used.
-
4:
– const double
Input
-
On entry:
, for , the values of the dependent variables at .
-
5:
– double
Output
-
On exit: the values of the derivatives
evaluated at given , for .
-
6:
– 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:
– Integer
Input
-
On entry: the number of left-hand boundary conditions (that is the number involving only).
Constraint:
.
-
5:
– Integer
Input
-
On entry: the number of coupled boundary conditions (that is the number involving both and ).
Constraint:
.
-
6:
– 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:
– Integer
Input
-
On entry:
, the number of equations.
-
2:
– double
Input
-
On entry: , the value of the continuation parameter. This is if continuation is not being used.
-
3:
– const double
Input
-
On entry: the value
, for .
-
4:
– const double
Input
-
On entry: the value
, for .
-
5:
– double
Output
-
On exit: the values
, for
. These must be ordered as follows:
-
(i)first, the conditions involving only (see numbeg);
-
(ii)next, the nummix coupled conditions involving both and (see nummix); and,
-
(iii)finally, the conditions involving only ().
-
6:
– 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:
– Nag_MeshSet
Input
-
On entry: indicates whether you wish to supply an initial mesh and approximate solution () or whether default values are to be used, ().
Constraint:
or .
-
8:
– Integer
Input
-
On entry:
mnp must be set to the maximum permitted number of points in the finite difference mesh.
Constraint:
.
-
9:
– Integer *
Input/Output
-
On entry: must be set to the number of points to be used in the initial mesh.
Constraint:
.
On exit: the number of points in the final mesh.
-
10:
– double
Input/Output
-
On entry: you must set and . If on entry a default equispaced mesh will be used, otherwise you must specify a mesh by setting
, for .
Constraints:
- if , ;
- if , .
On exit:
define the final mesh (with the returned value of
np) and
and
.
-
11:
– double
Input/Output
-
On entry: if
, then
y need not be set.
If
, then the array
y must contain an initial approximation to the solution such that
contains an approximation to
On exit: the approximate solution
satisfying
(5) on the final mesh, that is
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:
– double
Input
-
On entry: a positive absolute error tolerance. If
is the final mesh,
is the
th component of the approximate solution at
, and
is the
th component of the true solution of
(1) and
(2), then, except in extreme circumstances, it is expected that
Constraint:
.
-
13:
– double
Output
-
On exit:
, for , holds the largest estimated error (in magnitude) of the th component of the solution over all mesh points.
-
14:
– function, supplied by the user
External Function
-
jacobf evaluates the Jacobian
, for
and
, given
and
, for
.
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:
– Integer
Input
-
On entry:
, the number of equations.
-
2:
– double
Input
-
On entry: , the value of the independent variable.
-
3:
– double
Input
-
On entry: , the value of the continuation parameter. This is if continuation is not being used.
-
4:
– const double
Input
-
On entry:
, for , the values of the dependent variables at .
-
5:
– double
Output
-
On exit: must be set to the value of , evaluated at the point , for and .
-
6:
– 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:
– function, supplied by the user
External Function
-
jacobg evaluates the Jacobians
and
. 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:
– Integer
Input
-
On entry:
, the number of equations.
-
2:
– double
Input
-
On entry: , the value of the continuation parameter. This is if continuation is not being used.
-
3:
– const double
Input
-
On entry: the value
, for .
-
4:
– const double
Input
-
On entry: the value
, for .
-
5:
– double
Output
-
On exit: must be set to the value , for and .
-
6:
– double
Output
-
On exit: must be set to the value , for and .
-
7:
– 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:
– function, supplied by the user
External Function
-
jaceps evaluates the derivative
given
and
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:
– Integer
Input
-
On entry:
, the number of equations.
-
2:
– double
Input
-
On entry: , the value of the independent variable.
-
3:
– double
Input
-
On entry: , the value of the continuation parameter.
-
4:
– const double
Input
-
On entry: the solution values
, for , at the point .
-
5:
– double
Output
-
On exit: must contain the value at the point , for .
-
6:
– 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:
– function, supplied by the user
External Function
-
jacgep evaluates the derivatives
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:
– Integer
Input
-
On entry:
, the number of equations.
-
2:
– double
Input
-
On entry: , the value of the continuation parameter.
-
3:
– const double
Input
-
On entry: the value of
, for .
-
4:
– const double
Input
-
On entry: the value of
, for .
-
5:
– double
Output
-
On exit: must contain the value of , for .
-
6:
– 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:
– 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:
– 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, and . These arguments must not both be zero.
- NE_2_REAL_ARG_LE
-
On entry, and .
Constraint: .
- 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.
-
(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
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
).
- 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
).
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,
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, .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
On entry, .
Constraint: .
- NE_INT_RANGE_CONS
-
On entry, and .
Constraint: .
On entry, and .
Constraint: .
On entry, ,
and .
Constraint: .
- 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
.
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
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
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
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 , mesh point ,
but mesh point .
- NE_REAL_ARG_LE
-
On entry, .
Constraint: .
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.
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
with
and boundary conditions
to an accuracy specified by
. The continuation facility is used with the continuation parameter
introduced as in the differential equation above and with
initially. (The continuation facility is not needed for this problem and is used here for illustration.)
10.1
Program Text
10.2
Program Data
None.
10.3
Program Results