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_invlaplace_weeks (c06lb)

Purpose

nag_sum_invlaplace_weeks (c06lb) computes the inverse Laplace transform f(t) f(t)  of a user-supplied function F(s) F(s) , defined for complex ss. The function uses a modification of Weeks' method which is suitable when f(t) f(t)  has continuous derivatives of all orders. The function returns the coefficients of an expansion which approximates f(t) f(t)  and can be evaluated for given values of tt by subsequent calls of nag_sum_invlaplace_weeks_eval (c06lc).

Syntax

[sigma, b, m, acoef, errvec, ifail] = c06lb(f, sigma0, sigma, b, epstol, 'mmax', mmax)
[sigma, b, m, acoef, errvec, ifail] = nag_sum_invlaplace_weeks(f, sigma0, sigma, b, epstol, 'mmax', mmax)

Description

Given a function f (t) f (t)  of a real variable tt, its Laplace transform F(s) F(s)  is a function of a complex variable ss, defined by
F(s) = estf(t)dt,  Re(s) > σ0.
0
F (s) = 0 e-st f (t) dt ,   Re(s) > σ0 .
Then f(t) f(t)  is the inverse Laplace transform of F(s) F(s) . The value σ0 σ0  is referred to as the abscissa of convergence of the Laplace transform; it is the rightmost real part of the singularities of F(s) F(s) .
nag_sum_invlaplace_weeks (c06lb), along with its companion nag_sum_invlaplace_weeks_eval (c06lc), attempts to solve the following problem:
The method is a modification of Weeks' method (see Garbow et al. (1988a)), which approximates f(t) f(t)  by a truncated Laguerre expansion:
m1
(t) = eσtaie bt / 2 Li(bt),  σ > σ0,  b > 0
i = 0
f~ (t) = eσt i=0 m-1 ai e -bt / 2 Li (bt) ,   σ > σ0 ,   b > 0
where Li (x) Li (x)  is the Laguerre polynomial of degree ii. This function computes the coefficients ai ai  of the above Laguerre expansion; the expansion can then be evaluated for specified tt by calling nag_sum_invlaplace_weeks_eval (c06lc). You must supply the value of σ0 σ0 , and also suitable values for σσ and bb: see Section [Further Comments] for guidance.
The method is only suitable when f(t) f(t)  has continuous derivatives of all orders. For such functions the approximation (t) f~ (t)  is usually good and inexpensive. The function will fail with an error exit if the method is not suitable for the supplied function F(s) F(s) .
The function is designed to satisfy an accuracy criterion of the form:
|( f(t) (t) )/(eσt)| < εtol,   for all ​t
| f(t)- f~(t) e σt | < ε tol ,   for all ​t
where εtol εtol  is a user-supplied bound. The error measure on the left-hand side is referred to as the pseudo-relative error, or pseudo-error for short. Note that if σ > 0 σ>0  and tt is large, the absolute error in (t) f~ (t)  may be very large.
nag_sum_invlaplace_weeks (c06lb) is derived from the function MODUL1 in Garbow et al. (1988a).

References

Garbow B S, Giunta G, Lyness J N and Murli A (1988a) Software for an implementation of Weeks' method for the inverse laplace transform problem ACM Trans. Math. Software 14 163–170
Garbow B S, Giunta G, Lyness J N and Murli A (1988b) Algorithm 662: A Fortran software package for the numerical inversion of the Laplace transform based on Weeks' method ACM Trans. Math. Software 14 171–176

Parameters

Compulsory Input Parameters

1:     f – function handle or string containing name of m-file
f must return the value of the Laplace transform function F(s)F(s) for a given complex value of ss.
[result] = f(s)

Input Parameters

1:     s – complex scalar
The value of ss for which F(s)F(s) must be evaluated. The real part of s is greater than σ0σ0.

Output Parameters

