NAG FL Interface
d02hbf (bvp_​shoot_​genpar)

Settings help

FL Name Style:


FL Specification Language:


1 Purpose

d02hbf solves a two-point boundary value problem for a system of ordinary differential equations, using initial value techniques and Newton iteration; it generalizes subroutine d02haf to include the case where parameters other than boundary values are to be determined.

2 Specification

Fortran Interface
Subroutine d02hbf ( p, n1, pe, e, n, soln, m1, fcn, bc, range, w, sdw, ifail)
Integer, Intent (In) :: n1, n, m1, sdw
Integer, Intent (Inout) :: ifail
Real (Kind=nag_wp), Intent (In) :: pe(n1), e(n)
Real (Kind=nag_wp), Intent (Inout) :: p(n1)
Real (Kind=nag_wp), Intent (Out) :: soln(n,m1), w(n,sdw)
External :: fcn, bc, range
C Header Interface
#include <nag.h>
void  d02hbf_ (double p[], const Integer *n1, const double pe[], const double e[], const Integer *n, double soln[], const Integer *m1,
void (NAG_CALL *fcn)(const double *x, const double y[], double f[], const double p[]),
void (NAG_CALL *bc)(double g1[], double g2[], const double p[]),
void (NAG_CALL *range)(double *a, double *b, const double p[]),
double w[], const Integer *sdw, Integer *ifail)
The routine may be called by the names d02hbf or nagf_ode_bvp_shoot_genpar.

3 Description

d02hbf solves a two-point boundary value problem by determining the unknown parameters p1,p2,,pn1 of the problem. These parameters may be, but need not be, boundary values; they may include eigenvalue parameters in the coefficients of the differential equations, length of the range of integration, etc. The notation and methods used are similar to those of d02haf and you are advised to study this first. (The parameters p1,p2,,pn1 correspond precisely to the unknown boundary conditions in d02haf.) It is assumed that we have a system of n first-order ordinary differential equations of the form:
dyi dx =fi(x,y1,y2,,yn),  i=1,2,,n,  
and that the derivatives fi are evaluated by fcn. The system, including the boundary conditions given by bc and the range of integration given by range, involves the n1 unknown parameters p1,p2,,pn1 which are to be determined, and for which initial estimates must be supplied. The number of unknown parameters n1 must not exceed the number of equations n. If n1<n, we assume that (n-n1) equations of the system are not involved in the matching process. These are usually referred to as ‘driving equations’; they are independent of the parameters and of the solutions of the other n1 equations. In numbering the equations for fcn, the driving equations must be put first.
The estimated values of the parameters are corrected by a form of Newton iteration. The Newton correction on each iteration is calculated using a Jacobian matrix whose (i,j)th element depends on the derivative of the ith component of the solution, yi, with respect to the jth parameter, pj. This matrix is calculated by a simple numerical differentiation technique which requires n1 evaluations of the differential system.
If the argument ifail is set appropriately, the routine automatically prints messages to inform you of the flow of the calculation. These messages are discussed in detail in Section 9.
d02hbf is a simplified version of d02saf which is described in detail in Gladwell (1979).

4 References

Gladwell I (1979) The development of the boundary value codes in the ordinary differential equations chapter of the NAG Library 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

You are strongly recommended to read Sections 3 and 9 in conjunction with this section.
1: p(n1) Real (Kind=nag_wp) array Input/Output
On entry: an estimate for the ith argument, pi, for i=1,2,,n1.
On exit: the corrected value for the ith argument, unless an error has occurred, when it contains the last calculated value of the argument.
2: n1 Integer Input
On entry: n1, the number of arguments.
Constraint: 1n1n.
3: pe(n1) Real (Kind=nag_wp) array Input
On entry: the elements of pe must be given small positive values. The element pe(i) is used
  1. (i)in the convergence test on the ith argument in the Newton iteration, and
  2. (ii)in perturbing the ith argument when approximating the derivatives of the components of the solution with respect to this argument for use in the Newton iteration.
