NAG FL Interface
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 ( f, 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 :: f
C Header Interface
#include <nag.h>
void  c06laf_ (
void (NAG_CALL *f)(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)
The routine may be called by the names c06laf or nagf_sum_invlaplace_crump.

3 Description

Given a function Fp 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))
f tj e atj τ 12 F a - k=1 Re F a + kπi τ cos kπtj τ - Im F a + kπi τ sin kπtj τ ,  
for j=1,2,,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 f0.0 can be obtained when it has a finite value. τ is calculated as tfac × max0.01,tj , where tfac > 0.5 . You specify tfac and αb , an upper bound on the exponential order α of the inverse function ft . α has two alternative interpretations:
  1. (i)α is the smallest number such that
    ft m × expαt  
    for large t,
  2. (ii)α is the real part of the singularity of Fp with largest real part.
The method depends critically on the value of α. See Section 9 for further details. The routine calculates at least two different values of the argument a, such that a>αb , in an attempt to achieve the requested relative error and provide error estimates. The values of tj , for j=1,2,,n, must be supplied in monotonically increasing order. The routine calculates the values of the inverse function ftj 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 emSn transformation Math. Tables Aids Comput. 10 91–96

5 Arguments

1: f Subroutine, supplied by the user. External Procedure
f must evaluate the real and imaginary parts of the function Fp for a given value of p.
The specification of f is:
Fortran Interface
Subroutine f ( pr, pi, fr, fi)
Real (Kind=nag_wp), Intent (In) :: pr, pi
Real (Kind=nag_wp), Intent (Out) :: fr, fi
C Header Interface
void  f_ (const double *pr, const double *pi, double *fr, double *fi)
1: pr Real (Kind=nag_wp) Input
2: pi Real (Kind=nag_wp) Input
On entry: the real and imaginary parts of the argument p.
3: fr Real (Kind=nag_wp) Output
4: fi Real (Kind=nag_wp) Output
On exit: the real and imaginary parts of the value Fp.
f 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: f should not return floating-point 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: n Integer Input
On entry: n, the number of points at which the value of the inverse Laplace transform is required.
Constraint: n1.
3: tn Real (Kind=nag_wp) array Input
On entry: each tj must specify a point at which the inverse Laplace transform is required, for j=1,2,,n.
Constraint: 0.0 t1 < t2 << tn.
4: valinvn Real (Kind=nag_wp) array Output
On exit: an estimate of the value of the inverse Laplace transform at t=tj, for j=1,2,,n.
5: errestn Real (Kind=nag_wp) array Output
On exit: an estimate of the error in valinvj. This is usually an estimate of relative error but, if valinvj<relerr, errestj estimates the absolute error. errestj is unreliable when valinvj is small but slightly greater than relerr.
6: 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.0relerr<1.0. If relerr is set too small or to 0.0, the routine uses a value sufficiently larger than machine precision.
7: alphab Real (Kind=nag_wp) Input
On entry: αb, an upper bound for α (see Section 3). Usually, αb should be specified equal to, or slightly larger than, the value of α. If αb < α then the prescribed accuracy may not be achieved or completely incorrect results may be obtained. If αb is too large c06laf will be inefficient and convergence may not be achieved.
Note:  it is as important to specify αb correctly as it is to specify the correct function for inversion.
8: tfac Real (Kind=nag_wp) Input
On entry: tfac, a factor to be used in calculating the parameter τ. Larger values (e.g., 5.0) may be specified for difficult problems, but these may require very large values of mxterm.
Suggested value: tfac=0.8.
Constraint: tfac>0.5.
9: mxterm Integer Input
On entry: the maximum number of (complex) terms to be used in the evaluation of the Fourier series.
Suggested value: mxterm100, except for very simple problems.
Constraint: mxterm1.
10: nterms Integer Output
On exit: the number of (complex) terms actually used.
11: na Integer Output
On exit: the number of values of a used by the routine. See Section 9.
12: 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 alphab- see Section 9.
13: 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 alphab- see Section 9.
14: nfeval Integer Output
On exit: the number of calls to f made by the routine.
15: work 4×mxterm+2 Real (Kind=nag_wp) array Workspace
16: 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 c06laf may return useful information.
ifail=1
On entry, mxterm=value.
Constraint: mxterm1.
On entry, n=value.
Constraint: n1.
On entry, relerr=value.
Constraint: 0.0relerr<1.0.
On entry, tfac=value.
Constraint: tfac>0.5.
ifail=2
On entry, t1=value.
Constraint: 0.0t1<t2<<tn.
On entry, the elements of t are not strictly increasing.
ifail=3
On entry, tn is too large: tn=value. If necessary, scale the problem as described in Section 9.
ifail=4
Required accuracy cannot be obtained. Try increasing tfac, alphab, or both. tfac=value and alphab=value.
ifail=5
Convergence failure in epsilon algorithm: mxterm=value. Some values of valinvj may be calculated to the desired accuracy; this may be determined by examining the values of errestj . Try reducing the range of t or increasing mxterm. If ifail=5 still results, try reducing tfac.
ifail=6
All values of valinv have been calculated, but are not of the required accuracy; the values of errestj should be examined carefully. Try reducing the range of t, or increasing tfac, alphab or both.
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 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 ft . If there is a discontinuity at t=c then the method converges such that f c fc- + fc+ / 2 .
Apparent loss of accuracy, when ft 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 ft is estimated to be less than relerr then this routine switches to absolute error estimation. However, when ft is slightly larger than relerr the relative error estimates are likely to cause ifail=6. If this is found inconvenient it can sometimes be avoided by adding k/p to the function Fp , which shifts the inverse to k+ft .
Loss of accuracy may also occur for highly oscillatory functions.
More serious loss of accuracy can occur if α is unknown and is incorrectly estimated. See Section 9.

