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_sum_accelerate (c06ba)

Purpose

nag_sum_accelerate (c06ba) accelerates the convergence of a given convergent sequence to its limit.

Syntax

[ncall, result, abserr, work, ifail] = c06ba(seqn, ncall, work)
[ncall, result, abserr, work, ifail] = nag_sum_accelerate(seqn, ncall, work)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: lwork has been removed from the interface
.

Description

nag_sum_accelerate (c06ba) performs Shanks' transformation on a given sequence of real values by means of the Epsilon algorithm of Wynn (1956). A (possibly unreliable) estimate of the absolute error is also given.
The function must be called repetitively, once for each new term in the sequence.

References

Shanks D (1955) Nonlinear transformations of divergent and slowly convergent sequences J. Math. Phys. 34 1–42
Wynn P (1956) On a device for computing the em(Sn)em(Sn) transformation Math. Tables Aids Comput. 10 91–96

Parameters

Compulsory Input Parameters

1:     seqn – double scalar
The next term of the sequence to be considered.
2:     ncall – int64int32nag_int scalar
On the first call ncall must be set to 00. Thereafter ncall must not be changed between calls.
3:     work(lwork) – double array
lwork, the dimension of the array, must satisfy the constraint lwork7lwork7.
Used as workspace, but must not be changed between calls.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

lwork

Output Parameters

1:     ncall – int64int32nag_int scalar
The number of terms in the sequence that have been considered.
2:     result – double scalar
The estimate of the limit of the sequence. For the first two calls, result = seqnresult=seqn.
3:     abserr – double scalar
An estimate of the absolute error in result. For the first three calls, abserr is set to a large machine-dependent constant.
4:     work(lwork) – double array
5:     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
On entry,ncall < 0ncall<0.
  ifail = 2ifail=2
On entry,lwork < 7lwork<7.

Accuracy

The accuracy of the absolute error estimate abserr varies considerably with the type of sequence to which the function is applied. In general it is better when applied to oscillating sequences than to monotonic sequences where it may be a severe underestimate.

Further Comments

Timing

The time taken is approximately proportional to the final value of ncall.

Choice of lwork

For long sequences, a ‘window’ of the last nn values can be used instead of all the terms of the sequence. Tests on a variety of problems indicate that a suitable value is n = 50 n=50 ; this implies a value for lwork of 5656. You are advised to experiment with other values for your own specific problems.

Convergence

nag_sum_accelerate (c06ba) will induce convergence in some divergent sequences. See Shanks (1955) for more details.

Example

This example attempts to sum the infinite series
((1)n + 1)/(n2) = (π2)/12
n = 1
n=1 (-1) n+1 n2 = π212
by considering the sequence of partial sums
1 2 3 10
,,,,
n = 1 n = 1 n = 1 n = 1
n= 1 1 , n= 1 2 , n= 1 3 , , n= 1 10
function nag_sum_accelerate_example
% Set up initial values before calling NAG routine for the first time.
work = zeros(16, 1);
answer = pi^2/12.0;
ncall = 0;
sig = 1.0;
seqn = 0.0;

% Initialize arrays for plotting.
x = zeros(1,10);
seqnArr = zeros(1,10);
resArr = zeros(1,10);
errArr = zeros(1,10);
absArr = zeros(1,10);

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

% Output headings.
disp('No. of  Term   Estimate  Estimated    Actual');
disp(' terms  value  of limit  abs error     error');

% Loop over terms.
for i = 1:10;
    r = double(i);
    seqn = seqn + sig/(r^2);

    % NB - input and output ncall *must* be the same variable (not e.g.
    % ncallOut).  Ditto for work.
    [ncall, result, abserr, work, ifail] = ...
        nag_sum_accelerate(seqn, int64(ncall), work);
    if ifail ~= 0
        % Incorrect parameter values.  Print message and exit.
        error('Warning: nag_sum_accelerate returned with ifail = %1d ',ifail);
    end
    err = result-answer;
    sig = -sig;
    if i <= 3
        % First three calls of nag_sum_accelerate return no error estimate.
        fprintf('%4d %8.4f %8.4f      -      %11.2e\n', ...
            i, seqn, result, err);
    else
        fprintf('%4d %8.4f %8.4f %11.2e %11.2e\n', ...
            i, seqn, result, abserr, err);
    end

    % Accumulate results for plotting.
    x(i) = i;
    seqnArr(i) = seqn;
    resArr(i) = result;
    errArr(i) = err;
    if i > 3
        absArr(i) = abserr;
    end
end

% Plot results.
fig = figure('Number', 'off');
display_plot(x, resArr, errArr, absArr)

function display_plot(x, y1, y2, y3)
% 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.
% Use a log plot for both curves.
[haxes, hline1, hline2] = plotyy(x, y1, x, abs(y2), 'semilogy', 'semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [0.7 1.1]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.7 0.8 0.9 1.0 1.1]);
% At this point, we must specify these labels explicitly (don't know why).
set(haxes(1), 'YTickLabel', [0.7; 0.8; 0.9; 1.0; 1.1]);
set(haxes(2), 'YLim', [1.0e-10 1]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1.0e-10 1.0e-8 1.0e-6 1.0e-4 1.0e-2 1.0]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [1 10]);
    set(haxes(iaxis), 'FontSize', 13);
    set(haxes(iaxis), 'XTick', [1 2 3 4 5 6 7 8 9 10]);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Estimate sum of infinite series by sequence of partial sums', ...
    titleFmt{:});
