hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_ode_bvp_shoot_genpar_intern (d02ag)

Purpose

nag_ode_bvp_shoot_genpar_intern (d02ag) solves a two-point boundary value problem for a system of ordinary differential equations, using initial value techniques and Newton iteration; it generalizes nag_ode_bvp_shoot_bval (d02ha) to include the case where parameters other than boundary values are to be determined.

Syntax

[h, param, c, ifail] = d02ag(h, e, parerr, param, m1, aux, bcaux, raaux, prsol, 'n', n, 'n1', n1)
[h, param, c, ifail] = nag_ode_bvp_shoot_genpar_intern(h, e, parerr, param, m1, aux, bcaux, raaux, prsol, 'n', n, 'n1', n1)

Description

nag_ode_bvp_shoot_genpar_intern (d02ag) solves a two-point boundary value problem by determining the unknown parameters p1,p2,,pn1p1,p2,,pn1 of the problem. These parameters may be, but need not be, boundary values (as they are in nag_ode_bvp_shoot_bval (d02ha)); they may include eigenvalue parameters in the coefficients of the differential equations, length of the range of integration, etc. The notation and methods used are similar to those of nag_ode_bvp_shoot_bval (d02ha) and you are advised to study this first. (There the parameters p1,p2,,pn1p1,p2,,pn1 correspond to the unknown boundary conditions.) It is assumed that we have a system of nn first-order ordinary differential equations of the form
(dyi)/(dx) = fi(x,y1,y2,,yn),  i = 1,2,,n,
dyi dx =fi(x,y1,y2,,yn),  i=1,2,,n,
and that derivatives fifi are evaluated by aux. The system, including the boundary conditions given by bcaux, and the range of integration and matching point, rr, given by raaux, involves the n1n1 unknown parameters p1,p2,,pn1p1,p2,,pn1 which are to be determined, and for which initial estimates must be supplied. The number of unknown parameters n1n1 must not exceed the number of equations nn. If n1 < nn1<n, we assume that (nn1)(n-n1) equations of the system are not involved in the matching process. These are usually referred to as ‘driving equations’; they are independent of the parameters and of the solutions of the other n1n1 equations. In numbering the equations for aux, the driving equations must be put last.
The estimated values of the parameters are corrected by a form of Newton iteration. The Newton correction on each iteration is calculated using a matrix whose (i,j)(i,j)th element depends on the derivative of the iith component of the solution, yiyi, with respect to the jjth parameter, pjpj. This matrix is calculated by a simple numerical differentiation technique which requires n1n1 evaluations of the differential system.

References

None.

Parameters

You are strongly recommended to read Sections [Description] and [Further Comments] in conjunction with this section.

Compulsory Input Parameters

1:     h – double scalar
h must be set to an estimate of the step size, hh, needed for integration.
2:     e(n) – double array
e(i)ei must be set to a small quantity to control the iith solution component. The element e(i)ei is used:
(i) in the bound on the local error in the iith component of the solution yiyi during integration,
(ii) in the convergence test on the iith component of the solution yiyi at the matching point in the Newton iteration.
The elements e(i)ei should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.
3:     parerr(n1) – double array
n1, the dimension of the array, must satisfy the constraint n1nn1n.
parerr(i)parerri must be set to a small quantity to control the iith parameter component. The element parerr(i)parerri is used:
(i) in the convergence test on the iith parameter in the Newton iteration,
(ii) in perturbing the iith parameter when approximating the derivatives of the components of the solution with respect to the iith parameter, for use in the Newton iteration.
The elements parerr(i)parerri should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.
4:     param(n1) – double array
n1, the dimension of the array, must satisfy the constraint n1nn1n.
param(i)parami must be set to an estimate for the iith parameter, pipi, for i = 1,2,,n1i=1,2,,n1.
5:     m1 – int64int32nag_int scalar
Determines whether or not the final solution is computed as well as the parameter values.
m1 = 1m1=1
The final solution is not calculated;
m1 > 1m1>1
The final values of the solution at interval (length of range)/(m11)(m1-1) are calculated and stored sequentially in the array c starting with the values of yiyi evaluated at the first end point (see raaux) stored in c(1,i)c1i.
6:     aux – function handle or string containing name of m-file
aux must evaluate the functions fifi (i.e., the derivatives yiyi) for given values of its arguments, x,y1,,ynx,y1,,yn, p1,,pn1.p1,,pn1.
[f] = aux(y, x, param)

