Numerical Algorithms manuscript No.

(will be inserted by the editor)

Review of Inverse Laplace Transform Algorithms for Laplace-Space

Numerical Approaches

Kristopher L. Kuhlman
arXiv:1204.4754v3 [math.NA] 20 Jul 2012

the date of receipt and acceptance should be inserted later

Abstract A boundary element method (BEM) simulation is used to compare the efficiency of numerical in-

verse Laplace transform strategies, considering general requirements of Laplace-space numerical approaches.

The two-dimensional BEM solution is used to solve the Laplace-transformed diffusion equation, producing

a time-domain solution after a numerical Laplace transform inversion. Motivated by the needs of numer-

ical methods posed in Laplace-transformed space, we compare five inverse Laplace transform algorithms

and discuss implementation techniques to minimize the number of Laplace-space function evaluations. We

investigate the ability to calculate a sequence of time domain values using the fewest Laplace-space model

evaluations. We find Fourier-series based inversion algorithms work for common time behaviors, are the

most robust with respect to free parameters, and allow for straightforward image function evaluation re-use

across at least a log cycle of time.

Keywords numerical Laplace transform inversion · boundary element method · 2D diffusion · Helmholtz

equation · Laplace-space numerical methods · groundwater modeling

1 1 Introduction

2 Simulation methods that are posed in Laplace-transformed space, then numerically inverted back to the

3 time domain (i.e., Laplace-space methods), are a viable alternative to the more standard use of finite

Repository Performance Department, Sandia National Laboratories

4100 National Parks Highway, Carlsbad, New Mexico, USA

Tel: +1 575-234-0084

E-mail: klkuhlm@sandia.gov
4 differences in time. We use the the two-dimensional boundary element method (BEM) as an example of

5 this type of approach, to solve the Laplace-transformed diffusion equation (i.e., the Yukawa or modified

6 Helmholtz equation). We investigate five numerical inverse Laplace transform methods and implementa-

7 tion approaches, namely the methods of [Stehfest(1970)], [Schapery(1962)], [Weeks(1966)], [Talbot(1979)],

8 and [de Hoog et al(1982)]. Naively implemented Laplace-space simulations can be more computationally

9 expensive than using finite differences in time, but they have the advantage of allowing evaluation at

10 any time, without evolving from an initial condition, and image function calculations are trivially par-

11 allelized across Laplace parameters [Davies and Crann(2002)]. When Laplace-space numerical models are

12 used in parameter estimation, hundreds or thousands of forward simulations may be required – making

13 forward model efficiency critical. Although parameter estimation may be done directly in Laplace space

14 [Barnhart and Illangasekare(2012)], choosing an efficient inversion strategy is important in most applica-

15 tions.

16 The Laplace transform has a long history of use to derive analytical solutions to diffusion and wave

17 problems (e.g., see list of citations by [Duffy(2004), pp. 191-220]). Often the analytical inverse trans-

18 form is too difficult to find or evaluate in closed form. A researcher then resorts to approximate ana-

19 lytical methods (e.g., [Hantush(1960), Sternberg(1969)]) or numerical inversion (e.g., [Malama et al(2009),

20 Mishra and Neuman(2010)]). Numerical methods can similarly benefit from the Laplace transform, convert-

21 ing the time-dependence of a differential equation to parameter dependence. Laplace-space finite-element

22 approaches have seen application to groundwater flow and solute transport (e.g., [Sudicky and McLaren(1992),

23 Morales-Casique and Neuman(2009)]), and Laplace-space BEM has also been used in groundwater applica-

24 tions (e.g., [Kythe(1995), §10.3] or [Liggett and Liu(1982), §10.1]). The Laplace transform analytic element

25 method [Kuhlman and Neuman(2009)] is a transient extension of the analytic element method. These dif-

26 ferent Laplace-space approaches may have diverse spatial solution strategies, but they have a common

27 requirement of effective Laplace transform numerical inversion algorithms. We couple a BEM model in the

28 Laplace domain with a numerical Laplace transform inversion routine, but our conclusions should be valid

29 for both gridded and mesh-free Laplace-space numerical methods. Any Laplace-space numerical approach

30 begins with determination of optimal Laplace parameter values. Then each image function evaluation is

31 computed from the simulation. The final step involves approximating the time-domain solution from the

32 vector of image function values using the algorithm of choice.

2
33 [Bellman et al(1966)] was an early review book on numerical Laplace transform inversion for linear

34 and non-linear problems, but without the benefit of the many algorithms that have since been devel-

35 oped. [Davies and Martin(1979)] performed a thorough survey, assessing numerical Laplace transform in-

36 version algorithm accuracy for techniques available in 1979, using simple functions for their benchmarks.

37 [Duffy(1993)] reviewed the numerical inversion characteristics for more pathological time behaviors using

38 the Fourier series, Talbot, and Weeks inversion methods. The review book by [Cohen(2007)] summarizes

39 historical reviews and discusses commonly used inversion their variations. More details and examples can

40 be found in these reference regarding the convergence behavior of the five inversion algorithms discussed

41 here.

42 While these published numerical inverse Laplace transform algorithm reviews are thorough and useful,

43 they focus on computing a single time-domain solution as accurately as possible. These reviews did not try

44 to minimize Laplace-space function evaluations, since their functions were simple closed-form expressions,

45 not simulations. We investigate Laplace transform inversions algorithms that can compute a sequence of

46 time domain values using the fewest Laplace-space model evaluations possible, a desirable property for use

