NAG Library Routine Document
c06laf (invlaplace_crump)
1
Purpose
c06laf estimates values of the inverse Laplace transform of a given function using a Fourier series approximation. Real and imaginary parts of the function, and a bound on the exponential order of the inverse, are required.
2
Specification
Fortran Interface
Subroutine c06laf ( 
fun, n, t, valinv, errest, relerr, alphab, tfac, mxterm, nterms, na, alow, ahigh, nfeval, work, ifail) 
Integer, Intent (In)  ::  n, mxterm  Integer, Intent (Inout)  ::  ifail  Integer, Intent (Out)  ::  nterms, na, nfeval  Real (Kind=nag_wp), Intent (In)  ::  t(n), relerr, alphab, tfac  Real (Kind=nag_wp), Intent (Out)  ::  valinv(n), errest(n), alow, ahigh, work(4*mxterm+2)  External  ::  fun 

C Header Interface
#include nagmk26.h
void 
c06laf_ ( void (NAG_CALL *fun)(const double *pr, const double *pi, double *fr, double *fi), const Integer *n, const double t[], double valinv[], double errest[], const double *relerr, const double *alphab, const double *tfac, const Integer *mxterm, Integer *nterms, Integer *na, double *alow, double *ahigh, Integer *nfeval, double work[], Integer *ifail) 

3
Description
Given a function
$F\left(p\right)$ defined for complex values of
$p$,
c06laf estimates values of its inverse Laplace transform by Crump's method (see
Crump (1976)). (For a definition of the Laplace transform and its inverse, see the
C06 Chapter Introduction.)
Crump's method applies the epsilon algorithm (see
Wynn (1956)) to the summation in Durbin's Fourier series approximation (see
Durbin (1974))
for
$j=1,2,\dots ,n$, by choosing
$a$ such that a prescribed relative error should be achieved. The method is modified slightly if
$t=0.0$ so that an estimate of
$f\left(0.0\right)$ can be obtained when it has a finite value.
$\tau $ is calculated as
${t}_{\mathrm{fac}}\times \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(0.01,{t}_{j}\right)$, where
${t}_{\mathrm{fac}}>0.5$. You specify
${t}_{\mathrm{fac}}$ and
${\alpha}_{b}$, an upper bound on the exponential order
$\alpha $ of the inverse function
$f\left(t\right)$.
$\alpha $ has two alternative interpretations:
(i) 
$\alpha $ is the smallest number such that
for large $t$, 
(ii) 
$\alpha $ is the real part of the singularity of $F\left(p\right)$ with largest real part. 
The method depends critically on the value of
$\alpha $. See
Section 9 for further details. The routine calculates at least two different values of the argument
$a$, such that
$a>{\alpha}_{b}$, in an attempt to achieve the requested relative error and provide error estimates. The values of
${t}_{\mathit{j}}$, for
$\mathit{j}=1,2,\dots ,n$, must be supplied in monotonically increasing order. The routine calculates the values of the inverse function
$f\left({t}_{j}\right)$ in decreasing order of
$j$.
4
References
Crump K S (1976) Numerical inversion of Laplace transforms using a Fourier series approximation J. Assoc. Comput. Mach. 23 89–96
Durbin F (1974) Numerical inversion of Laplace transforms: An efficient improvement to Dubner and Abate's method Comput. J. 17 371–376
Wynn P (1956) On a device for computing the ${e}_{m}\left({S}_{n}\right)$ transformation Math. Tables Aids Comput. 10 91–96
5
Arguments
 1: $\mathbf{fun}$ – Subroutine, supplied by the user.External Procedure

fun must evaluate the real and imaginary parts of the function
$F\left(p\right)$ for a given value of
$p$.
The specification of
fun is:
Fortran Interface
Real (Kind=nag_wp), Intent (In)  ::  pr, pi  Real (Kind=nag_wp), Intent (Out)  ::  fr, fi 

