1. No category

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 Vajda2
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: [email protected]
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: [email protected]
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.
1
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)
ii)
iii)
the image function can be evaluated on the real axis, but not necessarily on the
complex plane;
the original is infinitely differentiable for all times t > 0,
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.
2
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. D’Amore 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 Gaver–Stehfest … 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 points—the 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 reality, the above language is somewhat misleading because it is a generalization of the
concepts used originally for a special class of transforms (rational polynomial functions) that
arise in linear circuit theory.
Practical problems. Duffy [3] argues that “The standard method used to test a numerical
inversion scheme is to employ a transform whose inverse is known exactly. Although this is
3
acceptable when developing the scheme, such guidance applied to transforms encountered is in
general not useful….” He introduces the concept of tests “representative of those found in
applied problems”. Since this is also the governing concept of our work, we will elaborate on it
in more detail.
Why do we need new sets of test problems?
Davies and Martin [2] were the first to recognize the importance of a standardized set of test
problems. Their selection has been somewhat problematic for several reasons.
(1) They selected image functions that could be evaluated without additional programming in
standard FORTRAN IV.
(2) Relying on the traditions of Laplace transform tables, they included image functions
whose original may not be infinitely many times differentiable for all positive times.
(3) Usabel [8] notes that only relatively short times were considered.
(4) Duffy [3] stated that all the inverses could be "found analytically", and hence were not
representative of real life practical problems.
Of course, problem (1) is related to the inherent limitations of the computational environment of
their time and should not be a constraint now.
In our opinion, problem (2) is more fundamental. The existence of singularities is a major
problem, because numerical algorithms based on calculus concepts will always fail in the
vicinity of such points. Since these large errors are essentially independent of the method, they
may mask the real differences occurring in problems in which the original is infinitely
differentiable for all t > 0, and hence the performance of different methods can not be rigorously
assessed.
The problem described in (3) is not inherent, since the time is “short” or “long” only in the
context of a particular image function and its parameters. Therefore, the problem is avoided if
one defines a large number of test problems with randomly selected numerical parameters,
affecting the “time constants” of the functions.
As far as problem (4) is considered, we disagree with the criticism of Duffy [3]. The
classification of a test problem as practical or not is controversial. First of all, a test has to have a
solution. If the solution can be found only by another numerical method, which itself is complex
enough to make it practically impossible to estimate the error, then the usefulness of the test is
questionable. Second, problems that are found in the engineering literature and have a well
established solution, are almost always taken from published tables. Therefore, the requirement
to find “more complicated” practical problems that have well-defined solutions is difficult to
meet. In fact, the tables of Laplace transforms contain practically all the known image-original
pairs ever encountered in practice. Therefore, an inversion method clearly should be considered
4
valuable if it can reproduce a great number of these pairs, particularly if the parameters in the
images are selected somewhat randomly.
We show that the emphasis on practical problems to be solved numerically may lead to major
difficulties. A good example is Test Problem No 1 of Duffy [3]. The image function is
F ( z) =
1 1
1 
− 2z

z ( z + 1)  2 z e − 1
(1)
and the inverse is required at t = 1, 2, 4, and 8. The performance of several inversion algorithms
was measured by comparing the results to f(t) values calculated independently from the
expression
f (t ) =
[
]
1
e 2  1 N sin nπt − tan −1 (nπ )
1
+ e −t  − 2  − ∑
2
n n 2π 2 + 1
 2 e − 1 π n =1
