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  ,  ,  , and can be considered in terms of typical applications  ,  .
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  ,  ,  ,  , in solving partial differential equations  , and in the pricing and hedging of financial derivatives  ,  ,  .
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  and  . The Gaver-Stehfest algorithm  is one of the most powerful algorithms for this purpose. The convergence of this algorithm has been examined in  . Unfortunately despite its theoretical advantages, in many practical applications, numerical inversion often encounters numerical accuracy problems      . 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  , C++ and MATLAB numerical class library  ,  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
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  .
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.
If is the random variable with the probability density function and the cumulative distribution function , this gives
2.4. Numerical Laplace Double Transformation Technique
We define the following double transformation technique for the Laplace transform of the inversion:
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
3) Take the Numerical Laplace Transform of , resulting in
4) Compare the functions and , and define the error-function as:
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  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:
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  ,
4. Gaver-Stehfest Algorithm of Inverse Laplace Transforms
The Gaver-Stehfest method  uses the summation:
where is the Laplace transform of . The coefficients depends only on the (necessarily even) number of expansion terms, , given by:
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  :
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:
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
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
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
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
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
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
Figure 5. Given . Compute . The determining function is .
Figure 6. Given . Compute . The determining function is , Bessel function of order 0.
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 , .
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
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
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  :
Applying this to the transformed integrand from the Equation (23) we get for all with . Therefore,
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  .
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:
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
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
where is the exact Laplace transform solution. Let the number of subintervals . From the Equation (26) we have
The integral is used, where . Now . This yields
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
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
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 .
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:
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 .
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.
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 firstname.lastname@example.org
 Abate, J., Chudhuru, G. and Whitt, W. (2000) An Introduction to Numerical Transform Inversion and Its Application to Probability Models. In: Grassmann, W., Ed., Computational Probability, Springer, Berlin, 257-323.
 Badescu, A., Breuer, L., Da Silva Soares, A., Latouche, G., Remiche, M.A. and Stanford, D.A. (2005) Risk Processes Analyzed as Fluid Queues. Scandinavian Actuarial Journal, 2005, 127-141.
 Deakin, A. and Davison, M. (2010) An Analytic Solution for a Vasicek Interest Rate Convertible Bond Model. Journal of Applied Mathematics, 2010, Article ID: 263451.
 Fusai, G. and Roncoroni, A. (2008) Implementing Models in Quantitative Finance: Methods and Cases. In: Avellaneda, M., Barone-Adesi, G., Broadie, M., Davis, M., Derman, E., Klüppelberg, C. and Schachermayer, W., Eds., Springer Finance, Springer-Verlag, Berlin, Heidelberg.
 Krougly, Z.L. and Jeffrey, D.J. (2009) Implementation and Application of Extended Precision in Matlab. Proceedings of the 11th WSEAS International Conference on Mathematical Methods and Computational Techniques in Electrical Engineering, Athens, 28-30 September 2009, 103-108.
 Krougly, Z.L., Jeffrey, D.J. and Tsarapkina, D. (2013) Software Implementation of Numerical Algorithms in Arbitrary Precision. 15th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC), Timisoara, 23-26 September 2013, 131-137.