C Header Interface
#include nagmk26.h
void 
fun (const double *pr, const double *pi, double *fr, double *fi) 

 1: $\mathbf{pr}$ – Real (Kind=nag_wp)Input
 2: $\mathbf{pi}$ – Real (Kind=nag_wp)Input

On entry: the real and imaginary parts of the argument $p$.
 3: $\mathbf{fr}$ – Real (Kind=nag_wp)Output
 4: $\mathbf{fi}$ – Real (Kind=nag_wp)Output

On exit: the real and imaginary parts of the value $F\left(p\right)$.
fun must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
c06laf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: fun should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
c06laf. If your code inadvertently
does return any NaNs or infinities,
c06laf is likely to produce unexpected results.
 2: $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of points at which the value of the inverse Laplace transform is required.
Constraint:
${\mathbf{n}}\ge 1$.
 3: $\mathbf{t}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: each ${\mathbf{t}}\left(\mathit{j}\right)$ must specify a point at which the inverse Laplace transform is required, for $\mathit{j}=1,2,\dots ,n$.
Constraint:
$0.0\le {\mathbf{t}}\left(1\right)<{\mathbf{t}}\left(2\right)<\cdots <{\mathbf{t}}\left(n\right)$.
 4: $\mathbf{valinv}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: an estimate of the value of the inverse Laplace transform at
$t={\mathbf{t}}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,n$.
 5: $\mathbf{errest}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: an estimate of the error in
${\mathbf{valinv}}\left(j\right)$. This is usually an estimate of relative error but, if
${\mathbf{valinv}}\left(j\right)<{\mathbf{relerr}}$,
${\mathbf{errest}}\left(j\right)$ estimates the absolute error.
${\mathbf{errest}}\left(j\right)$ is unreliable when
${\mathbf{valinv}}\left(j\right)$ is small but slightly greater than
relerr.
 6: $\mathbf{relerr}$ – Real (Kind=nag_wp)Input

On entry: the required relative error in the values of the inverse Laplace transform. If the absolute value of the inverse is less than
relerr, absolute accuracy is used instead.
relerr must be in the range
$0.0\le {\mathbf{relerr}}<1.0$. If
relerr is set too small or to
$0.0$, the routine uses a value sufficiently larger than
machine precision.
 7: $\mathbf{alphab}$ – Real (Kind=nag_wp)Input

On entry:
${\alpha}_{b}$, an upper bound for
$\alpha $ (see
Section 3). Usually,
${\alpha}_{b}$ should be specified equal to, or slightly larger than, the value of
$\alpha $. If
${\alpha}_{b}<\alpha $ then the prescribed accuracy may not be achieved or completely incorrect results may be obtained. If
${\alpha}_{b}$ is too large
c06laf will be inefficient and convergence may not be achieved.
Note: it is as important to specify ${\alpha}_{b}$ correctly as it is to specify the correct function for inversion.
 8: $\mathbf{tfac}$ – Real (Kind=nag_wp)Input

On entry:
${t}_{\mathrm{fac}}$, a factor to be used in calculating the parameter
$\tau $. Larger values (e.g.,
$5.0$) may be specified for difficult problems, but these may require very large values of
mxterm.
Suggested value:
${\mathbf{tfac}}=0.8$.
Constraint:
${\mathbf{tfac}}>0.5$.
 9: $\mathbf{mxterm}$ – IntegerInput

On entry: the maximum number of (complex) terms to be used in the evaluation of the Fourier series.
Suggested value:
${\mathbf{mxterm}}\ge 100$, except for very simple problems.
Constraint:
${\mathbf{mxterm}}\ge 1$.
 10: $\mathbf{nterms}$ – IntegerOutput

On exit: the number of (complex) terms actually used.
 11: $\mathbf{na}$ – IntegerOutput

On exit: the number of values of
$a$ used by the routine. See
Section 9.
 12: $\mathbf{alow}$ – Real (Kind=nag_wp)Output

On exit: the smallest value of
$a$ used in the algorithm. This may be used for checking the value of
${\mathbf{alphab}}$ see
Section 9.
 13: $\mathbf{ahigh}$ – Real (Kind=nag_wp)Output

