NAG Toolbox: nag_ode_bvp_shoot_genpar (d02hb)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{p}\left(\mathbf{n1}\right)$ – double array

An estimate for the $\mathit{i}$th argument, $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$.

2:     $\mathrm{pe}\left(\mathbf{n1}\right)$ – double array

The elements of pe must be given small positive values. The element $\mathbf{pe}\left(i\right)$ is used

(i)	in the convergence test on the $i$th argument in the Newton iteration, and
(ii)	in perturbing the $i$th argument when approximating the derivatives of the components of the solution with respect to this argument for use in the Newton iteration.

The elements $\mathbf{pe}\left(i\right)$ should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.

Constraint: $\mathbf{pe}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,\mathbf{n1}$.

3:     $\mathrm{e}\left(\mathbf{n}\right)$ – double array

The elements of e must be given positive values. The element $\mathbf{e}\left(i\right)$ is used in the bound on the local error in the $i$th component of the solution $y_i$ during integration.

The elements $\mathbf{e}\left(i\right)$ should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.

Constraint: $\mathbf{e}\left(\mathit{i}\right)>0.0$, for $\mathit{i}=1,2,\dots ,\mathbf{n}$.

4:     $\mathrm{m1}$ – int64int32nag_int scalar

A value which controls exit values.

$\mathbf{m1}=1$

The final solution is not calculated.

$\mathbf{m1}>1$

The final values of the solution at interval (length of range)/$\left(\mathbf{m1}-1\right)$ are calculated and stored sequentially in the array soln starting with the values of the solutions evaluated at the first end point (see range) stored in the first column of soln.

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

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

fcn must evaluate the functions $f_{\mathit{i}}$ (i.e., the derivatives $y_{\mathit{i}}^{\prime }$), for $\mathit{i}=1,2,\dots ,\mathit{n}$, at a general point $x$.

[f] = fcn(x, y, p)

Input Parameters

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

$x$, the value of the argument.

2:     $\mathrm{y}\left(:\right)$ – double array

$y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, the value of the argument.

3:     $\mathrm{p}\left(:\right)$ – double array

The current estimate of the argument $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$.

Output Parameters

1:     $\mathrm{f}\left(:\right)$ – double array

The value of $f_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$. The $f_i$ may depend upon the parameters $p_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,{\mathit{n}}_1$. If there are any driving equations (see Description) then these must be numbered first in the ordering of the components of f in fcn.

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

bc must place in g1 and g2 the boundary conditions at $a$ and $b$ respectively (see range).

[g1, g2] = bc(p)

Input Parameters

1:     $\mathrm{p}\left(:\right)$ – double array

An estimate of the argument $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$.

Output Parameters

1:     $\mathrm{g1}\left(:\right)$ – double array