Input Parameters

1:     y( : :) – double array
yiyi, for i = 1,2,,ni=1,2,,n, the value of the argument.
2:     x – double scalar
xx, the value of the argument.
3:     param( : :) – double array
pipi, for i = 1,2,,n1i=1,2,,n1, the value of the parameters.

Output Parameters

1:     f( : :) – double array
The value of fifi, for i = 1,2,,ni=1,2,,n.
7:     bcaux – function handle or string containing name of m-file
bcaux must evaluate the values of yiyi at the end points of the range given the values of p1,,pn1p1,,pn1.
[g0, g1] = bcaux(param)

Input Parameters

1:     param( : :) – double array
pipi, for i = 1,2,,ni=1,2,,n, the value of the parameters.

Output Parameters

1:     g0( : :) – double array
The values yiyi, for i = 1,2,,ni=1,2,,n, at the boundary point x0x0 (see raaux).
2:     g1( : :) – double array
The values yiyi, for i = 1,2,,ni=1,2,,n, at the boundary point x1x1 (see raaux).
8:     raaux – function handle or string containing name of m-file
raaux must evaluate the end points, x0x0 and x1x1, of the range and the matching point, rr, given the values p1,p2,,pn1p1,p2,,pn1.
[x0, x1, r] = raaux(param)

Input Parameters

1:     param( : :) – double array
pipi, for i = 1,2,,n1i=1,2,,n1, the value of the parameters.

Output Parameters

1:     x0 – double scalar
Must contain the left-hand end of the range, x0x0.
2:     x1 – double scalar
Must contain the right-hand end of the range x1x1.
3:     r – double scalar
Must contain the matching point, rr.
9:     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 pipi, for i = 1,2,,n1i=1,2,,n1, their errors, eiei, and the sum of squares of the errors at the matching point, rr.
prsol(param, res, n1, err)

Input Parameters

1:     param(n1) – double array
pipi, for i = 1,2,,n1i=1,2,,n1, the current value of the parameters.
2:     res – double scalar
The sum of squares of the errors in the parameters, i = 1n1ei2i=1n1ei2.
3:     n1 – int64int32nag_int scalar
n1n1, the number of parameters.
4:     err(n1) – double array
The errors in the parameters, eiei, for i = 1,2,,n1i=1,2,,n1.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the array e.
nn, the total number of differential equations.
2:     n1 – int64int32nag_int scalar
Default: The dimension of the arrays parerr, param. (An error is raised if these dimensions are not equal.)
n1n1, the number of parameters.
If n1 < nn1<n, the last nn1n-n1 differential equations (in aux) are driving equations (see Section [Description]).
Constraint: n1nn1n.

Input Parameters Omitted from the MATLAB Interface

mat copy wspace wspac1 wspac2

Output Parameters

1:     h – double scalar
The last step length used.
2:     param(n1) – double array
The corrected value for the iith parameter, unless an error has occurred, when it contains the last calculated value of the parameter (possibly perturbed by parerr(i) × (1 + |param(i)|)parerri×(1+|parami|) if the error occurred when calculating the approximate derivatives).
3:     c(m1,n) – double array
The solution when m1 > 1m1>1 (see m1).
If m1 = 1m1=1, the elements of c are not used.
4:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
This indicates that n1 > nn1>n on entry, that is the number of parameters is greater than the number of differential equations.
  ifail = 2ifail=2
As for ifail = 4ifail=4 except that the integration failed while calculating the matrix for use in the Newton iteration.
  ifail = 3ifail=3
The current matching point rr does not lie between the current end points x0x0 and x1x1. If the values x0x0, x1x1 and rr depend on the parameters pipi, this may occur at any time in the Newton iteration if care is not taken to avoid it when coding raaux.
  ifail = 4ifail=4
The step length for integration h has halved more than 1313 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 pipi. If, on failure, h has the sign of rx0r-x0 then failure has occurred whilst integrating from x0x0 to rr, otherwise it has occurred whilst integrating from x1x1 to rr.
  ifail = 5ifail=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)).
  ifail = 6ifail=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.
  ifail = 7ifail=7