47 in Laplace-space numerical methods. Using numerical Laplace transform inversion in a simulation approach,

48 rather than a time-marching method, allows the researcher to readily switch between fast and accurate by

49 changing the number of approximation terms in the inversion.

50 In the next section we define the mathematical formulation of the governing equation and Laplace

51 transform. In the third section we introduce the five inverse Laplace transform algorithms. In the final

52 section we compare results using five different inversion algorithms to invert the BEM modified Helmholtz

53 solution on the same domain with four different boundary conditions, leading to recommendations for

54 Laplace-space numerical approaches.

55 2 Governing Equation and Laplace Transform

56 The BEM model generally simulates substance flow (e.g., energy or groundwater), which can be related to

57 a potential φ (e.g., temperature or hydraulic head). The medium property α is diffusivity [L2 /T], the ratio

58 of the conductance in the substance flux and potential gradient relation (e.g., Fourier’s or Darcy’s law)

59 to the substance capacity per unit mass (e.g., heat capacity or storativity). The BEM (e.g., [Kythe(1995),

3
60 Liggett and Liu(1982), Brebbia et al(1984)]) can be used to solve the diffusion equation

1 ∂φ
∇2 φ = , (1)
α ∂t

61 where α is a real constant in space and time. We consider (1) in a domain subject to a combination of

62 Dirichlet φ (Γu (s), t) = fu (s, t) and Neumann n̂ · ∇φ (Γq (s), t) = fq (s, t) boundary conditions along the

63 perimeter of the 2D domain Γ = Γu ∪ Γq , where n̂ is the boundary unit normal, and s is a boundary

64 arc-length parameter. Without loss of generality, we only consider homogeneous initial conditions.

65 The Laplace transform is

Z ∞
L{f (t)} ≡ f¯(p) = f (t)e−pt dt, (2)
0

66 where p is the generally complex-valued Laplace parameter, and the over-bar denotes a transformed variable.

67 The transformed diffusion equation with zero initial conditions is the homogeneous Yukawa or modified

68 Helmholtz equation,

∇2 φ − q 2 φ = 0 , (3)

69 where q 2 = p/α. Equation 3 arises in several groundwater applications, including transient, leaky, and lin-

70 earized unsaturated flow [Bakker and Kuhlman(2011)]. The transformed boundary conditions are φ̄ (Γu (s)) =

71 fu (s)f¯t (p) and n̂ · ∇φ̄(Γq (s)) = fq (s)f¯t (p), where the temporal and spatial behaviors have been decomposed

72 as in separation of variables. Arbitrary time behavior can be developed through convolution in t (Duhamel’s

73 theorem), which is multiplication of image functions in Laplace space. Here, f¯t (p) represents the Laplace

74 transform of the time behavior applied to the boundary conditions. The Laplace transformation makes it

75 possible to solve transient diffusion (a parabolic equation) using the BEM, which is well-suited for elliptical-

76 type equations.

77 The inverse Laplace transform is defined as the Bromwich contour integral,

Z σ+i∞
1
L−1 f¯(p) = f (t) = f¯(p)ept dp,

(4)
2πi σ−i∞

78 where the abscissa of convergence σ > 0 is a real constant chosen to put the contour to the right of all

79 singularities in f¯(p). In Laplace-space numerical approaches, (3) is solved by a suitable numerical method,

80 therefore only samples of f¯(p) are available; this precludes an analytical inversion. Five numerical inverse

81 Laplace transform algorithms are discussed in the following section.

4
82 3 Numerical Inverse Laplace Transform Methods

83 Equation 4 is an integral equation for unknown f (t) given f¯(p); its numerical solution is broadly split into

84 two categories. Methods are either based on quadrature or functional expansion using analytically invertible

85 basis functions. [Davies(2005), Chap. 19] relates most major classes of inverse Laplace transform methods

86 using a unified theoretical foundation; we adopt a simplified form of their general notation. The Fourier

87 series and Talbot methods are quadrature-based, directly approximating (4). Weeks’ and Piessen’s methods

88 are f¯(p) expansions using complex-valued basis functions, while the Gaver-Stehfest and Schapery methods

89 use real-valued functions to accomplish this.

90 The numerical inverse Laplace transform is in general an ill-posed problem (e.g., [Al-Shuaibi(2001)]).

91 No single approach is optimal for all circumstances and temporal behaviors, leading to the diversity of

92 viable numerical approaches in the literature (e.g., [Cohen(2007)]).

93 3.1 Gaver-Stehfest Method

94 The Post-Widder formula [Widder(1941), Al-Shuaibi(2001)] is an approximation to (4) that only requires

95 f¯(p) for real p to represent (2) as an asymptotic Taylor series expansion. The formula requires high-order

96 analytic image function derivatives, and is impractical for numerical computation. Stehfest proposed a

97 discrete version of the Post-Widder formula using finite differences and Salzer summation [Stehfest(1970)],

N  
ln 2 X ¯ ln 2
f (t, N ) = Vk f k . (5)
t t
k=1

98 The Vk coefficients only depend on the number of expansion terms, N (which must be even), which are

min(k,N/2) N
j 2 (2j )!
Vk = (−1)k+N/2
X
. (6)
j =⌊(k+1)/2⌋
(N
2 − j )! j ! (j − 1)! (k − j )! (2j − k)!

99 These become very large and alternate in sign for increasing k. The sum (5) begins to suffer from cancellation

100 for N ≥ the number of decimal digits of precision (e.g., double precision = 16). For f¯t (p) that are non-