The value of $y_{\mathit{i}}\left(a\right)$, (where this may be a known value or a function of the parameters $p_{\mathit{j}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$ and $\mathit{j}=1,2,\dots ,{\mathit{n}}_1$).

2:     $\mathrm{g2}\left(:\right)$ – double array

The value of $y_{\mathit{i}}\left(b\right)$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, (where these may be known values or functions of the parameters $p_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,{\mathit{n}}_1$). If $\mathit{n}>{\mathit{n}}_1$, so that there are some driving equations, then the first $\mathit{n}-{\mathit{n}}_1$ values of g2 need not be set since they are never used.

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

range must evaluate the boundary points $a$ and $b$, each of which may depend on the arguments $p_1,p_2,\dots ,p_{{\mathit{n}}_1}$. The integrations in the shooting method are always from $a$ to $b$.

[a, b] = range(p)

Input Parameters

1:     $\mathrm{p}\left(:\right)$ – double array

The current estimate of the $\mathit{i}$th argument, $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$.

Output Parameters

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

$a$, one of the boundary points.

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

The second boundary point, $b$. Note that $\mathbf{b}>\mathbf{a}$ forces the direction of integration to be that of increasing $x$. If a and b are interchanged the direction of integration is reversed.

Optional Input Parameters

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

Default: the dimension of the arrays p, pe. (An error is raised if these dimensions are not equal.)

${\mathit{n}}_1$, the number of arguments.

Constraint: $1\le \mathbf{n1}\le \mathbf{n}$.

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

Default: the dimension of the array e.

$\mathit{n}$, the total number of differential equations.

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

Output Parameters

1:     $\mathrm{p}\left(\mathbf{n1}\right)$ – double array

The corrected value for the $i$th argument, unless an error has occurred, when it contains the last calculated value of the argument.

2:     $\mathrm{soln}\left(\mathbf{n},\mathbf{m1}\right)$ – double array

The solution when $\mathbf{m1}>1$.

3:     $\mathrm{w}\left(\mathbf{n},\mathit{sdw}\right)$ – double array

$\mathit{sdw}=3\mathbf{n}+14+\max \left(11,\mathbf{n}\right)$.

With $\mathbf{ifail}=\mathbf{2}$, $\mathbf{3}$, $\mathbf{4}$ or $\mathbf{5}$ (see Error Indicators and Warnings), $\mathbf{w}\left(\mathit{i},1\right)$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, contains the solution at the point $x$ when the error occurred. $\mathbf{w}\left(1,2\right)$ contains $x$.

4:     $\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$

One or more of the arguments n, n1, m1, sdw, e or pe is incorrectly set.

   $\mathbf{ifail}=2$

The step length for the integration became too short whilst calculating the residual (see Further Comments).

   $\mathbf{ifail}=3$

No initial step length could be chosen for the integration whilst calculating the residual.

   $\mathbf{ifail}=4$

As for $\mathbf{ifail}=\mathbf{2}$ but the error occurred when calculating the Jacobian.

   $\mathbf{ifail}=5$

As for $\mathbf{ifail}=\mathbf{3}$ but the error occurred when calculating the Jacobian.

   $\mathbf{ifail}=6$

The calculated Jacobian has an insignificant column. This can occur because a parameter $p_i$ is incorrectly entered when posing the problem.

   $\mathbf{ifail}=7$

The linear algebra function used (nag_lapack_dgesvd (f08kb)) has failed. This error exit should not occur and can be avoided by changing the initial estimates $p_i$.

   $\mathbf{ifail}=8$

The Newton iteration has failed to converge. This can indicate a poor initial choice of parameters $p_i$ or a very difficult problem. Consider varying the elements $\mathbf{pe}\left(i\right)$ if the residuals are small in the monitoring output. If the residuals are large, try varying the initial parameters $p_i$.

   $\mathbf{ifail}=9$

   $\mathbf{ifail}=10$

   $\mathbf{ifail}=11$

   $\mathbf{ifail}=12$

   $\mathbf{ifail}=13$

Indicates that a serious error has occurred. Check all array subscripts and function argument lists in the call to nag_ode_bvp_shoot_genpar (d02hb). Seek expert help.

   $\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: d02hb_example

function d02hb_example


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

% Set up initial values for first case.
n = int64(2);
n1 = int64(2);
p = [0.2;0];
pe = [1e-05;0.001];
e = [0.0001;0.0001];
m1 = int64(6);
marray = zeros(1,m1);

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

fprintf('\nCase 1 \n\n');

[pOut1, solnOut1, w1, ifail] = ...
  d02hb(...
         p, pe, e, m1, @fcn1, @bc1, @range1, 'n1', n1, 'n', n);

% Output results.
fprintf('Final parameters \n')
fprintf('    %1.3e   ',pOut1);
fprintf('\n\n Final solution \n')
fprintf('X-value     Components of solution \n')

[x0, x1] = range1(p);
h = (x1-x0)/double(m1-1);
for i = 1:m1;
  m = x0 + double(i-1)*h;
  fprintf('%7.3f', m);
  fprintf('%12.4f',solnOut1(1:n,i));
  fprintf('\n');
  % Save results for plotting.
  marray(i) = m;
end

fprintf('\nCase 2 \n\n');
% Set up initial values for second case.
n = int64(3);
n1 = n;
p = [32;6000;0.54];
pe = [1e-05;1e-4;1e-4];
e = [1e-2;1e-2;1e-2];
m1 = int64(6);
qarray = zeros(1,m1);

[pOut2, solnOut2, w2, ifail] = ...
  d02hb(...
         p, pe, e, int64(m1), @fcn2, @bc2, @range2, 'n1', n1, 'n', n);
% Output results.
fprintf('Final parameters \n')
fprintf('    %1.3e   ',pOut2);
fprintf('\n\n Final solution \n')
fprintf(' X-value        Components of solution \n')

[x0, x1] = range2(pOut2);
h = (x1-x0)/double(m1-1);
for i = 1:m1; 
  q = x0 + double(i-1)*h;
  fprintf('%6.0f', q);
  fprintf('%12.4f',solnOut2(1:n,i));
  fprintf('\n');
  % Save results for plotting.
  qarray(i) = q;
end

% Plot results.
fig1 = figure;
display_plot(marray, solnOut1, 'Derivative');
fig2 = figure;
display_plot(marray, solnOut2, 'Angle');


function [g1, g2] = bc1(p)
  % Evaluate boundary conditions.
  z=0.1;
  g1(1) = 0.1+p(1)*sqrt(z)*0.1+0.01*z;
  g1(2) = p(1)*0.05/sqrt(z) + 0.01;
  g2(1) = 1/6;
  g2(2) = p(2);

function [g1, g2] = bc2(p)
  % Evaluate boundary conditions.
  g1(1) = 0.0;
  g1(2) = 500.0;
  g1(3) = 0.5;
  g2(1) = 0.0;
  g2(2) = 450.0;
  g2(3) = p(3);

function f = fcn1(x, y, p)
  % Evaluate derivatives.
  f = zeros(2,1);
  f(1) = y(2);
  f(2) = (y(1)^3-y(2))/2/x;

function f = fcn2(x, y, p)
  % Evaluate derivatives.
  f = zeros(3,1);
  f(1) = tan(y(3));
  f(2) = -p(1)*tan(y(3))/y(2) - 0.00002*y(2)/cos(y(3));
  f(3) = -p(1)/y(2)^2;

function [x0, x1] = range1(p)
  % Evaluate boundary points.
  x0 = 0.1;
  x1 = 16.0;

function [x0, x1] = range2(p)
  % Evaluate boundary points.
  x0 = 0.0;
  x1 = p(2);

function display_plot(x, y, ylabelString)
  % Choose the type of plot based on the y axis label.
  if strncmp(ylabelString, 'Derivative', 8)
    % Plot the results as a linear graph (2 y axes).
    [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.1 0.18]);
    set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
    set(haxes(1), 'YTick', [0.1 0.12 0.14 0.16 0.18]);
    set(haxes(2), 'YLim', [0 0.02]);
    set(haxes(2), 'YMinorTick', 'on');
    set(haxes(2), 'YTick', [0:0.005:0.02]);
    for iaxis = 1:2
      % These properties must be the same for both sets of axes.
      set(haxes(iaxis), 'XLim', [0 16]);
      set(haxes(iaxis), 'XTick', [0:4:16]);
    end
    set(gca, 'box', 'off'); 
    % Set the title.
    title({['Parameterised Two-point BVP using'], ...
           ['Initial Value Techniques & Newton Iteration']}, ...
          'position',[8,0.17]);
    % Label the axes.
    xlabel('x');
    ylabel(haxes(1),'Solution');
    yh2 = ylabel(haxes(2),ylabelString);
    set(yh2,'position',[15,0.01]);
    % Add legend.
    legend('solution y(x)','derivative y''(x)','Location','South');
    % Set some features of the two lines.
    set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
    set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '--');
  else
    % Plot the results as a linear graph (2 y axes).
    [haxes, hline1, hline2] = plotyy(x, y(1,:), x, y(3,:));
    % Add a third curve,
    hold on
    hline3 = plot(x, y(2,:),'Color','Magenta');
    % Set the axis limits and the tick specifications to beautify the plot.
    set(haxes(1), 'YLim', [-1000 1000]);
    set(haxes(1), 'XMinorTick', 'on', 'YMinorTick', 'on');
    set(haxes(1), 'YTick', [-1000:500:1000]);
    set(haxes(2), 'YLim', [-1 1]);
    set(haxes(2), 'YMinorTick', 'on');
    set(haxes(2), 'YTick', [-1:0.5:1]);
    for iaxis = 1:2
      % These properties must be the same for both sets of axes.
      set(haxes(iaxis), 'XLim', [0 16]);
      set(haxes(iaxis), 'XTick', [0:4:16]);
    end
    set(gca, 'box', 'off');
    % Set the title.
    title('Parameterized projectile problem');
    % Label the axes.
    xlabel('x');
    ylabel(haxes(1),'Height and Velocity');
    ylabel(haxes(2),ylabelString);
    % Add legend.
    legend('height','velocity','angle','Location','South');
    % Set some features of the three lines.
    set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
    set(hline2, 'Linewidth', 0.5, 'Marker', '*', 'LineStyle', ':');
    set(hline3, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '--');
end

d02hb example results

d02hb example program results 

Case 1 

Final parameters 
    4.629e-02       3.494e-03   

 Final solution 
X-value     Components of solution 
  0.100      0.1025      0.0173
  3.280      0.1217      0.0042
  6.460      0.1338      0.0036
  9.640      0.1449      0.0034
 12.820      0.1557      0.0034
 16.000      0.1667      0.0035

Case 2 

Final parameters 
    3.239e+01       5.962e+03       -5.353e-01   

 Final solution 
 X-value        Components of solution 
     0      0.0000    500.0000      0.5000
  1192    529.6471    451.5554      0.3280
  2385    807.2140    420.3050      0.1231
  3577    820.3972    409.4254     -0.1032
  4769    556.1322    419.9890     -0.3296
  5962     -0.0003    450.0000     -0.5353