NUMERICAL METHODS FOR LAPLACE TRANSFORM INVERSION

261 Pages·2007·2.38 MB·English

NUMERICAL METHODS FOR
LAPLACE TRANSFORM INVERSION
Numerical Methods and Algorithms

VOLUME 5

Series Editor:

Claude Brezinski
Université des Sciences et Technologies de Lille, France
NUMERICAL METHODS FOR
LAPLACE TRANSFORM INVERSION

By

ALAN M. COHEN
Cardiff University
Library of Congress Control Number: 2006940349
ISBN-13: 978-0-387-28261-9 e-ISBN-13: 978-0-387-68855-8

Printed on acid-free paper.

AMS Subject Classifications: 44A10, 44-04, 65D30, 65D32, 65Bxx

 2007 Springer Science+Business Media, LLC
All rights reserved. This work may not be translated or copied in whole or in part without the written
permission of the publisher (Springer Science+Business Media, LLC, 233 Spring Street, New York, NY
10013, USA), except for brief excerpts in connection with reviews or scholarly analysis. Use in
connection with any form of information storage and retrieval, electronic adaptation, computer software,
or by similar or dissimilar methodology now known or hereafter developed is forbidden.
The use in this publication of trade names, trademarks, service marks, and similar terms, even if they are
not identified as such, is not to be taken as an expression of opinion as to whether or not they are subject
to proprietary rights.

Printed in the United States of America.

9 8 7 6 5 4 3 2 1