(2)
According to Duffy, N = 107 terms in the above sum yields an error of O(10 - 7).
Duffy selected this problem because “it did not contain any branch point singularities, and the
transform and its inverse were similar to those that arise in solutions of the wave equation over a
finite domain” [3]. In addition, according to Duffy “… first, unlike most transforms in the
literature, we have an infinite number of singularities to resolve rather than just a few.
Consequently, for any finite number of terms in the inversion, we badly misrepresent the
contribution from the poles far from the origin. Second, these poles lie along the imaginary axis
rather than the negative real axis and they represent sinusoidal solutions. Although the
contributions of poles far from the origin are not dominant, they do represent higher harmonics
which are necessary for accurate results.”
While most authors considered the above problem challenging because it was difficult to obtain
the f(t) function from Eq. (2), it has also been recognized that something may be wrong.
Weideman [12] wrote “though unrelated to our algorithm, we should point out that the very
small errors observed … should be viewed with circumspection because… f(t) is nondifferentiable at t = 2, 4, 6, 8; as may be deduced from the infinite series representation”.
D’Amore et al. [13] mentioned that the “infinite sum … was calculated numerically using the
NAG Library”, but did not present explicitly the numerical results.
Our major criticism with regard to the above test problem is that, employing some basic
properties of the Laplace transform, in the considered interval t ≤ 8 , the original can be
represented as the simple function
f (t ) =
(
)
(
)
(
)
(
)
1 −t
e + t − 1 + u (t − 2) e 2−t − 1 + u (t − 4) e 4−t − 1 + u (t − 6) e 6−t − 1
2
(3)
5
where u(t) is the Heaviside function. We suspect that: should the simple form (3) be included in
Duffy’s paper, the above Laplace pair would not be considered a “tough test problem” and would
not be included in any subsequent studies.
In general, when the original contains delayed functions, and therefore is suspect to being nondifferentiable (or even jumping) at discrete points, the user knows this fact a priori. In the case
of underground flow applications [14-16], this is due to abrupt production rate changes. In
compartmental modeling of pharmacokinetics [17-20] the phenomenon is related to
administering the drug, say, by intravenous injection at a specific time. The standard way to
handle this situation is to decompose the response into the sum of delayed functions, each of
which is infinitely differentiable (for positive times) in itself, and to obtain only the individual
differentiable constituents from the table of transforms, or by numerical inversion. Therefore, we
will exclude the problem just discussed from our investigation.
We also believe that the fact that an independent direct calculation of f(t) is “difficult”, does not
mean that the test problem is challenging or representative. As an example, we consider Duffy’s
Test Problem No 2. The image function is given as
F ( z) =
(100 − z )sinh (
)
z [z sinh ( z ) + z cosh ( z )]
z /2
(4)
and the f(t) function can be independently evaluated by the series
1 N 
1
f (t ) = − + ∑ 100 + 2
2 n =1 
bn
 2bn sin (bn / 2)e − bn t

2
 2 + bn cos(bn )
