Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_ode_bvp_coll_nlin_interp (d02ty)

## Purpose

nag_ode_bvp_coll_nlin_interp (d02ty) interpolates on the solution of a general two-point boundary value problem computed by nag_ode_bvp_coll_nlin_solve (d02tl).

## Syntax

[y, rcomm, ifail] = d02ty(x, neq, mmax, rcomm, icomm)
[y, rcomm, ifail] = nag_ode_bvp_coll_nlin_interp(x, neq, mmax, rcomm, icomm)

## Description

nag_ode_bvp_coll_nlin_interp (d02ty) and its associated functions (nag_ode_bvp_coll_nlin_solve (d02tl), nag_ode_bvp_coll_nlin_setup (d02tv), nag_ode_bvp_coll_nlin_contin (d02tx) and nag_ode_bvp_coll_nlin_diag (d02tz)) solve the two-point boundary value problem for a nonlinear mixed order system of ordinary differential equations
 $y1m1 x = f1 x,y1,y11,…,y1m1-1,y2,…,ynmn-1 y2m2 x = f2 x,y1,y11,…,y1m1-1,y2,…,ynmn-1 ⋮ ynmn x = fn x,y1,y11,…,y1m1-1,y2,…,ynmn-1$
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=\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
 $gizya=0, i=1,2,…,p,$
and the right boundary conditions at $b$ as
 $g-jzyb=0, j=1,2,…,q,$
where $y=\left({y}_{1},{y}_{2},\dots ,{y}_{n}\right)$ and
 $zyx = y1x, y11 x ,…, y1m1-1 x ,y2x,…, ynmn-1 x .$