Convergence has not been obtained after 1212 satisfactory iterations of the Newton method.
A further discussion of these errors and the steps which might be taken to correct them is given in Section [Further Comments].

Accuracy

If the process converges, the accuracy to which the unknown parameters are determined is usually close to that specified by you; and the solution, if requested, is usually determined to the accuracy specified.

Further Comments

The time taken by nag_ode_bvp_shoot_genpar_intern (d02ag) depends on the complexity of the system, and on the number of iterations required. In practice, integration of the differential equations is by far the most costly process involved.
There may be particular difficulty in integrating the differential equations in one direction (indicated by ifail = 2ifail=2 or 44). The value of rr should be adjusted to avoid such difficulties.
If the matching point rr is at one of the end points x0x0 or x1x1 and some of the parameters are used only to determine the boundary values at this point, then good initial estimates for these parameters are not required, since they are completely determined by the function (for example, see p2p2 in EX1 of Section [Example]).
Wherever they occur in the procedure, the error parameters contained in the arrays e and parerr are used in ‘mixed’ form; that is e(i)ei always occurs in expressions of the form e(i) × (1 + |yi|)ei×(1+|yi|), and parerr(i)parerri always occurs in expressions of the form parerr(i) × (1 + |pi|)parerri×(1+|pi|). Though not ideal for every application, it is expected that this mixture of absolute and relative error testing will be adequate for most purposes.
Note that convergence is not guaranteed. You are strongly advised to provide an output prsol, as shown in EX1 of Section [Example], in order to monitor the progress of the iteration. Failure of the Newton iteration to converge (see ifail = 6ifail=6 or 77) usually results from poor starting approximations to the parameters, though occasionally such failures occur because the elements of one or both of the arrays parerr or e are too small. (It should be possible to distinguish these cases by studying the output from prsol.) Poor starting approximations can also result in the failure described under ifail = 4ifail=4 and 55 in Section [Error Indicators and Warnings] (especially if these errors occur after some Newton iterations have been completed, that is, after two or more calls of prsol). More frequently, a singular matrix in the Newton method (monitored as ifail = 5ifail=5) occurs because the mathematical problem has been posed incorrectly. The case ifail = 4ifail=4 usually occurs because hh or rr has been poorly estimated, so these values should be checked first. If ifail = 2ifail=2 is monitored, the solution y1,y2,,yny1,y2,,yn is sensitive to perturbations in the parameters pipi. Reduce the size of one or more values parerr(i)parerri to reduce the perturbations. Since only one value pipi is perturbed at any time when forming the matrix, the perturbation which is too large can be located by studying the final output from prsol and the values of the parameters returned by nag_ode_bvp_shoot_genpar_intern (d02ag). If this change leads to other types of failure improve the initial values of pipi by other means.
The computing time for integrating the differential equations can sometimes depend critically on the quality of the initial estimates for the parameters pipi. If it seems that too much computing time is required and, in particular, if the values err(i)erri (available on each call of prsol) are much larger than the expected values of the solution at the matching point rr, then the coding of aux, bcaux and raaux should be checked for errors. If no errors can be found, an independent attempt should be made to improve the initial estimates for param(i)parami.
The function can be used to solve a very wide range of problems, for example:
(a) eigenvalue problems, including problems where the eigenvalue occurs in the boundary conditions;
(b) problems where the differential equations depend on some parameters which are to be determined so as to satisfy certain boundary conditions (see EX1 in Section [Example]);
(c) problems where one of the end points of the range of integration is to be determined as the point where a variable yiyi takes a particular value (see EX2 in Section [Example]);
(d) singular problems and problems on infinite ranges of integration where the values of the solution at x0x0 or x1x1 or both are determined by a power series or an asymptotic expansion (or a more complicated expression) and where some of the coefficients in the expression are to be determined (see EX1 in Section [Example]); and
(e) differential equations with certain terms defined by other independent (driving) differential equations.

Example

