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 (d02hb)

Purpose

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

Syntax

[p, soln, w, ifail] = d02hb(p, pe, e, m1, fcn, bc, range, 'n1', n1, 'n', n)
[p, soln, w, ifail] = nag_ode_bvp_shoot_genpar(p, pe, e, m1, fcn, bc, range, 'n1', n1, 'n', n)

Description

nag_ode_bvp_shoot_genpar (d02hb) 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; 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. (The parameters p1,p2,,pn1p1,p2,,pn1 correspond precisely to the unknown boundary conditions in nag_ode_bvp_shoot_bval (d02ha).) 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 the derivatives fifi are evaluated by fcn. The system, including the boundary conditions given by bc and the range of integration given by range, 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 fcn, the driving equations must be put first.
The estimated values of the parameters are corrected by a form of Newton iteration. The Newton correction on each iteration is calculated using a Jacobian 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.
If the parameter ifail is set appropriately, the function automatically prints messages to inform you of the flow of the calculation. These messages are discussed in detail in Section [Further Comments].
nag_ode_bvp_shoot_genpar (d02hb) is a simplified version of nag_ode_bvp_shoot_genpar_algeq (d02sa) which is described in detail in Gladwell (1979).

References

Gladwell I (1979) The development of the boundary value codes in the ordinary differential equations chapter of the NAG Library Codes for Boundary Value Problems in Ordinary Differential Equations. Lecture Notes in Computer Science (eds B Childs, M Scott, J W Daniel, E Denman and P Nelson) 76 Springer–Verlag

Parameters

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

Compulsory Input Parameters

1:     p(n1) – double array
n1, the dimension of the array, must satisfy the constraint 1n1n1n1n.
An estimate for the iith parameter, pipi, for i = 1,2,,n1i=1,2,,n1.
2:     pe(n1) – double array
n1, the dimension of the array, must satisfy the constraint 1n1n1n1n.
The elements of pe must be given small positive values. The element pe(i)pei is used
(i) in the convergence test on the iith parameter in the Newton iteration, and
(ii) in perturbing the iith parameter when approximating the derivatives of the components of the solution with respect to this parameter for use in the Newton iteration.
The elements pe(i)pei should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.
Constraint: pe(i) > 0.0pei>0.0, for i = 1,2,,n1i=1,2,,n1.
3:     e(n) – double array
n, the dimension of the array, must satisfy the constraint n2n2.
The elements of e must be given positive values. The element e(i)ei is used in the bound on the local error in the iith component of the solution yiyi during integration.
The elements e(i)ei should not be chosen too small. They should usually be several orders of magnitude larger than machine precision.
Constraint: e(i) > 0.0ei>0.0, for i = 1,2,,ni=1,2,,n.
4:     m1 – int64int32nag_int scalar
A value which controls exit 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 soln starting with the values of the solutions evaluated at the first end point (see range) stored in the first column of soln.
Constraint: m11m11.
5:     fcn – function handle or string containing name of m-file
fcn must evaluate the functions fifi (i.e., the derivatives yiyi), for i = 1,2,,ni=1,2,,n, at a general point xx.
[f] = fcn(x, y, p)

Input Parameters

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

Output Parameters

1:     f( : :) – double array
The value of fifi, for i = 1,2,,ni=1,2,,n. The fifi may depend upon the parameters pjpj, for j = 1,2,,n1j=1,2,,n1. If there are any driving equations (see Section [Description]) then these must be numbered first in the ordering of the components of f in fcn.
6:     bc – function handle or string containing name of m-file
bc must place in g1 and g2 the boundary conditions at aa and bb respectively (see range).
[g1, g2] = bc(p)

Input Parameters

1:     p( : :) – double array
An estimate of the parameter pipi, for i = 1,2,,n1i=1,2,,n1.

Output Parameters

1:     g1( : :) – double array
The value of yi(a)yi(a), (where this may be a known value or a function of the parameters pjpj, for i = 1,2,,ni=1,2,,n and j = 1,2,,n1j=1,2,,n1).
2:     g2( : :) – double array
The value of yi(b)yi(b), for i = 1,2,,ni=1,2,,n, (where these may be known values or functions of the parameters pjpj, for j = 1,2,,n1j=1,2,,n1). If n > n1n>n1, so that there are some driving equations, then the first nn1n-n1 values of g2 need not be set since they are never used.
7:     range – function handle or string containing name of m-file
range must evaluate the boundary points aa and bb, each of which may depend on the parameters p1,p2,,pn1p1,p2,,pn1. The integrations in the shooting method are always from aa to bb.
[a, b] = range(p)