First, nag_ode_bvp_coll_nlin_setup (d02tv) must be called to specify the initial mesh, error requirements and other details. Then, nag_ode_bvp_coll_nlin_solve (d02tl) can be used to solve the boundary value problem. After successful computation, nag_ode_bvp_coll_nlin_diag (d02tz) can be used to ascertain details about the final mesh and other details of the solution procedure, and nag_ode_bvp_coll_nlin_interp (d02ty) can be used to compute the approximate solution anywhere on the interval $\left[a,b\right]$ using interpolation.
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
Grossman C (1992) Enclosures of the solution of the Thomas–Fermi equation by monotone discretization J. Comput. Phys. 98 26–32
Keller H B (1992) Numerical Methods for Two-point Boundary-value Problems Dover, New York

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{x}$ – double scalar
$x$, the independent variable.
Constraint: $a\le {\mathbf{x}}\le b$, i.e., not outside the range of the original mesh specified in the initialization call to nag_ode_bvp_coll_nlin_setup (d02tv).
2:     $\mathrm{neq}$int64int32nag_int scalar
The number of differential equations.
Constraint: ${\mathbf{neq}}$ must be the same value as supplied to nag_ode_bvp_coll_nlin_setup (d02tv).
3:     $\mathrm{mmax}$int64int32nag_int scalar
The maximal order of the differential equations, $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({m}_{\mathit{i}}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$.
Constraint: ${\mathbf{mmax}}$ must contain the maximum value of the components of the argument m as supplied to nag_ode_bvp_coll_nlin_setup (d02tv).
4:     $\mathrm{rcomm}\left(*\right)$ – double array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument rcomm in the previous call to nag_ode_bvp_coll_nlin_solve (d02tl).
This must be the same array as supplied to nag_ode_bvp_coll_nlin_solve (d02tl) and must remain unchanged between calls.
5:     $\mathrm{icomm}\left(*\right)$int64int32nag_int array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument icomm in the previous call to nag_ode_bvp_coll_nlin_solve (d02tl).
This must be the same array as supplied to nag_ode_bvp_coll_nlin_solve (d02tl) and must remain unchanged between calls.

None.

### Output Parameters

1:     $\mathrm{y}\left({\mathbf{neq}},{\mathbf{mmax}}\right)$ – double array
${\mathbf{y}}\left(\mathit{i},\mathit{j}\right)$ contains an approximation to ${y}_{\mathit{i}}^{\left(\mathit{j}\right)}\left(x\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{neq}}$ and $\mathit{j}=0,1,\dots ,{m}_{\mathit{i}}-1$. The remaining elements of y (where ${m}_{i}<{\mathbf{mmax}}$) are initialized to $0.0$.
2:     $\mathrm{rcomm}\left(*\right)$ – double array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument rcomm in the previous call to nag_ode_bvp_coll_nlin_solve (d02tl).
Contains information about the solution for use on subsequent calls to associated functions.
3:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Note: nag_ode_bvp_coll_nlin_interp (d02ty) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:
${\mathbf{ifail}}=1$
Constraint: ${\mathbf{mmax}}=\mathrm{max}\phantom{\rule{0.25em}{0ex}}{\mathbf{m}}\left(i\right)$.
Constraint: ${\mathbf{neq}}={\mathbf{neq}}$ in nag_ode_bvp_coll_nlin_setup (d02tv).
Constraint: ${\mathbf{x}}\le _$.
The solver function did not produce any results suitable for interpolation.
The solver function does not appear to have been called.
x is too small.
${\mathbf{ifail}}=2$
The solver function did not converge to a suitable solution.
A converged intermediate solution has been used.
Interpolated values should be treated with caution.
The solver function did not satisfy the error requirements.
Interpolated values should be treated with caution.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

If nag_ode_bvp_coll_nlin_interp (d02ty) returns the value ${\mathbf{ifail}}={\mathbf{0}}$, the computed values of the solution components ${y}_{i}$ should be of similar accuracy to that specified by the argument tols of nag_ode_bvp_coll_nlin_setup (d02tv). Note that during the solution process the error in the derivatives ${y}_{i}^{\left(\mathit{j}\right)}$, for $\mathit{j}=1,2,\dots ,{m}_{i}-1$, has not been controlled and that the derivative values returned by nag_ode_bvp_coll_nlin_interp (d02ty) are computed via differentiation of the piecewise polynomial approximation to ${y}_{i}$. See also Further Comments.

If nag_ode_bvp_coll_nlin_interp (d02ty) returns the value ${\mathbf{ifail}}={\mathbf{2}}$ in this function and ${\mathbf{ifail}}={\mathbf{5}}$ in nag_ode_bvp_coll_nlin_solve (d02tl), then the accuracy of the interpolated values may be proportional to the quantity ermx as returned by nag_ode_bvp_coll_nlin_diag (d02tz).
If nag_ode_bvp_coll_nlin_solve (d02tl) returned a value for ifail other than ${\mathbf{ifail}}={\mathbf{0}}$, then nothing can be said regarding either the quality or accuracy of the values computed by nag_ode_bvp_coll_nlin_interp (d02ty).

## Example

The following example is used to illustrate that a system with singular coefficients can be treated without modification of the system definition. See also nag_ode_bvp_coll_nlin_solve (d02tl), nag_ode_bvp_coll_nlin_setup (d02tv), nag_ode_bvp_coll_nlin_contin (d02tx) and nag_ode_bvp_coll_nlin_diag (d02tz), for the illustration of other facilities.
Consider the Thomas–Fermi equation used in the investigation of potentials and charge densities of ionized atoms. See Grossman (1992), for example, and the references therein. The equation is
 $y′′=x-1/2y3/2$
with boundary conditions
 $y0= 1, ya= 0, a> 0.$
The coefficient ${x}^{-1/2}$ implies a singularity at the left-hand boundary $x=0$.
We use the initial approximation $y\left(x\right)=1-x/a$, which satisfies the boundary conditions, on a uniform mesh of six points. For illustration we choose $a=1$, as in Grossman (1992). Note that in ffun and fjac (see nag_ode_bvp_coll_nlin_solve (d02tl)) we have taken the precaution of setting the function value and Jacobian value to $0.0$ in case a value of $y$ becomes negative, although starting from our initial solution profile this proves unnecessary during the solution phase. Of course the true solution $y\left(x\right)$ is positive for all $x.
```function d02ty_example

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

global a; % For communication with local functions

% Initialize variables and arrays.
neq  = int64(1);
nlbc = int64(1);
nrbc = int64(1);
ncol = int64(4);
mmax = int64(2);
m    = int64([2]);
tols = [1.0e-05];

% Set values for problem-specific physical parameters.
a = 1.0;

% Set up the mesh.
nmesh  = int64(6);
mxmesh = int64(100);
ipmesh = zeros(mxmesh, 1, 'int64');
mesh   = zeros(mxmesh, 1);

% Set initial mesh with constant spacing.
dx = a/double(nmesh-1);
mesh(1:nmesh) = [0:dx:a];

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

% d02tv is a setup routine to be called prior to d02tk.
[work, iwork, ifail] = d02tv(...
m, nlbc, nrbc, ncol, tols, nmesh, mesh, ipmesh);
fprintf('\n Tolerance = %8.1e  A = %8.2f\n\n', tols(1), a);

% 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\n', ...
' Mesh points:\n'], ermx, iermx, ijermx);
for imesh = 1:nmesh
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, 2);
fprintf('\n\n Computed solution\n');
fprintf('       x     solution   derivative\n');
for imesh = 1:nmesh
% Call d02ty to perform interpolation on the solution.
xx = mesh(imesh);
[y, work, ifail] = d02ty(...
xx, neq, mmax, work, iwork);
fprintf(' %8.2f %10.5f %10.5f\n', xx, y(1,1), y(1,2));
xarray(imesh) = xx;
yarray(imesh, :) = y(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).

f = zeros(neq, 1);
if y(1,1) < 0
f(1,1) = 0;
fprintf(' In ffun\n');
else
f(1,1) = y(1,1)^1.5/sqrt(x);
end

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

dfdy = zeros(neq, neq, 2);
if y(1,1) < 0
dfdy(1,1,1) = 0;
fprintf(' In fjac\n');
else
dfdy(1,1,1) = 1.5*sqrt(y(1,1))/sqrt(x);
end

function [ga] = gafun(ya, neq, m, nlbc)
% Evaluate boundary conditions at left-hand end of range.

ga = zeros(nlbc, 1);
ga(1,1) = ya(1,1) - 1;

function [dgady] = gajac(ya, neq, m, nlbc)
% Evaluate Jacobians (partial derivatives of ga).

function [gb] = gbfun(yb, neq, m, nrbc)
% Evaluate boundary conditions at right-hand end of range.

gb = zeros(nrbc, 1);
gb(1,1) = yb(1,1);

function [dgbdy] = gbjac(yb, neq, m, nrbc)
% Evaluate Jacobians (partial derivatives of gb).

dgbdy = zeros(nrbc, neq, 2);
dgbdy(1,1,1) = 1;

function [y, dym] = guess(x, neq, m)
% Evaluate initial approximations to solution components and derivatives.

global a
y = zeros(neq, 2);
dym = zeros(neq, 1);
y(1,1) =  1 - x/a;
y(1,2) = -1/a;

dym(1) = 0;

function display_plot(x, y)
% Plot both curves.
[hline] = plot(x, y(:,1), x, y(:,2));
title({['Thomas-Fermi Equation for Determining'],...
['Effective Nuclear Charge in Heavy Atoms']});
% Label the axes.
xlabel('Relative Distance');
ylabel('Effective Nuclear Charge');
% Set some features of the three lines.
set(hline(1), 'Linewidth', 0.25, 'Marker', '+', 'LineStyle', '-');
set(hline(2), 'Linewidth', 0.25, 'Marker', 'x', 'LineStyle', '--');
```
```d02ty example results

Tolerance =  1.0e-05  A =     1.00

Used a mesh of 11 points
Maximum error =   3.09e-06 in interval 1 for component 1

Mesh points:
1(1) 0.0000   2(3) 0.1000   3(2) 0.2000   4(3) 0.3000
5(2) 0.4000   6(3) 0.5000   7(2) 0.6000   8(3) 0.7000
9(2) 0.8000  10(3) 0.9000  11(1) 1.0000

Computed solution
x     solution   derivative
0.00    1.00000   -1.84496
0.10    0.84944   -1.32330
0.20    0.72721   -1.13911
0.30    0.61927   -1.02776
0.40    0.52040   -0.95468
0.50    0.42754   -0.90583
0.60    0.33867   -0.87372
0.70    0.25239   -0.85369
0.80    0.16764   -0.84248
0.90    0.08368   -0.83756
1.00    0.00000   -0.83655
```