2
(
)
(5)
where the coefficient bn can be obtained as the subsequent solutions to a nonlinear
(transcendent) equation: x tan(x) = 1 . The combination of finding an (infinite) number of roots
numerically and summing an (infinite) number of terms is not a trivial task, and one can hardly
guarantee the accuracy of the independently calculated f(t) values. Therefore, the above pair is
not an ideal test problem.
Suggested Sets of Test Problems
Set A
The most extensive, methodical, and reliable tables of Laplace transforms, such as Roberts and
Kaufman [21] and Prudnikov, Brychov and Marichev [22] list pairs on hundreds of pages.
However desirable is to use such a complete table, that task should be postponed for better times.
For this work, we use Table 29.3 of Abramowitz and Stegun [23]. This excellent compilation
contains a relatively small number of entries (129), but it is very economical in the sense that few
special cases are listed, if the parametrization of a more general item allows to generate the
special one. Therefore, we select this table as our starting point.
6
From the 129 entries, 105 f(t) functions are infinitely differentiable for all t > 0. The nondifferentiable f(t) functions are easy to filter out, because they are represented by a figure, and/or
they contain the Heaviside function or the Dirac delta symbol. We have programmed the 105
pairs in Mathematica (version 4.1 [24]), using a notation as near to the original compilation as
possible. Particular care was exercised to use the same symbol for a similar type of parameter,
therefore we use the notation: real numbers: a, b, c, k; postive integers: nu, mu, n, for the
generating parameters. As in the original table, all parameters should be different, and b > a. In
the runs shown later, we generated the actual test problems selecting a = 3/5, b = 5/7, c = -9/7,
k = 9/11, nu = 3, mu = 4 and n = 5 . The times at which we require the calculation of the
numerical inverse are t = 1/10, 1, 10.
In the resulting 105 problems that will be referred to as Set A, the original is either known
“analytically” or can be evaluated by a finite number of calls to named functions plus some
arithmetic operations. Mathematica is capable to evaluate any of these expressions with any
required number of significant digits. In addition, the f(t) functions are infinitely many times
differentiable at any positive t (but not necessarily at t = 0).
We consider Set A as the most objective measure of the performance of a numerical inversion
algorithm, because i) this set has been created from a table representing the wide spectrum of
Laplace transform pairs, ii) no additional subjective consideration was involved in the selection
process, and iii) the f(t) values are known with any required accuracy.
Set B
The next set is a collection of test problems with inverse that can be calculated with controlled
precision, using an “infinite” summation whose coefficients are still explicitly available. Many
problems of heat conduction, diffusion, flow in porous media, electricity fall into this category.
A typical example would be test B-1:
F ( z) =
( )
sinh ( z )
sinh x z
(6)
where 0 ≤ x ≤ 1 . The corresponding f(t) is known to have the form:
∞
f (t ) = 2π ∑ (−1) n n e − n π t sin( nπx)
2
2
(7)
n =1
at least that is given on page 217 of the current textbook by Schiff [5]. The selection of x follows
our established practice, we use x = 2/5. The time values are the same as previously, t = 0.1, 1,
and 10. Several similar problems are collected in Set B. As will be described, the results present
some surprises.
7
Set C
Set C is in the tradition of Duffy [3], and includes his Test Problems No 2, 3, 4 and 5 as C-1
through C-4. The independent solutions of these problems can be generated using an arsenal of
numerical techniques, and some human interaction is necessary to accept the “right”
independently calculated value of the inverse. As it turns out, the available other methods could
not always provide an accurate enough f(t) and we added a question mark after the “known”
values to denote this fact.
The next problem, our test C-5, is taken from De Chant [31]. The Laplace transform is related to
the one-dimensional heat equation with an important boundary condition. The problem used to
have an “analytical” solution in the form of infinite series that seemed to be reasonable. As De
Chant suggests, there has been a typo in the published solution that had remained undiscovered
for half a century. We wished to see if a numerical inversion method can be used to tell which
solution is correct: the one suggested in a classical book or the one presented by De Chant.
Two additional problems are related to the Mittag-Leffler distribution. Test C-6 has a closedform solution [32], but for test C-7 only a series representation is available which is difficult to
calculate to high precision.
Set D
In this set we collected some important problems, somewhat on the “verge” of symbolic and
numerical mathematics.
( D-1: ) F ( z ) = exp(− z a ) where 0 < a < 1,
and
( D-2: ) F ( z ) = exp(− z b ) where b < 0
are taken from Mikusinski and Boehme [33]. (We use a = 2/5 and b = -2/5). (According to
Mikusinski and Bohme, no “analytical” representation is known for the f(t) function of D-2,
though the inverse itself exists.)
Problem ( D-3: ) F ( z ) = z ln ( z ) is interesting, because strictly speaking the F(z) function is not a
Laplace transform. Indeed, we know that a Laplace transform cannot have the property
F ( z ) → ∞ as Re( z ) → ∞ (see Schiff [5], page 22). Nevertheless, the F(z) is a transform, in
the Heaviside-Mikusinski operational sense, and the inverse would be
f (t ) = t −2
It is interesting to see what happens if we try to invert Problem D-3 using a particular inversion
algorithm.
Problem ( D-4: ) F ( z ) = z n + a should result in f (t ) = Γ(− n − a ) / t n + a in the HeavisideMikusinski operational sense ( we select n = 3 and a = 2/5 ), though it is not a Laplace transform
in the strict sense.
8
should result in the generalization of the
Problem ( D-5: ) F ( z ) = 1 /(1 + z a ) − n = (1 + z a ) n
Mittag-Leffler distribution. We select ( n = 3 and a = 2/5 ). This exponent is not included in the
original definition [31], so we do not know if f(t) exists or not.
Problem ( D-6: ) F ( z ) = exp(− z c ) with c = 3/2 is of great importance. We know from [32] that
there is no inverse in the Heaviside-Mikusinkski (or Laplace) sense. We wish to see whether an
algorithm provides a false solution or indicates somehow that there is no inverse.
Problem ( D-7: ) F ( z ) = Ei(− z ) with
if t < 1
0

f (t ) = - 0.5 if t = 1
- 1/t if t > 1

is included as a counterexample. Obviously, f(t) is non-differentiable in t = 1 so this is an
example of a transform that is outside the scope of this study. The only reason we include it is to
see how a particular algorithm “reacts” to such a problem.
Using the Test Sets
Of course no set of test problems is meaningful without demonstrating that it is realistic to
anticipate a reasonable performance of at least one available inversion algorithm on it, and that
the problems included provide useful information on the performance of a particular method.
The BigNumber-Stehfest algorithm has been selected here because it provides results that agree
with the “analytical” ones whenever the latter are available.
The BigNumber-Stehfest Algorithm
The BigNumber-Stehfest algorithm is a slightly modified version of the Stehfest [1] method, also
called Gaver-Stehfest method. The improvement is achieved using the following ideas:
1) The number of terms in the Stehfest summation (N) is successively doubled
2) The precision of the internal arithmetic calculations and the F(z) evaluations is increased
simultaneously with the number of terms in the algorithm (for N terms N decimal digit
precision is required)
3) A simple "convergence criterion" is applied: the digits being the same in the
approximants calculated with the number of terms N and 2N are considered already
"final" e.g. significant.
4) A “processed” form of the Laplace transform is used. The original Stehfest algorithm
requires, that the F(z) function evaluates to a real number at any positive z. Obviously,
this is not true for all entries in our Set A. As an example consider our test no A-89 (entry
9
 z 2 − 9 / 25 
 that will not evaluate
29.3.109 in Abramowitz and Stegun [22]): F ( z ) = ln
z2


