hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_ode_bvp_coll_nlin_setup (d02tv)

Purpose

nag_ode_bvp_coll_nlin_setup (d02tv) is a setup function which must be called prior to the first call of the nonlinear two-point boundary value solver nag_ode_bvp_coll_nlin (d02tk).

Syntax

[rwork, iwork, ifail] = d02tv(m, nlbc, nrbc, ncol, tols, nmesh, mesh, ipmesh, 'neq', neq, 'mxmesh', mxmesh, 'lrwork', lrwork, 'liwork', liwork)
[rwork, iwork, ifail] = nag_ode_bvp_coll_nlin_setup(m, nlbc, nrbc, ncol, tols, nmesh, mesh, ipmesh, 'neq', neq, 'mxmesh', mxmesh, 'lrwork', lrwork, 'liwork', liwork)

Description

nag_ode_bvp_coll_nlin_setup (d02tv) and its associated functions (nag_ode_bvp_coll_nlin (d02tk), nag_ode_bvp_coll_nlin_contin (d02tx), nag_ode_bvp_coll_nlin_interp (d02ty) and nag_ode_bvp_coll_nlin_diag (d02tz)) solve the two-point boundary value problem for a nonlinear system of ordinary differential equations
y1(m1) (x) = f1 (x,y1,y1(1),,y1(m11),y2,,yn(mn1))
y2(m2) (x) = f2 (x,y1,y1(1),,y1(m11),y2,,yn(mn1))
yn(mn) (x) = fn (x,y1,y1(1),,y1(m11),y2,,yn(mn1))
y1(m1) (x) = f1 (x,y1,y1(1),,y1(m1-1),y2,,yn(mn-1)) y2(m2) (x) = f2 (x,y1,y1(1),,y1(m1-1),y2,,yn(mn-1)) yn(mn) (x) = fn (x,y1,y1(1),,y1(m1-1),y2,,yn(mn-1))
over an interval [a,b][a,b] subject to pp ( > 0>0) nonlinear boundary conditions at aa and qq ( > 0>0) nonlinear boundary conditions at bb, where p + q = i = 1n mi p+q = i=1 n mi . Note that yi(m)(x)yi (m) (x) is the mmth derivative of the iith solution component. Hence yi(0)(x) = yi(x)yi (0) (x)=yi(x). The left boundary conditions at aa are defined as
gi(z(y(a))) = 0,  i = 1,2,,p,
gi(z(y(a)))=0,  i=1,2,,p,
and the right boundary conditions at bb as
gj(z(y(b))) = 0,  j = 1,2,,q,
g-j(z(y(b)))=0,  j=1,2,,q,
where y = (y1,y2,,yn)y=(y1,y2,,yn) and
z(y(x)) = (y1(x), y1(1) (x) ,, y1(m11) (x) ,y2(x),, yn(mn1) (x) ) .
z(y(x)) = (y1(x), y1(1) (x) ,, y1(m1-1) (x) ,y2(x),, yn(mn-1) (x) ) .
See Section [Further Comments] for information on how boundary value problems of a more general nature can be treated.
nag_ode_bvp_coll_nlin_setup (d02tv) is used to specify an initial mesh, error requirements and other details. nag_ode_bvp_coll_nlin (d02tk) is then used to solve the boundary value problem.
The solution function nag_ode_bvp_coll_nlin (d02tk) proceeds as follows. A modified Newton method is applied to the equations
yi(mi)(x)fi(x,z(y(x))) = 0,   i = 1,,n
yi (mi) (x)-fi(x,z(y(x)))= 0,   i= 1,,n
and the boundary conditions. To solve these equations numerically the components yiyi are approximated by piecewise polynomials vijvij using a monomial basis on the jjth mesh sub-interval. The coefficients of the polynomials vijvij form the unknowns to be computed. Collocation is applied at Gaussian points
vij(mi)(xjk)fi(xjk,z(v(xjk))) = 0,  i = 1,2,,n,
vij (mi) (xjk)-fi(xjk,z(v(xjk)))=0,  i=1,2,,n,
where xjkxjk is the kkth collocation point in the jjth mesh sub-interval. Continuity at the mesh points is imposed, that is
vij(xj + 1)vi,j + 1(xj + 1) = 0,  i = 1,2,,n,
vij(xj+1)-vi,j+1(xj+1)=0,  i=1,2,,n,
where xj + 1xj+1 is the right-hand end of the jjth mesh sub-interval. The linearized collocation equations and boundary conditions, together with the continuity conditions, form a system of linear algebraic equations which are solved using nag_matop_real_gen_blkdiag_lu (f01lh) and nag_linsys_real_blkdiag_fac_solve (f04lh). 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 nag_ode_bvp_coll_nlin (d02tk).
The solver attempts to satisfy the conditions
(yivi)/((1.0 + vi))tols(i),  i = 1,2,,n,
yi-vi (1.0+vi) tolsi,  i=1,2,,n,
(1)
where vivi is the approximate solution for the iith 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).

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