Input Parameters

1:     p( : :) – double array
The current estimate of the iith parameter, pipi, for i = 1,2,,n1i=1,2,,n1.

Output Parameters

1:     a – double scalar
aa, one of the boundary points.
2:     b – double scalar
The second boundary point, bb. Note that b > ab>a forces the direction of integration to be that of increasing xx. If a and b are interchanged the direction of integration is reversed.

Optional Input Parameters

1:     n1 – int64int32nag_int scalar
Default: The dimension of the arrays p, pe. (An error is raised if these dimensions are not equal.)
n1n1, the number of parameters.
Constraint: 1n1n1n1n.
2:     n – int64int32nag_int scalar
Default: The dimension of the array e.
nn, the total number of differential equations.
Constraint: n2n2.

Input Parameters Omitted from the MATLAB Interface

sdw

Output Parameters

1:     p(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.
2:     soln(n,m1) – double array
The solution when m1 > 1m1>1.
3:     w(n,sdw) – double array
sdw3n + 14 + max (11,n)sdw3n+14+max(11,n).
With ifail = 2ifail=2, 33, 44 or 55 (see Section [Error Indicators and Warnings]), w(i,1)wi1, for i = 1,2,,ni=1,2,,n, contains the solution at the point xx when the error occurred. w(1,2)w12 contains xx.
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
One or more of the parameters n, n1, m1, sdw, e or pe is incorrectly set.
  ifail = 2ifail=2
The step length for the integration became too short whilst calculating the residual (see Section [Further Comments]).
  ifail = 3ifail=3
No initial step length could be chosen for the integration whilst calculating the residual.
Note: ifail = 2ifail=2 or 33 can occur due to choosing too small a value for e or due to choosing the wrong direction of integration. Try varying e and interchanging aa and bb. These error exits can also occur for very poor initial choices of the parameters in the array p and, in extreme cases, because nag_ode_bvp_shoot_genpar (d02hb) cannot be used to solve the problem posed.
  ifail = 4ifail=4
As for ifail = 2ifail=2 but the error occurred when calculating the Jacobian.
  ifail = 5ifail=5
As for ifail = 3ifail=3 but the error occurred when calculating the Jacobian.
  ifail = 6ifail=6
The calculated Jacobian has an insignificant column. This can occur because a parameter pipi is incorrectly entered when posing the problem.
Note: ifail = 4ifail=4, 55 or 66 usually indicate a badly scaled problem. You may vary the size of pe. Otherwise the use of the more general nag_ode_bvp_shoot_genpar_algeq (d02sa) which affords more control over the calculations is advised.
  ifail = 7ifail=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 pipi.
  ifail = 8ifail=8
The Newton iteration has failed to converge. This can indicate a poor initial choice of parameters pipi or a very difficult problem. Consider varying the elements pe(i)pei if the residuals are small in the monitoring output. If the residuals are large, try varying the initial parameters pipi.
  ifail = 9ifail=9
  ifail = 10ifail=10
  ifail = 11ifail=11
  ifail = 12ifail=12
  ifail = 13ifail=13
Indicates that a serious error has occurred. Check all array subscripts and function parameter lists in the call to nag_ode_bvp_shoot_genpar (d02hb). Seek expert help.

Accuracy

If the process converges, the accuracy to which the unknown parameters are determined is usually close to that specified by you; the solution, if requested, may be determined to a required accuracy by varying e.

Further Comments

The time taken by nag_ode_bvp_shoot_genpar (d02hb) 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.
Wherever they occur in the function, the error parameters contained in the arrays e and pe 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 pe(i)pei always occurs in expressions of the form
pe(i) × (1 + |pi|).
pei×(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.
You may determine a suitable direction of integration aa to bb and suitable values for e(i)ei by integrations with nag_ode_ivp_rkts_range (d02pe). The best direction of integration is usually the direction of decreasing solutions. You are strongly recommended to set ifail to obtain self-explanatory error messages, and also monitoring information about the course of the computation. You may select the channel numbers on which this output is to appear by calls of nag_file_set_unit_error (x04aa) (for error messages) or nag_file_set_unit_advisory (x04ab) (for monitoring information) – see Section [Example] for an example. Otherwise the default channel numbers will be used. The monitoring information produced at each iteration includes the current parameter values, the residuals and 22-norms: a basic norm and a current norm. At each iteration the aim is to find parameter values which make the current norm less than the basic norm. Both these norms should tend to zero as should the residuals. (They would all be zero if the exact parameters were used as input.) For more details, in particular about the other monitoring information printed, you are advised to consult the specification of nag_ode_bvp_shoot_genpar_algeq (d02sa), and especially the description of the parameter monit there.
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 of the residuals printed by the monitoring function are much larger than the expected values of the solution at bb, then the coding of fcn, bc and range should be checked for errors. If no errors can be found, an independent attempt should be made to improve the initial estimates for pipi.
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 Example 2 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 Example 2 in Section [Example]);
(d) singular problems and problems on infinite ranges of integration where the values of the solution at aa or bb 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 Example 1 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 (d02hb), 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 a truncated power series expansion
y(x) = 1 / 10 + p1 × sqrt(x) / 10 + x / 100
y(x)=1/10+p1×x/10+x/100
near the origin where p1p1 is one of the parameters to be determined. We choose the interval as [0.1,16][0.1,16] and setting p2 = y(16)p2=y(16), we can determine all the boundary conditions. We take X1 = 16X1=16. We write y = y(1)y=y1, y = y(2)y=y2, and estimate PARAM(1) = 0.2PARAM(1)=0.2, PARAM(2) = 0.0PARAM(2)=0.0. Note the call to nag_file_set_unit_advisory (x04ab) before the call to nag_ode_bvp_shoot_genpar (d02hb).
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)
y = tanϕ v = -(p1sinϕ+0.00002v2) vcosϕ ϕ = -p1v2
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. We estimate p1 = PARAM(1) = 32p1=PARAM(1)=32, p2 = PARAM(2) = 6000p2=PARAM(2)=6000 and p3 = PARAM(3) = 0.54p3=PARAM(3)=0.54 (though this last estimate is not important).
function nag_ode_bvp_shoot_genpar_example
% Set up initial values for first case.
n = 2;
n1 = 2;
p = [0.2;0];
pe = [1e-05;0.001];
e = [0.0001;0.0001];
m1 = 6;
marray = zeros(1,m1);

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

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