1:     result – complex scalar
The result of the function.
2:     sigma0 – double scalar
The abscissa of convergence of the Laplace transform, σ0σ0.
3:     sigma – double scalar
The parameter σσ of the Laguerre expansion. If on entry sigmaσ0sigmaσ0, sigma is reset to σ0 + 0.7σ0+0.7.
4:     b – double scalar
The parameter bb of the Laguerre expansion. If on entry b < 2 (σσ0)b < 2 ( σ - σ0 ), b is reset to 2.5(σσ0)2.5(σ-σ0).
5:     epstol – double scalar
The required relative pseudo-accuracy, that is, an upper bound on |f(t)(t)| eσt| f (t) - f~ (t) | e-σt.

Optional Input Parameters

1:     mmax – int64int32nag_int scalar
An upper bound on the number of Laguerre expansion coefficients to be computed. The number of coefficients actually computed is always a power of 22, so mmax should be a power of 22; if mmax is not a power of 22 then the maximum number of coefficients calculated will be the largest power of 22 less than mmax.
Default: 10241024
Constraint: mmax8mmax8.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     sigma – double scalar
The value actually used for σσ, as just described.
2:     b – double scalar
The value actually used for bb, as just described.
3:     m – int64int32nag_int scalar
The number of Laguerre expansion coefficients actually computed. The number of calls to f is m / 2 + 2m/2+2.
4:     acoef(mmax) – double array
The first m elements contain the computed Laguerre expansion coefficients, aiai.
5:     errvec(88) – double array
An 88-component vector of diagnostic information.
errvec(1)errvec1
Overall estimate of the pseudo-error
|f(t)(t)| eσt = errvec(2) + errvec(3) + errvec(4) | f (t) - f~ (t) | e-σt =errvec2 + errvec3 + errvec4.
errvec(2)errvec2
Estimate of the discretization pseudo-error.
errvec(3)errvec3
Estimate of the truncation pseudo-error.
errvec(4)errvec4
Estimate of the condition pseudo-error on the basis of minimal noise levels in function values.
errvec(5)errvec5
KK, coefficient of a heuristic decay function for the expansion coefficients.
errvec(6)errvec6
RR, base of the decay function for the expansion coefficients.
errvec(7)errvec7
Logarithm of the largest expansion coefficient.
errvec(8)errvec8
Logarithm of the smallest nonzero expansion coefficient.
The values KK and RR returned in errvec(5)errvec5 and errvec(6)errvec6 define a decay function KRiKR-i constructed by the function for the purposes of error estimation. It satisfies
|ai| < K Ri ,   ​ i = 1, 2, , m .
|ai| < K R -i ,   ​ i= 1, 2, , m .
6:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_sum_invlaplace_weeks (c06lb) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

  ifail = 1ifail=1
On entry,mmax < 8mmax<8.
W ifail = 2ifail=2
The estimated pseudo-error bounds are slightly larger than epstol. Note, however, that the actual errors in the final results may be smaller than epstol as bounds independent of the value of tt are pessimistic.
W ifail = 3ifail=3
Computation was terminated early because the estimate of rounding error was greater than epstol. Increasing epstol may help.
  ifail = 4ifail=4
The decay rate of the coefficients is too small. Increasing mmax may help.
  ifail = 5ifail=5
The decay rate of the coefficients is too small. In addition the rounding error is such that the required accuracy cannot be obtained. Increasing mmax or epstol may help.
W ifail = 6ifail=6
The behaviour of the coefficients does not enable reasonable prediction of error bounds. Check the value of sigma0. In this case, errvec(i) errveci  is set to 1.0 -1.0 , for i = 1,2,,5i=1,2,,5.
When ifail3ifail3, changing sigma or b may help. If not, the method should be abandoned.

Accuracy

The error estimate returned in errvec(1) errvec1  has been found in practice to be a highly reliable bound on the pseudo-error |f(t)(t)| eσt |f(t)-f~(t)| e-σt .

Further Comments

The Role of σ0