101 oscillatory and continuous, N ≤ 18 is usually sufficient [Stehfest(1970)]. If programmed using arbitrary

102 precision (e.g. Mathematica or a multi-precision library [Bailey et al(2002), Johansson(2011)]), the method

103 can be made accurate for most cases [Abate and Valkó(2004)]. Unfortunately, p is explicitly a function of

104 t; for each new t, a new f¯(p) vector is needed. In Laplace-space numerical approaches, each vector element

105 is constructed using a simulation, therefore this can be a large penalty.

5
106 The method is quite easy to program; the Vj can be computed once and saved as constants. This method

107 has been popular due to its simplicity and adequacy for exponentially decaying f¯t (p).

108 3.2 Schapery’s Method

109 We can expand the deviation of f (t) from steady-state fs using exponential basis functions [Schapery(1962)],

N
ai e−pi t ,
X
f (t, N ) = fs + (7)
i=1

110 where ai is a vector of unknown constants. Applying (2) to (7) gives

N
fs X ai
f¯(pj , N ) = + j = 1, 2, . . . , M. (8)
pj pi + pj
i=1

111 The pj are selected (a geometric series is recommended [Liggett and Liu(1982)]) to cover the important

112 fluctuations in f¯(p). After setting pi = pj the ai coefficients can be determined as the solution to Pij ai =

f¯(pj ) − fs /pj . The symmetric matrix to decompose is

113

 
−1 −1 −1
 (2p1 ) (p1 + p2 ) . . . (p1 + pN ) 
 
 
 ( p + p ) −1 (2p2)−1 ... −1 
(p2 + pN ) 
 2 1
Pij =
 ,
.. .. .. .. 

 . . . . 

 
 
( p N + p 1 ) −1 ( p N + p 2 ) −1 ... (2pN ) −1

114 which only depends on pj ; it can be decomposed independently of f¯(p) and fs .

115 This method is not difficult to implement when existing matrix decomposition libraries are available, and

116 only requires real computation. The method has been used for inverting BEM results [Liggett and Liu(1982)],

117 but has two main drawbacks. First, in its formulation above, it requires a steady-state solution, but (7)

118 could be posed without fs . Secondly, no theory is presented for optimally picking pj ; some trial and error

119 is required [Liggett and Liu(1982)].

120 3.3 Möbius Transformation Methods

121 We can use the Möbius transformation to conformally map the half-plane right of σ to the unit disc, mak-

122 ing the Laplace domain more amenable to approximation using orthonormal polynomials (e.g., Chebyshev

123 [Piessens(1972)], [Lanczos(1988), §28] or Laguerre [Weeks(1966)], [Lyness and Giunta(1986)], [Lanczos(1988),

124 §30]). If σ was chosen properly, f¯(p) is guaranteed to be analytic inside the unit circle. The most-used inverse

6
125 Laplace transform method from this class is Weeks’ method, which uses a complex power series to expand

126 f¯(p) inside the unit circle. Upon inverse Laplace transformation, the power series becomes a Laguerre

127 polynomial series.

128 Weeks method is

N
f (t, N + 1) = e(κ−b/2)t
X
an Ln (bt), (9)
n=0

129 where Ln (z ) is an n-order Laguerre polynomial and κ and b are free parameters. Weeks suggested κ =

130 σ + 1/tmax and b = N/tmax , where tmax is the maximum transformed time. The parameters b and κ are

131 chosen to optimize convergence; some schemes are given [Weideman(1999)] for finding optimum parameter

132 values for a given f¯t (p), but search techniques require hundreds of f¯(p) evaluations. A more general form of

133 (9) can also be used, which allows for more general asymptotic behavior of the image function [Davies(2005),

134 §19.5]. Weeks assumed pf¯(p) is analytic at infinity. The Laplace transform of (9) is known, but to make it

135 easier to represent with polynomials, f¯(p) is mapped inside the unit circle via z = (p − κ − 2b)/(p − κ + 2b).

136 The coefficients an are determined by the midpoint rule,

M−
X1
1 h  i  
an = Ψ exp iθj− 1 exp −inθj− 1 (10)
2M 2 2
j =−M

137 where θj = jπ/M and the conformally-mapped image function is

 
b ¯ b b
Ψ (z ) = f κ− + . (11)
1−z 2 1−z

138 The argument of f¯(z ) in (11) is the inverse mapping of z 7→ p, it shows p does not functionally depend on

139 t, but Weeks’ rules-of-thumb for b and κ depend on tmax .

140 There are other related methods which use different orthonormal polynomials to represent f¯(p) inside

141 the unit circle. Chebyshev polynomials (known as Piessen’s method [Piessens(1972)]) can be used to expand

142 the f¯(z ) on the real interval [−1, 1]. The Weeks method is moderately easy to program, requiring the use

143 of Clenshaw recurrence formula to accurately implement Laguerre polynomials. Piessen’s method is similar

144 to implement, with a similar recurrence formula for Chebyshev polynomials.

145 3.4 Talbot Method

146 We can deform the Bromwich contour into a parabola around the negative real axis if f¯(p) is analytic

147 in the region between the Bromwich and the deformed Talbot contours [Talbot(1979)]. Numerically, f¯(p)

7
148 must not overflow as p → −∞ (e.g., in the BEM implementation, the Green’s function is the second-kind

149 modified Bessel function, which grows exponentially as p → −∞). Oscillatory f¯t (p) often have pairs of poles

150 near the imaginary p axis; these poles must remain to the left of the deformed contour.

151 The Talbot method makes the Bromwich contour integral converge rapidly, since p becomes large and

