Inversion of Noise-free Laplace Transforms: Towards a Standardized Set ...

lack> Below is a cache of http://structure.bu.edu/Resource/pdfs/valko.pdf. It's a snapshot of the page taken as our search engine crawled the Web.
The web site itself may have changed. You can check the current page or check for previous versions at the Internet Archive. Yahoo! is not affiliated with the authors of this page or responsible for its content.
Inversion of Noise-free Laplace Transforms: Towards a Standardized Set of Test Problems
1
Inverse Problems in Engineering
(Publisher: Taylor & Francis)
Volume 10, Number 5, Year 2002, pp. 467 - 483

http://taylorandfrancis.metapress.com/openurl.asp?genre=issue&issn=1068-2767&volume=10&issue=5




Inversion of Noise-free Laplace Transforms: Towards a
Standardized Set of Test Problems


Peter P. Valkó
1
and Sandor Vajda
2


1
Harold Vance Department of Petroleum Engineering, Texas A&M University
mail: 3116 TAMU, College Station, TX 77843-3116
phone: (USA)-(979)-862 2757 fax: (USA)-(979)-862 1272 email: p-valko@tamu.edu

2
Department of Biomedical Engineering, Boston University
mail: 44 Cummington Street, Boston, MA 02215
phone: (USA)-(617)- 353-4757 fax: (USA)-(617)-353-6766 email: vajda@bu.edu

Abstract
The numerical inversion of Laplace transform arises in many applications of science and
engineering whenever ordinary and partial differential equations or integral equations are solved.
The increasing number of available numerical methods and computer codes has generated a need
for well-documented sets of test problems. Using such sets, algorithm developers can evaluate
the relative merits and drawbacks of their suggested new methods, and end-users can make
judgments on the applicability of an individual method for a specific problem. Many areas in
science and engineering, lead to problems that share three important properties: i) the image
function can be evaluated for real arguments, but not necessarily for complex ones; ii) the
original is known to be infinitely differentiable for times t > 0, iii) the values of the image
function can be obtained with any prescribed accuracy. The published test sets do not properly
cover these applications, as many included problems are beyond of the specific class, while the
remaining ones fail to address some of the potential difficulties arising in practice. The goal of
this paper is to establish a common ground for problem classification, to list the requirements for
the above class of problems, and to provide a carefully selected test set by addressing the
deficiencies of the ones currently available. The findings of the paper are used to solve a problem
of practical importance in modeling underground flow. Accompanying the paper, WEB links are
provided to a list of more than 800 relevant publications (going back to 1795), to a Mathematica
program to generate and solve the suggested set of test problems, and to a user friendly Java
program to solve inversion problems for a restricted class of image functions.

2
Introduction
The basic idea of integral transform methods is that in the image domain equations are usually
simpler than the original ones, and in many cases it is relatively easy to obtain the image of the
solution. Once the image is known, it has to be inverted to obtain the solution of the original
equation. In the case of Laplace transform, the inversion can be carried out using detailed Tables
of Laplace transform pairs and employing some basic properties. A more recent alternative is the
use Computer Algebra Systems capable both to look up tables and to apply basic properties.
However, in many cases the inverse is not a named function, or cannot be represented by any
simple formula, and one needs a numerical approach. To satisfy this need, a variety of powerful
inversion algorithms have been developed and tested. Nevertheless, users often find it difficult
to locate the algorithm that is most appropriate for a specific application, and each engineering or
scientific field uses only one or a few methods, in spite of the large number of available
algorithms. For example, the Stehfest [1] algorithm is rarely mentioned in recent computing
literature, and yet it is essentially the only method used in applications related to underground
flow. The best way to orient the user among various methods is to provide a large set of test
problems, and to demonstrate how a specific algorithm works on each of them. In their seminal
paper, Davies and Martin [2] presented a standard set of test problems with reproducible results.
However, as discussed by Duffy [3], this test set reflects the inherent limitations of inversion
techniques in the seventies, and hence has become obsolete.

The two authors of the present paper are interested in two very different applications, one being
the solution of underground flow and the other the distribution of drugs in pharmacokinetic
systems. These disparate fields, as well as many other areas in science and engineering, lead to
problems that share the following important properties:

i)
the image function can be evaluated on the real axis, but not necessarily on the
complex plane;
ii)
the original is infinitely differentiable for all times t > 0,
iii)
the values of the image function can be obtained with arbitrary precision.

The major goal of this paper is to create a set of test problems that satisfy the above
requirements. This means that, on one hand, we focus on a narrower class of problems than some
of the previous publications. On the other hand, we attempt to be more exhaustive, and to
represent all the potential difficulties that can occur in this particular class of applications. We
first discuss a number of concepts that help to specify the above class of inverse problems.
3

Noise. Let F(z) denote the function to be inverted. We restrict consideration to problems in
which F(z) is well defined in some mathematical sense and can be evaluated at any required
precision, possibly using Computer Algebra Systems (CAS). This case is sometimes referred to
as the Laplace transform F(z) is explicitly given [4]. Thus, we exclude problems in which
(a) F(z) is inferred from measurements, and hence only a limited number of values are available
and noise corruption is inherent; or (b) the Laplace transform is employed in conjunction with
another numerical method (such as solution of a system of equations in boundary element
methods applied in Laplace space, or another numerical inversion of the Fourier or Hankel type),
and thus the precision of F(z) values is limited in the practical sense.

Region of convergence. DAmore et al. [4] distinguish between the narrower image function
(that is defined only in the region of convergence of the Laplace integral) and a broader, less
restrictive F(z). In actual applications, only this latter, broader F(z) is available, that is the
image is given in the form of a symbolic expression or set of calculation instructions, without an
explicit reference to the region of convergence [5]. In the algorithmic literature it is taken for
granted that the Laplace transform function (in the narrower sense) can be analytically continued,
that the continuation can be done in a unique way, and that the broader function available for
inversion is identical to this unique continuation. In reality, however, some of these conditions
are impossible to check.

Regularization versus multiple precision. In contrast to the case of a noisy F(z), inverting a
transform that can be evaluated with any required precision may not require regularization. There
are other ways to proceed, first of all using extended precision calculations. The idea of extended
precision is as old as the desire to calculate inverses. Gaver [7] used 22 digits in his calculations
to generate probability functions (for a recent account, see [8]). Improving the method, Stehfest
[1] was aware of the limits that the word-length of the computer architecture (of his times) posed
on the ultimate results. Usabel [9] found that "Using GaverStehfest inversion technique, the
fact of performing calculations with more significant digits is clearly rewarded generously with a
very fast rate of convergence " The reason of this is explained by Kitaev and Frolov [10]:
Stehfest found the optimal coefficients in the sense that it eliminates N-1 members of
the Post-Widder formula. The most systematic application of multi-precision arithmetics in
Laplace inversion can be found in several works of Abate and Whitt, starting with their seminal
paper [11].

Poles and branch points. The notion that knowing something about the poles and zeros of F(z)
is useful was already suggested by Davies and Martins [2]. In the introduction of his paper Duffy
[3] states that most transforms contain both poles and branch pointsthe poles give the normal
mode solutions, while the branch points yield transient solutions. The present authors think that
some caution is beneficial when making judgments from the "properties" of the image function.
In realit