NAG Toolbox: nag_ode_bvp_shoot_genpar_intern (d02ag)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

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

h must be set to an estimate of the step size, $h$, needed for integration.

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

$\mathbf{e}\left(i\right)$ must be set to a small quantity to control the $i$th solution component. The element $\mathbf{e}\left(i\right)$ is used:

(i)	in the bound on the local error in the $i$th component of the solution $y_i$ during integration,
(ii)	in the convergence test on the $i$th component of the solution $y_i$ at the matching point in the Newton iteration.

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.

3:     $\mathrm{parerr}\left(\mathbf{n1}\right)$ – double array

$\mathbf{parerr}\left(i\right)$ must be set to a small quantity to control the $i$th parameter component. The element $\mathbf{parerr}\left(i\right)$ is used:

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

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

4:     $\mathrm{param}\left(\mathbf{n1}\right)$ – double array

$\mathbf{param}\left(\mathit{i}\right)$ must be set to an estimate for the $\mathit{i}$th parameter, $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathbf{n1}$.

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

Determines whether or not the final solution is computed as well as the parameter 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 c starting with the values of $y_i$ evaluated at the first end point (see raaux) stored in $\mathbf{c}\left(1,i\right)$.

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

aux must evaluate the functions $f_i$ (i.e., the derivatives $y_i^{\prime }$) for given values of its arguments, $x,y_1,\dots ,y_{\mathit{n}}$, $p_1,\dots ,p_{{\mathit{n}}_1}\text{.}$

[f] = aux(y, x, param)

Input Parameters

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

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

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

$x$, the value of the argument.

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

$p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$, the value of the parameters.

Output Parameters

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

The value of $f_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$.

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

bcaux must evaluate the values of $y_i$ at the end points of the range given the values of $p_1,\dots ,p_{{\mathit{n}}_1}$.

[g0, g1] = bcaux(param)

Input Parameters

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

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

Output Parameters

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

The values $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, at the boundary point $x_0$ (see raaux).

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

The values $y_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,\mathit{n}$, at the boundary point $x_1$ (see raaux).

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

raaux must evaluate the end points, $x_0$ and $x_1$, of the range and the matching point, $r$, given the values $p_1,p_2,\dots ,p_{{\mathit{n}}_1}$.

[x0, x1, r] = raaux(param)

Input Parameters

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

$p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$, the value of the parameters.

Output Parameters

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

Must contain the left-hand end of the range, $x_0$.

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

Must contain the right-hand end of the range $x_1$.

3:     $\mathrm{r}$ – double scalar

Must contain the matching point, $r$.

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

prsol is called at each iteration of the Newton method and can be used to print the current values of the parameters $p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$, their errors, $e_i$, and the sum of squares of the errors at the matching point, $r$.

prsol(param, res, n1, err)

Input Parameters

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

$p_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$, the current value of the parameters.

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

The sum of squares of the errors in the arguments, ${\displaystyle \sum _{i=1}^{{\mathit{n}}_1}}{e}_i^2$.

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

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

4:     $\mathrm{err}\left(\mathbf{n1}\right)$ – double array

The errors in the parameters, $e_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathit{n}}_1$.

Optional Input Parameters

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

Default: the dimension of the array e.

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

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

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

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

If $\mathbf{n1}<\mathbf{n}$, the last $\mathbf{n}-\mathbf{n1}$ differential equations (in aux) are driving equations (see Description).

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

Output Parameters

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

The last step length used.

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

The corrected value for the $i$th parameter, unless an error has occurred, when it contains the last calculated value of the parameter (possibly perturbed by $\mathbf{parerr}\left(i\right)\times \left(1+\left|\mathbf{param}\left(i\right)\right|\right)$ if the error occurred when calculating the approximate derivatives).

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

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

If $\mathbf{m1}=1$, the elements of c are not used.

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$

This indicates that $\mathbf{n1}>\mathbf{n}$ on entry, that is the number of parameters is greater than the number of differential equations.

   $\mathbf{ifail}=2$

As for $\mathbf{ifail}=\mathbf{4}$ except that the integration failed while calculating the matrix for use in the Newton iteration.

   $\mathbf{ifail}=3$

