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) = \int_{0}^{\infty} e^{- s t} f (t) d t, 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:

given a function $F (s)$ , compute values of its inverse Laplace transform $f (t)$ for specified values of $t$ .

The method is a modification of Weeks' method (see Garbow et al. (1988a)), which approximates

f (t)

by a truncated Laguerre expansion:

\tilde{f} (t) = e^{σ t} \sum_{i = 0}^{m - 1} a_{i} e^{- b t / 2} L_{i} (b t), σ > σ_{0}, b > 0

where

L_{i} (x)

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

σ_{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

\tilde{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:

| \frac{f (t) - \tilde{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

\tilde{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 (

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 \leq σ_{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) - \tilde{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:

mmax \geq 8

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,

a_{i}

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) - \tilde{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

K R^{- i}

constructed by the routine for the purposes of error estimation. It satisfies

| a_{i} | < K R^{- i},   ​ i = 1, 2, \dots, m .

10: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

−1

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

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

ifail \neq 0

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

−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: $mmax \geq 8$ .

$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) - \tilde{f} (t) | e^{- σ t}

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.

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

\tilde{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^{σ_{0} t}

, so a natural measure of relative accuracy of the approximation

\tilde{f} (t)

is:

ε_{nat} (t) = (\tilde{f} (t) - f (t)) / e^{σ_{0} t} .

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

\tilde{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

\tilde{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 \geq σ - σ_{0}

. The default value used by the routine is

2.5 (σ - σ_{0})

. A more refined choice is to set

b / 2 \geq \min_{j} | σ - s_{j} |

where

s_{j}

are the singularities of

F (s)

10 Example

This example computes values of the inverse Laplace transform of the function

F (s) = \frac{3}{s^{2} - 9} .

The exact answer is

f (t) = \sinh 3 t .

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

c06lb: FL CL CPP AD PY MB

NAG FL Interfacec06lbf (invlaplace_​weeks)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 The Role of σ0

9.2 Choice of σ

9.3 Choice of b

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
c06lbf (invlaplace_weeks)

9.1 The Role of $σ_{0}$

9.2 Choice of $σ$

9.3 Choice of $b$