NAG Library Routine Document
D02TVF
1 Purpose
D02TVF is a setup routine which must be called prior to the first call of the nonlinear twopoint boundary value solver
D02TKF.
2 Specification
SUBROUTINE D02TVF ( 
NEQ, M, NLBC, NRBC, NCOL, TOLS, MXMESH, NMESH, MESH, IPMESH, RWORK, LRWORK, IWORK, LIWORK, IFAIL) 
INTEGER 
NEQ, M(NEQ), NLBC, NRBC, NCOL, MXMESH, NMESH, IPMESH(MXMESH), LRWORK, IWORK(LIWORK), LIWORK, IFAIL 
REAL (KIND=nag_wp) 
TOLS(NEQ), MESH(MXMESH), RWORK(LRWORK) 

3 Description
D02TVF and its associated routines (
D02TKF,
D02TXF,
D02TYF and
D02TZF) solve the twopoint boundary value problem for a nonlinear system of ordinary differential equations
over an interval
$\left[a,b\right]$ subject to
$p$ (
$\text{}>0$) nonlinear boundary conditions at
$a$ and
$q$ (
$\text{}>0$) nonlinear boundary conditions at
$b$, where
$p+q={\displaystyle \sum _{i=1}^{n}}{m}_{i}$. Note that
${y}_{i}^{\left(m\right)}\left(x\right)$ is the
$m$th derivative of the
$i$th solution component. Hence
${y}_{i}^{\left(0\right)}\left(x\right)={y}_{i}\left(x\right)$. The left boundary conditions at
$a$ are defined as
and the right boundary conditions at
$b$ as
where
$y=\left({y}_{1},{y}_{2},\dots ,{y}_{n}\right)$ and
See
Section 8 for information on how boundary value problems of a more general nature can be treated.
D02TVF is used to specify an initial mesh, error requirements and other details.
D02TKF is then used to solve the boundary value problem.
The solution routine
D02TKF proceeds as follows. A modified Newton method is applied to the equations
and the boundary conditions. To solve these equations numerically the components
${y}_{i}$ are approximated by piecewise polynomials
${v}_{ij}$ using a monomial basis on the
$j$th mesh subinterval. The coefficients of the polynomials
${v}_{ij}$ form the unknowns to be computed. Collocation is applied at Gaussian points
where
${x}_{jk}$ is the
$k$th collocation point in the
$j$th mesh subinterval. Continuity at the mesh points is imposed, that is
where
${x}_{j+1}$ is the righthand end of the
$j$th mesh subinterval. The linearized collocation equations and boundary conditions, together with the continuity conditions, form a system of linear algebraic equations which are solved using
F01LHF and
F04LHF. For use in the modified Newton method, an approximation to the solution on the initial mesh must be supplied via the procedure argument
GUESS of
D02TKF.
The solver attempts to satisfy the conditions
where
${v}_{i}$ is the approximate solution for the
$i$th solution component and
TOLS is supplied by you. The mesh is refined by trying to equidistribute the estimated error in the computed solution over all mesh subintervals, and an extrapolationlike test (doubling the number of mesh subintervals) is used to check for
(1).
The routines are based on modified versions of the codes COLSYS and COLNEW (see
Ascher et al. (1979) and
Ascher and Bader (1987)). A comprehensive treatment of the numerical solution of boundary value problems can be found in
Ascher et al. (1988) and
Keller (1992).
4 References
Ascher U M and Bader G (1987) A new basis implementation for a mixed order boundary value ODE solver SIAM J. Sci. Stat. Comput. 8 483–500
Ascher U M, Christiansen J and Russell R D (1979) A collocation solver for mixed order systems of boundary value problems Math. Comput. 33 659–679
Ascher U M, Mattheij R M M and Russell R D (1988) Numerical Solution of Boundary Value Problems for Ordinary Differential Equations Prentice–Hall
Gill P E, Murray W and Wright M H (1981) Practical Optimization Academic Press
Keller H B (1992) Numerical Methods for Twopoint Boundaryvalue Problems Dover, New York
Schwartz I B (1983) Estimating regions of existence of unstable periodic orbits using computerbased techniques SIAM J. Sci. Statist. Comput. 20(1) 106–120
5 Parameters
 1: NEQ – INTEGERInput
On entry: $n$, the number of ordinary differential equations to be solved.
Constraint:
${\mathbf{NEQ}}\ge 1$.
 2: M(NEQ) – INTEGER arrayInput