Parameters

Compulsory Input Parameters

1:     m(neq) – int64int32nag_int array
neq, the dimension of the array, must satisfy the constraint neq1neq1.
mimi, the order of the iith differential equation, for i = 1,2,,ni=1,2,,n.
Constraint: 1m(i)41mi4, for i = 1,2,,ni=1,2,,n.
2:     nlbc – int64int32nag_int scalar
pp, the number of left boundary conditions defined at the left-hand end, aa ( = mesh(1)=mesh1).
Constraint: nlbc1nlbc1.
3:     nrbc – int64int32nag_int scalar
qq, the number of right boundary conditions defined at the right-hand end, bb ( = mesh(nmesh)=meshnmesh).
Constraints:
  • nrbc1nrbc1;
  • nlbc + nrbc = i = 1nm(i)nlbc+nrbc=i=1nmi.
4:     ncol – int64int32nag_int scalar
The number of collocation points to be used in each mesh sub-interval.
Constraint: mmaxncol7mmaxncol7, where mmax = max (m(i))mmax=max(mi).
5:     tols(neq) – double array
neq, the dimension of the array, must satisfy the constraint neq1neq1.
The error requirement for the iith solution component.
Constraint: 100 × machine precision < tols(i) < 1.0100×machine precision<tolsi<1.0, for i = 1,2,,ni=1,2,,n.
6:     nmesh – int64int32nag_int scalar
The number of points to be used in the initial mesh of the solution process.
Constraint: nmesh6nmesh6.
7:     mesh(mxmesh) – double array
mxmesh, the dimension of the array, must satisfy the constraint mxmesh2 × nmesh1mxmesh2×nmesh-1.
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][a,b].
mesh(1)mesh1 must contain the left boundary point, aa, and mesh(nmesh)meshnmesh must contain the right boundary point, bb.
Constraint: mesh(i) < mesh(i + 1)meshi<meshi+1, for i = 1,2,,nmesh1i=1,2,,nmesh-1.
8:     ipmesh(mxmesh) – int64int32nag_int array
mxmesh, the dimension of the array, must satisfy the constraint mxmesh2 × nmesh1mxmesh2×nmesh-1.
ipmesh(i)ipmeshi specifies whether or not the initial mesh point defined in mesh(i)meshi, for i = 1,2,,nmeshi=1,2,,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) = 1ipmeshi=1
Indicates that mesh(i)meshi should be a fixed point in all meshes.
ipmesh(i) = 2ipmeshi=2
Indicates that mesh(i)meshi is not a fixed point.
Constraints:
  • ipmesh(1) = 1ipmesh1=1 and ipmesh(nmesh) = 1ipmeshnmesh=1, (i.e., the left and right boundary points, aa and bb, must be fixed points, in all meshes);
  • ipmesh(i) = 1ipmeshi=1 or 22, for i = 2,3,,nmesh1i=2,3,,nmesh-1.

Optional Input Parameters

