The Role of High Precision Arithmetic in Calculating Numerical Laplace and Inverse Laplace Transforms ()
1. Introduction
Laplace transforms play a key role in many applications of mathematics to the fields of engineering, physics, and finance, whenever probability density functions, or linear differential equations or integral equations are involved. Laplace transform techniques may simplify the task of solving systems of differential equations [1] , [2] , [3] , and can be considered in terms of typical applications [4] , [5] .
Numerical inversion of Laplace transform is crucial for many applications. Unfortunately, when considering interesting examples, it is often difficult to find an analytical expression for the inverse Laplace transform. Inverting the Laplace transform is a challenging task. This challenge faced in many application areas including the finding of various performance measures in queueing and related probability models [6] , [7] , [8] , [9] , in solving partial differential equations [10] , and in the pricing and hedging of financial derivatives [11] , [12] , [13] .
This paper investigate the complicated and very interesting relationship between numerical precision and the number of terms in one particular Laplace transform inversion algorithm, the Gaver-Stehfest algorithm, and illustrates this relationship using several carefully chosen numerical examples.
For numerous practical situations the inverse of Laplace transform is com- plicated and either doesn’t have a closed form, or has a solution which cannot be represented by any simple formula, performed even in symbolic software (Maple or Mathematica). An alternative is to use a numerical technique for inversion. One way to choose among various alternative methods is to provide a large set of test problems, and to demonstrate how a specific algorithm works on each of them.
Several algorithms have been proposed for numerical Laplace transforms inversion, see for instance the surveys in [4] and [14] . The Gaver-Stehfest algorithm [15] is one of the most powerful algorithms for this purpose. The convergence of this algorithm has been examined in [16] . Unfortunately despite its theoretical advantages, in many practical applications, numerical inversion often encounters numerical accuracy problems [14] [17] [18] [19] [20] . As such, small rounding errors in computation in standard double arithmetic may signifi- cantly corrupt the results, rendering these algorithms impractical to apply. Ex- tended precision allows to add additional significant figures, and produce results that converge to the solution. Laplace and inverse transforms for the test functions used in numerical calculations are presented in Table 1. These complicated functions are used to test the accuracy of the numerical Laplace transform and its inverse.
In general, lowercase letters used to denote the function
to be transformed, and the uppercase letter
to denote its Laplace transform, for example,
. If the closed form of
inversion is unknown, the original
compared with numerical solution
after double transformation. The results are illustrated in the plots and error estimations.
With the help of the arprec library [21] , C++ and MATLAB numerical class library [22] , [23] applied, to investigate the role that extended precision can play in accuracy of Laplace transform and inversions.
The remainder of the paper is organized as follows. In Section 2, a brief description of the underlying theory is given, to introduce numerical Laplace double transformation technique. Sections 3 and 4, apply the Gaver-Stehfest algorithm to the test functions with various degrees numerical accuracy. In
Table 1. Laplace and inverse transforms for the test functions used in numerical calcu- lations.
Sections 5 and 6, the stability and accuracy of the Laplace transform inversion and the role that the number of expansion terms and precision of the arithmetic play in the numerical results is described. Section 7 describes the algorithm and software implementation of the numerical direct Laplace transform. This section gives background material needed to provide the method, described in the next Section 8. In Section 8 the numerical double transformation technique to confirm agreement of the numerical inversion results is presented. In Section 9 compares the execution time for various arbitrary precision calculations. Concluding remarks are given in Section 10. The Appendix introduces C++ code used to implement numerical Laplace and inverse Laplace transform in arbitrary precision, and illustrates the corresponding graphical user interface with the help of several screenshots.
2. Numerical Laplace Transforms and Their Inverses
2.1. Laplace Transform
Let
be a function defined for
. Then the integral
(1)
is said to be the Laplace transform of
, provided the integral converges. The symbol
is the Laplace transformation operator, which act on the function
and generates a new function,
.
2.2. Inverse Laplace Transform
If
represents the Laplace transform of a function
, that is,
, then
is the inverse Laplace transform of
and
. The inverse Laplace transform
is uniquely determined in the sense that if
and
and
are conti- nuous functions, then
. This result is known as Lerch’s theorem [4] .
2.3. Numerical Laplace Transform and Inversion
The Laplace transform can be inverted either algebraically or numerically. The notation
used for the numerical approximation to
(numerical inversion of the Laplace transform
), and
for the numerical Laplace transform of
.
These results, together with the following well known properties which provide very useful numerical checks, may be applied to numerical algorithms and corresponding software.
(2)
(3)
(4)
If
is the random variable with the probability density function
and the cumulative distribution function
, this gives
(5)
2.4. Numerical Laplace Double Transformation Technique
We define the following double transformation technique for the Laplace transform of the inversion:
(6)
This definition will be used to estimate the accuracy of the Laplace transform inversion, when its closed form is unknown.
After applying the Laplace transform, the problem is said to be in the Laplace domain and it is denoted as a function of
not
. While calculations might be easier in the Laplace domain, leaving the solution in the Laplace domain is typically not useful. To transform the result back into the time-domain, inverse Laplace transforms are used. When the analytical answer is unknown, it is difficult to know whether or not the numerical inversion results are accurate. Moreover, it is hard to judge whether or not changes to the method improve or degrade the inversion estimate.
The following steps are used:
1) Begin with the Laplace domain function
.
2) Compute the numerical inversion using some set of parameters. In this case, we will control the precision level and the number of terms in the approximation. Setting the precision level to
, we get
(7)
3) Take the Numerical Laplace Transform of
, resulting in
(8)
4) Compare the functions
and
, and define the error-function as:
(9)
5) Repeat the process with some other precision level
.
6) Compare
and
. The precision level that provides lower errors is superior, and the difference between the error functions can provide a way of quantifying the accuracy improvement gained from increasing the precision level.
3. Challenging Examples of Laplace Transform and Their Inverses
This demonstration applies the Gaver-Stehfest algorithm [15] to determine the inverse Laplace transforms of the test functions to various degrees of numerical accuracy. The inverse functions and corresponding test functions are presented in Table 1.
The inverse Laplace transform of
is
, defined such that the following must hold:
(10)
Table 1 lists seven functions which have been used as tests of the numerical Laplace transform and inverse transform. Test 6 involves the Bessel function of order 0 [24] ,
(11)
4. Gaver-Stehfest Algorithm of Inverse Laplace Transforms
The Gaver-Stehfest method [15] uses the summation:
(12)
where
is the Laplace transform of
. The coefficients
depends only on the (necessarily even) number of expansion terms,
, given by:
(13)
For each function an error
is calculated as the measure for the accuracy of the numerical solution. Let
be the analytical solution defined for
. We define by
the numerical solution. Then
gives the root-mean-square deviation between the analytical and numerical solutions for the
values [19] :
(14)
The sum (12) doesn’t provide convergence due to roundoff errors for the large number of terms
, usually if
exceeds the number of decimal digits of precision
(e.g.
greater than 16 for standard double precision arithmetics). The software implementation of the numerical Laplace transform and the Laplace transform inversion are given in the Appendix.
5. Accuracy of the Numerical Laplace Transform Inversion as a Function of the Number of Expansion Terms and Precision
Numerical inversion of the Laplace transform is an unstable process, so all algorithms are applied in arbitrary precision. In Table 2 and in Figures 2-8 numerical results related to the test functions presented in Table 1 are reported. The numerical solutions are given for different precisions:
(double precision),
,
,
and
. First, as presented in Table 2, assigned
, and the number of terms
in the approximation equals to the digits of precision
. The measures (14) are used to calculate the errors
to the
values 0.5, 1, 1.5, ..., 35 for the functions 1 - 6, and to the
values 0.5, 1, 1.5, ..., 140 for the function 7. The reason to use this presentation is because
values require a large amount of space to realize the correct accuracy results.
Example 1. The first function explored is:
(15)
where
is the exact solution of the inversion. The results correspond to finding numerical inverse Laplace transform are plotted in Figure 1. For this simple function standard double arithmetic algorithms work well. In Table 2 the
Table 2. Calculation errors of the inverse Laplace transform in arbitrary precision, Values 5.66e−6
.
Figure 1. Given
. Compute
. The determining function is
.
calculations in double precision are at the order of
. In the plot, the inver-
sion
and the exact solution
visibly overlapp-
ing even for
. The precision from
to
gives the sequence of improvements. High accuracy is borne out by the errors at the order
and
respectively to the precision
and
.
Example 2. The next test function inverted is
(16)
and the results correspond to numerical inverse Laplace transform are shown graphically in Figure 2. The accuracy, at double precision, in comparison with the exact solution is very poor. As can be seen from Table 2 the calculations in double precision are at the order of
. The precision level
and
show improvements in accuracy over double precision. In the plot, the inversion for
and the exact solution are visibly overlapping. High accuracy is borne out in Table 2 by the errors at the order
and
respectively to the precision level
and
.
Example 3. The Figure 3 provides the inverse results for the function
(17)
Steady improvement of the answer is observed through
. By
the inversion is indistinguishable from the exact solution on the inter- val shown, but will eventually diverge from the exact solution for some higher values of
.
Figure 2. Given
. Compute
. The determining function is
.
Figure 3. Given
. Compute
. The determining function is
.
Example 4. Figure 4 shows the inverse results to the function
(18)
Similar to the previous two test functions, as
increases, the interval on which the inversion is more accurate gets longer. Since the function is diverging and oscillating, the inaccuracies are more visible than in the previous two figures.
Figure 4. Given
. Compute
. The determining function is
.
Example 5. Consider the inverse of the function
(19)
In Figure 5 the precision level
and
show improvements in accuracy over double precision, but higher
do not result in a visibly better inversion. In the plot, the numerical inversion for
and the exact solution are visibly overlapping.
Example 6. The inverse results to the function
(20)
is plotted in Figure 6. Varying precision levels
,
,
and
show successive improvements in accuracy over double precision. The numerical inversion for
and the exact solution are visibly overlapping.
Naively it would seem that continuing to increase the number of terms
would improve the accuracy of the approximation, since more terms seem to produce results closer to the exact solution. Consider the Bessel example with the precision level
, increasing the number of terms to
. The result is slightly better, close to
, that for
and
. But using the precision level
and the number of terms in the approximation
(Figure 7), the error will be at the order of
due to numerical error dominating the solutions that are obtained this way. The algorithm was unable to provide even an order of magnitude estimate.
Example 7. Consider the inverse of the function
(21)
Figure 5. Given
. Compute
. The determining function is
.
Figure 6. Given
. Compute
. The determining function is
, Bessel function of order 0.
In Figure 8 steady improvement is observed through
with the error at the order of
(13 decimal places accuracy, Table 2).
As can be seen from Table 2, all test functions have low level of accuracy for
,
and for
(except the last example). For all functions the algorithm gave reasonable answer (at least 3 decimal places of accuracy), increasing the precision level up to 128. Finally all functions have high level of
Figure 7. Given
. The precision level is
and the number of expansion terms is
. Compute
. The determining function is
. The algorithm was unable to provide even an order of magnitude estimate.
Figure 8. Given
. Compute
. The determining function is
.
accuracy (at least 16 decimal place accuracy) increasing the precision level up to 256 digits.
6. The Role of the Number of Expansion Terms and Precision in the Numerical Accuracy
In Table 2 the numerical solutions are given for the number of expansion terms equal to the precision,
. Now the accuracy of the numerical inversion is investigated, varying the number of expansion terms and precision. The Table 3 and Figure 9 give the error estimates of the numerical inverse using the Gaver- Stehfest implementation for the function
having the La- place transform
. The solutions are observed while increasing the number of expansion terms
. However, there is a limitation to adding additional terms.
Let
. Increasing the number of terms in the computation to
it appears that the numerical inversion becomes unstable and functions are dominated by numerical errors. If
and
, the numerical inver-
Table 3. Numerical errors of the inverse Laplace transform
as a function of the number of expansion terms
and precision
, Values 3.44e−2
.
Figure 9. Numerical errors
of the Laplace transform inversion
as a function of the number of expansion terms
and precision
. Estimated data presented in terms of
, here 0 denotes the algorithm has failed.
sion becomes unstable for
. If
, the numerical inversion becomes unstable for
. On the other hand there is only slight improvement if the precision exceeds the number of terms,
.
Table 3 clearly shows that there is no point in increasing the precision beyond a point warranted by the accuracy of the number of terms in Gaver-Stehfest algorithm. On the other hand, using a large number of terms will not be of benefit if the precision with which each term is calculated, is insufficient. For example,
,
is a useless result while
gives high accuracy. The example,
,
gives the same accuracy as
,
.
Similar conclusions may be drawn from the results reported in Table 4 and in Figure 10.
Figure 10. Numerical errors
of the Laplace transform inversion
as a function of the number of expansion terms
and precision
Estimated data presented in terms of
, here 0 denotes the algorithm has failed.
Table 4. Numerical errors of the inverse Laplace transform
.
The error estimates are given for the example
(22)
Let
. Increasing the number of terms in the computation to
we can see that the numerical inversion becomes unstable. If
and
, the numerical inversion failed for
. In case of
, the numerical inversion becomes unstable for
. There is only a slight improvement if the precision exceeds the number of terms
. Again as can be seen from Table 3, the number of terms and precision must be in a harmonious balance for good results to be obtained.
Evidently the two plots in Figure 9 and Figure 10, are very much alike. They show similar tracking boundary movements and illustrate whether the algorithm has succeeded in obtaining high-order accuracy or fails due to numerical instability.
7. Numerical Computation of the Direct Laplace Transform
The Laplace transform of a function
is defined by (1) on the interval
. The problem of an infinite upper limit of integration may be removed by the substitution
which replaces infinite by finite limits.
When
,
and when
,
. Then
(23)
The behaviour of the function to be transformed must be considered at the new limits, and the exponential function inside the integral requires special examination in terms of high accuracy.
Compute the Direct Laplace Transform by Composite Simpson’s Rule
For integration over the interval
, an even
is chosen such that the function is adequately smooth over each subinterval
where
for all
with
. In particular,
and
. Then, the composite Simpson’s Rule is given by [25] :
(24)
Applying this to the transformed integrand from the Equation (23) we get
for all
with
. Therefore,
(25)
The basic Simpson’s rule formula divides the interval
of integration into two pieces. To apply the composite Simpson’s rule, the interval
must
be divided into an even number of subintervals
. Then
.
Let
be a function with four continuous derivatives. Then, the composite Simpson’s rule converges to the true value of the integral with rapidity
at worst. The error committed by the composite Simpson’s rule is bounded in absolute value by [25] .
(26)
For numerical Laplace transform in the obtaining approximate integral by this rule, some modifications must be made. First, the term
in (25) must be addressed. This term addresses the behaviour of the function at infinity. If the Laplace transform exists, the
, meaning that the exponential dampening term outweighs the value of
at infinity. Therefore,
must either exist or oscillate between some finite bounds. As such, it may be concluded that the term
vanishes and may be dropped from the formulation.
Next, we need to examine the last term,
. Evaluating this term should require that the function be defined at
. This can however prove problematic since many applications of the Laplace transform result in
-domain functions that have singularities at
. As such, we change the domain of integration to be
as opposed to
, where
is the machine epsilon depending on the precision level used. For example, in double precision
,
.
Therefore we have:
(27)
where
is the number of subintervals,
, and
is the machine epsilon at the precision level.
Our C++ software implementation of the numerical Laplace transform is based on Equation (27). The following improvements were made to speed up the calculations. Notice that only the powers of
depend on
in the Equation (27). As such, the function evaluations,
, need not be evaluated every time a new
value is calculated. This is especially useful in the double transformation calculations because each evaluation of the function
is the numerical inversion of
at some point
. Depending on the precision level, this can be an extremely time consuming step. The implementation of the method is shown in Appendix.
Example 8. The numerical Laplace transform for the Bessel function
is plotted in Figure 11. We used a composite Simpson’s Rule calculation with up to 25,000 subintervals to ensure high accuracy. Comparing Laplace transform
Figure 11. Given
. Compute
. The determining function is
.
and the exact solution
for Bessel function
is shown in Table 5.
Example 9. Consider the following Laplace transform example
(28)
Numerical Laplace transform results are shown in Figure 12 and Table 6.
8. Validation of the Numerical Inversion Using Double Transformation Technique
Example 10. The first computation presented in this section is the theoretical error for Numerical Laplace transform of the function
(29)
where
is the exact Laplace transform solution. Let the number of subintervals
. From the Equation (26) we have
(30)
The integral
is used, where
. Now
. This yields
(31)
The integrand
has four continuous derivatives and
.
Because
, then for example if
,
Table 5. Comparison of the Laplace transform
with the exact solution
for Bessel function
. Values 5.26e−2
.
Table 6. Comparison of the Laplace transform
with the exact solution
. Values 5.62E−6
.
Figure 12. Given
. Compute
. The determining function is
.
Next we compare this theoretical error bound result with the numerical error obtaining by two step double transformation technique, on
(32)
Figure 13 displays the numerical error
in varying precision levels for some range of
beginning at very low values. From Figure 13 it is clear that increasing the decimal precision greatly increases the accuracy of the estimates for low values of
. Each time
is doubled,
and
are nearly halved. Examining lower values of
reveals that the difference between the various
is not as pronounced as for
or
. This suggests that increasing
does not have a very significant impact on the double transfor- mation accuracy for higher values of
. As such, depending on the needs of the application, one might stop increasing precision past
since the double transformation errors are not showing a worthwhile improvement for the large increase in computational cost.
The error at
for
subintervals. For the precision level 64 the error is at the order of
, and for the precision level 128 the accuracy is much better, at the order
.
For the precision levels 16 and 32 the accuracy are at the order
and
respectively. For the precision levels 16 and 32, the accuracy of the answer has deteriorated due to roundoff error. Note that this error of the inverse Laplace transform is ignored in the first step of the calculation as it is much smaller than suggested by the composite Simpson’s rule in the second step of the double transformation algorithm.
Figure 13. Numerical error
for double transformation
in varying precision levels.
Example 11. This example illustrates two steps of the numerical double transformation calculation
(33)
Figure 14 depicts
in varying precision levels for some range of
from the very low values.
Example 12. Numerical double transformation for Bessel function of order 0 is given in Figure 15 and Table 7. The plots in Figure 11 and Figure 15 are indistinguishable. High accuracy is borne out by comparison of the approximation with the exact solution.
Figure 14. Numerical error
for double transformation
in varying precision levels.
Table 7. Comparison of
for Bessel function in Table 1 with the exact solution
. Values 7.3e−3
.
Figure 15. Given
. Compute
.
Example 13. Numerical double transformation results for the function
are given in Figure 16 and in Table 8. The plots in Figure 12 and
Figure 16 are nearly identical. As can be seen from Table 6 and Table 8, the accuracy is very high.
Consider again double transformation results for the function
.
Let
and the precision level N = 32. The results illustrated in the four plots (Figure 17) correspond to different number of subintervals
in composite Simpson’s rule to approximate the Laplace transform. The error estimations are at the order
and
respectively with n = 32, 64, 128 and 512 subintervals.
9. Comparison of Running Time in Arbitrary Precision Calculations
The execution time for the test functions is shown in Table 1. The direct Laplace transform and inverse Laplace transform are computed using different precision levels
. The number of subintervals used in the direct Laplace transform calculations is
. Table 9 gives the relative CPU time of the numerical solutions. All times are relative runs performing in double precision (16 digits). Computer configuration used: AMD 8350 8-Core processor 4 GHz, 8 GB RAM, 240 GB SDD, Windows 10 Pro.
Let
be the precision level. We approximate the relative CPU time
for inverse and direct Laplace transform algorithms by the polynomial
of degree 6. There are seven coefficients and the polynomial is:
(34)
For inverse Laplace transform algorithm, the coefficients in the polynomial
of degree 6 are:
Figure 16. Given
. Compute
.
Figure 17. Given
. The precision level
. Compute the error estimations for different number of subintervals with n = 32, 64, 128 and 512 in composite Simpson’s rule to approximate the Laplace transform.
Table 8. Comparison of
with original Laplace transform
. Values 4.3e−07
.
Table 9. Relative CPU time for inverse and direct Laplace transform algorithms as a function of the precision level. All times are relative to run in double precision (16 digits).
For the direct Laplace transform, the coefficients in the approximation polynomial
of degree 6 are:
In Figure 18 we plotted relative CPU time for inverse Laplace transform and for direct Laplace transform as a function of the precision level
.
10. Conclusion
Laplace Transform applications often require high accuracy beyond IEEE double precision. Common situations involve calculations that are numerically unstable, and even double-double precision is not sufficient to reach the necessary accuracy. Roundoff and underflow/overflow errors that occur during the computations can cause severe stability problems. There are several numerical inverse Laplace transform methods, each successful in some fields. The problem is to find methods successful in stability, accuracy and computational efficiency. Overall, the presented double transformation approach provides an effective way to compare the effectiveness of numerical inversion methods. In order to validate and improve the inversion solution, the numerical direct Laplace transform are used for this inversion, and compared it with the original Laplace transform. We implemented the composite Simpson’s Rule for direct Laplace transform, and investigate the role that extended precision can play obtaining accurate transform inversions. The high precision approach seems to be an effective way to handle the challenging problems dealing with periodic functions. We demonstrated that the level of precision chosen must match algorithms properly. In the Gaver-Stehfest algorithm the balance is between the truncation error and roundoff error. High precision in an inaccurate algorithm yield little benefit, while a potentially highly accurate algorithm may be defeated by roundoff error if inadequate accuracy is used.
Acknowledgements
Matt Davison thanks to Natural Sciences and Engineering Research Council of
Figure 18. Relative CPU time for inverse Laplace transform (left) and for direct Laplace transform (right) as a function of the precision level
. All times are relative to run in double precision (16 digits).
Canada for research funding. We are also grateful to the anonymous referees for useful suggestions which improved the present work.
Appendix: C++ Software Implementation of Numerical Laplace and Inverse Laplace Transform in Arbitrary Precision
The Gaver-Stehfest algorithm of inverse Laplace transform was implemented in multiple precision as described in Section 4. Numerical direct Laplace transform was implemented in multiple precision by the composite Simpson’s Rule, as covered in details in Section 7. The calculations can be performed in C++ up to 1000 places of decimals/1000 significant digits.
The C++ code of the inverse Laplace transform and direct Laplace transform are given below.
1. C++ implementation of the Gaver-Stehfest algorithm for inverse Laplace transform
2. C++ implementation of the composite Simpson’s Rule for direct Laplace transform
3. Demonstration the accuracy of Laplace transform inversionsThis demonstration shows the role that extended precision can play in accuracy of Laplace transform inversions by six screenshots (Figures 19-21)corresponding to the following periodic function. Given
. Compute
. The determining function is
.
Figure 19. Given
. Compute
. The determining function is
. Two screenshots are given for 512 number of terms and the precision level 128 and 256. We failed to get the solution.
Figure 20. Given
. Compute
. The determining function is
. Two screenshots are given for 16 and 256 precision levels.
Figure 21. Given
. Compute
. The determining function is
. Two screenshots are given for 512 and 1000 precision levels.
Previously Figure 4 illustrated a good approximation for
values, up to 35, if the precision level
.
Let increase
up to 100. The results in Figure 19 leads to unstable solutions, for the number of expansion terms
, if the precision level much smaller than
, as illustrated for
and 256. In this case, using too many terms causes rounding error to overtake the numerical solution, and we essentially obtain noise.
Next (Figure 20 and Figure 21) illustrate the solutions correspond to the precision level
equals 16, 256, 512 and 1000. The same number of expansion terms
used as the precision level
. The accuracy is very poor for
is 16 and 256, at order at most roughly
. We see significant improvements in accuracy as precision level increased to 512 and 1000. They are at order at least roughly
and
accordingly.
Submit or recommend next manuscript to SCIRP and we will provide best service for you:
Accepting pre-submission inquiries through Email, Facebook, LinkedIn, Twitter, etc.
A wide selection of journals (inclusive of 9 subjects, more than 200 journals)
Providing 24-hour high-quality service
User-friendly online submission system
Fair and swift peer-review system
Efficient typesetting and proofreading procedure
Display of the result of downloads and visits, as well as the number of cited articles
Maximum dissemination of your research work
Submit your manuscript at: http://papersubmission.scirp.org/
Or contact am@scirp.org