152 negative, making the ept term in (4) very small. A one-parameter “fixed” Talbot method was implemented

153 [Abate and Valkó(2004)]; the Bromwich contour is parametrized as p(θ) = rθ(cot(θ) + i), where 0 ≤ θ ≤ π ,

154 and as a rule of thumb r = 2M/(5tmax ). The fixed Talbot method is

N −1
" #
r f¯(r ) rt
X n
tp(θk ) ¯
o
f (t, N ) = e + ℜ e f [p(θk )] [1 + iζ (θk )] , (12)
N 2
k=1

155 where ζ (θk ) = θk + [θk cot(θk ) − 1] cot(θk ) and θk = kπ/N [Abate and Valkó(2004)]. Although f¯(p) doesn’t

156 depend on t, the free parameter r depends on tmax .

157 Step change f¯t (p) for non-zero time become very large as p → −∞, since L [H (t − τ )] = e−τ p /p, where

158 H (t − τ ) is the Heaviside step function centered on time τ . This can lead to precision loss, and stability

159 or convergence issues with the underlying numerical model, although Mathematica’s arbitrary precision

160 capabilities have been used to get around this problem [Abate and Valkó(2004)].

161 The fixed Talbot method is very simple to program; [Abate and Valkó(2004)] provide a ten-line Math-

162 ematica implementation.

163 3.5 Fourier Series Method

164 We can manipulate (4) into a Fourier transform; first it is expanded into real and imaginary parts (p =

165 γ + iω ),
∞
eγt
Z
[cos(ωt) + i sin(ωt)] ℜ f¯(p) + iℑ f¯(p) i dω.
    
f ( t) =
2πi −∞

166 Multiplying out the terms, keeping only the real part due to f (t) symmetry, and halving the integration

167 range due to symmetry again, leaves

∞
eγt
Z
ℜ f¯(p) cos(ωt) − ℑ f¯(p) sin(ωt) dω.
   
f ( t) = (13)
π 0

168 When f (t) is real, (13) can be represented using the complex form or just its real or imaginary parts.

169 Although these three representations are equivalent, when evaluating (13) with the trapezoid rule, the full

8
170 complex form gives the smallest discretization error [Davies(2005)]. The trapezoid rule approximation to

171 (13) is essentially a discrete Fourier transform,

2N
eγt X ′
    
iπk iπt
f (t, 2N + 1) = ℜ f¯ γ0 + exp , (14)
T T T
k=0

172 where γ0 = σ − log(ǫ)/T , ǫ is the desired relative accuracy (typically 10−8 in double precision), T is a scaling

173 parameter (often 2tmax ), and the prime indicates the k = 0 summation term is halved. The p in (14) do

174 not depend on t, but the free parameter T depends on tmax .

175 The non-accelerated Fourier series inverse algorithm 14 is almost useless because it requires thousands

176 of f¯(p) evaluations [Antia(2002), §9.8]. Practical approaches accelerate the convergence of the sum in 14.

177 Although this is sometimes called a fast-Fourier transform (FFT) method (e.g., [Cohen(2007), Chap 4.4]),

178 rarely do the number of f¯(p) evaluations in an accelerated approach justify an FFT approach. The method

179 implemented uses non-linear double acceleration with Padé approximation and an analytic expression for

180 the remainder in the series [de Hoog et al(1982)]. Although there are several other ways to accelerate the

181 Fourier series approach [Cohen(2007)], this method is popular and straightforward. Non-linear acceleration

182 techniques drastically reduce the required number of function evaluations, but can lead to numerical dis-

183 persion [Kano et al(2005), Morales-Casique and Neuman(2009)]. For diffusion, dispersion associated with

184 non-linear acceleration is not noticeable. Schapery’s, Talbot’s, and Weeks’ methods are not accelerated in

185 a non-linear manner, and therefore may lead to less numerical dispersion, which may be more important

186 in wave systems with sharp fronts.

187 The creation of the Padé approximation [de Hoog et al(1982)] is relatively straightforward in program-

188 ming languages that facilitate matrix manipulations (e.g., modern Fortran, Matlab, or NumPy [Oliphant(2007)]).

189 There is no dependence on matrix decomposition routines.

190 3.6 Algorithm Properties Summary

191 Table 1 summarizes aspects of the five inverse methods. The third column indicates whether p is explicitly

192 a function of t, the fourth column indicates if the rules-of-thumb used for the optimum parameters depend

193 on tmax , and the fifth column indicates whether the transform requires complex p and f¯(p).

194 For all methods considered here, computational effort to compute f (t) from the vector of f¯(p) values

195 was insignificant compared to the effort required to compute the BEM solution used to fill the f¯(p) vector.

9
 Table 1 Algorithmic Summary

Method Limitations on f¯(p) and f (t) p(t)? p(tmax )? p

Stehfest no oscillations, no discontinuities in f (t) yes no real

Schapery smoothly varying f (t), fs exists no no real

Weeks none no yes complex

fixed Talbot no high-frequency f (t), f¯(p → −∞) exists no yes complex

Fourier series none no yes complex

196 This suggests a more complicated method, which allows re-use of f¯(p) across more values of t and converges

197 in less evaluations of f¯(p), would be efficient for Laplace-domain numerical methods. If existing libraries

198 or simulations only support real arguments, then the Stehfest, Schapery, or Piessen’s methods must be

199 used. Complex p methods will pay a slight penalty in computational overhead compared to real-only p

200 routines. Computing with arbitrary or higher-than-double precision (e.g., [Abate and Valkó(2004)]) will

