Inverting the Laplace Transform by a Direct Integration Method

Another method for inverting the Laplace transform is based on the trapezoidal approximation of the Bromwich integral. The resulting series may converge slowly and therefore a sequence acceleration method is used to improve the rate of convergence.

The method is discussed in, for example, the paper by Dean G. 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). The code given here is based on eq. (2) of this reference.

The code consists of two parts: The main function is direct.m which computes the terms in the series. It then calls the function wynneps.m, which accelerates the convergenge via Wynn's epsilon algorithm.

If you choose to download these codes, note that they should be saved as two separate text files, direct.m and wynneps.m.

Note that nonlinear sequence accelerators such as the epsilon algorithm are known to be unstable. No special provision was made to check for this nor to prevent this, so these codes should be used wisely. In particular, the codes should not be used if accuracy of better than 10(-8) or so is required, nor for large values of t.

Nonetheless, the code does remarkably well on a wide range of transforms, even ones that involve discontinuities such as step functions---see Example 2 below.

Examples of usage:

Example 1:

>> F = '1./(s.*sqrt(s+1))';
>> t = 2;
>> f = direct(F,100,t,10)

f   =  9.544997360895490e-01

>> fex = erf(sqrt(t))

fex =  9.544997361036416e-01
 
>> error = f-fex

         = -1.409261596307942e-11

Example 2:

>> F = '(exp(-s)+exp(-2*s))./(s.*(1-exp(-s)).^2)';
>> tt = 0.02:0.02:3;
>> ff = [];   
>> for t = tt;
f = direct(F,100,t,10); ff = [ff; f];
end
>> plot(tt,ff)

Last updated: February 5, 1998.