1:     neq – int64int32nag_int scalar
Default: The dimension of the arrays m, tols. (An error is raised if these dimensions are not equal.)
nn, the number of ordinary differential equations to be solved.
Constraint: neq1neq1.
2:     mxmesh – int64int32nag_int scalar
Default: The dimension of the arrays mesh, ipmesh. (An error is raised if these dimensions are not equal.)
The maximum number of mesh points to be used during the solution process.
Constraint: mxmesh2 × nmesh1mxmesh2×nmesh-1.
3:     lrwork – int64int32nag_int scalar
The dimension of the array rwork as declared in the (sub)program from which nag_ode_bvp_coll_nlin_setup (d02tv) is called.
Suggested value: lrwork = mxmesh × (109 × neq2 + 78 × neq + 7)lrwork=mxmesh×(109×neq2+78×neq+7), which will permit mxmesh mesh points for a system of neq differential equations regardless of their order or the number of collocation points used.
Default: mxmesh × (109 × neq2 + 78 × neq + 7)mxmesh×(109×neq2+78×neq+7)
Constraint: lrwork 50 + neq × (mmax × (1 + neq + max (nlbc,nrbc)) + 6)   kn × (kn + 6) m* × (kn + m*2) + mxmesh × ((m* + 3)(2m* + 3)3 + kn(kn + m* + 6))   + mxmesh / 2 lrwork 50 + neq × ( mmax × ( 1 + neq + max(nlbc,nrbc) ) + 6 ) - k n × ( kn + 6 ) - m* × ( kn + m* - 2 ) + mxmesh × ( ( m* + 3 ) ( 2m* + 3 ) - 3 + kn ( kn + m* + 6 ) ) + mxmesh/2 , where m* = i = 1nm(i)m*=i=1nmi and kn = ncol × neqkn=ncol×neq.
4:     liwork – int64int32nag_int scalar
The dimension of the array iwork as declared in the (sub)program from which nag_ode_bvp_coll_nlin_setup (d02tv) is called.
Suggested value: liwork = mxmesh × (11 × neq + 6)liwork=mxmesh×(11×neq+6), which will permit mxmesh mesh points for a system of neq differential equations regardless of their order or the number of collocation points used.
Default: mxmesh × (11 × neq + 6)mxmesh×(11×neq+6)
Constraint: liwork23 + 3 × neqkn + mxmesh × (m* + kn + 4)liwork23+3×neq-kn+mxmesh×(m*+kn+4).

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     rwork(lrwork) – double array
Contains information for use by nag_ode_bvp_coll_nlin (d02tk). This must be the same array as will be supplied to nag_ode_bvp_coll_nlin (d02tk). The contents of this array must remain unchanged between calls.
2:     iwork(liwork) – int64int32nag_int array
Contains information for use by nag_ode_bvp_coll_nlin (d02tk). This must be the same array as will be supplied to nag_ode_bvp_coll_nlin (d02tk). The contents of this array must remain unchanged between calls.
3:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,neq < 1neq<1,
orm(i) < 1mi<1, for some ii,
orm(i) > 4mi>4, for some ii,
ornmesh < 6nmesh<6,
orthe values of mesh are not strictly increasing,
oripmesh(i)ipmeshi is invalid for some ii,
ormxmesh < 2 × nmesh1mxmesh<2×nmesh-1,
orncol < mmaxncol<mmax, where mmax = max (m(i))mmax=max(mi),
orncol > 7ncol>7,
ornlbc < 1nlbc<1,
ornrbc < 1nrbc<1,
ora value of tols is invalid,
or nlbc + nrbc i = 1n m(i) nlbc+nrbc i=1 n mi ,
orlrwork or liwork is too small.

Accuracy

Not applicable.

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

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
y1 = f1(x,y1,y2)
y2 = f2(x,y1,y2)
y1=f1(x,y1,y2) y2=f2(x,y1,y2)
on [a,b][a,b] subject to the boundary conditions
g1(y1(a)) = 0
g2(y2(a),y2(b)) = 0.
g1(y1(a))=0 g2(y2(a),y2(b))=0.
By adjoining the trivial ordinary differential equation
r = 0,
r=0,
which implies r(a) = r(b)r(a)=r(b), and letting r(b) = y2(b)r(b)=y2(b), say, we have a new system
y1 = f1(x,y1,y2)
y2 = f2(x,y1,y2)
r = 0,
y1 = f1(x,y1,y2) y2 = f2(x,y1,y2) r = 0,
subject to the separated boundary conditions
g1(y1(a)) = 0
g2(y2(a),r(a)) = 0
y2(b)r(b) = 0.
g1(y1(a))=0 g2(y2(a),r(a))=0 y2(b)-r(b)=0.
There is an obvious overhead in adjoining an extra differential equation: the system to be solved is increased in size.

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][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]
y (3) =f(t,y,y (1) ,y (2) )  on  [a,c]
subject to the conditions
y(a) = A
y(b) = B
y(1)(c) = C
y(a)=A y(b)=B y(1)(c)=C
where a < b < ca<b<c. This can be transformed to the system
y1(3) = f(t,y1,y1(1),y1(2)) y2(3) = f(t,y2,y2(1),y2(2)) }
  on ​ [0,1]
y1 (3) = f(t,y1,y1 (1) ,y1 (2) ) y2 (3) = f(t,y2,y2 (1) ,y2 (2) ) }   on ​ [0,1]
where
y1 y   on   [a,b]
y2 y   on   [b,c],
y1 y   on   [a,b] y2 y   on   [b,c],
subject to the boundary conditions
y1(0) = A
y1(1) = B
y2(1)(1) = C
y2(0) = B(from ​y1(1) = y2(0))
y1(1)(1) = y2(1)(0)
y1(2)(1) = y2(2)(0).
y1(0) = A y1(1) = B y2 (1) (1) = C y2(0) = B(from ​ y1(1)=y2(0)) y1 (1) (1) = y2 (1) (0) y1 (2) (1) = y2 (2) (0).
In this instance two of the resulting boundary conditions are nonseparated but they may next be treated as described above.

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 nag_ode_bvp_coll_nlin_setup (d02tv) and its related functions. For example suppose we have the sixth-order equation
y(6) = y.
y (6) =-y.
Writing the variables y1 = yy1=y and y2 = y(4)y2=y (4)  we obtain the system
y1(4) = y2
y2(2) = y1
y1 (4) = y2 y2 (2) = -y1
which has maximal order four, or writing the variables y1 = yy1=y and y2 = y(3)y2=y (3)  we obtain the system
y1(3) = y2
y2(3) = y1
y1 (3) = y2 y2 (3) = -y1
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 y1y1 and y2y2. 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.

