The Role of High Precision Arithmetic in Calculating Numerical Laplace and Inverse Laplace Transforms

In order to find stable, accurate, and computationally efficient methods for performing the inverse Laplace transform, a new double transformation approach is proposed. To validate and improve the inversion solution obtained using the Gaver-Stehfest algorithm, direct Laplace transforms are taken of the numerically inverted transforms to compare with the original function. The numerical direct Laplace transform is implemented with a composite Simpson’s rule. Challenging numerical examples involving periodic and oscillatory functions, are investigated. The numerical examples illustrate the computational accuracy and efficiency of the direct Laplace transform and its inverse due to increasing the precision level and the number of terms included in the expansion. It is found that the number of expansion terms and the precision level selected must be in a harmonious balance in order for correct and stable results to be obtained.


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 complicated 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 significantly corrupt the results, rendering these algorithms impractical to apply.Extended 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  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

( )
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.

Laplace Transform
Let ( ) is said to be the Laplace transform of ( ) The symbol  is the Laplace transformation operator, which act on the func- tion ( ) f t and generates a new function, ( )

( )
C s represents the Laplace transform of a function ( ) f t is the inverse Laplace transform of ( ) ( ) f t and ( ) This result is known as Lerch's theorem [4].

Numerical Laplace Transform and Inversion
The Laplace transform can be inverted either algebraically or numerically.The notation ( ) If X is the random variable with the probability density function f and the cumulative distribution function F , this gives ( ) ( ) ( )

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 s not t .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 1 N , we get 3) Take the Numerical Laplace Transform of ( ) ( ) ( ) ( ) 5) Repeat the process with some other precision level 2 N .

6) Compare
( ) . 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.

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 ( ) C s is ( ) f t , defined such that the fol- lowing must hold:

Gaver-Stehfest Algorithm of Inverse Laplace Transforms
The Gaver-Stehfest method [15] uses the summation: where ( ) For each function an error E is calculated as the measure for the accuracy of the numerical solution.Let ( ) f t be the analytical solution defined for  the numerical solution.Then E gives the root-mean-square deviation between the analytical and numerical solutions for the t values [19]: The sum (12) doesn't provide convergence due to roundoff errors for the large number of terms L , usually if L exceeds the number of decimal digits of pre- cision N (e.g.L 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.

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  .First, as presented in Table 2, assigned L N = , and the number of terms L in the approximation equals to the digits of precision N .The measures ( 14) are used to calculate the errors E to the t values 0.5, 1, 1.5, ..., 35 for the functions 1 -6, and to the t values 0.5, 1, 1.5, ..., 140 for the function 7. The reason to use this presentation is because t values require a large amount of space to realize the correct accuracy results.
Example 1.The first function explored is: where f t e − = 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 determining function is ( ) calculations in double precision are at the order of 6 10 − .In the plot, the inver- sion and the exact solution ( ) Example 3. The Figure 3 provides the inverse results for the function Steady improvement of the answer is observed through 128 N = . By

N =
the inversion is indistinguishable from the exact solution on the interval shown, but will eventually diverge from the exact solution for some higher values of t .
. The determining function is ( ) ( ) Example 4. Figure 4 shows the inverse results to the function Similar to the previous two test functions, as N 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.
. The determining function is ( ) ( ) Example 5. Consider the inverse of the function In Figure 5 the precision level 32 N = and 64 N = show improvements in accuracy over double precision, but higher N do not result in a visibly better inversion.In the plot, the numerical inversion for 64 N ≥ and the exact solution are visibly overlapping.
Example 6.The inverse results to the function is plotted in Figure 6.Varying precision levels 32 N = , 64 N = , 64 N = and 128 N = show successive improvements in accuracy over double precision.The numerical inversion for 128 N ≥ and the exact solution are visibly overlapping.Naively it would seem that continuing to increase the number of terms L 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 32 N = , increasing the number of terms to 64 L = .The re- sult is slightly better, close to Figure 5.Given ( ) . The determining function is ( ) ( ) Figure 6.Given ( ) . The determining function is ( ) ( ) Bessel function of order 0.
In Figure 8 steady improvement is observed through 64 N = with the error at the order of 13 10 − (13 decimal places accuracy, Table 2).
As can be seen from . The determining function is ( ) ( ) The algorithm was unable to provide even an order of magnitude estimate.accuracy (at least 16 decimal place accuracy) increasing the precision level up to 256 digits.

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, L N = .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- as a function of the number of expansion terms L and precision N , Values 3.44e−2     .There is only a slight provement if the precision exceeds the number of terms N L > .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.