On entry: ${m}_{\mathit{i}}$, the order of the $\mathit{i}$th differential equation, for $\mathit{i}=1,2,\dots ,n$.
Constraint:
$1\le {\mathbf{M}}\left(\mathit{i}\right)\le 4$, for $\mathit{i}=1,2,\dots ,n$.
 3: NLBC – INTEGERInput
On entry: $p$, the number of left boundary conditions defined at the lefthand end, $a$ ($\text{}={\mathbf{MESH}}\left(1\right)$).
Constraint:
${\mathbf{NLBC}}\ge 1$.
 4: NRBC – INTEGERInput
On entry: $q$, the number of right boundary conditions defined at the righthand end, $b$ ($\text{}={\mathbf{MESH}}\left({\mathbf{NMESH}}\right)$).
Constraints:
 ${\mathbf{NRBC}}\ge 1$;
 ${\mathbf{NLBC}}+{\mathbf{NRBC}}={\displaystyle \sum _{i=1}^{n}}{\mathbf{M}}\left(i\right)$.
 5: NCOL – INTEGERInput
On entry: the number of collocation points to be used in each mesh subinterval.
Constraint:
${m}_{\mathrm{max}}\le {\mathbf{NCOL}}\le 7$, where ${m}_{\mathrm{max}}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{M}}\left(i\right)\right)$.
 6: TOLS(NEQ) – REAL (KIND=nag_wp) arrayInput
On entry: the error requirement for the $i$th solution component.
Constraint:
$100\times \mathit{machineprecision}{\mathbf{TOLS}}\left(\mathit{i}\right)1.0$, for $\mathit{i}=1,2,\dots ,n$.
 7: MXMESH – INTEGERInput
On entry: the maximum number of mesh points to be used during the solution process.
Constraint:
${\mathbf{MXMESH}}\ge 2\times {\mathbf{NMESH}}1$.
 8: NMESH – INTEGERInput
On entry: the number of points to be used in the initial mesh of the solution process.
Constraint:
${\mathbf{NMESH}}\ge 6$.
 9: MESH(MXMESH) – REAL (KIND=nag_wp) arrayInput
On entry: the positions of the initial
NMESH mesh points. The remaining elements of
MESH need not be set. You should try to place the mesh points in areas where you expect the solution to vary most rapidly. In the absence of any other information the points should be equally distributed on
$\left[a,b\right]$.
${\mathbf{MESH}}\left(1\right)$ must contain the left boundary point, $a$, and ${\mathbf{MESH}}\left({\mathbf{NMESH}}\right)$ must contain the right boundary point, $b$.
Constraint:
${\mathbf{MESH}}\left(\mathit{i}\right)<{\mathbf{MESH}}\left(\mathit{i}+1\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{NMESH}}1$.
 10: IPMESH(MXMESH) – INTEGER arrayInput
On entry:
${\mathbf{IPMESH}}\left(\mathit{i}\right)$ specifies whether or not the initial mesh point defined in
${\mathbf{MESH}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,{\mathbf{NMESH}}$, should be a fixed point in all meshes computed during the solution process. The remaining elements of
IPMESH need not be set.
 ${\mathbf{IPMESH}}\left(i\right)=1$
 Indicates that ${\mathbf{MESH}}\left(i\right)$ should be a fixed point in all meshes.
 ${\mathbf{IPMESH}}\left(i\right)=2$
 Indicates that ${\mathbf{MESH}}\left(i\right)$ is not a fixed point.
Constraints:
 ${\mathbf{IPMESH}}\left(1\right)=1$ and ${\mathbf{IPMESH}}\left({\mathbf{NMESH}}\right)=1$, (i.e., the left and right boundary points, $a$ and $b$, must be fixed points, in all meshes);
 ${\mathbf{IPMESH}}\left(\mathit{i}\right)=1$ or $2$, for $\mathit{i}=2,3,\dots ,{\mathbf{NMESH}}1$.
 11: RWORK(LRWORK) – REAL (KIND=nag_wp) arrayCommunication Array
On exit: contains information for use by
D02TKF. This
must be the same array as will be supplied to
D02TKF. The contents of this array
must remain unchanged between calls.
 12: LRWORK – INTEGERInput
On entry: the dimension of the array
RWORK as declared in the (sub)program from which D02TVF is called.
Suggested value:
${\mathbf{LRWORK}}={\mathbf{MXMESH}}\times \left(109\times {{\mathbf{NEQ}}}^{2}+78\times {\mathbf{NEQ}}+7\right)$, which will permit
MXMESH mesh points for a system of
NEQ differential equations regardless of their order or the number of collocation points used.
Constraint:
${\mathbf{LRWORK}}\ge 50+{\mathbf{NEQ}}\times \left({m}_{\mathrm{max}}\times \left(1+{\mathbf{NEQ}}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{NLBC}},{\mathbf{NRBC}}\right)\right)+6\right)\phantom{\rule{0ex}{0ex}}{k}_{n}\times \left({k}_{n}+\phantom{\rule{0ex}{0ex}}6\right){m}^{*}\times \left({k}_{n}+{m}^{*}2\right)+{\mathbf{MXMESH}}\times \left(\left({m}^{*}+3\right)\left(2{m}^{*}+3\right)3+{k}_{n}\left({k}_{n}+{m}^{*}+6\right)\right)\phantom{\rule{0ex}{0ex}}+{\mathbf{MXMESH}}/2$, where ${m}^{*}={\displaystyle \sum _{i=1}^{n}}{\mathbf{M}}\left(i\right)$ and ${k}_{n}={\mathbf{NCOL}}\times {\mathbf{NEQ}}$.
 13: IWORK(LIWORK) – INTEGER arrayCommunication Array