On exit: the largest value of
$a$ used in the algorithm. This may be used for checking the value of
${\mathbf{alphab}}$ see
Section 9.
 14: $\mathbf{nfeval}$ – IntegerOutput

On exit: the number of calls to
fun made by the routine.
 15: $\mathbf{work}\left(4\times {\mathbf{mxterm}}+2\right)$ – Real (Kind=nag_wp) arrayWorkspace

 16: $\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: c06laf 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{n}}<1$, 
or  ${\mathbf{mxterm}}<1$, 
or  ${\mathbf{relerr}}<0.0$, 
or  ${\mathbf{relerr}}\ge 1.0$, 
or  ${\mathbf{tfac}}\le 0.5$. 
 ${\mathbf{ifail}}=2$

On entry,  ${\mathbf{t}}\left(1\right)<0.0$, 
or  ${\mathbf{t}}\left(1\right),{\mathbf{t}}\left(2\right),\dots ,{\mathbf{t}}\left({\mathbf{n}}\right)$ are not in strictly increasing order. 
 ${\mathbf{ifail}}=3$

${\mathbf{t}}\left({\mathbf{n}}\right)$ is too large for this value of
alphab. If necessary, scale the problem as described in
Section 9.
 ${\mathbf{ifail}}=4$

The required accuracy cannot be obtained. It is possible that
alphab is less than
$\alpha $. Alternatively, the problem may be especially difficult. Try increasing
tfac,
alphab or both.
 ${\mathbf{ifail}}=5$

Convergence failure in the epsilon algorithm. Some values of
${\mathbf{valinv}}\left(j\right)$ may be calculated to the desired accuracy; this may be determined by examining the values of
${\mathbf{errest}}\left(j\right)$. Try reducing the range of
t or increasing
mxterm. If
${\mathbf{ifail}}={\mathbf{5}}$ still results, try reducing
tfac.
 ${\mathbf{ifail}}=6$