201 incur a much larger penalty than the change from real to complex double precision. Generally, complex

202 p methods have better convergence properties than real-only methods. Expansion of f¯(p) along the real p

203 axis is separation of non-orthogonal exponentials, while expansion along the imaginary p axis is separation

204 of oscillatory functions [Lanczos(1988), §29].

205 4 Numerical Comparison

206 Four test problems were solved using the BEM for values of p required by each algorithm’s rules of thumb.

207 The test problem domain is a 3 × 2 rectangle, with homogeneous initial conditions and specified potential

208 at two ends φ̄(x = 0) = −2f¯t (p), and φ̄(x = 3) = 2f¯t (p), and zero normal flux along the other sides

209 ∂ φ̄/∂y (y = {0, 2}) = 0. All plots show the solution computed at a point closer to the x = 0 boundary

210 (x = 1/3), midway between the insulated boundaries (y = 1).

211 The first problem computes f¯(p) using the optimum p at each t (like most inverse Laplace transform

212 surveys), according to the rules-of-thumb for each method. While this is most accurate, it is very inefficient

213 – especially when many values of t are required. In the following sections, all methods except Stehfest

214 use the same f¯(p) to invert all t. A method’s sensitivity to non-optimal free parameters is important in

215 practical use for Laplace-space numerical approaches. By inverting more than one time with the same set

216 of Laplace-space function evaluations, large gains in efficiency can be made. The t range used in the plots

10
 0.5
Fourier series, M=5
0 Talbot, M=5
0
Stehfest, N=6
Schapery, M=5 + SS
-0.5 Weeks, M=5
-0.5
potential

-1

x-flux
-1 -1.5
Fourier series, M=5
Talbot, M=5
Stehfest, N=6 -2
-1.5 Schapery, M=5 + SS
Weeks, M=5 -2.5
FD solution
-2 -3
0.01 0.1 1 0.01 0.1 1
time time

Fig. 1 Plots of potential and flux through time with five methods for f¯t (p) = 1/p, using optimum p at each t. 15 × 5 = 75

total f¯(p) evaluations are used by each method.

217 spans three orders of magnitude; it was chosen to show the evolution of potential and substance flux from

218 initial conditions to steady state.

219 4.1 Steady Boundary Conditions, Optimum p

220 The first problem has steady-state boundary conditions. The transient behavior is solely due to evolution

221 from the zero initial condition, f¯t (p) = L[H (t)] = 1/p; f¯(p) has a pole at the origin. All methods performed

222 equally well with this simple test problem, although the Fourier series method deviates from the finite

223 difference solutions at larger time. Figure 1 shows the inverted potential and flux using as few evaluations

224 of f¯(p) possible, without major deviations from the finite difference benchmark solution. Some trial and

225 error was needed to use the Schapery method (i.e., further optimization may be possible).

226 As shown in Figure 2, all the methods performed very well for nine f¯(p) terms per t but at least 135 f¯(p)

227 evaluations are needed total for each method. Schapery’s method does the worst in this case, but this may

228 be improved with further optimization of pj values. The finite-difference approach took at least an order

229 of magnitude less computational effort for the given accuracy. Making Laplace-space numerical methods

230 useful alternatives to traditional time-marching approaches, requires improvements to this inefficiency.

231 4.2 Steady Boundary Conditions, Same p

232 All methods had more difficulty obtaining accurate results for a wide t range using only one vector of f¯(p)

233 (no Stehfest method, since p explicitly depends on t). Only the last log-cycle of times is inverted accurately

234 when using nine f¯(p) (Figure 3). All the methods – except possibly Schapery’s – have a more difficult

11
 0.5
Fourier series, M=9
0 Talbot, M=9
0
Stehfest, N=10
Schapery, M=9 + SS
-0.5 Weeks, M=9
-0.5
potential

-1

x-flux
-1 -1.5
Fourier series, M=9
Talbot, M=9
Stehfest, N=10 -2
-1.5 Schapery, M=9 + SS
Weeks, M=9 -2.5
FD solution
-2 -3
0.01 0.1 1 0.01 0.1 1
time time

Fig. 2 Plots of potential and flux through time with five methods for f¯t (p) = 1/p, using optimum p at each t. 15 × 9 = 135

total f¯(p) evaluations are used by each method. Fourier series, Talbot, Stehfest, and Weeks curves are nearly coincident.

0.5
Fourier series, M=9
0 0 Talbot, M=9
Schapery, M=9 + SS
-0.5 Weeks, M=9
-0.5
-1
potential

x-flux

-1.5
-1
Fourier series, M=9 -2
Talbot, M=9
Schapery, M=9 + SS -2.5
-1.5 Weeks, M=9
FD solution -3
-2 -3.5
0.01 0.1 1 0.01 0.1 1
time time

Fig. 3 Plots of potential and flux through time with four methods for f¯t (p) = 1/p, same p used across all t. Nine total

f¯(p) evaluations are used by each method.

235 time with the flux at early time (especially the fixed Talbot method). The apparent success of Schapery’s

236 method can be attributed to the expansion of the deviation from steady-state, which in this case decays

237 exponentially with time.

238 Figure 4 shows that when increasing to 51 f¯(p) terms, most convergence problems disappear, except

239 at small times. Grouping t values by log-cycles and inverting them together using the same f¯(p) is more

240 economical than using the optimal p for each t and is still relatively accurate. The results shown in Figure 4

241 are nearly as accurate as those shown in Figure 2, but required 1/3 the f¯(p) model evaluations.

