NAG FL Interface
c06lbf (invlaplace_​weeks)

Settings help

FL Name Style:


FL Specification Language:


1 Purpose

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

2 Specification

Fortran Interface
Subroutine c06lbf ( f, sigma0, sigma, b, epstol, mmax, m, acoef, errvec, ifail)
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 <nag.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)
The routine may be called by the names c06lbf or nagf_sum_invlaplace_weeks.

3 Description

Given a function f (t) of a real variable t, its Laplace transform F(s) is a function of a complex variable s, defined by
F (s) = 0 e-st f (t) dt ,   Re(s) > σ0 .  
Then f(t) is the inverse Laplace transform of F(s) . 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 F(s) .
c06lbf, along with its companion c06lcf, attempts to solve the following problem:
The method is a modification of Weeks' method (see Garbow et al. (1988a)), which approximates f(t) 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 routine computes the coefficients ai of the above Laguerre expansion; the expansion can then be evaluated for specified t by calling c06lcf. You must supply the value of σ0 , and also suitable values for σ and b: see Section 9 for guidance.
The method is only suitable when f(t) has continuous derivatives of all orders. For such functions the approximation f~ (t) is usually good and inexpensive. The routine will fail with an error exit if the method is not suitable for the supplied function F(s) .
The routine is designed to satisfy an accuracy criterion of the form:
| f(t)- 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.
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: f Complex (Kind=nag_wp) Function, supplied by the user. External Procedure
f must return the value of the Laplace transform function F(s) for a given complex value of s.
The specification of f is:
Fortran Interface
Function f ( s)
Complex (Kind=nag_wp) :: f
Complex (Kind=nag_wp), Intent (In) :: s
C Header Interface
Complex  f (const Complex *s)
1: s Complex (Kind=nag_wp) Input
On entry: the value of s for which F(s) must be evaluated. The real part of s is greater than σ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 floating-point 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: sigma0 Real (Kind=nag_wp) Input
On entry: the abscissa of convergence of the Laplace transform, σ0.
3: sigma Real (Kind=nag_wp) Input/Output
On entry: the parameter σ of the Laguerre expansion. If on entry sigmaσ0, sigma is reset to σ0+0.7.
On exit: the value actually used for σ, as just described.
4: b Real (Kind=nag_wp) Input/Output
On entry: the parameter b of the Laguerre expansion. If on entry b < 2(σ-σ0), b is reset to 2.5(σ-σ0).
On exit: the value actually used for b, as just described.
5: epstol Real (Kind=nag_wp) Input
On entry: the required relative pseudo-accuracy, that is, an upper bound on |f(t)-f~(t)| e-σt.
6: mmax Integer Input
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: mmax=1024 is sufficient for all but a few exceptional cases.
Constraint: mmax8.
7: m Integer Output
On exit: the number of Laguerre expansion coefficients actually computed. The number of calls to f is m/2+2.
8: acoef(mmax) Real (Kind=nag_wp) array Output
On exit: the first m elements contain the computed Laguerre expansion coefficients, ai.
9: errvec(8) Real (Kind=nag_wp) array Output
On exit: an 8-component vector of diagnostic information.
errvec(1)
Overall estimate of the pseudo-error |f(t)-f~(t)| e-σt =errvec(2) + errvec(3) + errvec(4).
errvec(2)
Estimate of the discretization pseudo-error.
errvec(3)
Estimate of the truncation pseudo-error.
errvec(4)
Estimate of the condition pseudo-error on the basis of minimal noise levels in function values.
errvec(5)
K, coefficient of a heuristic decay function for the expansion coefficients.
errvec(6)
R, base of the decay function for the expansion coefficients.
errvec(7)
Logarithm of the largest expansion coefficient.
errvec(8)
Logarithm of the smallest nonzero expansion coefficient.
The values K and R returned in errvec(5) and errvec(6) define a decay function KR-i constructed by the routine for the purposes of error estimation. It satisfies
|ai| < K R -i ,   ​ i= 1, 2, , m .  
10: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of -1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value -1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value -1 is recommended since useful values can be provided in some output arguments even when ifail0 on exit. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
Note: in some cases c06lbf may return useful information.
ifail=1
On entry, mmax=value.
Constraint: mmax8.
ifail=2
The estimate of the pseudo-error is slightly larger than epstol. Pseudo-error estimate errvec(1)=value and epstol=value.
ifail=3
The round-off error level is larger than epstol. Increasing epstol may help. Pseudo-error estimate errvec(1)=value and epstol=value. Changing sigma or b may help. If not, the method should be abandoned.
ifail=4
The decay rate of the coefficients is too small. Increasing mmax may help. mmax=value. Pseudo-error estimate errvec(1)=value and epstol=value. Changing sigma or b may help. If not, the method should be abandoned.
ifail=5
The decay rate of the coefficients is too small and round-off error is such that the required accuracy cannot be obtained. Increasing mmax or epstol may help. mmax=value. Pseudo-error estimate errvec(1)=value and epstol=value. Changing sigma or b may help. If not, the method should be abandoned.
ifail=6
Error bounds cannot be predicted. Check sigma0. sigma0=value. Changing sigma or b may help. If not, the method should be abandoned.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

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

8 Parallelism and Performance

c06lbf is not threaded in any implementation.

9 Further Comments

9.1 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 routine, 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 F(s) .
The value of σ0 is also relevant in defining a natural accuracy criterion. For large t, f(t) 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 .  
c06lbf uses the supplied value of σ0 only in determining the values of σ and b (see Sections 9.2 and 9.3); thereafter it bases its computation entirely on σ and b.

9.2 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 c06lbf is σ0+0.7 . There is no theoretical justification for this.

9.3 Choice of b

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

10 Example

This example computes values of the inverse Laplace transform of the function
F(s) = 3 s2-9 .  
The exact answer is
f(t) = sinh3t .  
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)
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 1e−10 1e−08 1e−06 0.0001 0.01 1 100 10000 1e+06 1e+08 0 1 2 3 4 5 1e−10 1e−08 1e−06 0.0001 0.01 1 f(t) Pseudo Error t Example Program Inverse Laplace Transform of 3/(s2-9) f(t) pseudo error gnuplot_plot_1 gnuplot_plot_2