[pOut1, solnOut1, w1, ifail] = nag_ode_bvp_shoot_genpar(p, pe, e, int64(m1), @fcn1,...
    @bc1, @range1, 'n1', int64(n1), 'n', int64(n));
if (ifail ~= 0)
    % Problems in integration, or incorrect parameters.  Print
    % message and exit.
    if (ifail < 5 && ifail ~= 1)
        for i = 1:n
            fprintf('W(1,2) = %2.2f  ', w(1,2));
            fprintf('    W(.,1) = %2.2f',w(i,1));
            fprintf('\n');
        end
    end
    error('Warning: nag_ode_bvp_shoot_genpar returned with ifail = %d ',ifail);
else
    % 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; j = 1:n;
        m = x0 + double(i-1)*h;
        fprintf(' %6.3f        ', m);
        fprintf('%2.4f    ',solnOut1(j,i));
        fprintf('\n');
        % Save results for plotting.
        marray(i) = m;
    end
end

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

[pOut2, solnOut2, w2, ifail] = nag_ode_bvp_shoot_genpar(p, pe, e, int64(m1), @fcn2,...
    @bc2, @range2, 'n1', int64(n1), 'n', int64(n));
if (ifail ~= 0)
    % Problems in integration, or incorrect parameters.  Print
    % message and exit.
    if (ifail < 5 && ifail ~= 1)
        for i = 1:n
            fprintf('W(1,2) = %2.2f  ', w(1,2));
            fprintf('    W(.,1) = %2.2f\n',w(i,1));
        end
    end
    error('Warning: nag_ode_bvp_shoot_genpar returned with ifail = %d ',ifail);