Nearly all techniques for inversion of the Laplace transform require you to supply the value of σ0 σ0 , the convergence abscissa, or else an upper bound on σ0 σ0 . For this function, one of the reasons for having to supply σ0 σ0  is that the parameter σσ must be greater than σ0 σ0 ; otherwise the series for (t) f~(t)  will not converge.
If you do not know the value of σ0 σ0 , you must be prepared for significant preliminary effort, either in experimenting with the method and obtaining chaotic results, or in attempting to locate the rightmost singularity of F(s) F(s) .
The value of σ0 σ0  is also relevant in defining a natural accuracy criterion. For large tt, f(t) f(t)  is of uniform numerical order k eσ0t k e σ0t , so a natural measure of relative accuracy of the approximation (t) f~ (t)  is:
εnat (t) = ( (t) f (t) ) / eσ0t .
εnat (t) = ( f~ (t) - f (t) ) / e σ0t .
nag_sum_invlaplace_weeks (c06lb) uses the supplied value of σ0 σ0  only in determining the values of σσ and bb (see Sections [Choice of ] and [Choice of ]); thereafter it bases its computation entirely on σσ and bb.

Choice of σ

Even when the value of σ0 σ0  is known, choosing a value for σσ is not easy. Briefly, the series for (t) f~(t)  converges slowly when σσ0 σ-σ0  is small, and faster when σσ0 σ-σ0  is larger. However the natural accuracy measure satisfies
|εnat(t)| < εtol e (σσ0) t
| εnat (t) | < εtol e ( σ - σ0 ) t
and this degrades exponentially with tt, the exponential constant being σσ0 σ-σ0 .
Hence, if you require meaningful results over a large range of values of tt, you should choose σσ0 σ-σ0  small, in which case the series for (t) f~(t)  converges slowly; while for a smaller range of values of tt, you can allow σσ0 σ-σ0  to be larger and obtain faster convergence.
The default value for σσ used by nag_sum_invlaplace_weeks (c06lb) is σ0 + 0.7 σ0+0.7 . There is no theoretical justification for this.

Choice of b

The simplest advice for choosing bb is to set b / 2 σ σ0 b/2 σ - σ0 . The default value used by the function is 2.5 (σσ0) 2.5 ( σ - σ0 ) . A more refined choice is to set
b / 2min |σsj|
j
b/2 minj |σ-sj|
where sj sj  are the singularities of F(s) F(s) .

Example

This example computes values of the inverse Laplace transform of the function
F(s) = 3/(s29) .
F(s) = 3 s2-9 .
The exact answer is
f(t) = sinh3t .
f(t) = sinh3t .
The program first calls nag_sum_invlaplace_weeks (c06lb) to compute the coefficients of the Laguerre expansion, and then calls nag_sum_invlaplace_weeks_eval (c06lc) to evaluate the expansion at t = 0 t=0, 11, 22, 33, 44, 55.
function nag_sum_invlaplace_weeks_example
% Initialize variables and arrays.
sigma0 = 3;
epstol = 1e-5;
b = 0;
sigma = 0;
n = 5;
earray = zeros(1, n+1);
jarray = zeros(1, n+1);
farray = zeros(1, n+1);
parray = zeros(1, n+1);

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

[sigmaOut, bOut, m, acoef, errvec, ifail] = ...
    nag_sum_invlaplace_weeks(@f, sigma0, sigma, b, epstol);
if ifail ~= 0
    % Convergence problems.  Print message and exit.
    error('Warning: nag_sum_invlaplace_weeks returned with ifail = %1d ',ifail);
end