The current matching point $r$ does not lie between the current end points $x_0$ and $x_1$. If the values $x_0$, $x_1$ and $r$ depend on the parameters $p_i$, this may occur at any time in the Newton iteration if care is not taken to avoid it when coding raaux.

   $\mathbf{ifail}=4$

The step length for integration h has halved more than $13$ times (or too many steps were needed to reach the end of the range of integration) in attempting to control the local truncation error whilst integrating to obtain the solution corresponding to the current values $p_i$. If, on failure, h has the sign of $r-x_0$ then failure has occurred whilst integrating from $x_0$ to $r$, otherwise it has occurred whilst integrating from $x_1$ to $r$.

   $\mathbf{ifail}=5$

The matrix of the equations to be solved for corrections to the variable parameters in the Newton method is singular (as determined by nag_lapack_dgetrf (f07ad)).

   $\mathbf{ifail}=6$

A satisfactory correction to the parameters was not obtained on the last Newton iteration employed. A Newton iteration is deemed to be unsatisfactory if the sum of the squares of the residuals (which can be printed using prsol) has not been reduced after three iterations using a new Newton correction.

   $\mathbf{ifail}=7$

Convergence has not been obtained after $12$ satisfactory iterations of the Newton method.

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

function d02ag_example


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

% For communication with prsol.
global iprint;

% Set up initial values before calling NAG routine for the first time.
h = 0.1;
e = [0.0001; 0.0001];
parerr = [1e-05; 0.001];
param = [0.2; 0];
m1 = 6;
n = 2;
marray = zeros(1,m1);

% First example - parameterized two-point boundary-value problem.
fprintf('Case 1 \n\n');

% Set iprint to 1 to obtain output from prsol at each iteration.
iprint = 0;

[h, param, c, ifail] = ...
d02ag(...
       h, e, parerr, param, int64(m1), @aux1, @bcaux1, @raaux1, @prsol);

% Output results.
fprintf('Final parameters \n');
disp(sprintf('   %10.2e',param));
fprintf('\n Final solution \n');
fprintf('   X-value      Components of solution\n');

% Call the aux routine to evaluate the derivatives.
[x0, x1, r] = raaux1(param);

% h is steplength for which derivatives were evaluated.
h = (x1-x0)/double(m1-1);
for i = 1:m1;
  m = x0 + double(i-1)*h;
  fprintf('%8.2f    ',m);
  fprintf('%14.4e',c(i,1:n));
  fprintf('\n');
  % Store the m values so they can be plotted.
  marray(i) = m;
end

% Second example - find gravitational constant and range given
% projectile terminal velocity.
fprintf('\n\n Case 2 \n\n');

% Set iprint to 1 to obtain output from prsol at each iteration.
iprint = 0;

h = 10.0;
param = [32.0; 6000.0; 0.54];
parerr = [1.0e-5; 1.0e-4; 1.0e-4];
e = [1.0e-2; 1.0e-2; 1.0e-2];
n = 3;
m1 = 6;
qarray = zeros(1,m1);

[h, param, c1, ifail] = ...
d02ag( ...
       h, e, parerr, param, int64(m1), @aux2, @bcaux2, @raaux2, @prsol);

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

% Call the aux routine to evaluate the derivatives.
[x0, x1, r] = raaux2(param);

% h is steplength for which derivatives were evaluated.
h = (x1-x0)/double(m1-1);
for i = 1:m1;
  q = x0 + double(i-1)*h;
  fprintf('%8.0f    ',q);
  fprintf('%14.4e',c1(i,1:n));
  fprintf('\n');
  % Store the q values so they can be used for plotting.
  qarray(i) = q;
end

% Plot results.
fig1 = figure;
display_plot(marray, c, 'Solution and Derivative');
fig2 = figure;
display_plot(qarray, c1, 'Height and Velocity');


function [f] = aux1(y,x,param)
  % Evaluate the derivatives (case 1).
  f(1) = y(2);
  f(2) = (y(1)^3-y(2))/(2.0*x);

function [g0,g1] = bcaux1(param)
  % Evaluate the derivatives at the end-points (case 1).
  z = 0.1;
  g0(1) = 0.1 + param(1)*sqrt(z)*0.1 + 0.01*z;
  g0(2) = param(1)*0.05/sqrt(z) + 0.01;
  g1(1) = 1.0/6.0;
  g1(2) = param(2);

