Table Of Content

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é, Padé-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 Queue M/MN/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é 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 the variate t, where 0 ≤ t < ∞, then the expected value of est is the Moment Generating Function which is defined by (cid:90) ∞ M(s)= estφ(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 aberrationwhichpavedthewayfortheconstructionofmoderncameras. Petzval [167] wrote a two volume treatise on the Laplace transform and its application to ordinary linear differential equations. Because of this substantial contribution 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è 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- netictheory. Heaviside’stransformwasamultipleoftheLaplacetransformand, given a transform, he devised various rules for finding the original function but withoutmuchconcernforrigour. Ifweconsiderthesimpledifferentialequation d2y +y =1, t>0 dt2 withinitialconditionsy(0)=y(cid:48)(0)=0thenHeavisidewouldwritepy fordy/dt, p2y for d2y/dt2 and so on. Thus the given equation is equivalent to (p2+1)y =1, and the ‘operational solution’ is 1 y ≡ . p2+1 Expanding the right hand side in powers of 1/p we obtain 1 1 1 y ≡ − + −··· . p2 p4 p6 PREFACE ix (cid:82) (cid:82) Heaviside regarded 1/p as equivalent to t1dt, i.e. t, 1/p2 as ttdt = t2/2!, 0 0 etc., so that the solution of the given differential equation is t2 t4 t6 y = − + −··· , 2! 4! 6! which is readily identified with 1−cost, the correct solution. For a differential equation of the form (again using the notation py = dy/dt, etc.) (a pn+a pn−1+···+a p+a )y =1, 0 1 n−1 n satisfying (cid:175) dry(cid:175) (cid:175) =0, r =0,1,··· ,n−1 dtr(cid:175) t=0 Heaviside has the operational solution 1 y = , φ(p) where we denote the nth degree polynomial by φ(p). If all the roots p , r = r 1,··· ,n of the nth degree algebraic equation φ(p) = 0 are distinct Heaviside gave the formula (known as the ‘Expansion Theorem’) 1 (cid:88)n eprt y = + . (2) φ(0) p φ(cid:48)(p ) r r 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 transform to solve integral equations in the early part of the 20th century. Based on notes left by Bateman, Erdélyi [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 foundationforoperationalmethods. Fortheexamplegivenaboveherecognized that the solution of the second order equation could be expressed as (cid:90) 1 γ+i∞ dp y = ept , 2πi p(p2+1) γ−i∞ where γ >0. For the more general nth order equation we have (cid:90) 1 γ+i∞ dp y = ept , (3) 2πi pφ(p) γ−i∞ x PREFACE where all the roots of φ(p) = 0 lie to the left of (cid:60)p = γ. The integral can be replaced by (cid:90) 1 dp y = ept , (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 (cid:90) 1 ∞ =p e−pty(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 (cid:90) ∞ f¯(s)= e−stf(t)dt. (6) 0 AnotherresearcherwhomadesignificantcontributionstothetheoryofLaplace transformswasWidderandinhisbook[253]hegivesanexpositionofthetheory of the Laplace-Stieltjes transform (cid:90) ∞ f¯(s)= e−stdα(t), (7) 0 where the function α(t) is of bounded variation in the interval 0 ≤ t ≤ R and the improper integral is defined by (cid:90) (cid:90) ∞ R e−stdα(t)= lim e−stdα(t). 0 R→∞ 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. ManyothermathematiciansmadecontributionstothetheoryofLaplacetrans- 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.