All values of
${\mathbf{valinv}}\left(j\right)$ have been calculated but not all are to the requested accuracy; the values of
${\mathbf{errest}}\left(j\right)$ should be examined carefully. Try reducing the range of
$t$, or increasing
tfac,
alphab or both.
 ${\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.
7
Accuracy
The error estimates are often very close to the true error but, because the error control depends on an asymptotic formula, the required error may not always be met. There are two principal causes of this: Gibbs' phenomena, and zero or small values of the inverse Laplace transform.
Gibbs' phenomena (see the
C06 Chapter Introduction) are exhibited near
$t=0.0$ (due to the method) and around discontinuities in the inverse Laplace transform
$f\left(t\right)$. If there is a discontinuity at
$t=c$ then the method converges such that
$f\left(c\right)\to \left(f\left(c\right)+f\left(c+\right)\right)/2$.
Apparent loss of accuracy, when
$f\left(t\right)$ is small, may not be serious. Crump's method keeps control of relative error so that good approximations to small function values may appear to be very inaccurate. If
$\leftf\left(t\right)\right$ is estimated to be less than
relerr then this routine switches to absolute error estimation. However, when
$\leftf\left(t\right)\right$ is slightly larger than
relerr the relative error estimates are likely to cause
${\mathbf{ifail}}={\mathbf{6}}$. If this is found inconvenient it can sometimes be avoided by adding
$k/p$ to the function
$F\left(p\right)$, which shifts the inverse to
$k+f\left(t\right)$.
Loss of accuracy may also occur for highly oscillatory functions.
More serious loss of accuracy can occur if
$\alpha $ is unknown and is incorrectly estimated. See
Section 9.
8
Parallelism and Performance
c06laf is not threaded in any implementation.
The value of
$n$ is less important in general than the value of
nterms. Unless
fun is very inexpensive to compute, the timing is proportional to
${\mathbf{na}}\times {\mathbf{nterms}}$. For simple problems
${\mathbf{na}}=2$ but in difficult problems
na may be somewhat larger.
You are referred to the
C06 Chapter Introduction for advice on simplifying problems with particular difficulties, e.g., where the inverse is known to be a step function.
The method does not work well for large values of
$t$ when
$\alpha $ is positive. It is advisable, especially if
${\mathbf{ifail}}={\mathbf{3}}$ is obtained, to scale the problem if
$\left\alpha \right$ is much greater than
$1.0$. See the
C06 Chapter Introduction.
The range of values of $t$ specified for a particular call should not be greater than about $10$ units. This is because the method uses arguments based on the value ${\mathbf{t}}\left(n\right)$ and these tend to be less appropriate as $t$ becomes smaller. However, as the timing of the routine is not especially dependent on $n$, it is usually far more efficient to evaluate the inverse for ranges of $t$ than to make separate calls to the routine for each value of $t$.
The most important argument to specify correctly is
alphab, an upper bound for
$\alpha $. If, on entry,
alphab is sufficiently smaller than
$\alpha $ then completely incorrect results will be obtained with
${\mathbf{ifail}}={\mathbf{0}}$. Unless
$\alpha $ is known theoretically it is strongly advised that you should test any estimated value used. This may be done by specifying a single value of
$t$ (i.e
${\mathbf{t}}\left(n\right)$,
$n=1$) with two sets of suitable values of
tfac,
relerr and
mxterm, and examining the resulting values of
alow and
ahigh. The value of
${\mathbf{t}}\left(1\right)$ should be chosen very carefully and the following points should be borne in mind:
(i) 
${\mathbf{t}}\left(1\right)$ should be small but not too close to $0.0$ because of Gibbs' phenomenon (see Section 7), 
(ii) 
the larger the value of ${\mathbf{t}}\left(1\right)$, the smaller the range of values of $a$ that will be used in the algorithm, 
(iii) 
${\mathbf{t}}\left(1\right)$ should ideally not be chosen such that $f\left({\mathbf{t}}\left(1\right)\right)=0.0$ or a very small value. For suitable problems ${\mathbf{t}}\left(1\right)$ might be chosen as, say, $0.1$ or $1.0$ depending on these factors. The routine calculates alow from the formula

Additional values of
$a$ are computed by adding
$1/\tau $ to the previous value. As
$\tau ={\mathbf{tfac}}\times {\mathbf{t}}\left(n\right)$, it will be seen that large values of
tfac and
relerr will test for
$a$ close to
alphab. Small values of
tfac and
relerr will test for
$a$ large. If the result of both tests is
${\mathbf{ifail}}={\mathbf{0}}$, with comparable values for the inverse, then this gives some credibility to the chosen value of
alphab. You should note that this test could be more computationally expensive than the calculation of the inverse itself. The example program (see
Section 10) illustrates how such a test may be performed.
10
Example
This example estimates the inverse Laplace transform of the function
$F\left(p\right)=1/\left(p+1/2\right)$. The true inverse of
$F\left(p\right)$ is
$\mathrm{exp}\left(t/2\right)$. Two preliminary calls to the routine are made to verify that the chosen value of
alphab is suitable. For these tests the single value
${\mathbf{t}}\left(1\right)=1.0$ is used. To test values of
$a$ close to
alphab, the values
${\mathbf{tfac}}=5.0$ and
${\mathbf{relerr}}=0.01$ are chosen. To test larger
$a$, the values
${\mathbf{tfac}}=0.8$ and
${\mathbf{relerr}}=\text{1.0E\u22123}$ are used. Because the values of the computed inverse are similar and
${\mathbf{ifail}}={\mathbf{0}}$ in each case, these tests show that there is unlikely to be a singularity of
$F\left(p\right)$ in the region
$0.04\le \mathrm{Re}\left(p\right)\le 6.51$.
10.1
Program Text
Program Text (c06lafe.f90)
10.2
Program Data
Program Data (c06lafe.d)
10.3
Program Results
Program Results (c06lafe.r)