to a real value for z < 3/5. Our defined rule of processing is simple: Whenever a noninteger exponent or logarithm appears, insert an absolute value operator. For the example
 z 2 − 9 / 25 
 . Note that in all these cases
considered, this creates F ( z ) _ Processed = ln
2

z


we change only the symbolic expression representing the broader F(z), and do not change
the Laplace transform itself. In all cases, the suggested processing can be done without
contradicting
the
principle
of
analytic
continuation,
because
neither
2
2
 z − 9 / 25 
 z − 9 / 25 
 is the “right” continuation,
 or F ( z ) _ Processed = ln
F ( z ) = ln
2


z
z2




being F(z) not analytic at the real argument z = 3/5.
5) The Stehfest coefficients (being rational fractions) are stored in an integer representation,
resulting in the required precision. (The complete data file is about 40 Megabyte and it
allows to calculate an approximant with up to N = 212 = 4096 terms, in contrast to the
usual practice of using 16 terms. The calculation of the integer representation required
several days of computation on five PCs, running simultaneously.)
Understanding the Output
In the output of the test calculations, the symbolic expression for F(z) is denoted similarly, and
the known inverse is denoted by f_Known. In addition, we attempt to use the built-in procedure
of Mathematica to find symbolically the F from the known f and the f from the known F.
Whenever a meaningful result was produced, we incorporated it into the output as
F_From_fKnown and f_Mma_Orig, respectively. The “Orig” refers to the built-in procedure:
InverseLaplaceTransform[] of Mathematica 4.1. Another result for the symbolic inverse,
denoted by f_Mma_UG is obtained using the Mathematica add-on package Lz of U. Graf [24].
In the output we incorporated another form of F(z), denoted by F_Processed and defined in the
previous section. A typical part of the output is shown in Fig. 1.
10
t
f_Known
f_MmaOrig
f_MmaUG
f_niL
1.0e-1
-0.36010801296083318e-1
-0.36010801296083318e-1
NA
-0.36010801296083318e-1
1.0e0
-0.37093043648453541e0
-0.37093043648453541e0
NA
-0.37093043648453541e0
1.0e1
-0.40143127224491179e2
-0.40143127224491179e2
NA
Message: no evidence of convergence
Fig. 1 Excerpt from the output
11
As we see from the sample output, in the particular case:
•
The built-in Mathematica procedure was unable to create the symbolic Laplace transform from the known
original f(t).
•
The Processed F(z) contains the absolute value “operator” wrapped around the argument of the logarithm.
This is the actual function evaluated in the BigNumber-Stehfest algorithm.
•
The built-in procedure of Mathematica did find the inverse symbolically from the known F(z) (even if the
identity of the f_Known and f_MmaOrig expressions is not obvious at the first sight).
•
In this case U. Graf’s package did not provide a symbolic inverse of the known F(z).
•
The numerical inverse calculated by the BigNumber-Stehfest procedure (called niL) is exact to 17
significant digits for time t = 0.1 and t = 1.
•
For time t = 10, there is an error message. The procedure was unable to check the convergence (though the
result was already accurate.)
This run was made with maximum N = 28 = 256 Stehfest terms. Using N = 210, there was no
error message for this particular test problem, or any of the 105 entries included in Set A. The
complete output is available from our WEB site (see below).
Results for Set A
For all 105 test problems, the calculated numerical inverses are already accurate to 17 significant
digits when we allowed maximum N = 210 terms in the BigNumber-Stehfest algorithm. If time
t = 100 was also included (or time t = 0.01), the maximum allowed number of terms had to be
increased to 212, and the computations needed more than a day on a current PC. If time t = 1000
is also included, seven tests fail at that t-value. There is no doubt that these calculations could
also be done to 17 digit accuracy, but we currently do not have a data file of the Stehfest
coefficients for more than 4096 terms. Based on the results, it is safe to say that for the type of
problems considered, the BigNumber-Stehfest algorithm is a reliable method. (It is not clear,
however, whether it is the most economical one.)
Some additional information, not directly related to the issue of numerical inversion was also
obtained. In particular, it can be seen from the results that the Lz package of U. Graf [25]
provides the symbolic inverse correctly, whenever it comes up with an answer. The built-in
InverseLaplaceTransform function of Mathematica does not perform so well. Though the
current (4.1) version is better than the previous one, it still inserts unnecessary DiracDelta and
UnitStep terms, sometimes just superfluous, sometimes mathematically incorrect. This
inconvenience is relatively easy to filter out. For several cases, however, the answer seems to be
perfectly reasonable, but is in fact incorrect. For instance, in problem A_S 29.3.101 where the
image function is F(z) = ln(z)/(1+z2), the symbolic inverse provided by Mathematica 4.1
contains a surplus “EulerGamma” .
Results for Set B
All tests of this set could be solved using the BigNumber-Stehfest algorithm. We note that for
our test problem B-1a, reference [5] contains a typographical error (a negative sign is missing).
After correcting this error, the inversion algorithm reproduced perfectly the infinite series
12
solution, up to 17 digits. (We note that a similar typo was found regarding test B-2a also taken
from the table included in Schiff’s otherwise excellent treatise [5]).
Results for Set C
The results of the BigNumber-Stehfest algorithm on these problems are excellent as far as we
understand that in some of the cases the most reliable way to calculate the result is the numerical
Laplace inversion and hence some of the results cannot be corroborated independently. As a byproduct of our calculations, we were able to show that De Chant‘s [31] solution is correct, and
Lin [32] was right in correcting a previous typo in the literature.
Results for Set D
The BigNumber-Stehfest algorithm solves the inversion problem for the Heaviside-Mikusinski
image functions, whenever the original exists and is (infinitely) differentiable (D-1 to D-4). In
case of D-5 we had no prior information whether the inverse exists or not, but the results indicate
the inverse does exists. We have no reason to doubt that the results are exact up to 17 significant
digits. For test D-6, that is F ( z ) = exp(− z 3 / 2 ) we got the message of “no evidence of
convergence” for all times, in accordance with the fact, that the inverse does not exist, neither in
the strict Laplace sense nor in the broader Heaviside-Mikusinski sense.
In the case of test D-7 the message “no evidence of convergence” appears at t = 1/10 and t = 1,
but at t = 10 the algorithm “converges” to the right answer f(t) = -1/10, or at least the first 17
digits are already correct at N = 512 terms causing the algorithm to stop. This convergence is,
however, apparent. In fact, at very large number of terms in the Stehfest series the result would
destabilize. Intuitively: the existence of one non-differentiable point at t = τ does ultimately
prevent the BigNumber-Stehfest algorithm to converge at any t, but the effect may be masked, if
t − τ is large enough.
13
Availability and Additional Information
The WEB site http://pumpjack.tamu.edu/~valko/Nil serves as a starting point for additional information.
•
It provides the output of the runs described in this paper (in HTML format.)
•
A list of more than 800 relevant publications going back to 1795 are also made available.
•
Interested readers can download the original Mathematica programs containing all the test Sets and their
solutions with the BigNumber-Stehfest algorithm.
•
A user friendly Java program can be launched directly from the WEB site. This “non-professional” version
of the BigNumber-Stehfest method, programmed in Java, can be run on virtually any current computer and
platform. The limitations of the Java version are the following:
i)
It can use maximum 1024 terms,
ii) The class of F(z) images is limited by the availability of the functions (abs, sqrt, exp, ln, sin, cos, tan,
arcsin, arccos, arctan, sinh, cosh, tanh, arcsinh, arccosh, arctanh .
iii) The Java version uses the Java WEB Start technology and needs a one time installation of the launcher.
(For details see our WEB site.)
iv) Execution is relatively slow, because the multiple precision calculation of the elementary functions is
accomplished with a relatively unsophisticated (but reliable) Taylor expansion algorithm programmed
by Jason Tiscione [34].
It is remarkable, that Java provides the possibility of multiple precision arithmetic in contrast to, for instance, fortran
90/95. The name of the Java class is BigNumber (hence the name of our algorithm.) We found that for not too large
number of digits (in the range of 32 –128), the speed of Java is competitive to Mathematica.
In addition to the BigNumber-Stehfest algorithm, the test problems in Set A were solved using four other inversion
methods due to Durbin [26], Crump [27], Weeks [28] and Piessens [29]. All methods have been programmed in
Mathematica using the so-called default values of the parameters included in a Mathematica “Add-On” package
prepared by A. Mallet [30]. The results of these runs are also available on our WEB site. In the paper we did not
comment on the results, because it is not fair to use the methods without selecting the necessary additional
parameters (first of all, the estimate of the abscissa of convergence) carefully, and differently from image to image
and time point to time point.
Application Example
To illustrate the practical importance of the findings of this paper we take a problem from
underground flow. An injection well in the center of a cylindrical reservoir is considered, where
the permeability of the formation increases away from the well, according to a power law with
constant exponent n>0 (“damaged well”). The mathematical model describing the pressure as
function of location and time is similar to the variable coefficient heat equation and can be
written as:
1 ∂  n +1 ∂p  ∂p
r
=
∂r  ∂t
κ 2 r ∂r 
with boundary conditions
14
1  n +1 ∂p 
 = −1 (constant wellbore rate)