For this function two examples are presented. There is a single example program for nag_ode_bvp_shoot_genpar_intern (d02ag), with a main program and the code to solve the two example problems given in Example 1 (EX1) and Example 2 (EX2).
Example 1 (EX1)
This example finds the solution of the differential equation
y = (y3y)/(2x)
y=y3-y 2x
on the range 0x160x16, with boundary conditions y(0) = 0.1y(0)=0.1 and y(16) = 1 / 6y(16)=1/6.
We cannot use the differential equation at x = 0x=0 because it is singular, so we take the truncated series expansion
y(x) = (1/10) + p1(sqrt(x))/10 + x/100
y(x)=110+p1x10+x100
near the origin (which is correct to the number of terms given in this case). Here p1p1 is one of the parameters to be determined. We choose the range as [0.1,16][0.1,16] and setting p2 = y(16)p2=y(16), we can determine all the boundary conditions. We take the matching point to be 1616, the end of the range, and so a good initial guess for p2p2 is not necessary. We write y = y(1)y=y1, y = y(2)y=y2, and estimate p1 = param(1) = 0.2p1=param1=0.2, p2 = param(2) = 0.0p2=param2=0.0.
Example 2 (EX2)
This example finds the gravitational constant p1p1 and the range p2p2 over which a projectile must be fired to hit the target with a given velocity. The differential equations are
y = tanφ
v = ((p1sinφ + 0.00002v2))/(vcosφ)
φ = (p1)/(v2)k
y=tanϕ v= -(p1sinϕ+0.00002v2) vcosϕ ϕ=-p1v2k
on the range 0 < x < p20<x<p2 with boundary conditions
y = 0, v = 500, φ = 0.5 at x = 0
y = 0, v = 450, φ = p3 at x = p2.
y=0, v=500, ϕ=0.5 at x=0 y=0, v=450, ϕ=p3 at x=p2.
We write y = y(1)y=y1, v = y(2)v=y2, φ = y(3)ϕ=y3, and we take the matching point r = p2r=p2. We estimate p1 = param(1) = 32p1=param1=32, p2 = param(2) = 6000p2=param2=6000 and p3 = param(3) = 0.54p3=param3=0.54 (though this estimate is not important).
function nag_ode_bvp_shoot_genpar_intern_example
% 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);

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

% 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;

[hOut, paramOut, c, ifail] = ...
    nag_ode_bvp_shoot_genpar_intern(h, e, parerr, param, int64(m1), @aux1, @bcaux1, ...
    @raaux1, @prsol);
if (ifail ~= 0)
    % Problems in integration.  Print message and exit.
    error('Warning: nag_ode_bvp_shoot_genpar_intern returned with ifail = %1d ',ifail);
else
    % Output results.
    fprintf('Final parameters \n');
    disp(sprintf('   %2.6e   ',paramOut));
    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; j = 1:n;
        m = x0 + double(i-1)*h;
        fprintf('%8.2f   ',m);
        disp(sprintf('   %2.4e  ',c(i,j)));
        % Store the m values so they can be plotted.
        marray(i) = m;
    end
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);

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

if ifail ~= 0
    % Problems in integration.  Print message and exit.
    error('Warning: nag_ode_bvp_shoot_genpar_intern returned with ifail = %1d ',ifail);
else
    % Output results.
    fprintf('Final parameters\n');
    disp(sprintf('   %10.6e   ',paramOut1));
    fprintf('\nFinal solution\n');
    fprintf(' X-value     Components of solution\n');

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

    % h is steplength for which derivatives were evaluated.
    h = (x1-x0)/double(m1-1);
    for i = 1:m1; j = 1:n;
        q = x0 + double(i-1)*h;
        fprintf('%6.0f   ',q);
        disp(sprintf('%14.4e  ',c1(i,j)));
        % Store the q values so they can be used for plotting.
        qarray(i) = q;
    end
end
% Plot results.
fig = figure('Number', 'off');
display_plot(marray, c, 'Solution and Derivative');

fig = figure('Number', 'off');
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: ');
    for k=1:int32(n) % can't use int64 in loop range
        fprintf('%15.6e ', param(k));
    end
    fprintf('\nResiduals: ');
    for k=1:int32(n) % can't use int64 in loop range
        fprintf('%13.6e ', err(k));
    end
    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.0;