function [] = prsol(param, resid, n, err)
  % For communication with main routine.
  global iprint;
  if iprint == 1
    fprintf('Current parameters: ');
    fprintf('%15.6e ', param(1:n));
    fprintf('\nResiduals: ');
    fprintf('%13.6e ', err(1:n));
    fprintf('\nSum of residuals squared = %13.6e\n', resid);
  end

function [x0, x1, r] = raaux1(param)
  % Evaluate the end-points (case 1).
  x0 = 0.1;
  x1 = 16;
  r  = 16;

function [f] = aux2(y,x,param)
  % Evaluate the derivatives (case 2).
  c = cos(y(3));
  s = sin(y(3));
  f(1) = s/c;
  f(2) = -(param(1)*s+0.00002*y(2)*y(2))/(y(2)*c);
  f(3) = -param(1)/(y(2)*y(2));
  
function [g0,g1] = bcaux2(param)
  % Evaluate the derivatives at the end-points (case 2).
  g0(1) = 0;
  g0(2) = 500;
  g0(3) = 0.5;
  g1(1) = 0;
  g1(2) = 450;
  g1(3) = param(3);

function [x0, x1, r] = raaux2(param)
  % Evaluate the end-points (case 2).
  x0 = double(0);
  x1 = param(2);
  r  = param(2);

function display_plot(data, index, ylabelString)
  % Choose the type of plot based on y axis label.
  if strncmp(ylabelString, 'Solution', 8)
    % Plot the results as a linear graph.
    plot(data,index(:,1),'-+',data,index(:,2),'--x');
    % Set the axis limits and the tick specifications to beautify the plot.
    set(gca, 'XLim', [0 16]);
    set(gca, 'XTick', [0 4 8 12 16]);
    % Set the title.
    title('Parameterized Two-point Boundary-value Problem');
    % Label the axes.
    xlabel('x');
    ylabel(ylabelString);
    % Add a legend.
    legend('solution y(x)','derivative y''(x)','Location','Best');
  else
    % Plot the results as a linear graph (2 y axes).
    [haxes, hline1, hline2] = plotyy(data,index(:,1),data,index(:,3));
    hold on
    % Add a third curve.
    hline3 = plot(data,index(:,2));
    % Set the axis limits and the tick specifications to beautify the plot.
    set(haxes(1), 'YLim', [0 850]);
    set(haxes(1), 'YMinorTick', 'on');
    set(haxes(1), 'YTick', [0 200 400 600 800]);
    set(haxes(2), 'YLim', [-1 1]);
    set(haxes(2), 'YMinorTick', 'on');
    set(haxes(2), 'YTick', [-1 -0.5 0 0.5 1]);
    for iaxis = 1:2
        % These properties must be the same for both sets of axes.
        set(haxes(iaxis), 'XLim', [0 6000]);
        set(haxes(iaxis), 'XTick', [0 1000 2000 3000 4000 5000 6000]);
    end
    set(gca, 'box', 'off'); % no ticks on opposite axes.
    % Set the title.
    title(['Range from Projectile Terminal Velocity']);
    % Label the axes.
    xlabel('x');
    ylabel(haxes(1),ylabelString);
    ylabel(haxes(2),'Angle');
    % Add a legend.
    legend('angle','height','velocity','Location','Best');
    % 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', '--', ...
        'Color', 'Magenta');
  end

d02ag example results

Case 1 

Final parameters 
     4.64e-02     3.49e-03

 Final solution 
   X-value      Components of solution
    0.10        1.0247e-01    1.7341e-02
    3.28        1.2170e-01    4.1796e-03
    6.46        1.3382e-01    3.5764e-03
    9.64        1.4488e-01    3.4178e-03
   12.82        1.5572e-01    3.4142e-03
   16.00        1.6667e-01    3.4943e-03


 Case 2 

Final parameters
     3.24e+01     5.96e+03    -5.35e-01

Final solution
   X-value      Components of solution
       0        0.0000e+00    5.0000e+02    5.0000e-01
    1193        5.2981e+02    4.5156e+02    3.2807e-01
    2385        8.0765e+02    4.2030e+02    1.2315e-01
    3578        8.2080e+02    4.0944e+02   -1.0316e-01
    4771        5.5625e+02    4.2001e+02   -3.2958e-01
    5963        0.0000e+00    4.5000e+02   -5.3523e-01