The elements pe(i) should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.
Constraint: pe(i)>0.0, for i=1,2,,n1.
4: e(n) Real (Kind=nag_wp) array Input
On entry: the elements of e must be given positive values. The element e(i) is used in the bound on the local error in the ith component of the solution yi during integration.
The elements e(i) should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.
Constraint: e(i)>0.0, for i=1,2,,n.
5: n Integer Input
On entry: n, the total number of differential equations.
Constraint: nn1.
6: soln(n,m1) Real (Kind=nag_wp) array Output
On exit: the solution when m1>1.
7: m1 Integer Input
On entry: a value which controls exit values.
m1=1
The final solution is not calculated.
m1>1
The final values of the solution at interval (length of range)/(m1-1) are calculated and stored sequentially in the array soln starting with the values of the solutions evaluated at the first end point (see range) stored in the first column of soln.
Constraint: m11.
8: fcn Subroutine, supplied by the user. External Procedure
fcn must evaluate the functions fi (i.e., the derivatives yi), for i=1,2,,n, at a general point x.
The specification of fcn is:
Fortran Interface
Subroutine fcn ( x, y, f, p)
Real (Kind=nag_wp), Intent (In) :: x, y(*), p(*)
Real (Kind=nag_wp), Intent (Out) :: f(*)
C Header Interface
void  fcn (const double *x, const double y[], double f[], const double p[])
In the description of the arguments of d02hbf below, n and n1 denote the numerical values of n and n1 in the call of d02hbf.
1: x Real (Kind=nag_wp) Input
On entry: x, the value of the argument.
2: y(*) Real (Kind=nag_wp) array Input
On entry: yi, for i=1,2,,n, the value of the argument.
3: f(*) Real (Kind=nag_wp) array Output
On exit: the value of fi, for i=1,2,,n. The fi may depend upon the parameters pj, for j=1,2,,n1. If there are any driving equations (see Section 3) then these must be numbered first in the ordering of the components of f in fcn.
4: p(*) Real (Kind=nag_wp) array Input
On entry: the current estimate of the argument pi, for i=1,2,,n1.
fcn must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02hbf is called. Arguments denoted as Input must not be changed by this procedure.
Note: fcn should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02hbf. If your code inadvertently does return any NaNs or infinities, d02hbf is likely to produce unexpected results.
9: bc Subroutine, supplied by the user. External Procedure
bc must place in g1 and g2 the boundary conditions at a and b respectively (see range).
The specification of bc is:
Fortran Interface
Subroutine bc ( g1, g2, p)
Real (Kind=nag_wp), Intent (In) :: p(*)
Real (Kind=nag_wp), Intent (Out) :: g1(*), g2(*)
C Header Interface
void  bc (double g1[], double g2[], const double p[])
In the description of the arguments of d02hbf below, n and n1 denote the numerical values of n and n1 in the call of d02hbf.
1: g1(*) Real (Kind=nag_wp) array Output
On exit: the value of yi(a), (where this may be a known value or a function of the parameters pj, for i=1,2,,n and j=1,2,,n1).
2: g2(*) Real (Kind=nag_wp) array Output
On exit: the value of yi(b), for i=1,2,,n, (where these may be known values or functions of the parameters pj, for j=1,2,,n1). If n>n1, so that there are some driving equations, the first n-n1 values of g2 need not be set since they are never used.
3: p(*) Real (Kind=nag_wp) array Input
On entry: an estimate of the argument pi, for i=1,2,,n1.
bc must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02hbf is called. Arguments denoted as Input must not be changed by this procedure.
Note: bc should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02hbf. If your code inadvertently does return any NaNs or infinities, d02hbf is likely to produce unexpected results.
10: range Subroutine, supplied by the user. External Procedure
range must evaluate the boundary points a and b, each of which may depend on the arguments p1,p2,,pn1. The integrations in the shooting method are always from a to b.
The specification of range is:
Fortran Interface
Subroutine range ( a, b, p)
Real (Kind=nag_wp), Intent (In) :: p(*)
Real (Kind=nag_wp), Intent (Out) :: a, b
C Header Interface
void  range (double *a, double *b, const double p[])
In the description of the arguments of d02hbf below, n1 denotes the actual value of n1 in the call of d02hbf.
1: a Real (Kind=nag_wp) Output
On exit: a, one of the boundary points.
2: b Real (Kind=nag_wp) Output
On exit: the second boundary point, b. Note that b>a forces the direction of integration to be that of increasing x. If a and b are interchanged the direction of integration is reversed.
3: p(*) Real (Kind=nag_wp) array Input
On entry: the current estimate of the ith argument, pi, for i=1,2,,n1.
range must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d02hbf is called. Arguments denoted as Input must not be changed by this procedure.
Note: range should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d02hbf. If your code inadvertently does return any NaNs or infinities, d02hbf is likely to produce unexpected results.
11: w(n,sdw) Real (Kind=nag_wp) array Output
Used mainly for workspace.
On exit: with ifail=2, 3, 4 or 5 (see Section 6), w(i,1), for i=1,2,,n, contains the solution at the point x when the error occurred. w(1,2) contains x.
12: sdw Integer Input
On entry: the second dimension of the array w as declared in the (sub)program from which d02hbf is called.
Constraint: sdw3n+14+max(11,n).
13: ifail Integer Input/Output
This routine uses an ifail input value codification that differs from the normal case to distinguish between errors and warnings (see Section 4 in the Introduction to the NAG Library FL Interface).
On entry: ifail must be set to one of the values below to set behaviour on detection of an error; these values have no effect when no error is detected. The behaviour relate to whether or not program execution is halted and whether or not messages are printed when an error or warning is detected.
ifail Execution Error Printing Warning Printed
000 halted No No
001 continue No No
010 halted Yes No
011 continue Yes No
100 halted No Yes
101 continue No Yes
110 halted Yes Yes
111 continue Yes Yes
For environments where it might be inappropriate to halt program execution when an error is detected, the value 1, 11, 101 or 111 is recommended. If the printing of messages is undesirable, then the value 1 is recommended. Otherwise, the recommended value is 110. When the value 1, 11, 101 or 111 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or −1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, m1=value.
Constraint: m11.
On entry, n1=value.
Constraint: n11.
On entry, n=value and n1=value.
Constraint: nn1.
On entry, sdw=value and n=value.
Constraint: sdw3×n+max(25,n+14); that is, value.
On entry a negative or zero convergence test tolerance has been set.
On entry a negative or zero local error tolerance has been set.
ifail=2
In the integration with initial or final parameters, the step size was reduced too far for the integration to proceed. Either this routine is not a suitable method for solving the problem, or the initial choice of parameters is very poor.
The step length for the integration became too short to proceed when calculating the residual.
ifail=3
An initial step-length could be found for integration to proceed with the current parameters.
In the integration with initial or final parameters, a suitable initial step could not be found. Either this routine is not suitable for solving the problem, or the initial choice of parameters is very poor.
ifail=4
The step-length required to calculate the Jacobian to sufficient accuracy is too small
ifail=5
An initial step-length could be found for Jacobian calculation to proceed with the current parameters.
ifail=6
The Jacobian has an insignificant column. Make sure that the solution vector depends on all the parameters.
ifail=7
An internal singular value decomposition has failed.
This error can be avoided by changing the initial parameter estimates.
ifail=8
The Newton iteration has failed to converge.
This can indicate a poor initial choice of parameters or a very difficult problem.
Consider varying elements of the parameter convergence control if the residuals are small; otherwise vary initial parameter estimates.
ifail=9
Internal error in Newton method. Please contact NAG.
ifail=10
Internal error in calculating Jacobian. Please contact NAG.
ifail=11
Internal error in calculating residual. Please contact NAG.
ifail=12
Internal error in calculating residual. Please contact NAG.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

