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)

▸▿ Contents

    1  Purpose
    2  Syntax
    3  Description
    4  References
▸▿ 5  Parameters
    6  Error Indicators and Warnings
    7  Accuracy
▸▿ 8  Further Comments
    9  Example

Purpose

nag_sum_invlaplace_weeks (c06lb) computes the inverse Laplace transform ft  of a user-supplied function Fs , defined for complex s. The function uses a modification of Weeks' method which is suitable when ft  has continuous derivatives of all orders. The function returns the coefficients of an expansion which approximates ft  and can be evaluated for given values of t 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  of a real variable t, its Laplace transform Fs  is a function of a complex variable s, defined by
F s = ∫0∞ e-st f t dt ,   Res > σ0 .  
Then ft  is the inverse Laplace transform of Fs . The value σ0  is referred to as the abscissa of convergence of the Laplace transform; it is the rightmost real part of the singularities of Fs .
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 ft  by a truncated Laguerre expansion:
f~ t = eσt ∑ i=0 m-1 ai e -bt / 2 Li bt ,   σ > σ0 ,   b > 0  
where Li x  is the Laguerre polynomial of degree i. This function computes the coefficients ai  of the above Laguerre expansion; the expansion can then be evaluated for specified t by calling nag_sum_invlaplace_weeks_eval (c06lc). You must supply the value of σ0 , and also suitable values for σ and b: see Further Comments for guidance.
The method is only suitable when ft  has continuous derivatives of all orders. For such functions the approximation 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 Fs .
The function is designed to satisfy an accuracy criterion of the form:
ft- f~t e σt < ε tol ,   for all ​t  
where ε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  and t is large, the absolute error in 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 Fs for a given complex value of s.
[result] = f(s)

Input Parameters

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

Output Parameters

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

Optional Input Parameters

1:     mmax – int64int32nag_int scalar
Suggested value: mmax=1024 is sufficient for all but a few exceptional cases.
Default: 1024
An upper bound on the number of Laguerre expansion coefficients to be computed. The number of coefficients actually computed is always a power of 2, so mmax should be a power of 2; if mmax is not a power of 2 then the maximum number of coefficients calculated will be the largest power of 2 less than mmax.
Constraint: mmax≥8.

Output Parameters

1:     sigma – double scalar
The value actually used for σ, as just described.
2:     b – double scalar
The value actually used for b, 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+2.
4:     acoefmmax – double array
The first m elements contain the computed Laguerre expansion coefficients, ai.
5:     errvec8 – double array
An 8-component vector of diagnostic information.
errvec1
Overall estimate of the pseudo-error
f t - f~ t e-σt =errvec2 + errvec3 + errvec4.
errvec2
Estimate of the discretization pseudo-error.
errvec3
Estimate of the truncation pseudo-error.
errvec4
Estimate of the condition pseudo-error on the basis of minimal noise levels in function values.
errvec5
K, coefficient of a heuristic decay function for the expansion coefficients.
errvec6
R, base of the decay function for the expansion coefficients.
errvec7
Logarithm of the largest expansion coefficient.
errvec8
Logarithm of the smallest nonzero expansion coefficient.
The values K and R returned in errvec5 and errvec6 define a decay function KR-i constructed by the function for the purposes of error estimation. It satisfies
ai < K R -i ,   ​ i= 1, 2, …, m .  
6:     ifail – int64int32nag_int scalar
ifail=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=1
On entry,mmax<8.
W  ifail=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 t are pessimistic.
W  ifail=3
Computation was terminated early because the estimate of rounding error was greater than epstol. Increasing epstol may help.
   ifail=4
The decay rate of the coefficients is too small. Increasing mmax may help.
   ifail=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=6
The behaviour of the coefficients does not enable reasonable prediction of error bounds. Check the value of sigma0. In this case, errveci  is set to -1.0 , for i=1,2,…,5.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   ifail=-999
Dynamic memory allocation failed.
When ifail≥3, changing sigma or b may help. If not, the method should be abandoned.

Accuracy

The error estimate returned in errvec1  has been found in practice to be a highly reliable bound on the pseudo-error ft-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 , the convergence abscissa, or else an upper bound on σ0 . For this function, one of the reasons for having to supply σ0  is that the argument σ must be greater than σ0 ; otherwise the series for f~t  will not converge.
If you do not know the value of σ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 Fs .
The value of σ0  is also relevant in defining a natural accuracy criterion. For large t, ft  is of uniform numerical order k e σ0t , so a natural measure of relative accuracy of the approximation f~ t  is:
εnat t = f~ t - f t / e σ0t .  
nag_sum_invlaplace_weeks (c06lb) uses the supplied value of σ0  only in determining the values of σ and b (see Choice of and Choice of ); thereafter it bases its computation entirely on σ and b.

Choice of σ

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

Choice of b

The simplest advice for choosing b is to set b/2 ≥ σ - σ0 . The default value used by the function is 2.5 σ - σ0 . A more refined choice is to set
b/2 ≥ minj σ-sj  
where sj  are the singularities of Fs .

Example

This example computes values of the inverse Laplace transform of the function
Fs = 3 s2-9 .  
The exact answer is
ft = sinh⁡3t .  
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, 1, 2, 3, 4, 5.

Open in the MATLAB editor: c06lb_example

function c06lb_example


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

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

[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.
fig1 = figure;
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)
% 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), '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]);
end
set(gca, 'box', 'off'); % no ticks on opposite axes.
% Set the title.
title('Inverse Laplace Transform of 3/(s^2-9)');
% Label the x axis.
xlabel('t');
% Label the left & right y axes.
ylabel(haxes(1),'f(t)');
ylabel(haxes(2),'Pseudo Error');
% 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)
% 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), '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]);
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)');
% Label the x axis.
xlabel('t');
% Label the left & right y axes.
ylabel(haxes(1),'f(t)');
yh=ylabel(haxes(2),'Pseudo Error');
set(yh, 'position',[4.7,5e-10]);
% 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 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
c06lb_fig1.png
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–2015