% Prepare to output results.
disp(['No. of coefficients returned by nag_sum_invlaplace_weeks = ',num2str(m)]);
disp('  ');
disp('      Computed       Exact         Pseudo');
disp('t       f(t)          f(t)          error');
% Evaluate inverse transform for different values of t.  We use nag_sum_invlaplace_weeks_eval,
% which calculates the transform from the coefficients given by nag_sum_invlaplace_weeks.
for j = 0:5
    t = double(j);

    [finv, ifail] = nag_sum_invlaplace_weeks_eval(t, sigmaOut, bOut, acoef, errvec, 'm', m);
    if ifail ~= 0
        % Approximation is too large or too small.  Print message and exit.
        error('Warning: nag_sum_invlaplace_weeks_eval returned with ifail = %1d ',ifail);
    end
    exact = sinh(3.0*t);
    pserr = abs(finv-exact)/exp(sigmaOut*t);
    fprintf('%d   %10.4d   %11.4d     %8.4d\n', t, finv, exact, pserr);
    % Create arrays to be used for plotting.
    jarray(j+1) = t;
    farray(j+1) = finv;
    parray(j+1) = pserr;
end

% Plot results.
fig = figure('Number', 'off');
display_plot(jarray, farray, parray)


function [f] = f(s)
% Evaluate the Laplace transform function.
f=3.0/(s^2-9.0);
if isreal(f)
    f=complex(f);
end
function plot(jarray, farray, parray)
% 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(jarray, farray, jarray, parray,...
    'semilogy','semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [1.0e-10 1.0e10]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [1.0e-10 1.0e-5 1.0 1.0e5 1.0e10]);
set(haxes(2), 'YLim', [1.0e-10 1.0e-8]);
%set(haxes(2), 'FontSize', 11);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-9 1e-8]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 5]);
    set(haxes(iaxis), 'XTick', [0 1 2 3 4 5]);
    set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)', titleFmt{:});
% Label the x axis.
xlabel('t', labFmt{:});
% Label the left & right y axes.
ylabel(haxes(1),'f(t)', labFmt{:});
ylabel(haxes(2),'Pseudo Error', labFmt{:});
% Label the curves.
legend('f(t)','Pseudo Error','location','North')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '-');
function display_plot(jarray, farray, parray)
% 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(jarray, farray, jarray, parray,...
    'semilogy','semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [1.0e-10 1.0e10]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [1.0e-10 1.0e-5 1.0 1.0e5 1.0e10]);
set(haxes(2), 'YLim', [1.0e-10 1.0e-8]);
%set(haxes(2), 'FontSize', 11);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-9 1e-8]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 5]);
    set(haxes(iaxis), 'XTick', [0 1 2 3 4 5]);
    set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)', titleFmt{:});
% Label the x axis.
xlabel('t', labFmt{:});
% Label the left & right y axes.
ylabel(haxes(1),'f(t)', labFmt{:});
ylabel(haxes(2),'Pseudo Error', labFmt{:});
% Label the curves.
legend('f(t)','Pseudo Error','location','North')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '-');
 
nag_sum_invlaplace_weeks example program results

No. of coefficients returned by nag_sum_invlaplace_weeks = 64
  
      Computed       Exact         Pseudo
t       f(t)          f(t)          error
0   1.5129e-09          0000     1.5129e-09
1   1.0018e+01    1.0018e+01     1.7394e-09
2   2.0171e+02    2.0171e+02     1.2471e-10
3   4.0515e+03    4.0515e+03     9.7722e-10
4   8.1377e+04    8.1377e+04     3.0221e-10
5   1.6345e+06    1.6345e+06     1.6991e-09

function c06lb_example
% Initialize variables and arrays.
sigma0 = 3;
epstol = 1e-5;
b = 0;
sigma = 0;
n = 5;
earray = zeros(1, n+1);
jarray = zeros(1, n+1);
farray = zeros(1, n+1);
parray = zeros(1, n+1);

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

[sigmaOut, bOut, m, acoef, errvec, ifail] = ...
    c06lb(@f, sigma0, sigma, b, epstol);
if ifail ~= 0
    % Convergence problems.  Print message and exit.
    error('Warning: c06lb returned with ifail = %1d ',ifail);
end