12
 0.5
Fourier series, M=51
0 0 Talbot, M=51
Schapery, M=51 + SS
-0.5 Weeks, M=51
-0.5
-1
potential

x-flux
-1.5
-1
Fourier series, M=51 -2
Talbot, M=51
Schapery, M=51 + SS -2.5
-1.5 Weeks, M=51
FD solution -3
-2 -3.5
0.01 0.1 1 0.01 0.1 1
time time

Fig. 4 Plots of potential and flux through time with four methods for f¯t (p) = 1/p, same p used across all t. 51 total f¯(p)

evaluations used by each method. Schapery and Weeks curves are nearly coincident.

1.5 3
Fourier series, M=19 Fourier series, M=19
1 Talbot, M=19 2 Talbot, M=19
Schapery, M=19 + SS Schapery, M=19 + SS
Weeks, M=19 Weeks, M=19
0.5 1
potential

0 0
x-flux

-0.5 -1

-1 -2

-1.5 -3

-2 -4
0.01 0.1 1 0.01 0.1 1
time time

Fig. 5 Plots of potential and flux through time with four methods for f¯t (p) = p/(p2 + 16), same p used across all t. 19

total f¯(p) evaluations are used by each method.

242 4.3 Sinusoidal Boundary Conditions, Same p

p
243 This problem uses temporally sinusoidal boundary conditions, f¯t (p) = L(cos 4t) = p2 +16 . This boundary

244 condition violates some assumptions of the inverse transform algorithms (i.e., no steady-state solution,

245 oscillatory in time), but the behavior is still relatively simple and smooth, with singularities at p = ±4i.

246 Figure 5 shows the Schapery method fails since there is no fs , but the other methods do well for 19

247 terms across one t log cycle. Figure 6 shows all methods besides Schapery do well for 51 terms, across
pj
248 at least two t log cycles. A modified version of (8) substituting p2j +16
for fs /pj could extend Schapery’s

249 approach to this case, but this solution was not considered here because of its problem specificity.

13
 1.5 3
Fourier series, M=51 Fourier series, M=51
Talbot, M=51 Talbot, M=51
1 Schapery, M=51 + SS 2 Schapery, M=51 + SS
Weeks, M=51 Weeks, M=51
0.5 1
potential

x-flux
0
0
-1
-0.5
-2
-1
-3
-1.5
0.01 0.1 1 0.01 0.1 1
time time

Fig. 6 Plots of potential and flux through time with four methods for f¯t (p) = p/(p2 + 16), same p used across all t. 51 total

f¯(p) evaluations are used by each method. Fourier series, Talbot, and Weeks curves are nearly coincident for potential.

0.5 3
Fourier series, M=19
2 Talbot, M=19
0 Schapery, M=19 + SS
Weeks, M=19
1
-0.5
potential

0
x-flux

-1 -1
Fourier series, M=19
Talbot, M=19 -2
-1.5 Schapery, M=19 + SS
Weeks, M=19 -3
FD solution
-2 -4
0.01 0.1 1 0.01 0.1 1
time time

Fig. 7 Plots of potential and flux through time with four methods for f¯t (p) = exp(−0.08p)/p, same p used across all t. 19

total f¯(p) evaluations are used by each method. Weeks’ solution is undefined for t < 0.08.

250 4.4 Step-Change Boundary Condition for τ > 0, same p

251 Finally, the same domain was simulated but with step-change boundary conditions at τ = 0.08, or f¯t (p) =

252 L(H (t − 0.08)) = e−0.08p /p, with singularities at the origin and p = −∞. This function, and those derived

253 from it (e.g., a pulse or a square wave) are difficult functions to invert accurately, because f (t) is discontin-

254 uous. Figures 7 and 8 show the Talbot method does not work for t < τ in double precision, since f¯t (p) grows

255 exponentially as p → −∞. The Weeks and Schapery methods do worse than the Fourier series approach

256 (even with N = 51), but their parameters can be optimized further to improve these methods.

257 Although this step boundary condition could be implemented more accurately by shifting the results

258 from the first example by t = 0.08, other step-derived time behaviors involving a pulse or square wave

259 cannot be simplified in this way.

14
 0.5 1

0 0

-0.5 -1
potential

x-flux
-1 -2
Fourier series, M=51 Fourier series, M=51
Talbot, M=51 Talbot, M=51
-1.5 Schapery, M=51 + SS -3 Schapery, M=51 + SS
Weeks, M=51 Weeks, M=51
FD solution
-2 -4
0.01 0.1 1 0.01 0.1 1
time time

Fig. 8 Plots of potential and flux through time with four methods for f¯t (p) = exp(−0.08p)/p, same p used across all t. 51

total f¯(p) evaluations are used by each method. Weeks’ solution is undefined for t < 0.08.

Table 2 Numerical summary

Method Number of Terms Free Parameters Implementation

Stehfest N ≤ decimal precision none easiest

Schapery depends on choice of pj pj via trial & error moderate

Weeks p → i∞ slowly as N grows κ & b (very sensitive to b) moderate

2M
fixed Talbot p → −∞ quickly as N grows r= 5tmax
(automatic) easy

Fourier series p → i∞ slowly as N grows T = 2tmax (automatic) most difficult

260 4.5 Numerical Results Summary

261 Table 2 summarizes results from numerical testing with these four simple boundary condition time behav-

262 iors. The second column indicates what limit there is on the number of terms in the approximation and

263 therefore the accuracy of the method. The size of p required by the Weeks and Fourier series methods grow