Fixed Points and Singularities

The solver function nag_ode_bvp_coll_nlin (d02tk) 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.

Numerical Jacobians

The solver function nag_ode_bvp_coll_nlin (d02tk) requires an external procedure fjac to evaluate the partial derivatives of fifi with respect to the elements of z(y)z(y) ( = (y1,y11,,y1(m11),y2,,yn(mn1))=(y1,y11,,y1 (m1-1) ,y2,,yn (mn-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 [Example] in (d02tz) for an example using numerical approximations to the Jacobian. There central differences are used and each fifi is assumed to depend on all the components of zz. This requires two evaluations of the system of differential equations for each component of zz. The perturbation used depends on the size of each component of zz 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 zz which appear in the definitions of the fifi. 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).

Example

The following example is used to illustrate the treatment of nonseparated boundary conditions. See also nag_ode_bvp_coll_nlin (d02tk), nag_ode_bvp_coll_nlin_contin (d02tx), nag_ode_bvp_coll_nlin_interp (d02ty) and nag_ode_bvp_coll_nlin_diag (d02tz), 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
y1 = μβ(x)y1y3
y2 = β(x)y1y3y2 / λ
y3 = y2 / λy3 / η
y1=μ-β(x)y1y3 y2=β(x)y1y3-y2/λ y3=y2/λ-y3/η
subject to the periodic boundary conditions
yi(0) = yi(1) ,   i = 1,2,3.
yi(0)=yi(1) ,   i= 1,2,3.
Here y1,y2y1,y2 and y3y3 are respectively the proportions of susceptibles, infectives and latents to the whole population. λλ ( = 0.0279=0.0279 years) is the latent period, ηη ( = 0.01=0.01 years) is the infectious period and μμ ( = 0.02=0.02) is the population birth rate. β(x) = β0(1.0 + cos2πx)β(x)=β0(1.0+cos2πx) is the contact rate where β0 = 1575.0β0=1575.0.
The nonseparated boundary conditions are treated as described in Section [Further Comments] by adjoining the trivial differential equations
y4 = 0
y5 = 0
y6 = 0
y4=0 y5=0 y6=0
that is y4,y5y4,y5 and y6y6 are constants. The boundary conditions of the augmented system can then be posed in the separated form
y1(0)y4(0) = 0
y2(0)y5(0) = 0
y3(0)y6(0) = 0
y1(1)y4(1) = 0
y2(1)y5(1) = 0
y3(1)y6(1) = 0.
y1(0)-y4(0)=0 y2(0)-y5(0)=0 y3(0)-y6(0)=0 y1(1)-y4(1)=0 y2(1)-y5(1)=0 y3(1)-y6(1)=0.
This is a relatively easy problem and an (arbitrary) initial guess of 11 for each component suffices, even though two components of the solution are much smaller than 11.
function nag_ode_bvp_coll_nlin_setup_example
global beta0 eta lambda mu % For communication with local functions

% Initialize variables and arrays.
neq = 6;
nlbc = 3;
nrbc = 3;
ncol = 5;
mmax = 1;

% Set up the mesh.
nmesh = 11;
mxmesh = 2*nmesh;
ipmesh = zeros(mxmesh, 1);
mesh = zeros(mxmesh, 1);

% Set location of mesh boundaries, then calculate initial spacing.
mesh(1) = 0.0;
mesh(nmesh) = 1.0;
mstep = (mesh(nmesh) - mesh(1))/double(nmesh-1);
for i = 2:nmesh-1
    mesh(i) = mstep*double(i-1);
    ipmesh(i) = 2;
end

% Specify mesh end points as fixed.
ipmesh(1)     = 1;
ipmesh(nmesh) = 1;

m = [1; 1; 1; 1; 1; 1];
tols = [1e-05; 1e-05; 1e-05; 1e-05; 1e-05; 1e-05];

fprintf('nag_ode_bvp_coll_nlin_setup example program results\n\n');

% nag_ode_bvp_coll_nlin_setup is a setup routine to be called prior to nag_ode_bvp_coll_nlin.
[work, iwork, ifail] = nag_ode_bvp_coll_nlin_setup(int64(m), int64(nlbc), int64(nrbc), ...
    int64(ncol), tols, int64(nmesh), mesh, int64(ipmesh));
if ifail ~= 0
    % Unsuccessful call.  Print message and exit.
    error('Warning: nag_ode_bvp_coll_nlin_setup returned with ifail = %1d ',ifail);
end

% Set values for problem-specific physical parameters.
eta = 0.01;
mu = 0.02;
lambda = 0.0279;
beta0 = 1575.0;

% Call nag_ode_bvp_coll_nlin to solve BVP for this set of parameters.
[work, iwork, ifail] = nag_ode_bvp_coll_nlin(@ffun, @fjac, ...
    @gafun, @gbfun, @gajac, @gbjac,...
    @guess, work, iwork);
if ifail ~= 0
    % Unsuccessful call.  Print message and exit.
    error('Warning: nag_ode_bvp_coll_nlin returned with ifail = %1d ',ifail);
end

% Call nag_ode_bvp_coll_nlin_diag to extract mesh from solution.
[nmesh, mesh, ipmesh, ermx, iermx, ijermx, ifail] = nag_ode_bvp_coll_nlin_diag( ...
    int64(mxmesh), work, iwork);

% Output mesh results.
fprintf(' Used a mesh of %d points\n', nmesh);
fprintf([' Maximum error = %10.2e in interval %d for component %d\n\n', ...
    ' Mesh points:\n'], ermx, iermx, ijermx);
for imesh = 1:int32(nmesh) % can't use int64 in loop range.
    fprintf( '%4d(%d) %6.4f', imesh, ipmesh(imesh), mesh(imesh));
    if mod(imesh, 4) == 0
        fprintf('\n');
    end
end

% Output solution, and store it for plotting.
xarray = zeros(nmesh, 1);
yarray = zeros(nmesh, 3);
fprintf('\n\n Computed solution at mesh points\n');
fprintf('    x       y1         y2         y3\n');
for imesh = 1:int32(nmesh) % can't use int64 in loop range.
    % Call nag_ode_bvp_coll_nlin_interp to perform interpolation on the solution.
    [y, work, ifail] = nag_ode_bvp_coll_nlin_interp(mesh(imesh), int64(neq), int64(mmax), ...
        work, iwork);
    fprintf(' %6.3f  %8.2e  %8.2e  %8.2e\n', mesh(imesh), ...
        y(1,1), y(2,1), y(3,1));
    xarray(imesh) = mesh(imesh);
    for jcomp = 1:3
        yarray(imesh, jcomp) = y(jcomp,1);
    end
end
% Plot results.
fig = figure('Number', 'off');
display_plot(xarray, yarray);


function [f] = ffun(x, y, neq, m)
% Evaluate derivative functions (rhs of system of ODEs).

global beta0 eta lambda mu % For communication with main routine.
beta = beta0*(1.0 + cos(2.0*pi*x));
f = zeros(neq, 1);
f(1) = mu - beta*y(1,1)*y(3,1);
f(2) = beta*y(1,1)*y(3,1) - y(2,1)/lambda;
f(3) = y(2,1)/lambda - y(3,1)/eta;
f(4) = 0.0;
f(5) = 0.0;
f(6) = 0.0;
function [dfdy] = fjac(x, y, neq, m)
% Evaluate Jacobians (partial derivatives of f).

global beta0 eta lambda mu % For communication with main routine.
beta = beta0*(1.0 + cos(2.0*pi*x));
dfdy = zeros(neq, neq, 1);
dfdy(1,1,1) = -beta*y(3,1);
dfdy(1,3,1) = -beta*y(1,1);
dfdy(2,1,1) =  beta*y(3,1);
dfdy(2,2,1) = -1.0/lambda;
dfdy(2,3,1) =  beta*y(1,1);
dfdy(3,2,1) =  1.0/lambda;
dfdy(3,3,1) = -1.0/eta;
function [ga] = gafun(ya, neq, m, nlbc)
% Evaluate boundary conditions at left-hand end of range.

ga = zeros(nlbc, 1);
ga(1) = ya(1) - ya(4);
ga(2) = ya(2) - ya(5);
ga(3) = ya(3) - ya(6);
function [dgady] = gajac(ya, neq, m, nlbc)
% Evaluate Jacobians (partial derivatives of ga).

dgady = zeros(nlbc, neq, 1);
dgady(1,1,1) =  1.0;
dgady(1,4,1) = -1.0;
dgady(2,2,1) =  1.0;
dgady(2,5,1) = -1.0;
dgady(3,3,1) =  1.0;
dgady(3,6,1) = -1.0;
function [gb] = gbfun(yb, neq, m, nrbc)
% Evaluate boundary conditions at right-hand end of range.

gb = zeros(nrbc, 1);
gb(1) = yb(1) - yb(4);
gb(2) = yb(2) - yb(5);
gb(3) = yb(3) - yb(6);
function [dgbdy] = gbjac(yb, neq, m, nrbc)
% Evaluate Jacobians (partial derivatives of gb).

dgbdy = zeros(nrbc, neq, 1);
dgbdy(1,1,1) =  1.0;
dgbdy(1,4,1) = -1.0;
dgbdy(2,2,1) =  1.0;
dgbdy(2,5,1) = -1.0;
dgbdy(3,3,1) =  1.0;
dgbdy(3,6,1) = -1.0;
function [y, dym] = guess(x, neq, m)
% Evaluate initial approximations to solution components and derivatives.

y = zeros(neq, 1);
dym = zeros(neq, 1);
y(1) = 1.0;
y(2) = 1.0;
y(3) = 1.0;
y(4) = y(1);
y(5) = y(2);
y(6) = y(3);
for i = 1:int32(neq) % can't use int64 in loop range.
    dym(i) = 0.0;
end
function display_plot(x, y)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot two of the curves, then add the other one.
[haxes, hline1, hline2] = plotyy(x, y(:,1), x, y(:,2), 'plot', ...
    'semilogy');
% We want the third curve to be plotted on the right-hand y-axis.
axes(haxes(2));
hold on
hline3 = plot(x, y(:,3));
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [0.06 0.08]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.06 0.065 0.07 0.075 0.08]);
set(haxes(2), 'YLim', [1e-10 1]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-8 1e-6 1e-4 1e-2 1]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 1]);
    set(haxes(iaxis), 'XTick', [0 0.2 0.4 0.6 0.8 1]);
    set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Add title.