% Prepare to output results.
disp(['No. of coefficients returned by c06lb = ',num2str(m)]);
disp('  ');
disp('      Computed       Exact         Pseudo');
disp('t       f(t)          f(t)          error');
% Evaluate inverse transform for different values of t.  We use c06lc,
% which calculates the transform from the coefficients given by c06lb.
for j = 0:5
    t = double(j);

    [finv, ifail] = c06lc(t, sigmaOut, bOut, acoef, errvec, 'm', m);
    if ifail ~= 0
        % Approximation is too large or too small.  Print message and exit.
        error('Warning: c06lc returned with ifail = %1d ',ifail);
    end
    exact = sinh(3.0*t);
    pserr = abs(finv-exact)/exp(sigmaOut*t);
    fprintf('%d   %10.4d   %11.4d     %8.4d\n', t, finv, exact, pserr);
    % Create arrays to be used for plotting.
    jarray(j+1) = t;
    farray(j+1) = finv;
    parray(j+1) = pserr;
end

% Plot results.
fig = figure('Number', 'off');
display_plot(jarray, farray, parray)


function [f] = f(s)
% Evaluate the Laplace transform function.
f=3.0/(s^2-9.0);
if isreal(f)
    f=complex(f);
end
function plot(jarray, farray, parray)
% 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(jarray, farray, jarray, parray,...
    'semilogy','semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [1.0e-10 1.0e10]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [1.0e-10 1.0e-5 1.0 1.0e5 1.0e10]);
set(haxes(2), 'YLim', [1.0e-10 1.0e-8]);
%set(haxes(2), 'FontSize', 11);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-9 1e-8]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 5]);
    set(haxes(iaxis), 'XTick', [0 1 2 3 4 5]);
    set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)', titleFmt{:});
% Label the x axis.
xlabel('t', labFmt{:});
% Label the left & right y axes.
ylabel(haxes(1),'f(t)', labFmt{:});
ylabel(haxes(2),'Pseudo Error', labFmt{:});
% Label the curves.
legend('f(t)','Pseudo Error','location','North')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '-');
function display_plot(jarray, farray, parray)
% 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(jarray, farray, jarray, parray,...
    'semilogy','semilogy');
% Set the axis limits and the tick specifications to beautify the plot.
set(haxes(1), 'YLim', [1.0e-10 1.0e10]);
set(haxes(1), 'YMinorTick', 'on');
set(haxes(1), 'YTick', [1.0e-10 1.0e-5 1.0 1.0e5 1.0e10]);
set(haxes(2), 'YLim', [1.0e-10 1.0e-8]);
%set(haxes(2), 'FontSize', 11);
set(haxes(2), 'YMinorTick', 'on');
set(haxes(2), 'YTick', [1e-10 1e-9 1e-8]);
for iaxis = 1:2
    % These properties must be the same for both sets of axes.
    set(haxes(iaxis), 'XLim', [0 5]);
    set(haxes(iaxis), 'XTick', [0 1 2 3 4 5]);
    set(haxes(iaxis), 'FontSize', 13);
end
set(gca, 'box', 'off'); % so ticks aren't shown on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)', titleFmt{:});
% Label the x axis.
xlabel('t', labFmt{:});
% Label the left & right y axes.
ylabel(haxes(1),'f(t)', labFmt{:});
ylabel(haxes(2),'Pseudo Error', labFmt{:});
% Label the curves.
legend('f(t)','Pseudo Error','location','North')
% Set some features of the three lines.
set(hline1, 'Linewidth', 0.5, 'Marker', '+', 'LineStyle', '-');
set(hline2, 'Linewidth', 0.5, 'Marker', 'x', 'LineStyle', '-');
 
c06lb example program results

No. of coefficients returned by c06lb = 64
  
      Computed       Exact         Pseudo
t       f(t)          f(t)          error
0   1.5129e-09          0000     1.5129e-09
1   1.0018e+01    1.0018e+01     1.7394e-09
2   2.0171e+02    2.0171e+02     1.2471e-10
3   4.0515e+03    4.0515e+03     9.7722e-10
4   8.1377e+04    8.1377e+04     3.0221e-10
5   1.6345e+06    1.6345e+06     1.6991e-09


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