264 much slower than those required by the fixed Talbot method. The third column indicates what parame-

265 ters are needed to be tuned by the implementer to increase convergence of the method, and whether a

266 good choice is critical to the success of the method – an automatic method should not require searching

267 or optimizing parameters to obtain a robust solution. We define robustness as the ability of a solution to

268 remain useful, even when not at optimality. We prioritize a solution that is good enough and stable over

269 one that is excellent but catastrophically sensitive to parameter choice. The fourth column indicates the

270 ease of implementation in modern Fortran, Matlab, or NumPy. The methods could also be implemented in

271 a variable-precision environment like Mathematica or mpmath [Johansson(2011)], but this would further

272 require the BEM model be implemented in such an environment.

15
273 The modest success of the Schapery method is a bit surprising, given its simplicity and use of real p.

274 The results of the previous section were the product of many iterations of trial and error, this effort was

275 not included in the implementation effort. A better rule or parametrization of pj might make this method

276 more widely useful.

277 The sensitivity of Weeks’ method to the parameter choices was also surprising; similarly, the method

278 could have been improved after some optimization [Weideman(1999)], but Weeks’ rule of thumb was used for

279 the parameters. One of the noted advantages of Weeks’ method is the need to only compute optimal p once,

280 then any time can accurately be inverted [Kano et al(2005), Weideman(1999), Duffy(1993)]. When using the

281 simple rules-of-thumb for the the free parameters, this was not found to be the case. The generalized form

282 of Weeks’ method can include information about behavior of f¯(p) → ∞ (related to behavior as t → 0), but

283 this requires problem-specific knowledge.

284 The Fourier series method is more robust with respect to non-optimal p values, even though [Duffy(1993)]

285 cites this as a reason to use Weeks’ method over the Fourier series approach.

286 5 Conclusions

287 Laplace-space numerical approaches to solve the diffusion equation have several viable alternative inverse

288 Laplace transform algorithms to choose from. Historically, most Laplace-space solutions to the diffusion

289 equation have used real-only methods (i.e., Gaver-Stehfest or Schapery). More robust methods require

290 complex arithmetic and f¯(p) evaluations, but have the benefits of:

291 1. handling a broader class of time behaviors (Fourier series method);

292 2. still being relatively simple to implement (fixed Talbot method);

293 3. only utilizing double-precision complex data types, which are handled natively by modern Fortran,

294 Matlab, or NumPy, and by common extensions in C++ (Fourier series and Weeks’ methods).

295 Several practical recommendations are made regarding Laplace-space numerical modeling:

296 1. If many observations are needed across several time log cycles, large gains in efficiency can come from

297 inverting groups of times with a single f¯(p) vector (e.g., grouped by log cycle). This complicates the

298 implementation, but leads to much faster simulations.

16
299 2. If calculating f¯(p) is very expensive, and some numerical dispersion is allowable (not solving a wave

300 problem with sharp fronts), then the Fourier series method approach is most economical, and is auto-

301 matic and robust regarding free-parameter selection.

302 3. If only a single f¯t (p) is needed, then it may be worthwhile to optimize free parameters needed by Weeks’

303 or Piessen’s methods, or incorporate information about asymptotic f¯(p) behavior. Selection of optimum

304 b is far from automatic, and the Weeks method is not robust for non-optimal free parameters values.

305 4. If implementation time is a large factor, the fixed Talbot is quite simple to code and was automatic (no

306 need to select optimum parameters). The fixed Talbot may not work for non-zero step-time behavior

307 without extended precision.

308 5. If complex-valued function evaluations are not feasible (e.g., only real matrix or special function libraries

309 are available), the Schapery or Piessen’s methods are capable of using the same p values to invert different

310 times, which the Gaver-Stehfest method cannot.

311 6. When appropriate, the strategy used by Schapery to expand the deviation from a reference state could

312 be incorporated as a strategy to improve other algorithms.

Acknowledgements The author thanks Professors Cho Lik Chan and Barry Ganapol from the Aerospace and Mechanical

Engineering Department at the University of Arizona for initial inspiration and direction on this manuscript.

Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly

owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Admin-

istration under contract DE-AC04-94AL85000.

References

Abate and Valkó(2004). Abate J, Valkó P (2004) Multi-precision Laplace transform inversion. International Journal for

Numerical Methods in Engineering 60:979–993, DOI 10.1002/nme.995

Al-Shuaibi(2001). Al-Shuaibi A (2001) Inversion of the Laplace transform via Post-Widder formula. Integral Transforms

and Special Functions 11(3):225–232, DOI 10.1080/10652460108819314

Antia(2002). Antia H (2002) Numerical Methods for Scientists and Engineers, 2nd edn. Birkhäuser-Verlag

Bailey et al(2002). Bailey DH, Hida Y, Li XS, Thompson B (2002) ARPREC: An arbitrary precision computation package.

Tech. Rep. LBNL-53651, Lawrence Berkley National Lab

Bakker and Kuhlman(2011). Bakker M, Kuhlman K (2011) Computational issues and applications of line-elements to

model subsurface flow governed by the modified Helmholtz equation. Advances in Water Resources DOI 10.1016/j.

advwatres.2011.02.008

17
Barnhart and Illangasekare(2012). Barnhart KS, Illangarsekare TH (2012) Automatic transport model data assimilation

in Laplace space. Water Resources Research 48:W01510, DOI 10.1029/2011WR010955

Bellman et al(1966). Bellman R, Kalaba RE, Lockett JA (1966) Numerical inversion of the Laplace transform: Applications