8 Parallelism and Performance

c06laf is not threaded in any implementation.

9 Further Comments

9.1 Timing

The value of n is less important in general than the value of nterms. Unless f is very inexpensive to compute, the timing is proportional to na×nterms . For simple problems na=2 but in difficult problems na may be somewhat larger.

9.2 Precautions

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 α is positive. It is advisable, especially if ifail=3 is obtained, to scale the problem if α 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 tn 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 α. If, on entry, alphab is sufficiently smaller than α then completely incorrect results will be obtained with ifail=0. Unless α 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 tn , 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 t1 should be chosen very carefully and the following points should be borne in mind:
  1. (i) t1 should be small but not too close to 0.0 because of Gibbs' phenomenon (see Section 7),
  2. (ii)the larger the value of t1 , the smaller the range of values of a that will be used in the algorithm,
  3. (iii) t1 should ideally not be chosen such that ft1=0.0 or a very small value. For suitable problems t1 might be chosen as, say, 0.1 or 1.0 depending on these factors. The routine calculates alow from the formula
    alow = alphab - ln0.1×relerr 2×τ .  
Additional values of a are computed by adding 1/τ to the previous value. As τ = tfac × tn , 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 ifail=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 Fp=1/p+1/2 . The true inverse of Fp is exp-t/2 . Two preliminary calls to the routine are made to verify that the chosen value of alphab is suitable. For these tests the single value t1=1.0 is used. To test values of a close to alphab, the values tfac=5.0 and relerr=0.01 are chosen. To test larger a, the values tfac=0.8 and relerr=1.0E−3 are used. Because the values of the computed inverse are similar and ifail=0 in each case, these tests show that there is unlikely to be a singularity of Fp in the region -0.04 Rep 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)
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 1 1.5 2 2.5 3 3.5 4 4.5 5 1e−06 1e−05 0.0001 0.001 f(t) Error t Example Program Inverse Laplace Transform of 1/(p+0.5) gnuplot_plot_1 result gnuplot_plot_2 exp(t/2) gnuplot_plot_3 relative error gnuplot_plot_4 estimated error