else
    % 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; j = 1:n;
        q = x0 + double(i-1)*h;
        fprintf('%6.0f     ', q);
        fprintf('%10.4f  ',solnOut2(j,i));
        fprintf('\n');
        % Save results for plotting.
        qarray(i) = q;
    end
end

% Plot results.
fig = figure('Number', 'off');
display_plot(marray, solnOut1, 'Derivative');

fig = figure('Number', 'off');
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)
% 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 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.01 0.015 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 8 12 16]);
        set(haxes(iaxis), 'FontSize', 13);
    end
    set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
    % Set the title.
    title(['Parameterised Two-point BVP using Initial Value ', ...
        'Techniques & Newton Iteration'], titleFmt{:});
    % Label the axes.
    xlabel('x', labFmt{:});
    ylabel(haxes(1),'Solution', labFmt{:});
    ylabel(haxes(2),ylabelString, labFmt{:});
    % Add legend.
    legend('solution y(x)','derivative y''(x)','Location','Best');
    % 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,:));
    % 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 0 500 1000]);
    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 16]);
        set(haxes(iaxis), 'XTick', [0 4 8 12 16]);
        set(haxes(iaxis), 'FontSize', 13);
    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),'Height and Velocity', labFmt{:});
    ylabel(haxes(2),ylabelString, labFmt{:});
    % Add legend.
    legend('height','velocity','angle','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 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  

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

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

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

[pOut1, solnOut1, w1, ifail] = d02hb(p, pe, e, int64(m1), @fcn1,...
    @bc1, @range1, 'n1', int64(n1), 'n', int64(n));
if (ifail ~= 0)
    % Problems in integration, or incorrect parameters.  Print
    % message and exit.
    if (ifail < 5 && ifail ~= 1)
        for i = 1:n
            fprintf('W(1,2) = %2.2f  ', w(1,2));
            fprintf('    W(.,1) = %2.2f',w(i,1));
            fprintf('\n');
        end
    end
    error('Warning: d02hb returned with ifail = %d ',ifail);
else
    % 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; j = 1:n;
        m = x0 + double(i-1)*h;
        fprintf(' %6.3f        ', m);
        fprintf('%2.4f    ',solnOut1(j,i));
        fprintf('\n');
        % Save results for plotting.
        marray(i) = m;
    end
end

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

[pOut2, solnOut2, w2, ifail] = d02hb(p, pe, e, int64(m1), @fcn2,...
    @bc2, @range2, 'n1', int64(n1), 'n', int64(n));
if (ifail ~= 0)
    % Problems in integration, or incorrect parameters.  Print
    % message and exit.
    if (ifail < 5 && ifail ~= 1)
        for i = 1:n
            fprintf('W(1,2) = %2.2f  ', w(1,2));
            fprintf('    W(.,1) = %2.2f\n',w(i,1));
        end
    end
    error('Warning: d02hb returned with ifail = %d ',ifail);
else
    % 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; j = 1:n;
        q = x0 + double(i-1)*h;
        fprintf('%6.0f     ', q);
        fprintf('%10.4f  ',solnOut2(j,i));
        fprintf('\n');
        % Save results for plotting.
        qarray(i) = q;
    end
end

% Plot results.
fig = figure('Number', 'off');
display_plot(marray, solnOut1, 'Derivative');

fig = figure('Number', 'off');
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)
% 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 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.01 0.015 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 8 12 16]);
        set(haxes(iaxis), 'FontSize', 13);
    end
    set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
    % Set the title.
    title(['Parameterised Two-point BVP using Initial Value ', ...
        'Techniques & Newton Iteration'], titleFmt{:});
    % Label the axes.
    xlabel('x', labFmt{:});
    ylabel(haxes(1),'Solution', labFmt{:});
    ylabel(haxes(2),ylabelString, labFmt{:});
    % Add legend.
    legend('solution y(x)','derivative y''(x)','Location','Best');
    % 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,:));
    % 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 0 500 1000]);
    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 16]);
        set(haxes(iaxis), 'XTick', [0 4 8 12 16]);
        set(haxes(iaxis), 'FontSize', 13);
    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),'Height and Velocity', labFmt{:});
    ylabel(haxes(2),ylabelString, labFmt{:});
    % Add legend.
    legend('height','velocity','angle','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
 
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  


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