Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

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\left(t\right)$ of a user-supplied function $F\left(s\right)$, defined for complex $s$. The function uses a modification of Weeks' method which is suitable when $f\left(t\right)$ has continuous derivatives of all orders. The function returns the coefficients of an expansion which approximates $f\left(t\right)$ 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\left(t\right)$ of a real variable $t$, its Laplace transform $F\left(s\right)$ is a function of a complex variable $s$, defined by
 $F s = ∫0∞ e-st f t dt , Res > σ0 .$
Then $f\left(t\right)$ is the inverse Laplace transform of $F\left(s\right)$. The value ${\sigma }_{0}$ is referred to as the abscissa of convergence of the Laplace transform; it is the rightmost real part of the singularities of $F\left(s\right)$.
nag_sum_invlaplace_weeks (c06lb), along with its companion nag_sum_invlaplace_weeks_eval (c06lc), attempts to solve the following problem:
• given a function $F\left(s\right)$, compute values of its inverse Laplace transform $f\left(t\right)$ for specified values of $t$.
The method is a modification of Weeks' method (see Garbow et al. (1988a)), which approximates $f\left(t\right)$ by a truncated Laguerre expansion:
 $f~ t = eσt ∑ i=0 m-1 ai e -bt / 2 Li bt , σ > σ0 , b > 0$
where ${L}_{i}\left(x\right)$ is the Laguerre polynomial of degree $i$. This function computes the coefficients ${a}_{i}$ 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 ${\sigma }_{0}$, and also suitable values for $\sigma$ and $b$: see Further Comments for guidance.
The method is only suitable when $f\left(t\right)$ has continuous derivatives of all orders. For such functions the approximation $\stackrel{~}{f}\left(t\right)$ is usually good and inexpensive. The function will fail with an error exit if the method is not suitable for the supplied function $F\left(s\right)$.
The function is designed to satisfy an accuracy criterion of the form:
 $ft- f~t e σt < ε tol , for all ​t$