title('Model of Spread of Measles', titleFmt{:});
% Label the axes.
xlabel('Time', labFmt{:});
ylabel(haxes(1), 'Susceptibles', labFmt{:});
ylabel(haxes(2), 'Infectives and Latents', labFmt{:});
% Add a legend.
legend('Infectives','Latents','Susceptibles','Location','Best')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.25, 'Marker', '+', 'Line', '-');
set(hline2, 'Linewidth', 0.25, 'Marker', 'x', 'Line', '--');
set(hline3, 'Linewidth', 0.25, 'Marker', '*', 'Line', ':');
 
nag_ode_bvp_coll_nlin_setup example program results

 Used a mesh of 21 points
 Maximum error =   1.42e-08 in interval 5 for component 1

 Mesh points:
   1(1) 0.0000   2(3) 0.0500   3(2) 0.1000   4(3) 0.1500
   5(2) 0.2000   6(3) 0.2500   7(2) 0.3000   8(3) 0.3500
   9(2) 0.4000  10(3) 0.4500  11(2) 0.5000  12(3) 0.5500
  13(2) 0.6000  14(3) 0.6500  15(2) 0.7000  16(3) 0.7500
  17(2) 0.8000  18(3) 0.8500  19(2) 0.9000  20(3) 0.9500
  21(1) 1.0000

 Computed solution at mesh points
    x       y1         y2         y3
  0.000  7.52e-02  1.80e-05  4.98e-06
  0.050  7.61e-02  7.89e-05  2.19e-05
  0.100  7.66e-02  3.15e-04  8.92e-05
  0.150  7.58e-02  1.01e-03  2.98e-04
  0.200  7.26e-02  2.25e-03  7.13e-04
  0.250  6.78e-02  3.11e-03  1.08e-03
  0.300  6.41e-02  2.56e-03  9.84e-04
  0.350  6.29e-02  1.29e-03  5.50e-04
  0.400  6.33e-02  4.14e-04  1.97e-04
  0.450  6.43e-02  9.12e-05  4.78e-05
  0.500  6.53e-02  1.59e-05  8.81e-06
  0.550  6.63e-02  2.77e-06  1.51e-06
  0.600  6.73e-02  6.28e-07  3.13e-07
  0.650  6.83e-02  2.19e-07  9.64e-08
  0.700  6.93e-02  1.24e-07  4.87e-08
  0.750  7.03e-02  1.16e-07  4.09e-08
  0.800  7.13e-02  1.70e-07  5.51e-08
  0.850  7.23e-02  3.70e-07  1.13e-07
  0.900  7.33e-02  1.11e-06  3.22e-07
  0.950  7.43e-02  4.20e-06  1.18e-06
  1.000  7.52e-02  1.80e-05  4.98e-06

