void	d02tvc (Integer neq, const Integer m[], Integer nlbc, Integer nrbc, Integer ncol, const double tols[], Integer mxmesh, Integer nmesh, const double mesh[], const Integer ipmesh[], double rcomm[], Integer lrcomm, Integer icomm[], Integer licomm, NagError *fail)

The function may be called by the names: d02tvc or nag_ode_bvp_coll_nlin_setup.

3 Description

d02tvc and its associated functions (d02tlc, d02txc, d02tyc and d02tzc) solve the two-point boundary value problem for a nonlinear system of ordinary differential equations

\begin{array}{rcl} y_{1}^{(m_{1})} (x) & = & f_{1} (x, y_{1}, y_{1}^{(1)}, \dots, y_{1}^{(m_{1} - 1)}, y_{2}, \dots, y_{n}^{(m_{n} - 1)}) \\ y_{2}^{(m_{2})} (x) & = & f_{2} (x, y_{1}, y_{1}^{(1)}, \dots, y_{1}^{(m_{1} - 1)}, y_{2}, \dots, y_{n}^{(m_{n} - 1)}) \\ ⋮ \\ y_{n}^{(m_{n})} (x) & = & f_{n} (x, y_{1}, y_{1}^{(1)}, \dots, y_{1}^{(m_{1} - 1)}, y_{2}, \dots, y_{n}^{(m_{n} - 1)}) \end{array}

over an interval

[a, b]

subject to

p

(

> 0

) nonlinear boundary conditions at

a

and

q

(

> 0

) nonlinear boundary conditions at

b

, where

p + q = \sum_{i = 1}^{n} m_{i}

. Note that

y_{i}^{(m)} (x)

is the

m

th derivative of the

i

th solution component. Hence

y_{i}^{(0)} (x) = y_{i} (x)

. The left boundary conditions at

a

are defined as

g_{i} (z (y (a))) = 0, i = 1, 2, \dots, p,

and the right boundary conditions at

b

{\bar{g}}_{j} (z (y (b))) = 0, j = 1, 2, \dots, q,

where

y = (y_{1}, y_{2}, \dots, y_{n})

and

z (y (x)) = (y_{1} (x), y_{1}^{(1)} (x), \dots, y_{1}^{(m_{1} - 1)} (x), y_{2} (x), \dots, y_{n}^{(m_{n} - 1)} (x)) .

See Section 9 for information on how boundary value problems of a more general nature can be treated.

d02tvc is used to specify an initial mesh, error requirements and other details. d02tlc is then used to solve the boundary value problem.

The solution function d02tlc proceeds as follows. A modified Newton method is applied to the equations

y_{i}^{(m_{i})} (x) - f_{i} (x, z (y (x))) = 0, i = 1, \dots, n

and the boundary conditions. To solve these equations numerically the components

y_{i}

are approximated by piecewise polynomials

v_{i j}

using a monomial basis on the

j

th mesh sub-interval. The coefficients of the polynomials

v_{i j}

form the unknowns to be computed. Collocation is applied at Gaussian points

v_{i j}^{(m_{i})} (x_{j k}) - f_{i} (x_{j k}, z (v (x_{j k}))) = 0, i = 1, 2, \dots, n,

where

x_{j k}

is the

k

th collocation point in the

j

th mesh sub-interval. Continuity at the mesh points is imposed, that is

v_{i j} (x_{j + 1}) - v_{i, j + 1} (x_{j + 1}) = 0, i = 1, 2, \dots, n,

where

x_{j + 1}

is the right-hand end of the

j

th mesh sub-interval. The linearized collocation equations and boundary conditions, together with the continuity conditions, form a system of linear algebraic equations, an almost block diagonal system which is solved using special linear solvers. To start the modified Newton process, an approximation to the solution on the initial mesh must be supplied via the procedure argument guess of d02tlc.

The solver attempts to satisfy the conditions

\frac{‖ y_{i} - v_{i} ‖}{(1.0 + ‖ v_{i} ‖)} \leq tols [i - 1], i = 1, 2, \dots, n,

(1)

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 sub-intervals, and an extrapolation-like test (doubling the number of mesh sub-intervals) is used to check for (1).

The functions 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 Two-point Boundary-value Problems Dover, New York