r = 16.0;
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.0;
g0(2) = 500.0;
g0(3) = 0.5;
g1(1) = 0.0;
g1(2) = 450.0;
g1(3) = param(3);
function [x0, x1, r] = raaux2(param)
% Evaluate the end-points (case 2).
x0 = double(0.0);
x1 = param(2);
r = param(2);
function display_plot(data, index, ylabelString)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% 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', titleFmt{:});
    % Label the axes.
    xlabel('x', labFmt{:});
    ylabel(ylabelString, labFmt{:});
    % 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), 'FontSize', 13);
        set(haxes(iaxis), 'XTick', [0 1000 2000 3000 4000 5000 6000]);
    end
    set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
    % Set the title.
    title(['Find Gravitational Constant & Range from Projectile ', ...
        'Terminal Velocity'], titleFmt{:});
    % Label the axes.
    xlabel('x', labFmt{:});
    ylabel(haxes(1),ylabelString, labFmt{:});
    ylabel(haxes(2),'Angle', labFmt{:});
    % 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', '--');
end
 
nag_ode_bvp_shoot_genpar_intern example program results 

Case 1 

Final parameters 
   4.642688e-02      3.494287e-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.237291e+01      5.963172e+03      -5.352306e-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  

function d02ag_example
% 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);

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

% 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;

[hOut, paramOut, c, ifail] = ...
    d02ag(h, e, parerr, param, int64(m1), @aux1, @bcaux1, ...
    @raaux1, @prsol);
if (ifail ~= 0)
    % Problems in integration.  Print message and exit.
    error('Warning: d02ag returned with ifail = %1d ',ifail);
else
    % Output results.
    fprintf('Final parameters \n');
    disp(sprintf('   %2.6e   ',paramOut));
    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; j = 1:n;
        m = x0 + double(i-1)*h;
        fprintf('%8.2f   ',m);
        disp(sprintf('   %2.4e  ',c(i,j)));
        % Store the m values so they can be plotted.
        marray(i) = m;
    end
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);

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

if ifail ~= 0
    % Problems in integration.  Print message and exit.
    error('Warning: d02ag returned with ifail = %1d',ifail);
else
    % Output results.
    fprintf('Final parameters\n');
    disp(sprintf('   %10.6e   ',paramOut1));
    fprintf('\nFinal solution\n');
    fprintf(' X-value     Components of solution\n');

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

    % h is steplength for which derivatives were evaluated.
    h = (x1-x0)/double(m1-1);
    for i = 1:m1; j = 1:n;
        q = x0 + double(i-1)*h;
        fprintf('%6.0f   ',q);
        disp(sprintf('%14.4e  ',c1(i,j)));
        % Store the q values so they can be used for plotting.
        qarray(i) = q;
    end
end
% Plot results.
fig = figure('Number', 'off');
display_plot(marray, c, 'Solution and Derivative');

fig = figure('Number', 'off');
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: ');
    for k=1:int32(n) % can't use int64 in loop range
        fprintf('%15.6e ', param(k));
    end
    fprintf('\nResiduals: ');
    for k=1:int32(n) % can't use int64 in loop range
        fprintf('%13.6e ', err(k));
    end
    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.0;
r = 16.0;
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.0;
g0(2) = 500.0;
g0(3) = 0.5;
g1(1) = 0.0;
g1(2) = 450.0;
g1(3) = param(3);
function [x0, x1, r] = raaux2(param)
% Evaluate the end-points (case 2).
x0 = double(0.0);
x1 = param(2);
r = param(2);
function display_plot(data, index, ylabelString)
% Formatting for title and axis labels.
titleFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 14};
labFmt = {'FontName', 'Helvetica', 'FontWeight', 'Bold', 'FontSize', 13};
set(gca, 'FontSize', 13); % for legend, axis tick labels, etc.
% 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', titleFmt{:});
    % Label the axes.
    xlabel('x', labFmt{:});
    ylabel(ylabelString, labFmt{:});
    % 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), 'FontSize', 13);
        set(haxes(iaxis), 'XTick', [0 1000 2000 3000 4000 5000 6000]);
    end
    set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
    % Set the title.
    title(['Find Gravitational Constant & Range from Projectile ', ...
        'Terminal Velocity'], titleFmt{:});
    % Label the axes.
    xlabel('x', labFmt{:});
    ylabel(haxes(1),ylabelString, labFmt{:});
    ylabel(haxes(2),'Angle', labFmt{:});
    % 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', '--');
end
 
d02ag example program results 

Case 1 

Final parameters 
   4.642688e-02      3.494287e-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.237291e+01      5.963172e+03      -5.352306e-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  


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013