If the process converges, the accuracy to which the unknown parameters are determined is usually close to that specified by you; the solution, if requested, may be determined to a required accuracy by varying e.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
d02hbf is not thread safe and should not be called from a multithreaded user program. Please see Section 1 in FL Interface Multithreading for more information on thread safety.
d02hbf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

The time taken by d02hbf depends on the complexity of the system, and on the number of iterations required. In practice, integration of the differential equations is by far the most costly process involved.
Wherever they occur in the routine, the error arguments contained in the arrays e and pe are used in ‘mixed’ form; that is e(i) always occurs in expressions of the form
e(i)×(1+|yi|)  
and pe(i) always occurs in expressions of the form
pe(i)×(1+|pi|).  
Though not ideal for every application, it is expected that this mixture of absolute and relative error testing will be adequate for most purposes.
You may determine a suitable direction of integration a to b and suitable values for e(i) by integrations with d02pef. The best direction of integration is usually the direction of decreasing solutions.
You are strongly recommended to set ifail to obtain self-explanatory error messages, and also monitoring information about the course of the computation. You may select the unit numbers on which this output is to appear by calls of x04aaf (for error messages) or x04abf (for monitoring information). Otherwise the default unit numbers will be used, as specified in the Users' Note. The monitoring information produced at each iteration includes the current parameter values, the residuals and 2-norms: a basic norm and a current norm. At each iteration the aim is to find parameter values which make the current norm less than the basic norm. Both these norms should tend to zero as should the residuals. (They would all be zero if the exact parameters were used as input.) For more details, in particular about the other monitoring information printed, you are advised to consult the specification of d02saf, and especially the description of the argument monit there.
The computing time for integrating the differential equations can sometimes depend critically on the quality of the initial estimates for the parameters pi. If it seems that too much computing time is required and, in particular, if the values of the residuals printed by the monitoring routine are much larger than the expected values of the solution at b, then the coding of fcn, bc and range should be checked for errors. If no errors can be found, an independent attempt should be made to improve the initial estimates for pi.
The subroutine can be used to solve a very wide range of problems, for example:
  1. (a)eigenvalue problems, including problems where the eigenvalue occurs in the boundary conditions;
  2. (b)problems where the differential equations depend on some parameters which are to be determined so as to satisfy certain boundary conditions;
  3. (c)problems where one of the end points of the range of integration is to be determined as the point where a variable yi takes a particular value;
  4. (d)singular problems and problems on infinite ranges of integration where the values of the solution at a or b or both are determined by a power series or an asymptotic expansion (or a more complicated expression) and where some of the coefficients in the expression are to be determined; and
  5. (e)differential equations with certain terms defined by other independent (driving) differential equations.