where ${\epsilon }_{\mathit{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 $\sigma >0$ and $t$ is large, the absolute error in $\stackrel{~}{f}\left(t\right)$ 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:     $\mathrm{f}$ – function handle or string containing name of m-file
f must return the value of the Laplace transform function $F\left(s\right)$ for a given complex value of $s$.
[result] = f(s)

Input Parameters

1:     $\mathrm{s}$ – complex scalar
The value of $s$ for which $F\left(s\right)$ must be evaluated. The real part of s is greater than ${\sigma }_{0}$.

Output Parameters

1:     $\mathrm{result}$ – complex scalar
The value of the Laplace transform function $F\left(s\right)$ for the given complex value, $s$.
2:     $\mathrm{sigma0}$ – double scalar
The abscissa of convergence of the Laplace transform, ${\sigma }_{0}$.
3:     $\mathrm{sigma}$ – double scalar
The parameter $\sigma$ of the Laguerre expansion. If on entry ${\mathbf{sigma}}\le {\sigma }_{0}$, sigma is reset to ${\sigma }_{0}+0.7$.
4:     $\mathrm{b}$ – double scalar
The parameter $b$ of the Laguerre expansion. If on entry ${\mathbf{b}}<2\left(\sigma -{\sigma }_{0}\right)$, b is reset to $2.5\left(\sigma -{\sigma }_{0}\right)$.
5:     $\mathrm{epstol}$ – double scalar
The required relative pseudo-accuracy, that is, an upper bound on $\left|f\left(t\right)-\stackrel{~}{f}\left(t\right)\right|{e}^{-\sigma t}$.

### Optional Input Parameters

1:     $\mathrm{mmax}$int64int32nag_int scalar
Suggested value: ${\mathbf{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: ${\mathbf{mmax}}\ge 8$.

### Output Parameters

1:     $\mathrm{sigma}$ – double scalar
The value actually used for $\sigma$, as just described.
2:     $\mathrm{b}$ – double scalar
The value actually used for $b$, as just described.
3:     $\mathrm{m}$int64int32nag_int scalar
The number of Laguerre expansion coefficients actually computed. The number of calls to f is ${\mathbf{m}}/2+2$.
4:     $\mathrm{acoef}\left({\mathbf{mmax}}\right)$ – double array
The first m elements contain the computed Laguerre expansion coefficients, ${a}_{i}$.
5:     $\mathrm{errvec}\left(8\right)$ – double array
An $8$-component vector of diagnostic information.
${\mathbf{errvec}}\left(1\right)$
Overall estimate of the pseudo-error
$\left|f\left(t\right)-\stackrel{~}{f}\left(t\right)\right|{e}^{-\sigma t}={\mathbf{errvec}}\left(2\right)+{\mathbf{errvec}}\left(3\right)+{\mathbf{errvec}}\left(4\right)$.
${\mathbf{errvec}}\left(2\right)$
Estimate of the discretization pseudo-error.
${\mathbf{errvec}}\left(3\right)$
Estimate of the truncation pseudo-error.
${\mathbf{errvec}}\left(4\right)$
Estimate of the condition pseudo-error on the basis of minimal noise levels in function values.
${\mathbf{errvec}}\left(5\right)$
$K$, coefficient of a heuristic decay function for the expansion coefficients.
${\mathbf{errvec}}\left(6\right)$
$R$, base of the decay function for the expansion coefficients.
${\mathbf{errvec}}\left(7\right)$
Logarithm of the largest expansion coefficient.
${\mathbf{errvec}}\left(8\right)$
Logarithm of the smallest nonzero expansion coefficient.
The values $K$ and $R$ returned in ${\mathbf{errvec}}\left(5\right)$ and ${\mathbf{errvec}}\left(6\right)$ define a decay function $K{R}^{-i}$ constructed by the function for the purposes of error estimation. It satisfies
 $ai < K R -i , ​ i= 1, 2, …, m .$
6:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{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.

${\mathbf{ifail}}=1$
 On entry, ${\mathbf{mmax}}<8$.
W  ${\mathbf{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  ${\mathbf{ifail}}=3$
Computation was terminated early because the estimate of rounding error was greater than epstol. Increasing epstol may help.
${\mathbf{ifail}}=4$
The decay rate of the coefficients is too small. Increasing mmax may help.
${\mathbf{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  ${\mathbf{ifail}}=6$
The behaviour of the coefficients does not enable reasonable prediction of error bounds. Check the value of sigma0. In this case, ${\mathbf{errvec}}\left(\mathit{i}\right)$ is set to $-1.0$, for $\mathit{i}=1,2,\dots ,5$.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
When ${\mathbf{ifail}}\ge {\mathbf{3}}$, changing sigma or b may help. If not, the method should be abandoned.

## Accuracy

The error estimate returned in ${\mathbf{errvec}}\left(1\right)$ has been found in practice to be a highly reliable bound on the pseudo-error $\left|f\left(t\right)-\stackrel{~}{f}\left(t\right)\right|{e}^{-\sigma t}$.

### The Role of σ0

Nearly all techniques for inversion of the Laplace transform require you to supply the value of ${\sigma }_{0}$, the convergence abscissa, or else an upper bound on ${\sigma }_{0}$. For this function, one of the reasons for having to supply ${\sigma }_{0}$ is that the argument $\sigma$ must be greater than ${\sigma }_{0}$; otherwise the series for $\stackrel{~}{f}\left(t\right)$ will not converge.
If you do not know the value of ${\sigma }_{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\left(s\right)$.
The value of ${\sigma }_{0}$ is also relevant in defining a natural accuracy criterion. For large $t$, $f\left(t\right)$ is of uniform numerical order $k{e}^{{\sigma }_{0}t}$, so a natural measure of relative accuracy of the approximation $\stackrel{~}{f}\left(t\right)$ is:
 $εnat t = f~ t - f t / e σ0t .$
nag_sum_invlaplace_weeks (c06lb) uses the supplied value of ${\sigma }_{0}$ only in determining the values of $\sigma$ and $b$ (see Choice of and Choice of ); thereafter it bases its computation entirely on $\sigma$ and $b$.

### Choice of σ

Even when the value of ${\sigma }_{0}$ is known, choosing a value for $\sigma$ is not easy. Briefly, the series for $\stackrel{~}{f}\left(t\right)$ converges slowly when $\sigma -{\sigma }_{0}$ is small, and faster when $\sigma -{\sigma }_{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 $\sigma -{\sigma }_{0}$.
Hence, if you require meaningful results over a large range of values of $t$, you should choose $\sigma -{\sigma }_{0}$ small, in which case the series for $\stackrel{~}{f}\left(t\right)$ converges slowly; while for a smaller range of values of $t$, you can allow $\sigma -{\sigma }_{0}$ to be larger and obtain faster convergence.
The default value for $\sigma$ used by nag_sum_invlaplace_weeks (c06lb) is ${\sigma }_{0}+0.7$. There is no theoretical justification for this.

### Choice of b

The simplest advice for choosing $b$ is to set $b/2\ge \sigma -{\sigma }_{0}$. The default value used by the function is $2.5\left(\sigma -{\sigma }_{0}\right)$. A more refined choice is to set
 $b/2 ≥ minj σ-sj$
where ${s}_{j}$ are the singularities of $F\left(s\right)$.

## Example

This example computes values of the inverse Laplace transform of the function
 $Fs = 3 s2-9 .$
 $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$.
```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
``` 