NAG Toolbox: nag_ode_bvp_coll_nth (d02ja)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{n}$ – int64int32nag_int scalar

$n$, the order of the differential equation.

Constraint: $\mathbf{n}\ge 1$.

2:     $\mathrm{cf}$ – function handle or string containing name of m-file

cf defines the differential equation (see Description). It must return the value of a function $f_j\left(x\right)$ at a given point $x$, where, for $1\le j\le n+1$, $f_j\left(x\right)$ is the coefficient of $y^{\left(j-1\right)}\left(x\right)$ in the equation, and $f_0\left(x\right)$ is the right-hand side.

[result] = cf(j, x)

Input Parameters

1:     $\mathrm{j}$ – int64int32nag_int scalar

The index of the function $f_j$ to be evaluated.

2:     $\mathrm{x}$ – double scalar

The point at which $f_j$ is to be evaluated.

Output Parameters

1:     $\mathrm{result}$ – double scalar

The value of $f_j\left(x\right)$ at the given point $x$.

3:     $\mathrm{bc}$ – function handle or string containing name of m-file

bc defines the boundary conditions, each of which has the form $y^{\left(k-1\right)}\left(x_1\right)=s_k$ or $y^{\left(k-1\right)}\left(x_0\right)=s_k$. The boundary conditions may be specified in any order.

[j, rhs] = bc(ii)

Input Parameters

1:     $\mathrm{ii}$ – int64int32nag_int scalar

The index of the boundary condition to be defined.

Output Parameters

1:     $\mathrm{j}$ – int64int32nag_int scalar

Must be set to $-k$ if the boundary condition is $y^{\left(k-1\right)}\left(x_0\right)=s_k$, and to $+k$ if it is $y^{\left(k-1\right)}\left(x_1\right)=s_k$.

j must not be set to the same value $k$ for two different values of ii.

2:     $\mathrm{rhs}$ – double scalar

Must be set to the value $s_k$.

4:     $\mathrm{x0}$ – double scalar

5:     $\mathrm{x1}$ – double scalar

The left- and right-hand boundaries, $x_0$ and $x_1$, respectively.

Constraint: $\mathbf{x1}>\mathbf{x0}$.

6:     $\mathrm{k1}$ – int64int32nag_int scalar

The number of coefficients to be returned in the Chebyshev series representation of the solution (hence the degree of the polynomial approximation is $\mathbf{k1}-1$).

Constraint: $\mathbf{k1}\ge \mathbf{n}+1$.

7:     $\mathrm{kp}$ – int64int32nag_int scalar

The number of collocation points to be used.

Constraint: $\mathbf{kp}\ge \mathbf{k1}-\mathbf{n}$.

Optional Input Parameters

None.

Output Parameters

1:     $\mathrm{c}\left(\mathbf{k1}\right)$ – double array

The computed Chebyshev coefficients; that is, the computed solution is:

$$\underset{i=1}{\overset{\mathbf{k1}}{{\sum }^{\prime }}}\mathbf{c}\left(i\right)T_{i-1}\left(x\right)$$

where $T_i\left(x\right)$ is the $i$th Chebyshev polynomial of the first kind, and ${\sum }^{\prime }$ denotes that the first coefficient, $\mathbf{c}\left(1\right)$, is halved.

2:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

   $\mathbf{ifail}=1$

On entry,	$\mathbf{n}<1$,
or	$\mathbf{x0}\ge \mathbf{x1}$,
or	$\mathbf{k1}<\mathbf{n}+1$,
or	$\mathbf{kp}<\mathbf{k1}-\mathbf{n}$.

   $\mathbf{ifail}=2$

On entry,

$\mathit{lw}<2\times \left(\mathbf{kp}+\mathbf{n}\right)\times \left(\mathbf{k1}+1\right)+7\times \mathbf{k1}$ (insufficient workspace).

   $\mathbf{ifail}=3$

Either the boundary conditions are not linearly independent (that is, in bc the variable j is set to the same value $k$ for two different values of ii), or the rank of the matrix of equations for the coefficients is less than the number of unknowns. Increasing kp may overcome this latter problem.

   $\mathbf{ifail}=4$

The least squares function nag_linsys_real_gen_lsqsol (f04am) has failed to correct the first approximate solution (see nag_linsys_real_gen_lsqsol (f04am)).

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: d02ja_example

function d02ja_example


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

% Set up initial values.
n     = int64(2);
k1max = 8;
kpmax = 15;
lw    = 2*(kpmax+n)*(k1max+1)+7*k1max;
x0    = -1;
x1    = 1;
c     = zeros(1,kpmax);
p     = zeros(1,kpmax);

fprintf('  KP   K1   Chebyshev coefficients\n');
for kp = int64(10:5:kpmax)
  for k1 = int64(4:2:k1max)

        [c, ifail] = d02ja(n, @cf, @bc, x0, x1, k1, kp);
        % Output results.
        fprintf('%4d ',kp, k1);

        for kind = 1:k1
            fprintf('%8.4f  ',c(kind));
            if mod(kind, 6) == 0 && kind ~= k1
                fprintf('\n          ');
            end
        end
        fprintf('\n');
    end
end

% Now prepare to evaluate and plot the last solution.
fprintf('\n');
k1     = int64(8);
k1m1   = int64(k1-1);
m      = 9;
ia1    = int64(1);
xarray = zeros(k1,1);
yarray = zeros(k1,1);

fprintf('Last computed solution evaluated at %1d equally spaced points\n\n', m);
fprintf('    X         Y\n');
for i = 1:m
    x = (x0*double(m-i)+x1*double(i-1))/double(m-1);

    % Calling e02ak to evaluate the polynomial from its Chebyshev
    % representation.
    [y, ifail] = e02ak(k1m1, x0, x1, c, ia1, x);

    % Output results, and save them for plotting.
    fprintf('%8.4f  %8.4f  \n',x,y);
    xarray(i) = x;
    yarray(i) = y;
end
fig1 = figure;
display_plot(xarray, yarray)


function [j, rhs] = bc(i)
% Define the boundary conditions.
rhs = 0;
if (i == 1)
    j = int64(1);
else
    j = int64(-1);
end

function result = cf(j, x)
% Define the differential equation to be solved.
if (j == 2)
    result = 0;
else
    result = 1;
end

function display_plot(x, y)
% Plot results.
plot(x, y, '-+');
% Add title.
title({'Two-point Boundary-value Problem for ODE',...
    'by Chebyshev-series using Collocation and Least-squares'})
% Label the axes.
xlabel('x');
ylabel('y');

d02ja example results

  KP   K1   Chebyshev coefficients
  10    4  -0.6108   -0.0000    0.3054    0.0000  
  10    6  -0.8316   -0.0000    0.4246    0.0000   -0.0088   -0.0000  
  10    8  -0.8325   -0.0000    0.4253    0.0000   -0.0092    0.0000  
            0.0001   -0.0000  
  15    4  -0.6174   -0.0000    0.3087    0.0000  
  15    6  -0.8316   -0.0000    0.4246    0.0000   -0.0088   -0.0000  
  15    8  -0.8325   -0.0000    0.4253    0.0000   -0.0092   -0.0000  
            0.0001    0.0000  

Last computed solution evaluated at 9 equally spaced points

    X         Y
 -1.0000    0.0000  
 -0.7500   -0.3542  
 -0.5000   -0.6242  
 -0.2500   -0.7933  
  0.0000   -0.8508  
  0.2500   -0.7933  
  0.5000   -0.6242  
  0.7500   -0.3542  
  1.0000    0.0000