% Label the x axis.
xlabel('Number of terms in sequence', labFmt{:});
% Label the left and right y axes.
%set(haxes(1), 'YLim', [0 1]);
ylabel(haxes(1),'Result', labFmt{:});
ylabel(haxes(2),'abs(Error)', labFmt{:});
% Add the third curve (it's still a log plot because the last one was).
axes(haxes(2));
hold on;
hline3 = plot(x(4:10), y3(4:10));
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', 'o');
set(hline2, 'Linewidth', 0.5, 'Marker', 's');
set(hline3, 'Linewidth', 0.5, 'Marker', 'd');
% Add legend.
legend('Estimate of limit', 'Actual error', 'Estimated error', ...
    'Location', 'NorthEast');
 
nag_sum_accelerate example program results

No. of  Term   Estimate  Estimated    Actual
 terms  value  of limit  abs error     error
   1   1.0000   1.0000      -         1.78e-01
   2   0.7500   0.7500      -        -7.25e-02
   3   0.8611   0.8269      -         4.46e-03
   4   0.7986   0.8211    2.56e-01   -1.36e-03
   5   0.8386   0.8226    7.84e-02    1.23e-04
   6   0.8108   0.8224    5.97e-03   -3.26e-05
   7   0.8312   0.8225    1.52e-03    3.50e-06
   8   0.8156   0.8225    1.60e-04   -8.51e-07
   9   0.8280   0.8225    3.70e-05    1.01e-07
  10   0.8180   0.8225    4.48e-06   -2.32e-08

function c06ba_example
% Set up initial values before calling NAG routine for the first time.
work = zeros(16, 1);
answer = pi^2/12.0;
ncall = 0;
sig = 1.0;
seqn = 0.0;

% Initialize arrays for plotting.
x = zeros(1,10);
seqnArr = zeros(1,10);
resArr = zeros(1,10);
errArr = zeros(1,10);
absArr = zeros(1,10);

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

% Output headings.
disp('No. of  Term   Estimate  Estimated    Actual');
disp(' terms  value  of limit  abs error     error');

% Loop over terms.
for i = 1:10;
    r = double(i);
    seqn = seqn + sig/(r^2);

    % NB - input and output ncall *must* be the same variable (not e.g.
    % ncallOut).  Ditto for work.
    [ncall, result, abserr, work, ifail] = ...
        c06ba(seqn, int64(ncall), work);
    if ifail ~= 0
        % Incorrect parameter values.  Print message and exit.
        error('Warning: c06ba returned with ifail = %1d ',ifail);
    end
    err = result-answer;
    sig = -sig;
    if i <= 3
        % First three calls of c06ba return no error estimate.
        fprintf('%4d %8.4f %8.4f      -      %11.2e\n', ...
            i, seqn, result, err);
    else
        fprintf('%4d %8.4f %8.4f %11.2e %11.2e\n', ...
            i, seqn, result, abserr, err);
    end

    % Accumulate results for plotting.
    x(i) = i;
    seqnArr(i) = seqn;
    resArr(i) = result;
    errArr(i) = err;
    if i > 3
        absArr(i) = abserr;
    end
end

% Plot results.
fig = figure('Number', 'off');
display_plot(x, resArr, errArr, absArr)

function display_plot(x, y1, y2, y3)
% 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.
% Use a log plot for both curves.
[haxes, hline1, hline2] = plotyy(x, y1, x, abs(y2), 'semilogy', 'semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [0.7 1.1]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [0.7 0.8 0.9 1.0 1.1]);
% At this point, we must specify these labels explicitly (don't know why).
set(haxes(1), 'YTickLabel', [0.7; 0.8; 0.9; 1.0; 1.1]);
set(haxes(2), 'YLim', [1.0e-10 1]);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1.0e-10 1.0e-8 1.0e-6 1.0e-4 1.0e-2 1.0]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [1 10]);
    set(haxes(iaxis), 'FontSize', 13);
    set(haxes(iaxis), 'XTick', [1 2 3 4 5 6 7 8 9 10]);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Estimate sum of infinite series by sequence of partial sums', ...
    titleFmt{:});
% Label the x axis.
xlabel('Number of terms in sequence', labFmt{:});
% Label the left and right y axes.
%set(haxes(1), 'YLim', [0 1]);
ylabel(haxes(1),'Result', labFmt{:});
ylabel(haxes(2),'abs(Error)', labFmt{:});
% Add the third curve (it's still a log plot because the last one was).
axes(haxes(2));
hold on;
hline3 = plot(x(4:10), y3(4:10));
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', 'o');
set(hline2, 'Linewidth', 0.5, 'Marker', 's');
set(hline3, 'Linewidth', 0.5, 'Marker', 'd');
% Add legend.
legend('Estimate of limit', 'Actual error', 'Estimated error', ...
    'Location', 'NorthEast');
 
c06ba example program results

No. of  Term   Estimate  Estimated    Actual
 terms  value  of limit  abs error     error
   1   1.0000   1.0000      -         1.78e-01
   2   0.7500   0.7500      -        -7.25e-02
   3   0.8611   0.8269      -         4.46e-03
   4   0.7986   0.8211    2.56e-01   -1.36e-03
   5   0.8386   0.8226    7.84e-02    1.23e-04
   6   0.8108   0.8224    5.97e-03   -3.26e-05
   7   0.8312   0.8225    1.52e-03    3.50e-06
   8   0.8156   0.8225    1.60e-04   -8.51e-07
   9   0.8280   0.8225    3.70e-05    1.01e-07
  10   0.8180   0.8225    4.48e-06   -2.32e-08


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