d02tv:: Ordinary Differential EquationsIntegrators for Stiff Ordinary Differential Systems (NAG Toolbox)

nag_ode_bvp_coll_nlin_setup (d02tv) and its associated functions (nag_ode_bvp_coll_nlin_solve (d02tl), 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

\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 Further Comments for information on how boundary value problems of a more general nature can be treated.

The solution function nag_ode_bvp_coll_nlin_solve (d02tl) 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 nag_ode_bvp_coll_nlin_solve (d02tl).

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

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.

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.

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 .

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 one, then you should make those derivatives new variables.

Fixed Points and Singularities

Numerical Jacobians

Example

function d02tv_example


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

% Global constants communicated to local functions ffun and fjac
global  beta0 eta lambda mu
% Set values for problem-specific physical parameters.
eta    = 0.01;
mu     = 0.02;
lambda = 0.0279;
beta0  = 1575.0;

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

% Set up the mesh.
nmesh  = int64(7);
mxmesh = int64(8*nmesh);
ipmesh = zeros(mxmesh, 1, 'int64');
mesh   = zeros(mxmesh, 1);

% Set location of mesh boundaries, then calculate initial spacing.
mstep           = 1/double(nmesh-1);
mesh(1:nmesh)   = 0:mstep:1;
ipmesh(1:nmesh) = int64(2);

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

m(1:neq)      = int64(1);
tols(1:neq)   = 1e-6;

% d02tv is a setup routine to be called prior to d02tk.
[work, iwork, ifail] = d02tv(...
                             m, nlbc, nrbc, ncol, tols, nmesh, mesh, ipmesh);

% Call d02tk to solve BVP for this set of parameters.
[work, iwork, ifail] = d02tk(...
                             @ffun, @fjac, @gafun, @gbfun, @gajac, @gbjac,...
                             @guess, work, iwork);

% Call d02tz to extract mesh from solution.
[nmesh, mesh, ipmesh, ermx, iermx, ijermx, ifail] = ...
d02tz(...
       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',...
        ermx, iermx, ijermx);

% Output solution, and store it for plotting.
yarray           = zeros(nmesh, 3);
xarray(1:nmesh)  = mesh(1:nmesh);
for i = 1:nmesh
  % Call d02ty to perform interpolation on the solution.
  [y, work, ifail] = d02ty(mesh(i), neq, mmax, work, iwork);
  yarray(i, 1:3)   = y(1:3,1);
end

% Plot results.
fig1 = figure;
display_plot(xarray, yarray);


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

  beta = beta0*(1 + cos(2*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;

function [dfdy] = fjac(x, y, neq, m)
  % Evaluate Jacobians (partial derivatives of f).
  global beta0 eta lambda mu

  beta = beta0*(1 + cos(2*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/lambda;
  dfdy(2,3,1) =  beta*y(1,1);
  dfdy(3,2,1) =  1/lambda;
  dfdy(3,3,1) = -1/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;
  dgady(1,4,1) = -1;
  dgady(2,2,1) =  1;
  dgady(2,5,1) = -1;
  dgady(3,3,1) =  1;
  dgady(3,6,1) = -1;

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;
  dgbdy(1,4,1) = -1;
  dgbdy(2,2,1) =  1;
  dgbdy(2,5,1) = -1;
  dgbdy(3,3,1) =  1;
  dgbdy(3,6,1) = -1;

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:3) = 1;
  y(4:6) = y(1:3);

function display_plot(x, y)
  % Plot two of the curves, then add the other one.
  [haxes, hline1, hline2] = plotyy(x, y(:,1), x, y(:,2), 'plot', 'semilogy');
  % 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(2), 'YLim', [1e-10 1]);
  set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
  set(haxes(2), 'YMinorTick', 'on');
  set(haxes(1), 'YTick', [0.06  0.065 0.07 0.075 0.08]);
  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:1]);
  end
  set(gca, 'box', 'off');
  % Add title.
  title('Model of Spread of Measles');
  % Label the axes.
  xlabel('Time');
  ylabel(haxes(1), 'Susceptibles');
  ylabel(haxes(2), 'Infectives and Latents');
  % Add a legend.
  legend('Infectives','Latents','Susceptibles','Location','NorthEast')
  % Set some features of the three lines.
  set(hline1, 'Linewidth', 0.25, 'Marker', '+', 'LineStyle', '-');
  set(hline2, 'Linewidth', 0.25, 'Marker', 'x', 'LineStyle', '--', ...
      'Color', 'Magenta');
  set(hline3, 'Linewidth', 0.25, 'Marker', '*', 'LineStyle', ':');

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

NAG Toolbox: nag_ode_bvp_coll_nlin_setup (d02tv)

▸▿ Contents

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Nonseparated Boundary Conditions

Multipoint Boundary Value Problems

High Order Systems

Fixed Points and Singularities

Numerical Jacobians

Example