Numerical Inversion of the Laplace Transform
“Shall I refuse my dinner because I do not fully understand…digestion?”The numerical inversion of Laplace transform arises in many applications of science and engineering whenever ordinary and partial differential equations or integral equations are to be solved. The basic idea of any integral transform method is that in the image domain equations are usually simpler than the original ones. Therefore, it is relatively easy to obtain the image of the solution. Once the image is known, one has to invert it to go back to the original domain, that is to obtain the solution of the original equation.
In the case of Laplace transform, for this purpose one can use detailed Tables of Laplace transform pairs and some tricks known as the basic properties of Laplace transform. A more recent alternative is to use Computer Algebra Systems (CAS) – to do the table lookup (and the tricks.)
Unfortunately, it is not always easy to find the inverse (that is the original of the image) using the tables or CAS. One possible reason is that the inverse is not a named function or can not be represented by any simple “formula”. In general, one would need a numerical algorithm for the inversion problem.
Many such algorithms have been suggested and this WEB page is devoted to them
Numerical Inversion of the Laplace Transform
The List
Go to The List [PDF] of more than 1500 publications maintained together with Vladimír Vojta
Go to the page: Sets of Test Problems
where you can find supporting material to the paper (published in the journal Inverse Problems in Engineering) written together with S. Vajda of Boston.
GWR paper
Valkó, P.P. and Abate, J. : Comparison of Sequence Accelerators for the Gaver Method of Numerical Laplace Transform Inversion, Computers and Mathematics with Application, Vol. 48 (Iss.3-40) 2004 pp. 629-636. [PDF]
Mathematica source code of the GWR algorithm: http://library.wolfram.com/infocenter/MathSource/4738/
FT paper
Abate, J. and Valkó, P.P.: Multi-precision Laplace transform inversion, International Journal for Numerical Methods in Engineering, Vol. 60 (Iss. 5-7) 2004 pp 979–993. [PDF]
Mathematica source code of the FT algorithm:
http://library.wolfram.com/infocenter/MathSource/5026/
What is the Gaver-Stehfest algorithm?
It is a relatively simple numerical method to invert the Lapalace transform. If you want to read the original paper, you will find its bibliographic info in The List. (You may easily search for Stehfest.) If you want to find more the publications dealing with the algorithm, you may want to check out letter “S” in the third column.
What is the Gaver-Wynn-Rho (GWR) algorithm?
The widely-used Stehfest algorithm applies a certain linear convergence acceleration scheme to a series (originally introduced by Gaver) that approximates the inverse.
In contrast, the GWR algorithm uses another (nonlinear) acceleration scheme to the same series. The acceleration scheme itself was first suggested by Wynn. Its use in numerical Laplace transform inversion is suggested in the GWR paper above.
What is the Fixed-Talbot (FT) algorithm?
The Talbot algorithm approximates the complex path integral called “Bromwich integral”, using a special path. The Fixed Talbot version uses a slightly modified path as described in the FT paper above.
Some links of interest:
Luisa D’Amore W. Whitt Abate and Whitt Urs E. GrafWeideman J.H. Matthews H. Weber A. Mallet (MathSource) K. Bryan Kano and Brio
Common in the GWR and FT Mathematica codes:
The Mathematica code controls internally the precision of the function evaluations and this way it avoids the usual pitfalls common to other published codes.
To put it simply: if the inverse exists you can calculate it accurately.
You can figure out how many significant digits are in the result or you can set a target for the required number of significant digits and reach it by increasing the number of terms: M.
In contrast, with codes written in traditional programming languages and using machine precision arithmetics the number of terms should be kept very low in any method, otherwise round off errors start to dominate.