function d02tv_example
global beta0 eta lambda mu % For communication with local functions

% Initialize variables and arrays.
neq = 6;
nlbc = 3;
nrbc = 3;
ncol = 5;
mmax = 1;

% Set up the mesh.
nmesh = 11;
mxmesh = 2*nmesh;
ipmesh = zeros(mxmesh, 1);
mesh = zeros(mxmesh, 1);

% Set location of mesh boundaries, then calculate initial spacing.
mesh(1) = 0.0;
mesh(nmesh) = 1.0;
mstep = (mesh(nmesh) - mesh(1))/double(nmesh-1);
for i = 2:nmesh-1
    mesh(i) = mstep*double(i-1);
    ipmesh(i) = 2;
end

% Specify mesh end points as fixed.
ipmesh(1)     = 1;
ipmesh(nmesh) = 1;

m = [1; 1; 1; 1; 1; 1];
tols = [1e-05; 1e-05; 1e-05; 1e-05; 1e-05; 1e-05];

fprintf('d02tv example program results\n\n');

% d02tv is a setup routine to be called prior to d02tk.
[work, iwork, ifail] = d02tv(int64(m), int64(nlbc), int64(nrbc), ...
    int64(ncol), tols, int64(nmesh), mesh, int64(ipmesh));
if ifail ~= 0
    % Unsuccessful call.  Print message and exit.
    error('Warning: d02tv returned with ifail = %1d ',ifail);