Schwartz I B (1983) Estimating regions of existence of unstable periodic orbits using computer-based techniques SIAM J. Sci. Statist. Comput. 20(1) 106–120

5 Arguments

1: $neq$ – Integer Input

On entry:

n

, the number of ordinary differential equations to be solved.

Constraint:

neq \geq 1

2: $m [neq]$ – const Integer Input

On entry:

m [i - 1]

must contain

m_{i}

, the order of the

i

th differential equation, for

i = 1, 2, \dots, n

Constraint:

1 \leq m [i - 1] \leq 4

, for

i = 1, 2, \dots, n

3: $nlbc$ – Integer Input

On entry:

p

, the number of left boundary conditions defined at the left-hand end,

a

(

= mesh [0]

Constraint:

nlbc \geq 1

4: $nrbc$ – Integer Input

On entry:

q

, the number of right boundary conditions defined at the right-hand end,

b

(

= mesh [nmesh - 1]

Constraints:

$nrbc \geq 1$ ;
$nlbc + nrbc = \sum_{i = 1}^{n} m [i - 1]$ .

5: $ncol$ – Integer Input

On entry: the number of collocation points to be used in each mesh sub-interval.

Constraint:

m_{\max} \leq ncol \leq 7

, where

m_{\max} = \max (m [i - 1])

6: $tols [neq]$ – const double Input

On entry:

tols [i - 1]

must contain the error requirement for the

i

th solution component.

Constraint:

100 \times machine precision < tols [i - 1] < 1.0

, for

i = 1, 2, \dots, n

7: $mxmesh$ – Integer Input

On entry: the maximum number of mesh points to be used during the solution process.

Constraint:

mxmesh \geq 2 \times nmesh - 1

8: $nmesh$ – Integer Input

On entry: the number of points to be used in the initial mesh of the solution process.

Constraint:

nmesh \geq 6

9: $mesh [mxmesh]$ – const double Input

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

[a, b]

mesh [0]

must contain the left boundary point,

a

, and

mesh [nmesh - 1]

must contain the right boundary point,

b

Constraint:

mesh [i - 1] < mesh [i]

, for

i = 1, 2, \dots, nmesh - 1

10: $ipmesh [mxmesh]$ – const Integer Input

On entry:

ipmesh [i - 1]

specifies whether or not the initial mesh point defined in

mesh [i - 1]

, for

i = 1, 2, \dots, nmesh

, should be a fixed point in all meshes computed during the solution process. The remaining elements of ipmesh need not be set.

$ipmesh [i - 1] = 1$: Indicates that $mesh [i - 1]$ should be a fixed point in all meshes.
$ipmesh [i - 1] = 2$: Indicates that $mesh [i - 1]$ is not a fixed point.

Constraints:

$ipmesh [0] = 1$ and $ipmesh [nmesh - 1] = 1$ , (i.e., the left and right boundary points, $a$ and $b$ , must be fixed points, in all meshes);
$ipmesh [i - 1] = 1$ or $2$ , for $i = 2, 3, \dots, nmesh - 1$ .

11: $rcomm [lrcomm]$ – double Communication Array

On exit: contains information for use by d02tlc. This must be the same array as will be supplied to d02tlc. The contents of this array must remain unchanged between calls.

12: $lrcomm$ – Integer Input

On entry: the dimension of the array rcomm. If

lrcomm = 0

, a communication array size query is requested. In this case there is an immediate return with communication array dimensions stored in icomm;

icomm [0]

contains the required dimension of rcomm, while

icomm [1]

contains the required dimension of icomm.

Constraint:

lrcomm = 0

, or

lrcomm \geq 51 + neq + mxmesh \times (2 + m^{*} + k_{n}) - k_{n} + mxmesh / 2

, where

m^{*} = \sum_{i = 1}^{n} m [i - 1]

and

k_{n} = ncol \times neq

13: $icomm [licomm]$ – Integer Communication Array

On exit: contains information for use by d02tlc. This must be the same array as will be supplied to d02tlc. The contents of this array must remain unchanged between calls. If

lrcomm = 0

, a communication array size query is requested. In this case, on immediate return,

icomm [0]

will contain the required dimension for rcomm while

icomm [1]

will contain the required dimension for icomm.

14: $licomm$ – Integer Input

On entry: the dimension of the array icomm. If

lrcomm = 0

, a communication array size query is requested. In this case icomm need only be of dimension

2

in order to hold the required communication array dimensions for the given problem and algorithmic parameters.

Constraints:

if $lrcomm = 0$ , $licomm \geq 2$ ;
otherwise $licomm \geq 23 + neq + mxmesh$ .

15: $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_ALLOC_FAIL: Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_BAD_PARAM: On entry, argument $⟨ value ⟩$ had an illegal value.
NE_INT: On entry, $ipmesh [0]$ or $ipmesh [nmesh - 1]$ does not equal $1$ .

On entry, $ipmesh [i] \neq 1$ or $2$ for some $i$ .

On entry, $licomm = ⟨ value ⟩$ .
Constraint: $licomm \geq ⟨ value ⟩$ .

On entry, $lrcomm = ⟨ value ⟩$ .
Constraint: $lrcomm = 0$ or $lrcomm \geq ⟨ value ⟩$ .

On entry, $m [⟨ value ⟩] = ⟨ value ⟩$ .
Constraint: $1 \leq m [i - 1] \leq 4$ for all $i$ .

On entry, $neq = ⟨ value ⟩$ .
Constraint: $neq \geq 1$ .

On entry, $nmesh = ⟨ value ⟩$ .
Constraint: $nmesh \geq 6$ .
NE_INT_2: On entry, $mxmesh = ⟨ value ⟩$ and $nmesh = ⟨ value ⟩$ .
Constraint: $mxmesh \geq 2 \times nmesh - 1$ .

On entry, $ncol = ⟨ value ⟩$ and $\max (m [i]) = ⟨ value ⟩$ .
Constraint: $\max (m [i]) \leq ncol \leq 7$ .

On entry, $nlbc = ⟨ value ⟩$ and $nrbc = ⟨ value ⟩$ .
Constraint: $nlbc \geq 1$ and $nrbc \geq 1$ .
NE_INT_3: On entry, $nlbc = ⟨ value ⟩$ , $nrbc = ⟨ value ⟩$ and $sum (m [i]) = ⟨ value ⟩$ .
Constraint: $nlbc + nrbc = sum (m [i])$ .
NE_INTERNAL_ERROR: 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.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_NO_LICENCE: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
NE_NOT_STRICTLY_INCREASING: On entry, the elements of mesh are not strictly increasing.
NE_REAL: On entry, $tols [⟨ value ⟩] = ⟨ value ⟩$ .
Constraint: $tols [i - 1] > ⟨ value ⟩$ for all $i$ .

7 Accuracy

Not applicable.

8 Parallelism and Performance

d02tvc is not threaded in any implementation.

9 Further Comments

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 argument 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.

9.1 Nonseparated Boundary Conditions

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

\begin{matrix} y_{1}^{'} = f_{1} (x, y_{1}, y_{2}) \\ y_{2}^{'} = f_{2} (x, y_{1}, y_{2}) \end{matrix}

[a, b]

subject to the boundary conditions

\begin{array}{r} g_{1} (y_{1} (a)) = 0 \\ g_{2} (y_{2} (a), y_{2} (b)) = 0 . \end{array}

By adjoining the trivial ordinary differential equation

r^{'} = 0,

which implies

r (a) = r (b)

, and letting

r (b) = y_{2} (b)

, say, we have a new system

\begin{array}{rcl} y_{1}^{'} & = & f_{1} (x, y_{1}, y_{2}) \\ y_{2}^{'} & = & f_{2} (x, y_{1}, y_{2}) \\ r^{'} & = & 0, \end{array}

subject to the separated boundary conditions

\begin{array}{r} g_{1} (y_{1} (a)) = 0 \\ g_{2} (y_{2} (a), r (a)) = 0 \\ y_{2} (b) - r (b) = 0 . \end{array}

There is an obvious overhead in adjoining an extra differential equation: the system to be solved is increased in size.

9.2 Multipoint Boundary Value Problems

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 sub-interval defined by the multipoint conditions can be transformed onto the interval

[0, 1]

, 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

y^{(3)} = f (t, y, y^{(1)}, y^{(2)}) on [a, c]

subject to the conditions

\begin{array}{r} y (a) = A \\ y (b) = B \\ y^{(1)} (c) = C \end{array}

where

a < b < c

. This can be transformed to the system

\begin{array}{rcl} y_{1}^{(3)} & = & f (t, y_{1}, y_{1}^{(1)}, y_{1}^{(2)}) \\ y_{2}^{(3)} & = & f (t, y_{2}, y_{2}^{(1)}, y_{2}^{(2)}) \end{array}} on ​ [0, 1]

where

\begin{array}{rcl} y_{1} & \equiv & y on [a, b] \\ y_{2} & \equiv & y on [b, c], \end{array}

subject to the boundary conditions

\begin{array}{l} y_{1} (0) & = & A \\ y_{1} (1) & = & B \\ y_{2}^{(1)} (1) & = & C \\ y_{2} (0) & = & B (from ​ y_{1} (1) = y_{2} (0)) \\ y_{1}^{(1)} (1) & = & y_{2}^{(1)} (0) \\ y_{1}^{(2)} (1) & = & y_{2}^{(2)} (0) . \end{array}

In this instance two of the resulting boundary conditions are nonseparated but they may next be treated as described above.

9.3 High Order Systems

Systems of ordinary differential equations containing derivatives of order greater than four can always be reduced to systems of order suitable for treatment by d02tvc and its related functions. For example suppose we have the sixth-order equation

y^{(6)} = - y .

Writing the variables

y_{1} = y

and

y_{2} = y^{(4)}

we obtain the system

\begin{array}{rcl} y_{1}^{(4)} & = & y_{2} \\ y_{2}^{(2)} & = & - y_{1} \end{array}

which has maximal order four, or writing the variables

y_{1} = y

and

y_{2} = y^{(3)}

we obtain the system

\begin{array}{rcl} y_{1}^{(3)} & = & y_{2} \\ y_{2}^{(3)} & = & - y_{1} \end{array}

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

1

, then you should make those derivatives new variables.

9.4 Fixed Points and Singularities

The solver function d02tlc employs collocation at Gaussian points in each sub-interval 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.

9.5 Numerical Jacobians

The solver function d02tlc requires an external function fjac to evaluate the partial derivatives of

f_{i}

with respect to the elements of

z (y)

(

= (y_{1}, y_{1}^{1}, \dots, y_{1}^{(m_{1} - 1)}, y_{2}, \dots, y_{n}^{(m_{n} - 1)})

). 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 10 in d02tzc 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).

10 Example

The following example is used to illustrate the treatment of nonseparated boundary conditions. See also d02tlc, d02txc, d02tyc and d02tzc, 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

\begin{array}{l} y_{1}^{'} = μ - β (x) y_{1} y_{3} \\ y_{2}^{'} = β (x) y_{1} y_{3} - y_{2} / λ \\ y_{3}^{'} = y_{2} / λ - y_{3} / η \end{array}

subject to the periodic boundary conditions

y_{i} (0) = y_{i} (1), i = 1, 2, 3 .

Here

y_{1}, y_{2}

and

y_{3}

are respectively the proportions of susceptibles, infectives and latents to the whole population.

λ

(

= 0.0279

years) is the latent period,

η

(

= 0.01

years) is the infectious period and

μ

(

= 0.02

) is the population birth rate.

β (x) = β_{0} (1.0 + \cos 2 π x)

is the contact rate where

β_{0} = 1575.0

The nonseparated boundary conditions are treated as described in Section 9 by adjoining the trivial differential equations

\begin{array}{r} y_{4}^{'} = 0 \\ y_{5}^{'} = 0 \\ y_{6}^{'} = 0 \end{array}

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

\begin{array}{r} y_{1} (0) - y_{4} (0) = 0 \\ y_{2} (0) - y_{5} (0) = 0 \\ y_{3} (0) - y_{6} (0) = 0 \\ y_{1} (1) - y_{4} (1) = 0 \\ y_{2} (1) - y_{5} (1) = 0 \\ y_{3} (1) - y_{6} (1) = 0 . \end{array}

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

d02tv: FL CL CPP AD

NAG CL Interfaced02tvc (bvp_​coll_​nlin_​setup)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Nonseparated Boundary Conditions

9.2 Multipoint Boundary Value Problems

9.3 High Order Systems

9.4 Fixed Points and Singularities

9.5 Numerical Jacobians

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG CL Interface
d02tvc (bvp_coll_nlin_setup)