# NAG FL Interfacec06lbf (invlaplace_​weeks)

## ▸▿ Contents

Settings help

FL Name Style:

FL Specification Language:

## 1Purpose

c06lbf computes the inverse Laplace transform $f\left(t\right)$ of a user-supplied 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.

## 2Specification

Fortran Interface
 Subroutine c06lbf ( f, b, mmax, m,
 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
#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.

## 3Description

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 , Re(s) > σ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)$.
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:
 $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 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:
 $| f(t)- 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.
c06lbf is derived from the subroutine MODUL1 in Garbow et al. (1988a).
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

## 5Arguments

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
 Function f ( s)
 Complex (Kind=nag_wp) :: f Complex (Kind=nag_wp), Intent (In) :: s
 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 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: $\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 pseudo-accuracy, that is, an upper bound on $|f\left(t\right)-\stackrel{~}{f}\left(t\right)|{e}^{-\sigma t}$.
6: $\mathbf{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: ${\mathbf{mmax}}=1024$ is sufficient for all but a few exceptional cases.
Constraint: ${\mathbf{mmax}}\ge 8$.
7: $\mathbf{m}$Integer Output
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) array Output
On exit: the first m elements contain the computed Laguerre expansion coefficients, ${a}_{i}$.
9: $\mathbf{errvec}\left(8\right)$Real (Kind=nag_wp) array Output
On exit: an $8$-component vector of diagnostic information.
${\mathbf{errvec}}\left(1\right)$
Overall estimate of the pseudo-error $|f\left(t\right)-\stackrel{~}{f}\left(t\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 routine for the purposes of error estimation. It satisfies
 $|ai| < K R -i , ​ i= 1, 2, …, m .$
10: $\mathbf{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 ${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit. When the value $-\mathbf{1}$ or $\mathbf{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).

## 6Error 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).
Errors or warnings detected by the routine:
Note: in some cases c06lbf may return useful information.
${\mathbf{ifail}}=1$
On entry, ${\mathbf{mmax}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{mmax}}\ge 8$.
${\mathbf{ifail}}=2$
The estimate of the pseudo-error is slightly larger than epstol. Pseudo-error estimate ${\mathbf{errvec}}\left(1\right)=⟨\mathit{\text{value}}⟩$ and ${\mathbf{epstol}}=⟨\mathit{\text{value}}⟩$.
${\mathbf{ifail}}=3$
The round-off error level is larger than epstol. Increasing epstol may help. Pseudo-error estimate ${\mathbf{errvec}}\left(1\right)=⟨\mathit{\text{value}}⟩$ and ${\mathbf{epstol}}=⟨\mathit{\text{value}}⟩$. Changing sigma or b may help. If not, the method should be abandoned.
${\mathbf{ifail}}=4$
The decay rate of the coefficients is too small. Increasing mmax may help. ${\mathbf{mmax}}=⟨\mathit{\text{value}}⟩$. Pseudo-error estimate ${\mathbf{errvec}}\left(1\right)=⟨\mathit{\text{value}}⟩$ and ${\mathbf{epstol}}=⟨\mathit{\text{value}}⟩$. Changing sigma or b may help. If not, the method should be abandoned.
${\mathbf{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. ${\mathbf{mmax}}=⟨\mathit{\text{value}}⟩$. Pseudo-error estimate ${\mathbf{errvec}}\left(1\right)=⟨\mathit{\text{value}}⟩$ and ${\mathbf{epstol}}=⟨\mathit{\text{value}}⟩$. Changing sigma or b may help. If not, the method should be abandoned.
${\mathbf{ifail}}=6$
Error bounds cannot be predicted. Check sigma0. ${\mathbf{sigma0}}=⟨\mathit{\text{value}}⟩$. Changing sigma or b may help. If not, the method should be abandoned.
${\mathbf{ifail}}=-99$
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{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.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

## 7Accuracy

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 $|f\left(t\right)-\stackrel{~}{f}\left(t\right)|{e}^{-\sigma t}$.

## 8Parallelism and Performance

c06lbf is not threaded in any implementation.

### 9.1The Role of ${\mathbit{\sigma }}_{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 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:
 $εnat (t) = ( f~ (t) - f (t) ) / e σ0t .$
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$.

### 9.2Choice of $\mathbit{\sigma }$

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 c06lbf is ${\sigma }_{0}+0.7$. There is no theoretical justification for this.

### 9.3Choice of $\mathbit{b}$

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
 $b/2 ≥ minj |σ-sj|$
where ${s}_{j}$ are the singularities of $F\left(s\right)$.

## 10Example

This example computes values of the inverse Laplace transform of the function
 $F(s) = 3 s2-9 .$
 $f(t) = sinh⁡3t .$
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.1Program Text

Program Text (c06lbfe.f90)

### 10.2Program Data

Program Data (c06lbfe.d)

### 10.3Program Results

Program Results (c06lbfe.r)