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

Description

nag_ode_bvp_coll_nth_comp (d02tg) calculates an approximate solution of a linear or linearized system of ordinary differential equations as a Chebyshev series. Suppose there are

n

differential equations for

n

variables

y_{1}, y_{2}, \dots, y_{n}

, over the range

(x_{0}, x_{1})

. Let the

i

th equation be

\sum_{j = 1}^{m_{i} + 1} \sum_{k = 1}^{n} f_{k j}^{i} (x) y_{k}^{(j - 1)} (x) = r^{i} (x)

where

y_{k}^{(j)} (x) = \frac{d^{j} y_{k} (x)}{d x^{j}}

. coeff evaluates the coefficients

f_{k j}^{i} (x)

and the right-hand side

r^{i} (x)

for each

i

1 \leq i \leq n

, at any point

x

. The boundary conditions may be applied either at the end points or at intermediate points; they are written in the same form as the differential equations, and specified by bdyc. For example the

j

th boundary condition out of those associated with the

i

th differential equation takes the form

\sum_{j = 1}^{l_{i} + 1} \sum_{k = 1}^{n} f_{k j}^{i j} (x^{i j}) y_{k}^{(j - 1)} (x^{i j}) = r^{i j} (x^{i j}),

where

x^{i j}

lies between

x_{0}

and

x_{1}

. It is assumed in this function that certain of the boundary conditions are associated with each differential equation. This is for your convenience; the grouping does not affect the results.

The degree of the polynomial solution must be the same for all variables. You specify the degree required,

k_{1} - 1

, and the number of collocation points,

k_{p}

, in the range. The function sets up a system of linear equations for the Chebyshev coefficients, with

n

equations for each collocation point and one for each boundary condition. The collocation points are chosen at the extrema of a shifted Chebyshev polynomial of degree

k_{p} - 1

. The boundary conditions are satisfied exactly, and the remaining equations are solved by a least squares method. The result produced is a set of Chebyshev coefficients for the

n

functions

y_{1}, y_{2}, \dots, y_{n}

, with the range normalized to

[- 1, 1]

nag_fit_1dcheb_eval2 (e02ak) can be used to evaluate the components of the solution at any point on the range

[x_{0}, x_{1}]

(see Example for an example). nag_fit_1dcheb_deriv (e02ah) and nag_fit_1dcheb_integ (e02aj) may be used to obtain Chebyshev series representations of derivatives and integrals (respectively) of the components of the solution.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

If the number of matching points

k_{p}

is equal to the number of coefficients

k_{1}

minus the average number of boundary conditions

\frac{1}{n} \sum_{i = 1}^{n} l_{i}

, then the least squares solution reduces to simple solution of linear equations and true collocation results. The accuracy of the solution may be checked by repeating the calculation with different values of

k_{1}

. If the Chebyshev coefficients decrease rapidly, the size of the last two or three gives an indication of the error. If they do not decrease rapidly, it may be desirable to use a different method. Note that the Chebyshev coefficients are calculated for the range normalized to

[- 1, 1]

A system of nonlinear differential equations must be linearized before using the function. The calculation is repeated iteratively. On each iteration the linearized equation is used. In the example in Example, the

y

variables are to be determined at the current iteration whilst the

z

variables correspond to the solution determined at the previous iteration, (or the initial approximation on the first iteration). For a starting approximation, we may take, say, a linear function, and set up the appropriate Chebyshev coefficients before starting the iteration. For example, if

y_{1} = a x + b

in the range

(x_{0}, x_{1})

, we set

B

, the array of coefficients,

$B (1, 1) = a \times (x_{0} + x_{1}) + 2 \times b$ ,
$B (1, 2) = a \times (x_{1} - x_{0}) / 2$ ,
and the remainder of the entries to zero.

Example

Suppose an approximate solution is

z_{1}

z_{2}

such that

y_{1} \sim z_{1}

y_{2} \sim z_{2}

: then the first equation gives, on linearizing,