end

% Set values for problem-specific physical parameters.
eta = 0.01;
mu = 0.02;
lambda = 0.0279;
beta0 = 1575.0;

% Call d02tk to solve BVP for this set of parameters.
[work, iwork, ifail] = d02tk(@ffun, @fjac, ...
    @gafun, @gbfun, @gajac, @gbjac,...
    @guess, work, iwork);
if ifail ~= 0
    % Unsuccessful call.  Print message and exit.
    error('Warning: d02tk returned with ifail = %1d ',ifail);
end

% Call d02tz to extract mesh from solution.
[nmesh, mesh, ipmesh, ermx, iermx, ijermx, ifail] = d02tz( ...
    int64(mxmesh), work, iwork);

% Output mesh results.
fprintf(' Used a mesh of %d points\n', nmesh);
fprintf([' Maximum error = %10.2e in interval %d for component %d\n\n', ...
    ' Mesh points:\n'], ermx, iermx, ijermx);
for imesh = 1:int32(nmesh) % can't use int64 in loop range.
    fprintf( '%4d(%d) %6.4f', imesh, ipmesh(imesh), mesh(imesh));
    if mod(imesh, 4) == 0
        fprintf('\n');
    end
end

% Output solution, and store it for plotting.
xarray = zeros(nmesh, 1);
yarray = zeros(nmesh, 3);
fprintf('\n\n Computed solution at mesh points\n');
fprintf('    x       y1         y2         y3\n');
for imesh = 1:int32(nmesh) % can't use int64 in loop range.
    % Call d02ty to perform interpolation on the solution.
    [y, work, ifail] = d02ty(mesh(imesh), int64(neq), int64(mmax), ...
        work, iwork);
    fprintf(' %6.3f  %8.2e  %8.2e  %8.2e\n', mesh(imesh), ...
        y(1,1), y(2,1), y(3,1));
    xarray(imesh) = mesh(imesh);
    for jcomp = 1:3
        yarray(imesh, jcomp) = y(jcomp,1);
    end
end
% Plot results.
fig = figure('Number', 'off');
display_plot(xarray, yarray);


function [f] = ffun(x, y, neq, m)
% Evaluate derivative functions (rhs of system of ODEs).

global beta0 eta lambda mu % For communication with main routine.
beta = beta0*(1.0 + cos(2.0*pi*x));
f = zeros(neq, 1);
f(1) = mu - beta*y(1,1)*y(3,1);
f(2) = beta*y(1,1)*y(3,1) - y(2,1)/lambda;
f(3) = y(2,1)/lambda - y(3,1)/eta;
f(4) = 0.0;
f(5) = 0.0;
f(6) = 0.0;
function [dfdy] = fjac(x, y, neq, m)
% Evaluate Jacobians (partial derivatives of f).

global beta0 eta lambda mu % For communication with main routine.
beta = beta0*(1.0 + cos(2.0*pi*x));
dfdy = zeros(neq, neq, 1);
dfdy(1,1,1) = -beta*y(3,1);
dfdy(1,3,1) = -beta*y(1,1);
dfdy(2,1,1) =  beta*y(3,1);
dfdy(2,2,1) = -1.0/lambda;
dfdy(2,3,1) =  beta*y(1,1);
dfdy(3,2,1) =  1.0/lambda;
dfdy(3,3,1) = -1.0/eta;
function [ga] = gafun(ya, neq, m, nlbc)
% Evaluate boundary conditions at left-hand end of range.

ga = zeros(nlbc, 1);
ga(1) = ya(1) - ya(4);
ga(2) = ya(2) - ya(5);
ga(3) = ya(3) - ya(6);
function [dgady] = gajac(ya, neq, m, nlbc)
% Evaluate Jacobians (partial derivatives of ga).

