Inverting the Laplace Transform by the Weeks Method

One of the most efficient methods for the numerical inversion of the Laplace transform is the Weeks method. It is based on a Laguerre expansion and bilinear transformations, and it can be implemented with the FFT. It is necessary that the transform F(s) can be computed at arbitrary points in the complex plane.

For a description of the method, see for example the paper by Garbow, Giunta, Lyness and Murli, "Software for an implementation of Weeks' method for the inverse of the Laplace Transform," published in ACM TOMS Vol. 14, pp. 163--170 (1988) or the author's paper Algorithms for Parameter Selection in the Weeks Method for Inverting the Laplace Transform, SIAM Journal on Scientific Computing, Vol. 21, pp. 111-128 (1999).

Please Note: The codes below were developed in MATLAB 5. Some may not work in MATLAB 6. For an update to MATLAB 6, please see the bottom of the page.

The author has written two versions of the Weeks method:

In addition, both of these functions call the following two functions The first function computes the Week's coefficients, a_n, and the second function evaluates a Laguerre series with Clenshaw's algorithm.

Example of the usage of weeks.m and weekse.m:

We compute the inverse of the transform F(s) = 1/sqrt(s^2+1), at various t. (The exact inverse is f(t) = J_0(t), which enables us to compute the actual error.) We first choose arbitrary parameters (sigma,b) = (1,1); below we will consider the choice of optimal parameters. 

>> F = '1./sqrt(s.^2+1)';
>> t = [0.1 1 10];
>> N = 16; sig = 1; b = 1; 
>> f = weeks(F,t,N,sig,b)'

f =

   9.9750e-01   7.6520e-01  -2.4501e-01

>> fex = bessel(0,t)

fex =

   9.9750e-01   7.6520e-01  -2.4594e-01

>> abserr = abs(f-fex)

abserr =

   1.4197e-07   2.5928e-08   9.2465e-04
If an error estimate is also required, rather use weekse.m: 

>> [f,est] = weekse(F,t,N,sig,b);
>> f'

ans =

   9.9750e-01   7.6520e-01  -2.4501e-01

>> est'

ans =

   4.3389e-07   1.0672e-06   8.6475e-03
Recall the actual errors computed above: 

abserr =

   1.4197e-07   2.5928e-08   9.2465e-04

Of course, by increasing N one can compute f(t) to high accuracy: 


>> N = 64;
>> f = weeks(F,t,N,sig,b)'

f =

   9.9750e-01   7.6520e-01  -2.4594e-01

>> abserr = abs(f-fex)

abserr =

   1.1102e-16            0   2.7273e-13
The Weeks method contains two free parameters sigma and b. In the example above the arbitrary choice (sigma,b) = (1,1) was made. The optimal choice of these parameters requires some trial-and-error, or the use of the algorithms described in the author's manuscript referenced above. These algorithms consist of two main functions:

The first function, wpar1.m, implements Algorithm 1. It is limited to transforms in Class A, and the user is expected to provide an optimal relationship b = b(sigma). For details, see the author's paper referenced above.

The function wpar1.m calls the following function

which computes an estimate for the total error, E, on the optimal curve b = b(sigma). Likewise, wpar2.m calls the functions the first of which computes an estimate of the truncation error, T, and the latter of which computes an estimate of the total error, E.

Example of usage of wpar1.m:

We continue the example above, but now we estimate the optimal (sigma,b) via Algorithm 1: 

>> F = '1./sqrt(s.^2+1)';
>> N = 16; t = 1;
>> bopt = 'sqrt(s^2+1)';
>> [so,bo] = wpar1(F,t,N,0,10,50,bopt)

so =

   4.3039e+00


bo =

   4.4186e+00

>> f = weeks(F,t,N,so,bo)  % use the computed sigma and b to compute f

f =

   7.6520e-01

>> fex = bessel(0,t)

fex =

   7.6520e-01

>> abserr = abs(f-fex)

abserr =

   5.5511e-16
With the arbitrary choice (sigma,b) = (1,1) made above, the Weeks method yielded an accuracy of 2.6e-8 with N = 16. With the parameters (sigma,b) = (4.3, 4.4) computed by Algorithm 1 the method manages an accuracy of 5.6e-16 which is essentially full precision.

Example of usage of wpar2.m:

We compute the inverse of the transform F(s) = exp(-4*sqrt(s)) at t = 1. We first use an arbitrary choice (sigma,b) = (0.7,1.75): 


>> F = 'exp(-4*sqrt(s))';
>> N = 32; t = 1;
>> sig = 0.7; b = 1.75;
>> f = weeks(F,t,N,sig,b)

f =

   2.0675e-02

>> fex = 2*exp(-4/t)/sqrt(pi*t^3)

fex =

   2.0667e-02

>> abserr = abs(f-fex)

abserr =

   8.1746e-06
Now we estimate the optimal (sigma, b) via Algorithm2, i.e., wpar2.m. (Note that Algorithm 1 cannot be applied to this transform, since it is not in class A.) 

>> [so,bo] = wpar2(F,t,N,0,30,30,50)

so =

   5.7678e+00


bo =

   1.8057e+01

>> f = weeks(F,t,N,so,bo)

f =

   2.0667e-02

>> abserr = abs(f-fex)

abserr =

   1.6578e-09
Observe the increase in accuracy.

The transforms considered above are, of course, toy problems. But in the author's paper Algorithms 1 and 2 have been applied to realistic problems taken from actual applications. (These examples were obtained from the paper by Dean Duffy, "Numerical Inversion of the Laplace Transform: Comparison of Three New Methods on Characteristic Problems from Applications," ACM TOMS Vol. 19 pp. 335--359 (1993).)

One of these applications is a transform that describes shock waves in diatomic chains. Algorithm 1 produces good estimates for the optimal (sigma,b) as may be observed from a color contour plot. The plot shows the actual relative error as a function of (sigma,b), along with the optimal hyperbola b = b(sigma) (white line), and the estimated value of (sigma,b) (white cross). Another example is taken from the theory of Timoshenko beams; see another color contour plot .

Update to MATLAB 6: weeks_update.zip Thanks to Patrick Kano, University of Arizona.

Last updated: May 27, 2004.