Page 469 - A First Course In Stochastic Models

464 APPENDICES

that the discretization error can be limited to 10 −8 by choosing a = 19.1 and to
10 −13 by choosing a = 28.3. However, the constant a should not be chosen unnec-
essarily large. The risk of losing significant digits when calculating the infinite series
in (F.1) increases when the constant a gets too large. A useful method of summa-
tion for the infinite series in (F.1) is the classical Euler summation method. This
method proves to be quite effective in accelerating the convergence of the alternat-
ing infinite series in (F.1). Also, the method decreases the risk of losing significant
k
digits in the calculations. In Euler summation the infinite series
∞ (−1) a k in
k=0
(F.1) is approximated by the Euler sum
m
m −m
E(m, n) = ( )2 S n+k
k
k=0
for appropriately chosen values of m and n, where
j
k
S j = (−1) a k .
k=0
Numerical experience shows that the Euler sum E(m, n) computes the infinite

∞ k −13
series k=0 (−1) a k in (F.1) usually with an error of 10 or less when n = 38
and m = 11 are taken (this requires the computation of only 50 terms). For more
details the interested reader is referred to Abate and Whitt (1992). The Abate–Whitt
algorithm gives excellent results for functions f (t) that are sufficiently smooth
(say, twice continuously differentiable). However, the inversion algorithm performs
less satisfactorily for points at which the function f (t) or its derivative is not
differentiable.

Inversion algorithm of Den Iseger
Another simple algorithm to invert Laplace transforms was given in Den Iseger
(2002). In general this algorithm outperforms the Abate–Whitt algorithm in sta-
bility and accuracy. The strength of the Den Iseger algorithm is the fact that in
essence it boils down to an application of the discrete FFT algorithm. The Den
Iseger algorithm has the additional advantage of inverting the Laplace transform
simultaneously at several points. Suppose you wish to calculate f (t) for a number
of points in the interval [0, t 0 ]. For appropriately chosen values of > 0 and
M > 1 with (M − 1) = t 0 , the algorithm calculates the function values f (ℓ )
for ℓ = 0, 1, . . . , M − 1. The algorithm is based on the representation

n 1
e bℓ b + iλ j + iπt
f (ℓ ) ≈ α j Re f ∗ cos(πℓ(t + 1)) dt (F.2)
−1
j=1
for appropriately chosen values of b and n, where the abscissae λ j and the weights
α j for j = 1, . . . , n are given numbers that depend only on n. The error in (F.2)
converges very fast to zero as n gets larger. For practical purposes it suffices to