to Biology, Economics, Engineering, and Physics. Elsevier

Brebbia et al(1984). Brebbia C, Telles J, Wrobel L (1984) Boundary Element Techniques: Theory and Practice in Engi-

neering. Springer-Verlag

Cohen(2007). Cohen AM (2007) Numerical Methods for Laplace Transform Inversion. Springer

Davies and Crann(2002). Davies A, Crann D (2002) Parallel Laplace transform methods for boundary element solutions

to diffusion-type problems. Electronic Journal of Boundary Elements BETEQ 2001(2):231–238

Davies(2005). Davies B (2005) Integral Transforms and their Applications, 3rd edn. Springer

Davies and Martin(1979). Davies B, Martin B (1979) Numerical inversion of the Laplace transform: a survey and compar-

ison of methods. Journal of Computational Physics 33:1–32, DOI 10.1016/0021-9991(79)90025-1

de Hoog et al(1982). de Hoog F, Knight J, Stokes A (1982) An improved method for numerical inversion of Laplace

transforms. SIAM Journal of Scientific and Statistical Computing 3:357–366, DOI 10.1137/0903022

Duffy(1993). Duffy DG (1993) On the numerical inversion of Laplace transforms: Comparison of three new methods on

characteristic problems from applications. ACM Transactions on Mathematical Software 19(3):333–359, DOI 10.1145/

155743.155788

Duffy(2004). Duffy DG (2004) Transform methods for solving partial differential equations. CRC Press

Hantush(1960). Hantush M (1960) Modification of the theory of leaky aquifers. Journal of Geophysical Research

65(11):3713–3725, DOI 10.1029/JZ065i011p03713

Johansson(2011). Johansson F (2011) mpmath: a Python library for arbitrary-precision floating-point arithmetic (version

0.17). http://code.google.com/p/mpmath/

Kano et al(2005). Kano PO, Brio M, Moloney JV (2005) Application of the Weeks method for the numerical inversion of

the Laplace transform to the matrix exponential. Communications in Mathematical Sciences 3(3):335–372

Kuhlman and Neuman(2009). Kuhlman KL, Neuman SP (2009) Laplace-transform analytic-element method for transient

porous-media flow. Journal of Engineering Mathematics 64(2):113–130, DOI 10.1007/s10665-008-9251-1

Kythe(1995). Kythe PK (1995) An Introduction to Boundary Element Methods. CRC Press

Lanczos(1988). Lanczos C (1988) Applied Analysis. Dover

Liggett and Liu(1982). Liggett JA, Liu PL (1982) The Boundary Integral Equation Method for Porous Media Flow. Unwin

Hyman

Lyness and Giunta(1986). Lyness J, Giunta G (1986) A modification of the Weeks method for numerical inversion of the

Laplace transform. Mathematics of Computation 47(175):313–322, DOI 10.2307/2008097

Malama et al(2009). Malama B, Kuhlman K, Revil A (2009) Theory of transient streaming potentials associated with axial-

symmetric flow in unconfined aquifers. Geophysical Journal International 179(2):990–1003, DOI 10.1111/j.1365-246X.

2009.04336.x

18
Mishra and Neuman(2010). Mishra PK, Neuman SP (2010) Improved forward and inverse analyses of saturated-

unsaturated flow toward a well in a compressible unconfined aquifer. Water Resources Research 46(7):W07508, DOI

10.1029/2009WR008899

Morales-Casique and Neuman(2009). Morales-Casique E, Neuman SP (2009) Laplace-transform finite element solution of

nonlocal and localized stochastic moment equations of transport. Communications in Computational Physics 6(1):131–

161

Oliphant(2007). Oliphant TE (2007) Python for scientific computing. Computing in Science and Engineering 9(3):10–20,

DOI 10.1109/MCSE.2007.58

Piessens(1972). Piessens R (1972) A new numerical method for the inversion of the Laplace transform. Journal of the

Institute of Mathematics and its Applications 10:185–192, DOI 10.1093/imamat/10.2.185

Schapery(1962). Schapery R (1962) Approximate methods of transform inversion for visco-elastic stress analysis. In: Pro-

ceedings of the Fourth US National Congress on Applied Mechanics, vol 2, pp 1075–1085

Stehfest(1970). Stehfest H (1970) Algorithm 368: numerical inversion of Laplace transforms. Communications of the ACM

13(1):47–49, DOI 10.1145/361953.361969

Sternberg(1969). Sternberg Y (1969) Flow to wells in the presence of radial discontinuities. Ground Water 7(6):17–20,

DOI 10.1111/j.1745-6584.1969.tb01666.x

Sudicky and McLaren(1992). Sudicky E, McLaren R (1992) The Laplace transform Galerkin technique for large-scale

simulation of mass transport in discretely fractured porous formations. Water Resources Research 28(2):499–514, DOI

10.1029/91WR02560

Talbot(1979). Talbot A (1979) The accurate numerical inversion of Laplace transforms. IMA Journal of Applied Mathe-

matics 23(1):97, DOI 10.1093/imamat/23.1.97

Weeks(1966). Weeks W (1966) Numerical inversion of Laplace transforms using Laguerre functions. Journal of the ACM

13(3):419–429, DOI 10.1145/321341.321351

Weideman(1999). Weideman J (1999) Algorithms for parameter selection in the Weeks method for inverting the Laplace

transform. SIAM Journal of Scientific Computing 21(1):111–128, DOI 10.1137/S1064827596312432

Widder(1941). Widder D (1941) The Laplace Transform. Princeton University Press

19