r
∂r  r =1
κ2 
∂p
outer:
= 0 (no flow outer boundary)
∂r r =R
and initial condition:
p (t , r ) t =0 = 0 .
inner:
In the above model r is the dimensionless radius (1 at the wellbore and R at the outer boundary),
t is the dimensionless time and p is the dimensionless pressure. The parameter κ is not
independent from n and R, the choice
2( R n + 2 − 1)
κ2 =
(n + 2)( R 2 − 1)
ensures that the areal average of the dimensionless permeability is unity (see Fig. 2).
Using the Laplace transform technique, the image function of the pressure at any location, r can
be obtained as

 1− n2 
 1− n2 
−n / 2


F (s) = r
 AI β  a1 r  + BI − β  a1 r 





where I β (x) is the modified Bessel function of the first kind and of order β.
The coefficients A and B are functions of the Laplace variable and are determined from the two
boundary conditions:
2b3κ 2
2b4κ 2
A=
B=−
and
s (b2 b3 − b1b4 )
s (b2 b3 − b1b4 )
In the above formulae the following notation is used:
2
n
α=
, β=
, γ = α − 2β
n−2
n−2
a 0 = κ s , a1 = − a0α , a 2 = a1 R
1−
n
2
n
2
, a3 = nR , a 4 = a 0 R
b1 = a 0 I α (a1 ) − nI β (a1 ) + a 0 I β +1 (a1 )
b 2 = a 0 I −α (a1 ) − nI − β (a1 ) + a 0 I γ (a1 )
b 3 = a 4 I α (a 2 ) − a3 I β (a 2 ) + a 4 I β +1 (a 2 )
b 4 = a 4 I −α (a 2 ) − a3 I − β (a 2 ) + a 4 I γ (a 2 )
(It is interesting to note that for n → 0 and κ → 1 the solution tends to the known solution of the
constant coefficient heat equation given e.g. in [14].)
We consider R = 5000 and n = 0.5 (κ = 7.5212…) and wish to calculate the wellbore pressure (
e.g. p at r = 1) at times t = 1 (when the flow comes from a region commensurate to the
wellbore) and t = 108 (when the effect of injection reaches the outer boundary). The particular
permeability distribution is shown in Fig. 2. Results are summarized in Table 1.
15
10
1
rn
r=R
κ2
0.1
r =1
0.01
0
1000
2000
3000
4000
5000
r
Fig. 2 Permeability normalized by its areal average value ( R = 5000 and n = 0.5 )
Time, t
BigNumber-Stehfest
Stehfest [1]
(standard
double-precision
N = 8 or 16)
Durbin*
[26,30]
(N=2000;
a=0; e=2)
Crump*
[27,30]
(N=30;
a=0; e=2)
Weeks*
[28,30]
(N=100;
a=0)
Piessens
[29,30]
(N=10;
b=1)
1
7.9020786793717900
NA
NA
NA
NA
NA
108
118.69900631000067
118.699017…
1.7… 102
118.699021…
118.2…
118.8…
*Needs evaluation of F(z) at complex z
Table 1 pw ( R = 5000 and n = 0.5 ) calculated from various methods
In Table 1, “NA” indicates that the particular method was unable to provide even an order of
magnitude estimate. As seen from the Table, the inversion at t = 1 turned out to be a difficult
problem, while at t =108 almost all the methods provided reasonable results. Because of the tests
carried out with Test Set A, we can be confident that the results shown in the BigNumberStehfest column are accurate to 17 digits. We are not aware of any other way to calculate those
results with the precision indicated.
16
To illustrate a case, when the BigNumber-Stehfest algorithm does not work, we consider a
related problem. Now we assume, that the well had been injecting with constant rate till time
t = 1, and then it was shut in (“fall-off test”). The Laplace transform of the wellbore pressure is
obtained using the first shifting theorem of Laplace transform: F2(s) = (1-e-s)F(s). We wish to
determine the inverse of this transform at the wellbore (r = 1) at time t = 2, that is when unit
time has elapsed after the shut-in. Table 2 shows the intermediate results from the BigNumberStehfest algorithm, before the message “no evidence of convergence” pops up:
Number of terms*, N
Approximant
8
3.0…
16
2.0…
32
2.95..
64
2.962…
128
2.9612..
256
6.3×1057
* also number of decimal places used in internal calculations
Table 2 Intermediate results of the BigNumber-Stehfest algorithm for F2(s)
We see, that in spite of the apparent convergence at moderate number of terms, the algorithm
ultimately “blows up” (similarly to our Test Problem D-7, for times t = 1/10 and 1). The reason
is, that the f2(t) function is not differentiable at time t = 1 (when the well was shut-in). The
engineer, however, knows this, and does not try to invert F2(s). Rather, she or he applies
superposition to obtain the correct result, taking the difference of two values of the inverse of
F(s): one at time t = 2 and the other at time t = 1:
pw = 10.863251997232379 - 7.902078679371900 = 2.961173317860588;
where the second number comes directly from Table 1, and the first number can be obtained
running the algorithm on F(s) for t = 2.
Therefore, it would be unnecessary and misleading to include F2(s) as a standard test problem
for the class of problems considered in this paper.
Conclusions
In this paper we constructed several sets of test problems for the numerical inversion of the
Laplace transform. We focused on the case when the image function i) can be evaluated (at least
for positive real arguments) with any required precision, and ii) the original is known to be
infinitely many times differentiable. Set A is most likely to reveal any weakness of a suggested
17
algorithm, due to the fact, that it has been created from a representative table of transform pairs
without any additional subjective consideration. In particular, we tested the BigNumber-Stehfest
algorithm and concluded that it gives the required result to 17 significant digits for Set A.
Additional Sets B and C show that in contrary to common perception, the numerical inversion of
the Laplace transform might be the most reliable method to calculate the f(t) function, even if a
more traditional representation is available. Of particular interest is Set D, showing that the
BigNumber-Stehfest algorithm works also for generalized transforms in the sense of the
Heaviside-Mikusinski operator calculus. In addition, the algorithm does not converge to a false
value if the inverse does not exist. Though Sets B, C and D helped us to clarify several issues, it
is clear from the results that Set A is basically sufficient to evaluate a particular numerical
inversion algorithm, as far as the class of problems is restricted in the above sense of conditions
i) and ii). Therefore, we suggest Set A as a standard test set. Detailed results, programs (both in
Mathematica and Java) and additional information are provided on the accompanying WEB site.
The application example taken from underground flow illustrates the practical significance of
defining the right set of test problems.
Acknowledgement
The authors would like to thank Joe Abate (Ridgewood, NJ), Vladimír Vojta (Prague, CZ), Urs
E. Graf, (Biel-Bienne, CH), Benedek Valkó (Budapest, H) for insightful discussions, and J.
Tiscione, D. Smith, D. Bailey (US) for their contributions related to multi-precision calculations.
References
1
Stehfest, H.: Algorithm 368: Numerical inversion of Laplace transforms, Comm. ACM 13 (1) pp 4749 and 624, 1970.
2
Davies, B. and Martin, B.: Numerical inversion of the Laplace transform: A survey and comparison
of methods. J. Comput. 33, 1 (Oct.), pp 1–32, 1979.
3
Duffy D. G.: On the numerical inversion of Laplace transforms: comparison of three new methods on
characteristic problems from applications, ACM Trans. Math. Soft. 19, p333-359, 1993.
4
D'Amore, L., Murli, A. and Rizzardi, M.: An extension of the Henrici formula for Laplace transform
onversion, Inverse Problems 16 pp1441-1456, 2000.
5
Schiff, J. L.: The Laplace Transform. Theory and Applications., Springer, New York, 1999.
6
D'Amore, L., Murli, A. and Rizzardi, M.: Recent developments related to the numerical inversion of a
Laplace Transform function: the real inversion problem, Proceedings of the Intern. Conf. on Inverse
Problems in Engineering: theory and practice, D. Delanunay, Y. Jarny, K.A. Woodbury (Eds.),
Engineering Foundation Publisher, 1996.
7
Gaver, D.P.,Jr.: Observing stochastic processes, and approximate transform inversion, Operations
Research, 14 (3) pp 444-459, 1965.
8
Whiteley, K. S. and Garriga, A.: Derivation of continuous molecular weight distributions from the
generating function, Computational and Theoretical , Polymer Science, 11 (4) pp 319-324, 2001.
18
9
Usábel, M.: Calculating multivariate ruin probabilities via Gaver-Stehfest inversion technique
insurance: Mathematics and Economics, 25, pp133-142,1999
10 Frolov, G.A. and Kitaev, M.Y.: Improvement of accuracy in numerical methods for inverting Laplace
transforms based on the Post-Widder formula, Computers & Mathematics with Applications, 36, (5)
pp 23-34, 1998.
11 Abate J. and Whitt, W.: The Fourier-series method for inverting transforms of probability
distributions, Queueing Systems, 10, 5 – 88, 1992.
12 Weideman, J. A. C.: Algorithms for parameter selection in the Weeks method for inverting the
Laplace transform, Siam J. Sci. Comput., 21, 111-128, 1999.
13 D'Amore, L, Laccetti, G. and Murli, A: An implementation of a Fourier series method for the
numerical inversion of the Laplace transform; ACM Trans. Math. Softw. 25, 3 pp 279 – 305. 1999.
14 Raghavan R. Well test analysis, Prentice-Hall, Englewood-Cliffs, 1993.
15 Motz, L.H.: Drawdowns for leaky-aquifer flow with storage in finite-width sink
Irrigation and Drainage Engineering-Asce, 120 (7) pp 820-827, 1994.
Journal of
16 Gieske, A.: Flow to a well in a water table aquifer: Speeding up the computation of type curves for
low beta values, Ground Water, 37 (4) pp171-174, 1999.
17 Moench, A.F.: Double-porosity models for a fissured groundwater reservoir with fracture skin, Water
Resources Research 20 (4/5) pp 831-846, 1984.
18 Schalla, M. and Weiss, M.: Pharmacokinetic curve fitting using numerical inverse Laplace
transformation. Eur. J. Pharm. Sci. 7(4), 305-309, 1999.
19 Popovic J. : Derivation of Laplace transform for the general disposition deconvolution equation in
drug metabolism kinetics, Experimental & Toxicologic Pathology. 51(4-5) pp409-411, 1999.
20 Fukumura, K., Yamaoka, K., Higashimori, M. and Nakagawa, T.: Analysis Program Based on Finite
Element Method, MULTI(FEM), for Evaluation of Dose-Dependent Local Disposition of Drug in
Liver, Pharm. Sci. 88 (5), 538-543, 1999.
21 Roberts, G. E. and Kaufman, H. Table of Laplace Transforms, Saunders Co., Philadelphia, 1966.
22 Prudnikov, A.P. Brychov, Y.A. and Marichev, O.I. Integrals and Series: Vols 4 and 5, Gordon and
Breach Science Publishers, New York, 1992.
23 Abramowitz, M and Stegun, I. A. Handbook of Mathematical Functions, Dover, N.Y. 1970.
24 Mathematica version 4.1, 2000, http://www.wolfram.com
25 Urs E. Graf: http://www.hta-bi.bfh.ch/~gru/LaplaceAndzTransforms_Pack_Info/e.htm
26 Durbin, F.: Numerical inversion of Laplace transforms: an efficient improvement to Dubner and
Abate's method, Computer Journal 17 (4) 371-376, 1974.
19
27 Crump, K. S.: Numerical inversion of Laplace transforms using a Fourier series approximation,
Journal of the Association for Computing Machinery, 23 (1) 89-96, 1976.
28 Weeks, W.T., Numerical inversion of Laplace transforms using Laguerre functions J. ACM., v. 13,
419-426, 1966.
29 Piessens, R.: Gaussian quadrature formulas for the numerical integration of Bromwich's integral and
the inversion of the Laplace transform', J. Eng. Math., v. 5, 1-9, 1971.
30 Mallet, A.: Numerical Inversion of Laplace Transform, 2000, Mathematica “Add On”
http://www.mathsource.com/Content/Enhancements/Algebraic/0210-968
31 De Chant, L. J.: Analytical Solutions for Diffusive Finite Reservoir Problems Using a Modified
Orthogonal Expansion Method, Math Comput Modelling, 28, pp 73-86, 1998.
32 Lin, G. D.: On the Mittag-Leffler distributions, Journal of Statistical Planning and Inference
(Elsevier) , vol 74 , pp1-9, 1998.
33 Mikusinski, J. and Boehme T.K.: Operational calculus, vol 2 (2nd ed.) Pergamon, Oxford 1987.
34 Jason Tiscione’s WEB site http://www.geocities.com/SiliconValley/2548/
20