On exit: contains information for use by
D02TKF. This
must be the same array as will be supplied to
D02TKF. The contents of this array
must remain unchanged between calls.
 14: LIWORK – INTEGERInput
On entry: the dimension of the array
IWORK as declared in the (sub)program from which D02TVF is called.
Suggested value:
${\mathbf{LIWORK}}={\mathbf{MXMESH}}\times \left(11\times {\mathbf{NEQ}}+6\right)$, which will permit
MXMESH mesh points for a system of
NEQ differential equations regardless of their order or the number of collocation points used.
Constraint:
${\mathbf{LIWORK}}\ge 23+3\times {\mathbf{NEQ}}{k}_{n}+{\mathbf{MXMESH}}\times \left({m}^{*}+{k}_{n}+4\right)$.
 15: IFAIL – INTEGERInput/Output

On entry:
IFAIL must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of IFAIL on exit.
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
If on entry
${\mathbf{IFAIL}}={\mathbf{0}}$ or
${{\mathbf{1}}}$, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
 ${\mathbf{IFAIL}}=1$
On entry,  ${\mathbf{NEQ}}<1$, 
or  ${\mathbf{M}}\left(i\right)<1$, for some $i$, 
or  ${\mathbf{M}}\left(i\right)>4$, for some $i$, 
or  ${\mathbf{NMESH}}<6$, 
or  the values of MESH are not strictly increasing, 
or  ${\mathbf{IPMESH}}\left(i\right)$ is invalid for some $i$, 
or  ${\mathbf{MXMESH}}<2\times {\mathbf{NMESH}}1$, 
or  ${\mathbf{NCOL}}<{m}_{\mathrm{max}}$, where ${m}_{\mathrm{max}}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{M}}\left(i\right)\right)$, 
or  ${\mathbf{NCOL}}>7$, 
or  ${\mathbf{NLBC}}<1$, 
or  ${\mathbf{NRBC}}<1$, 
or  a value of TOLS is invalid, 
or  ${\mathbf{NLBC}}+{\mathbf{NRBC}}\ne {\displaystyle \sum _{i=1}^{n}}{\mathbf{M}}\left(i\right)$, 
or  LRWORK or LIWORK is too small. 
7 Accuracy
Not applicable.
For problems where sharp changes of behaviour are expected over short intervals it may be advisable to:
– 
use a large value for NCOL; 
– 
cluster the initial mesh points where sharp changes in behaviour are expected; 
– 
maintain fixed points in the mesh using the parameter IPMESH to ensure that the remeshing process does not inadvertently remove mesh points from areas of known interest before they are detected automatically by the algorithm. 
A boundary value problem with nonseparated boundary conditions can be treated by transformation to an equivalent problem with separated conditions. As a simple example consider the system
on
$\left[a,b\right]$ subject to the boundary conditions
By adjoining the trivial ordinary differential equation
which implies
$r\left(a\right)=r\left(b\right)$, and letting
$r\left(b\right)={y}_{2}\left(b\right)$, say, we have a new system
subject to the separated boundary conditions
There is an obvious overhead in adjoining an extra differential equation: the system to be solved is increased in size.
Multipoint boundary value problems, that is problems where conditions are specified at more than two points, can also be transformed to an equivalent problem with two boundary points. Each subinterval defined by the multipoint conditions can be transformed onto the interval
$\left[0,1\right]$, say, leading to a larger set of differential equations. The boundary conditions of the transformed system consist of the original boundary conditions and the conditions imposed by the requirement that the solution components be continuous at the interior break points. For example, consider the equation
subject to the conditions
where
$a<b<c$. This can be transformed to the system
where
subject to the boundary conditions
In this instance two of the resulting boundary conditions are nonseparated but they may next be treated as described above.
Systems of ordinary differential equations containing derivatives of order greater than four can always be reduced to systems of order suitable for treatment by D02TVF and its related routines. For example suppose we have the sixthorder equation
Writing the variables
${y}_{1}=y$ and
${y}_{2}={y}^{\left(4\right)}$ we obtain the system
which has maximal order four, or writing the variables
${y}_{1}=y$ and
${y}_{2}={y}^{\left(3\right)}$ we obtain the system
which has maximal order three. The best choice of reduction by choosing new variables will depend on the structure and physical meaning of the system. Note that you will control the error in each of the variables
${y}_{1}$ and
${y}_{2}$. Indeed, if you wish to control the error in certain derivatives of the solution of an equation of order greater than one, then you should make those derivatives new variables.
The solver routine
D02TKF employs collocation at Gaussian points in each subinterval of the mesh. Hence the coefficients of the differential equations are not evaluated at the mesh points. Thus, fixed points should be specified in the mesh where either the coefficients are singular, or the solution has less smoothness, or where the differential equations should not be evaluated. Singular coefficients at boundary points often arise when physical symmetry is used to reduce partial differential equations to ordinary differential equations. These do not pose a direct numerical problem for using this code but they can severely impact its convergence.
The solver routine
D02TKF requires an external procedure
FJAC to evaluate the partial derivatives of
${f}_{i}$ with respect to the elements of
$z\left(y\right)$ (
$\text{}=\left({y}_{1},{y}_{1}^{1},\dots ,{y}_{1}^{\left({m}_{1}1\right)},{y}_{2},\dots ,{y}_{n}^{\left({m}_{n}1\right)}\right)$). In cases where the partial derivatives are difficult to evaluate, numerical approximations can be used. However, this approach might have a negative impact on the convergence of the modified Newton method. You could consider the use of symbolic mathematic packages and/or automatic differentiation packages if available to you.
See
Section 9 in D02TZF for an example using numerical approximations to the Jacobian. There central differences are used and each
${f}_{i}$ is assumed to depend on all the components of
$z$. This requires two evaluations of the system of differential equations for each component of
$z$. The perturbation used depends on the size of each component of
$z$ and a minimum quantity dependent on the
machine precision. The cost of this approach could be reduced by employing an alternative difference scheme and/or by only perturbing the components of
$z$ which appear in the definitions of the
${f}_{i}$. A discussion on the choice of perturbation factors for use in finite difference approximations to partial derivatives can be found in
Gill et al. (1981).
9 Example
The following example is used to illustrate the treatment of nonseparated boundary conditions. See also
D02TKF,
D02TXF,
D02TYF and
D02TZF, for the illustration of other facilities.
The following equations model of the spread of measles. See
Schwartz (1983). Under certain assumptions the dynamics of the model can be expressed as
subject to the periodic boundary conditions
Here
${y}_{1},{y}_{2}$ and
${y}_{3}$ are respectively the proportions of susceptibles, infectives and latents to the whole population.
$\lambda $ (
$\text{}=0.0279$ years) is the latent period,
$\eta $ (
$\text{}=0.01$ years) is the infectious period and
$\mu $ (
$\text{}=0.02$) is the population birth rate.
$\beta \left(x\right)={\beta}_{0}\left(1.0+\mathrm{cos}2\pi x\right)$ is the contact rate where
${\beta}_{0}=1575.0$.
The nonseparated boundary conditions are treated as described in
Section 8 by adjoining the trivial differential equations
that is
${y}_{4},{y}_{5}$ and
${y}_{6}$ are constants. The boundary conditions of the augmented system can then be posed in the separated form
This is a relatively easy problem and an (arbitrary) initial guess of
$1$ for each component suffices, even though two components of the solution are much smaller than
$1$.
9.1 Program Text
Program Text (d02tvfe.f90)
9.2 Program Data
Program Data (d02tvfe.d)
9.3 Program Results
Program Results (d02tvfe.r)