NAG Library Routine Document
c06lbf (invlaplace_weeks)
1
Purpose
c06lbf computes the inverse Laplace transform
$f\left(t\right)$ of a usersupplied function
$F\left(s\right)$, defined for complex
$s$. The routine uses a modification of Weeks' method which is suitable when
$f\left(t\right)$ has continuous derivatives of all orders. The routine 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
c06lcf.
2
Specification
Fortran Interface
Integer, Intent (In)  ::  mmax  Integer, Intent (Inout)  ::  ifail  Integer, Intent (Out)  ::  m  Real (Kind=nag_wp), Intent (In)  ::  sigma0, epstol  Real (Kind=nag_wp), Intent (Inout)  ::  sigma, b  Real (Kind=nag_wp), Intent (Out)  ::  acoef(mmax), errvec(8)  Complex (Kind=nag_wp), External  ::  f 

C Header Interface
#include nagmk26.h
void 
c06lbf_ ( Complex (NAG_CALL *f)(const Complex *s), const double *sigma0, double *sigma, double *b, const double *epstol, const Integer *mmax, Integer *m, double acoef[], double errvec[], Integer *ifail) 

3
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
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)$.
c06lbf, along with its companion
c06lcf, 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:
where
${L}_{i}\left(x\right)$ is the Laguerre polynomial of degree
$i$. This routine computes the coefficients
${a}_{i}$ of the above Laguerre expansion; the expansion can then be evaluated for specified
$t$ by calling
c06lcf. You must supply the value of
${\sigma}_{0}$, and also suitable values for
$\sigma $ and
$b$: see
Section 9 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 routine will fail with an error exit if the method is not suitable for the supplied function $F\left(s\right)$.
The routine is designed to satisfy an accuracy criterion of the form:
where
${\epsilon}_{\mathit{tol}}$ is a usersupplied bound. The error measure on the lefthand side is referred to as the
pseudorelative
error, or
pseudoerror for short. Note that if
$\sigma >0$ and
$t$ is large, the absolute error in
$\stackrel{~}{f}\left(t\right)$ may be very large.
c06lbf is derived from the subroutine MODUL1 in
Garbow et al. (1988a).
4
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
5
Arguments
 1: $\mathbf{f}$ – Complex (Kind=nag_wp) Function, supplied by the user.External Procedure

f must return the value of the Laplace transform function
$F\left(s\right)$ for a given complex value of
$s$.
The specification of
f is:
Fortran Interface
Complex (Kind=nag_wp)  ::  f  Complex (Kind=nag_wp), Intent (In)  ::  s 

C Header Interface
#include nagmk26.h
Complex 
f (const Complex *s) 

 1: $\mathbf{s}$ – Complex (Kind=nag_wp)Input

On entry: the value of
$s$ for which
$F\left(s\right)$ must be evaluated. The real part of
s is greater than
${\sigma}_{0}$.
f must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
c06lbf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: f should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
c06lbf. If your code inadvertently
does return any NaNs or infinities,
c06lbf is likely to produce unexpected results.
 2: $\mathbf{sigma0}$ – Real (Kind=nag_wp)Input

On entry: the abscissa of convergence of the Laplace transform, ${\sigma}_{0}$.
 3: $\mathbf{sigma}$ – Real (Kind=nag_wp)Input/Output

On entry: the parameter
$\sigma $ of the Laguerre expansion. If on entry
${\mathbf{sigma}}\le {\sigma}_{0}$,
sigma is reset to
${\sigma}_{0}+0.7$.
On exit: the value actually used for $\sigma $, as just described.
 4: $\mathbf{b}$ – Real (Kind=nag_wp)Input/Output

On entry: 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)$.
On exit: the value actually used for $b$, as just described.
 5: $\mathbf{epstol}$ – Real (Kind=nag_wp)Input

On entry: the required relative pseudoaccuracy, that is, an upper bound on $\leftf\left(t\right)\stackrel{~}{f}\left(t\right)\right{e}^{\sigma t}$.
 6: $\mathbf{mmax}$ – IntegerInput

On entry: 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.
Suggested value:
${\mathbf{mmax}}=1024$ is sufficient for all but a few exceptional cases.
Constraint:
${\mathbf{mmax}}\ge 8$.
 7: $\mathbf{m}$ – IntegerOutput

On exit: the number of Laguerre expansion coefficients actually computed. The number of calls to
f is
${\mathbf{m}}/2+2$.
 8: $\mathbf{acoef}\left({\mathbf{mmax}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the first
m elements contain the computed Laguerre expansion coefficients,
${a}_{i}$.
 9: $\mathbf{errvec}\left(8\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: an
$8$component vector of diagnostic information.
 ${\mathbf{errvec}}\left(1\right)$
 Overall estimate of the pseudoerror $\leftf\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 pseudoerror.
 ${\mathbf{errvec}}\left(3\right)$
 Estimate of the truncation pseudoerror.
 ${\mathbf{errvec}}\left(4\right)$
 Estimate of the condition pseudoerror 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 routine for the purposes of error estimation. It satisfies
 10: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, because for this routine the values of the output arguments may be useful even if
${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit, the recommended value is
$1$.
When the value $\mathbf{1}\text{ or}1$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Note: c06lbf may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

On entry,  ${\mathbf{mmax}}<8$. 
 ${\mathbf{ifail}}=2$

The estimated pseudoerror 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.
 ${\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.
 ${\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$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
When
${\mathbf{ifail}}\ge {\mathbf{3}}$, changing
sigma or
b may help. If not, the method should be abandoned.
7
Accuracy
The error estimate returned in ${\mathbf{errvec}}\left(1\right)$ has been found in practice to be a highly reliable bound on the pseudoerror $\leftf\left(t\right)\stackrel{~}{f}\left(t\right)\right{e}^{\sigma t}$.
8
Parallelism and Performance
c06lbf is not threaded in any implementation.
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 routine, 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:
c06lbf uses the supplied value of
${\sigma}_{0}$ only in determining the values of
$\sigma $ and
$b$ (see
Sections 9.2 and
9.3); thereafter it bases its computation entirely on
$\sigma $ and
$b$.
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
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 c06lbf is ${\sigma}_{0}+0.7$. There is no theoretical justification for this.
The simplest advice for choosing
$b$ is to set
$b/2\ge \sigma {\sigma}_{0}$. The default value used by the routine is
$2.5\left(\sigma {\sigma}_{0}\right)$. A more refined choice is to set
where
${s}_{j}$ are the singularities of
$F\left(s\right)$.
10
Example
This example computes values of the inverse Laplace transform of the function
The exact answer is
The program first calls
c06lbf to compute the coefficients of the Laguerre expansion, and then calls
c06lcf to evaluate the expansion at
$t=0$,
$1$,
$2$,
$3$,
$4$,
$5$.
10.1
Program Text
Program Text (c06lbfe.f90)
10.2
Program Data
Program Data (c06lbfe.d)
10.3
Program Results
Program Results (c06lbfe.r)