Numerical Computation of the Direct Laplace Transform
The Laplace transform of a function ( ) The problem of an infinite upper limit of integration may be removed by the substitution ( ) = which replaces infinite by finite limits.
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 [ ] , a b , an even n is chosen such that the function is adequately smooth over each subinterval 1 , In particular, 0 x a = and n x b = .Then, the composite Simpson's Rule is given by [25]: Applying this to the transformed integrand from the Equation ( 23) we get The basic Simpson's rule formula divides the interval [ ] Let f be a function with four continuous derivatives.Then, the composite Simpson's rule converges to the true value of the integral with rapidity Evaluating this term should require that the function be defined at 0 t = .This can however prove problematic since many applications of the Laplace transform result in t -domain functions that have singularities at 0 t = .As such, we change the do- main of integration to be ( ) 0,1 −  as opposed to ( ) 0,1 , where  is the ma- chine epsilon depending on the precision level used.For example, in double precision ( ) Therefore we have: ( where n is the number of subintervals, . The determining function is ( ) and the exact solution Example 9. Consider the following Laplace transform example , and Numerical Laplace transform results are shown in Figure 12 and Table 6.

Validation of the Numerical Inversion Using Double
Example 10.The first computation presented in this section is the theoretical error for Numerical Laplace transform of the function The integral ( ) ( ) This yields ( ) ( ) The integrand ( ) has four continuous derivatives and ( ) ( ) ( )( )( ) Because 0 1 u ≤ ≤ , then for example if 0.7 s = , ( ) ( )   ( ) with the exact solution ( ) Next we compare this theoretical error bound result with the numerical error obtaining by two step double transformation technique, on ( )

( )
N s ε is not as pronounced as for 1 s = or 2 s = .This suggests that increasing N does not have a very significant impact on the double transfor- mation accuracy for higher values of t .As such, depending on the needs of the application, one might stop increasing precision past 64 N = since the double transformation errors are not showing a worthwhile improvement for the large increase in computational cost.
The error at ( ) subintervals.For the precision level 64 the error is at the order of 17 10 − , and for the precision level 128 the accuracy is much better, at the order 30 10 − .For the precision levels 16 and 32 the accuracy are at the order 7 10 − and 12 10 − 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.( )  ( ) for Bessel function in Table 1 with the exact solution ( )

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.Canada for research funding.We are also grateful to the anonymous referees for useful suggestions which improved the present work.
results are illustrated in the plots and error estimations.

Figure 4 .
Figure 4. Given ( ) ( ) 2 2 2 1 1 s C s s the number of terms in the approximation 128 L = (Figure7), the error will be at the order of 39 10 due to numeri- cal 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 7.Given ( )

Figure 9 .
Figure 9. Numerical errors E of the Laplace transform inversion ( ) ( ) ( ) { } 1 2 1 0.2 1 f t s − − = + +   as a function of the number of expansion terms L and precision N .Estimated data presented in terms of

Figure 10 .
Figure 10.Numerical errors E of the Laplace transform inversion( )

, a b
of integration into two pieces.To apply the composite Simpson's rule, the interval [ ]

4 n
− at worst.The error committed by the composite Simpson's rule is bounded in absolute value by[25].
transform in the obtaining approximate integral by this rule, some modifications must be made.First, the term .This term addresses the behaviour of the function at infinity.If the Laplace transform exists, the that the ex- ponential dampening term outweighs the value of ( ) f t at infinity.Therefore, ( ) lim t f t →∞ must either exist or oscillate between some finite bounds.As such, it may be concluded that the term need to examine the last term, Figure 11.Given

Figure 13
Figure 13 displays the numerical error

Figure 17
Figure 17.Given ( ) In Figure18we plotted relative CPU time for inverse Laplace transform and for direct Laplace transform as a function of the precision level N .

Figure 18 .
Figure 18.Relative CPU time for inverse Laplace transform (left) and for direct Laplace transform (right) as a function of the precision level N .All times are relative to run in double precision (16 digits).

Table 1 .
Laplace and inverse transforms for the test functions used in numerical calculations.
Function ( )C s

Table 2 .
the Calculation errors of the inverse Laplace transform in arbitrary precision, Values 5.66e−6 Function( )

Table 2
, all test functions have low level of accuracy for 16 N = ,

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 al-Similar conclusions may be drawn from the results reported in Table4 and in gorithm.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,

Table 5 .
Comparison of the Laplace transform ( )

Table 6 .
Comparison of the Laplace transform ( )

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 poly-