2 y_{1}^{'} + (z_{2}^{2} - 1) y_{1} + (2 z_{1} z_{2} + 1) y_{2} = 2 z_{1} z_{2}^{2} .

The starting approximation is taken to be

z_{1} = 0

z_{2} = 3

. In the program below, the array

B

is used to hold the coefficients of the previous iterate (or of the starting approximation). We iterate until the Chebyshev coefficients converge to five figures. nag_fit_1dcheb_eval2 (e02ak) is used to calculate the solution from its Chebyshev coefficients.

function d02tg_example


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

global b x0 x1 ndegree

% Initialize variables and arrays.
n       = 2;
nderiv  = int64([1 2]);
nbc     = nderiv;
x0      = -1;
x1      = 1;
ncoeff  = int64(9);
ncolloc = int64(15);
ndegree = ncoeff-1;

% Iterate to covergence within tolerance of 1.0e0-5.
iter = 0;
emax = 1.0;
b = zeros(ncoeff, n);
b(1,2) = 6.0;
while emax >= 1.0e-05
  iter = iter + 1;

  [c, ifail] = d02tg(...
                     nderiv, nbc, x0, x1, ncoeff, ncolloc, ...
                     @coeff, @bdyc);
  % Output coefficients.
  fprintf('\n Iteration %2d. Chebyshev coefficients are:\n',iter);
  for j = 1:n
    fprintf('%8.4f',c(1:ncoeff,j));
    fprintf('\n');
  end
  
  emax = max(abs(c - b));
  b = c;
end

% Use these coefficients to evaluate the solution.
npt = 9;
fprintf('\n Solution evaluated at %d equally spaced points.\n\n', npt);
fprintf('%7s%16s\n','x','y');

% Prepare to store results for plotting.
dx = (x1-x0)/(double(npt-1));
xarray = [x0:dx:x1];
yarray = zeros(npt, n);

for ipt = 1:npt
  xx = xarray(ipt);
  fprintf('%10.4f', xx);
  for jeqn = 1:n
    [yarray(ipt, jeqn), ifail] = ...
    e02ak(...
           ndegree, x0, x1, c(:,jeqn), int64(1), xx);
    fprintf('%10.4f', yarray(ipt, jeqn));
  end
  fprintf('\n');
end

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


function [x, a, rhs] = bdyc(ii, jj, a, ia, ia1, rhs)
  % Evaluate coefficient functions and rhs of boundary conditions.

  x = -1;
  a(ii,jj) = 1;
  if (ii == 2 && jj == 1)
    rhs = 3;
  end

function [a, rhs] = coeff(x, ii, a, ia, ia1, rhs)
  % Evaluate coefficient functions and rhs of differential equations.

  global b x0 x1 ndegree

  if (ii <= 1)
    [z1, ifail] = e02ak(ndegree, x0, x1, b(:,1), int64(1), x);
    [z2, ifail] = e02ak(ndegree, x0, x1, b(:,2), int64(1), x);
    a(1,1) = z2*z2 - 1;
    a(1,2) = 2;
    a(2,1) = 2*z1*z2 + 1;
    rhs = 2*z1*z2*z2;
  else
    a(1,2) = -1;
    a(2,3) =  2;
  end

function display_plot(x, y)
  % Plot both curves.
  [haxes, hline1, hline2] = plotyy(x, y(:,1), x, y(:,2));
  % Set the axis limits and the tick specifications to beautify the plot.
  set(haxes(1), 'YLim', [-0.5 0.1]);
  set(haxes(1), 'YMinorTick', 'on');
  set(haxes(1), 'YTick', [-0.5:0.1:0]);
  set(haxes(2), 'YLim', [2.6 3.2]);
  set(haxes(2), 'YMinorTick', 'on');
  set(haxes(2), 'YTick', [2.6:0.1:3.2]);
  for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [-1 1]);
  end
  % Add title.
  title('Solution by Chebyshev Collocation');
  % Label the axes.
  xlabel('x');
  ylabel(haxes(1), 'y_1');
  ylabel(haxes(2), 'y_2');
  % Add a legend.
  legend('y_1','y_2','Location','Best')
  % Set some features of the three lines.
  set(hline1, 'Linewidth', 0.25, 'Marker', '+', 'LineStyle', '-');
  set(hline2, 'Linewidth', 0.25, 'Marker', 'x', 'LineStyle', '--');

d02tg example results


 Iteration  1. Chebyshev coefficients are:
 -0.5659 -0.1162  0.0906 -0.0468  0.0196 -0.0069  0.0021 -0.0006  0.0001
  5.7083 -0.1642 -0.0087  0.0059 -0.0025  0.0009 -0.0003  0.0001 -0.0000

 Iteration  2. Chebyshev coefficients are:
 -0.6338 -0.1599  0.0831 -0.0445  0.0193 -0.0071  0.0023 -0.0006  0.0001
  5.6881 -0.1792 -0.0144  0.0053 -0.0023  0.0008 -0.0003  0.0001 -0.0000

 Iteration  3. Chebyshev coefficients are:
 -0.6344 -0.1604  0.0828 -0.0446  0.0193 -0.0071  0.0023 -0.0006  0.0001
  5.6880 -0.1793 -0.0145  0.0053 -0.0023  0.0008 -0.0003  0.0001 -0.0000

 Iteration  4. Chebyshev coefficients are:
 -0.6344 -0.1604  0.0828 -0.0446  0.0193 -0.0071  0.0023 -0.0006  0.0001
  5.6880 -0.1793 -0.0145  0.0053 -0.0023  0.0008 -0.0003  0.0001 -0.0000

 Solution evaluated at 9 equally spaced points.

      x               y
   -1.0000    0.0000    3.0000
   -0.7500   -0.2372    2.9827
   -0.5000   -0.3266    2.9466
   -0.2500   -0.3640    2.9032
    0.0000   -0.3828    2.8564
    0.2500   -0.3951    2.8077
    0.5000   -0.4055    2.7577
    0.7500   -0.4154    2.7064
    1.0000   -0.4255    2.6538

On entry,	$n < 1$ ,
or	$m (i) < 1$ for some $i$ ,
or	$l (i) < 0$ for some $i$ ,
or	$x0 \geq x1$ ,
or	$k1 < 1 + m (i)$ for some $i$ ,
or	$n \times kp < n \times k1 + \sum_{i = 1}^{n} l (i)$ ,
or	$ldc < k1$ .

On entry,	lw is too small (see Arguments),
or	$liw < n \times k1$ .

NAG Toolbox: nag_ode_bvp_coll_nth_comp (d02tg)

▸▿ Contents

Purpose

Syntax