10 Example

For this routine two examples are presented. There is a single example program for d02hbf, with a main program and the code to solve the two example problems given in Example 1 (EX1) and Example 2 (EX2).
Example 1 (EX1)
This example finds the solution of the differential equation
y=(y3-y)/2x  
on the range 0x16, with boundary conditions y(0)=0.1 and y(16)=1/6. We cannot use the differential equation at x=0 because it is singular, so we take a truncated power series expansion
y(x)=1/10+p1×x/10+x/100  
near the origin where p1 is one of the parameters to be determined. We choose the interval as [0.1,16] and setting p2=y(16), we can determine all the boundary conditions. We take X1=16. We write y=y(1), y=y(2), and estimate PARAM(1)=0.2, PARAM(2)=0.0. Note the call to x04abf before the call to d02hbf.
Example 2 (EX2)
This example finds the gravitational constant p1 and the range p2 over which a projectile must be fired to hit the target with a given velocity.
The differential equations are
y = tanϕ v = -(p1sinϕ+0.00002v2) vcosϕ ϕ = -p1v2  
on the range 0<x<p2, with boundary conditions
y=0, v=500, ϕ=0.5 at  x=0, y=0, v=450, ϕ=p3 at  x=p2.  
We write y=y(1), v=y(2), ϕ=y(3). We estimate p1=PARAM(1)=32, p2=PARAM(2)=6000 and p3=PARAM(3)=0.54 (though this last estimate is not important).

10.1 Program Text

Program Text (d02hbfe.f90)

10.2 Program Data

Program Data (d02hbfe.d)

10.3 Program Results

Program Results (d02hbfe.r)
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0 2 4 6 8 10 12 14 16 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 Solution Derivative x Example Program 1 Parameterised Two-point Boundary-value Problem using Initial Value Techniques and Newton Iteration solution y(x) derivative y'(x) param(2) gnuplot_plot_1 gnuplot_plot_2
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 0 100 200 300 400 500 600 700 800 900 0 1000 2000 3000 4000 5000 6000 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 Height and Velocity Angle x Example Program 2 Find Gravitational Constant and Range Given Projectile Terminal Velocity height velocity angle gnuplot_plot_1 gnuplot_plot_2 gnuplot_plot_3