dgady = zeros(nlbc, neq, 1);
dgady(1,1,1) =  1.0;
dgady(1,4,1) = -1.0;
dgady(2,2,1) =  1.0;
dgady(2,5,1) = -1.0;
dgady(3,3,1) =  1.0;
dgady(3,6,1) = -1.0;
function [gb] = gbfun(yb, neq, m, nrbc)
% Evaluate boundary conditions at right-hand end of range.

gb = zeros(nrbc, 1);
gb(1) = yb(1) - yb(4);
gb(2) = yb(2) - yb(5);
gb(3) = yb(3) - yb(6);
function [dgbdy] = gbjac(yb, neq, m, nrbc)
% Evaluate Jacobians (partial derivatives of gb).

dgbdy = zeros(nrbc, neq, 1);
dgbdy(1,1,1) =  1.0;
dgbdy(1,4,1) = -1.0;
dgbdy(2,2,1) =  1.0;
dgbdy(2,5,1) = -1.0;
dgbdy(3,3,1) =  1.0;
dgbdy(3,6,1) = -1.0;
function [y, dym] = guess(x, neq, m)
% Evaluate initial approximations to solution components and derivatives.

y = zeros(neq, 1);
dym = zeros(neq, 1);
y(1) = 1.0;
y(2) = 1.0;
y(3) = 1.0;
y(4) = y(1);
y(5) = y(2);
y(6) = y(3);
for i = 1:int32(neq) % can't use int64 in loop range.
    dym(i) = 0.0;
end
function display_plot(x, y)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% Plot two of the curves, then add the other one.
[haxes, hline1, hline2] = plotyy(x, y(:,1), x, y(:,2), 'plot', ...
    'semilogy');
% We want the third curve to be plotted on the right-hand y-axis.
axes(haxes(2));
hold on
hline3 = plot(x, y(:,3));
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [0.06 0.08]);
set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.06 0.065 0.07 0.075 0.08]);
set(haxes(2), 'YLim', [1e-10 1]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-8 1e-6 1e-4 1e-2 1]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 1]);
    set(haxes(iaxis), 'XTick', [0 0.2 0.4 0.6 0.8 1]);
    set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Add title.
title('Model of Spread of Measles', titleFmt{:});
% Label the axes.
xlabel('Time', labFmt{:});
ylabel(haxes(1), 'Susceptibles', labFmt{:});
ylabel(haxes(2), 'Infectives and Latents', labFmt{:});
% Add a legend.
legend('Infectives','Latents','Susceptibles','Location','Best')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.25, 'Marker', '+', 'Line', '-');
set(hline2, 'Linewidth', 0.25, 'Marker', 'x', 'Line', '--');
set(hline3, 'Linewidth', 0.25, 'Marker', '*', 'Line', ':');
 
d02tv example program results

 Used a mesh of 21 points
 Maximum error =   1.42e-08 in interval 5 for component 1

 Mesh points:
   1(1) 0.0000   2(3) 0.0500   3(2) 0.1000   4(3) 0.1500
   5(2) 0.2000   6(3) 0.2500   7(2) 0.3000   8(3) 0.3500
   9(2) 0.4000  10(3) 0.4500  11(2) 0.5000  12(3) 0.5500
  13(2) 0.6000  14(3) 0.6500  15(2) 0.7000  16(3) 0.7500
  17(2) 0.8000  18(3) 0.8500  19(2) 0.9000  20(3) 0.9500
  21(1) 1.0000

 Computed solution at mesh points
    x       y1         y2         y3
  0.000  7.52e-02  1.80e-05  4.98e-06
  0.050  7.61e-02  7.89e-05  2.19e-05
  0.100  7.66e-02  3.15e-04  8.92e-05
  0.150  7.58e-02  1.01e-03  2.98e-04
  0.200  7.26e-02  2.25e-03  7.13e-04
  0.250  6.78e-02  3.11e-03  1.08e-03
  0.300  6.41e-02  2.56e-03  9.84e-04
  0.350  6.29e-02  1.29e-03  5.50e-04
  0.400  6.33e-02  4.14e-04  1.97e-04
  0.450  6.43e-02  9.12e-05  4.78e-05
  0.500  6.53e-02  1.59e-05  8.81e-06
  0.550  6.63e-02  2.77e-06  1.51e-06
  0.600  6.73e-02  6.28e-07  3.13e-07
  0.650  6.83e-02  2.19e-07  9.64e-08
  0.700  6.93e-02  1.24e-07  4.87e-08
  0.750  7.03e-02  1.16e-07  4.09e-08
  0.800  7.13e-02  1.70e-07  5.51e-08
  0.850  7.23e-02  3.70e-07  1.13e-07
  0.900  7.33e-02  1.11e-06  3.22e-07
  0.950  7.43e-02  4.20e-06  1.18e-06
  1.000  7.52e-02  1.80e-05  4.98e-06


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013