springer.com
Contents
Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii
Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii
1 Basic Results 1
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Transforms of Elementary Functions . . . . . . . . . . . . . . . . 2
1.2.1 Elementary Properties of Transforms . . . . . . . . . . . . 3
1.3 Transforms of Derivatives and Integrals . . . . . . . . . . . . . . 5
1.4 Inverse Transforms . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Convolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 The Laplace Transforms of some Special Functions . . . . . . . 11
1.7 Difference Equations and Delay Differential Equations . . . . . . 14
1.7.1 z-Transforms . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.8 Multidimensional Laplace Transforms . . . . . . . . . . . . . . . 18
2 Inversion Formulae and Practical Results 23
2.1 The Uniqueness Property . . . . . . . . . . . . . . . . . . . . . . 23
2.2 The Bromwich Inversion Theorem . . . . . . . . . . . . . . . . . 26
2.3 The Post-Widder Inversion Formula . . . . . . . . . . . . . . . . 37
2.4 Initial and Final Value Theorems . . . . . . . . . . . . . . . . . . 39
2.5 Series and Asymptotic Expansions . . . . . . . . . . . . . . . . . 42
2.6 Parseval’s Formulae . . . . . . . . . . . . . . . . . . . . . . . . . 43
3 The Method of Series Expansion 45
3.1 Expansion as a Power Series . . . . . . . . . . . . . . . . . . . . . 45
3.1.1 An alternative treatment of series expansions . . . . . . . 49
3.2 Expansion in terms of Orthogonal Polynomials . . . . . . . . . . 49
3.2.1 Legendre Polynomials . . . . . . . . . . . . . . . . . . . . 50
3.2.2 Chebyshev Polynomials . . . . . . . . . . . . . . . . . . . 52
3.2.3 Laguerre Polynomials . . . . . . . . . . . . . . . . . . . . 55
3.2.4 The method of Weeks . . . . . . . . . . . . . . . . . . . . 58
3.3 Multi-dimensional Laplace transform inversion . . . . . . . . . . 66
vi CONTENTS
4 Quadrature Methods 71
4.1 Interpolation and Gaussian type Formulae . . . . . . . . . . . . 71
4.2 Evaluation of Trigonometric Integrals . . . . . . . . . . . . . . . 75
4.3 Extrapolation Methods . . . . . . . . . . . . . . . . . . . . . . . . 77
4.3.1 The P -transformation of Levin . . . . . . . . . . . . . . . 77
4.3.2 The Sidi mW-Transformation for the Bromwich integral . 78
4.4 Methods using the Fast Fourier Transform (FFT) . . . . . . . . . 81
4.5 Hartley Transforms . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.6 Dahlquist’s “Multigrid” extension of FFT . . . . . . . . . . . . . 95
4.7 Inversion of two-dimensional transforms . . . . . . . . . . . . . . 100
5 Rational Approximation Methods 103
5.1 The Laplace Transform is Rational . . . . . . . . . . . . . . . . . 103
5.2 The least squares approach to rational Approximation . . . . . . 106
5.2.1 Sidi’s Window Function . . . . . . . . . . . . . . . . . . . 108
5.2.2 The Cohen-Levin Window Function . . . . . . . . . . . . 109
5.3 Pad´e, Pad´e-type and Continued Fraction Approximations . . . . 111
5.3.1 Prony’s method and z-transforms . . . . . . . . . . . . . . 116
5.3.2 The Method of Grundy . . . . . . . . . . . . . . . . . . . 118
5.4 Multidimensional Laplace Transforms . . . . . . . . . . . . . . . 119
6 The Method of Talbot 121
6.1 Early Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.2 A more general formulation . . . . . . . . . . . . . . . . . . . . . 123
6.3 Choice of Parameters . . . . . . . . . . . . . . . . . . . . . . . . . 125
6.4 Additional Practicalities . . . . . . . . . . . . . . . . . . . . . . . 129
6.5 Subsequent development of Talbot’s method . . . . . . . . . . . . 130
6.5.1 Piessens’ method . . . . . . . . . . . . . . . . . . . . . . . 130
6.5.2 The Modification of Murli and Rizzardi . . . . . . . . . . 132
6.5.3 Modifications of Evans et al . . . . . . . . . . . . . . . . . 133
6.5.4 The Parallel Talbot Algorithm . . . . . . . . . . . . . . . 137
6.6 Multi-precision Computation . . . . . . . . . . . . . . . . . . . . 138
7 Methods based on the Post-Widder Inversion Formula 141
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
7.2 Methods akin to Post-Widder . . . . . . . . . . . . . . . . . . . . 143
7.3 Inversion of Two-dimensional Transforms . . . . . . . . . . . . . 146
8 The Method of Regularization 147
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
8.2 Fredholm equations of the first kind — theoretical
considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
8.3 The method of Regularization . . . . . . . . . . . . . . . . . . . . 150
8.4 Application to Laplace Transforms . . . . . . . . . . . . . . . . . 151
CONTENTS vii
9 Survey Results 157
9.1 Cost’s Survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
9.2 The Survey by Davies and Martin . . . . . . . . . . . . . . . . . 158
9.3 Later Surveys . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
9.3.1 Narayanan and Beskos . . . . . . . . . . . . . . . . . . . . 160
9.3.2 Duffy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
9.3.3 D’Amore, Laccetti and Murli . . . . . . . . . . . . . . . . 161
9.3.4 Cohen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
9.4 Test Transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
10 Applications 169
10.1 Application 1. Transient solution for the Batch Service
N
Queue M/M /1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
10.2 Application 2. Heat Conduction in a Rod. . . . . . . . . . . . . . 178
10.3 Application 3. Laser Anemometry . . . . . . . . . . . . . . . . . 181
10.4 Application 4. Miscellaneous Quadratures. . . . . . . . . . . . . . 188
10.5 Application 5. Asian Options . . . . . . . . . . . . . . . . . . . . 192
11 Appendix 197
11.1 Table of Laplace Transforms . . . . . . . . . . . . . . . . . . . . . 198
11.1.1 Table of z-Transforms . . . . . . . . . . . . . . . . . . . . 203
11.2 The Fast Fourier Transform (FFT) . . . . . . . . . . . . . . . . . 204
11.2.1 Fast Hartley Transforms (FHT) . . . . . . . . . . . . . . . 206
11.3 Quadrature Rules . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
11.4 Extrapolation Techniques . . . . . . . . . . . . . . . . . . . . . . 212
11.5 Pad´e Approximation . . . . . . . . . . . . . . . . . . . . . . . . . 220
11.5.1 Continued Fractions. Thiele’s method . . . . . . . . . . . 223
11.6 The method of Steepest Descent . . . . . . . . . . . . . . . . . . 226
11.7 Gerschgorin’s theorems and the Companion Matrix . . . . . . . . 227
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 231
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 249
Preface
The Laplace transform, as its name implies, can be traced back to the work of
the Marquis Pierre-Simon de Laplace (1749-1827). Strange as it may seem no
reference is made to Laplace transforms in Rouse Ball’s “A Short Account of
the History of Mathematics”. Rouse Ball does refer to Laplace’s contribution
to Probability Theory and his use of the generating function. Nowadays it is
well-known that if φ(t) is the probability density in the distribution function of
st
the variate t, where 0 ≤ t < ∞, then the expected value of e is the Moment
Generating Function which is defined by
∫ ∞
st
M(s) = e φ(t)dt. (1)
0
The term on the right hand side of (1) is, if we replace s by −s, the quantity
that we now call the Laplace transform of the function φ(t).
One of the earliest workers in the field of Laplace transforms was J.M. Petzval
(1807-1891) although he is best remembered for his work on optical lenses and
aberration which paved the way for the construction of modern cameras. Petzval
[167] wrote a two volume treatise on the Laplace transform and its application
to ordinary linear differential equations. Because of this substantial contribu-
tion the Laplace transform might well have been called the Petzval transform
had not one of his students fallen out with him and accused him of plagiarising
Laplace’s work. Although the allegations were untrue it influenced Boole and
Poincar`e to call the transformation the Laplace transform.
The full potential of the Laplace transform was not realised until Oliver Heavi-
side (1850-1925) used his operational calculus to solve problems in electromag-
netic theory. Heaviside’s transform was a multiple of the Laplace transform and,
given a transform, he devised various rules for finding the original function but
without much concern for rigour. If we consider the simple differential equation
2
d y
+ y = 1, t > 0
2
dt
′
with initial conditions y(0) = y (0) = 0 then Heaviside would write py for dy/dt,
2 2 2
p y for d y/dt and so on. Thus the given equation is equivalent to
2
(p + 1)y = 1,
and the ‘operational solution’ is
1
y ≡ .
2
p + 1
Expanding the right hand side in powers of 1/p we obtain
1 1 1
y ≡ − + − · · · .
2 4 6
p p p
PREFACE ix
∫ ∫
t 2 t 2
Heaviside regarded 1/p as equivalent to 1dt, i.e. t, 1/p as tdt = t /2!,
0 0
etc., so that the solution of the given differential equation is
2 4 6
t t t
y = − + − · · · ,
2! 4! 6!
which is readily identified with 1 − cos t, the correct solution.
For a differential equation of the form (again using the notation py = dy/dt,
etc.)
n n−1
(a0p + a1p + · · · + an−1p + an)y = 1,
satisfying
∣
dry ∣
∣
= 0, r = 0, 1, · · · , n − 1
r ∣
dt
t=0
Heaviside has the operational solution
1
y = ,
φ(p)
where we denote the nth degree polynomial by φ(p). If all the roots pr, r =
1, · · · , n of the nth degree algebraic equation φ(p) = 0 are distinct Heaviside
gave the formula (known as the ‘Expansion Theorem’)
n
∑ prt
1 e
y = + . (2)
′
φ(0) prφ (pr)
r=0
Compare this to (1.23). Carslaw and Jaeger [31] give examples of Heaviside’s
approach to solving partial differential equations where his approach is very
haphazard. Curiously, in his obituaries, there is no reference to his pioneering
work in the Operational Calculus.
Bateman (1882-1944) seems to have been the first to apply the Laplace trans-
form to solve integral equations in the early part of the 20th century. Based
on notes left by Bateman, Erd´elyi [78] compiled a table of integral transforms
which contains many Laplace transforms.
Bromwich (1875-1929), by resorting to the theory of functions of a complex
variable helped to justify Heaviside’s methods to some extent and lay a firm
foundation for operational methods. For the example given above he recognized
that the solution of the second order equation could be expressed as
∫
γ+i∞
1 dp
pt
y = e ,
2
2πi p(p + 1)
γ−i∞
where γ > 0. For the more general nth order equation we have
∫
γ+i∞
1 dp
pt
y = e , (3)
2πi pφ(p)
γ−i∞
x PREFACE
where all the roots of φ(p) = 0 lie to the left of ℜp = γ. The integral can be
replaced by
∫
1 dp
pt
y = e , (4)
2πi pφ(p)
C
where C is any circular path with centre the origin which contains all the zeros
of φ(p) = 0 inside its circumference. Applying the Cauchy residue theorem to
this integral leads to the expansion theorem (2). Again no specific mention of
his contribution to the Operational Calculus was made in Bromwich’s obituary
(Hardy [111]). However, Jeffreys [118] gives an exposition of his methods in his
book.
Starting in the 1920’s considerable effort was put into research on transforms.
In particular, Carson [32] and van der Pol made significant contributions to
the study of Heaviside transforms. Thus Carson established, for the differential
equation considered above, that
∫
∞
1
−pt
= p e y(t)dt. (5)
φ(p)
0
Van der Pol gave a simpler proof of Carson’s formula and showed how it could
be extended to deal with non-zero initial conditions. Doetsch in his substantial
contributions to transform theory preferred to use the definition that is now
familiar as the Laplace transform
∫
∞
¯ −st
f(s) = e f(t)dt. (6)
0
Another researcher who made significant contributions to the theory of Laplace
transforms was Widder and in his book [253] he gives an exposition of the theory
of the Laplace-Stieltjes transform
∫
∞
¯ −st
f(s) = e dα(t), (7)
0
where the function α(t) is of bounded variation in the interval 0 ≤ t ≤ R and
the improper integral is defined by
∫ ∫
∞ R
−st −st
e dα(t) = lim e dα(t).
R→∞
0 0
This, of course, reduces to the standard definition of the Laplace transform
when dα(t) = f(t)dt and f(t) is a continuous function. An advantage of the
Laplace-Stieltjes approach is that it enables us to deal in a sensible manner
with functions of t which have discontinuities. In particular, we note that the
inversion theorem takes the form (2.14) at a point of discontinuity.
Many other mathematicians made contributions to the theory of Laplace trans-
forms and we shall just recall Tricomi [233] who used expansions in terms of
Laguerre polynomials in order to facilitate the inversion of Laplace transforms.