AM  Vol.8 No.4 , April 2017
The Role of High Precision Arithmetic in Calculating Numerical Laplace and Inverse Laplace Transforms
ABSTRACT
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.

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 f ( t ) to be transformed, and the uppercase letter C ( s ) to denote its Laplace transform, for example, L { f ( t ) } = C ( s ) . If the closed form of C ( s ) inversion is unknown, the original C ( s ) compared with numerical solution C ˜ ˜ ( s ) 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 f ( t ) be a function defined for t 0 . Then the integral

L { f ( t ) } = 0 e s t f ( t ) d t (1)

is said to be the Laplace transform of f ( t ) , provided the integral converges. The symbol L is the Laplace transformation operator, which act on the function f ( t ) and generates a new function, C ( s ) = L { f ( t ) } .

2.2. Inverse Laplace Transform

If C ( s ) represents the Laplace transform of a function f ( t ) , that is, L { f ( t ) } = C ( s ) , then f ( t ) is the inverse Laplace transform of C ( s ) and f ( t ) = L 1 { C ( s ) } . The inverse Laplace transform L 1 { C ( s ) } is uniquely determined in the sense that if C ( s ) = G ( s ) and f ( t ) and g ( t ) are conti- nuous functions, then f ( t ) = g ( t ) . 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 f ˜ ( t ) used for the numerical approximation to f ( t ) (numerical inversion of the Laplace transform C ( s ) ), and C ˜ ( s ) for the numerical Laplace transform of f ( t ) .

These results, together with the following well known properties which provide very useful numerical checks, may be applied to numerical algorithms and corresponding software.

Integralproperty , 0 e s t f ( t ) d t = C ( 0 ) (2)

Initialvaluetheorem , lim t 0 f ( t ) = lim s s C ( s ) (3)

Finalvaluetheorem , lim t f ( t ) = lim s 0 s C ( s ) (4)

If X is the random variable with the probability density function f and the cumulative distribution function F , this gives

C ( 0 ) = 0 e s t d F ( t ) = 0 e s t f ( t ) d t = 1 (5)

2.4. Numerical Laplace Double Transformation Technique

We define the following double transformation technique for the Laplace transform of the inversion:

C ˜ ˜ ( s ) = L { L 1 { C ( s ) } } (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 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 C ( s ) .

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 N 1 , we get

f ^ N 1 ( t ) = L N 1 1 { C ( s ) } (7)

3) Take the Numerical Laplace Transform of f ^ N 1 ( t ) , resulting in

L { f ^ N 1 ( t ) } = C ˜ ˜ N 1 ( s ) (8)

4) Compare the functions C ( s ) and C ˜ ˜ N 1 ( s ) , and define the error-function as:

ε N 1 ( s ) = | C ( s ) C ˜ ˜ N 1 ( s ) | (9)

5) Repeat the process with some other precision level N 2 .

6) Compare ε N 1 ( s ) and ε N 2 ( s ) . 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 C ( s ) is f ( t ) , defined such that the following must hold:

C ( s ) = 0 e s t f ( t ) d t (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] ,

J 0 ( t ) = n = 0 ( 1 ) n t 2 n 2 2 n ( n ! ) n , and J 0 ( s ) = 1 1 + s 2 (11)

4. Gaver-Stehfest Algorithm of Inverse Laplace Transforms

The Gaver-Stehfest method [15] uses the summation:

f ( t ) ln ( 2 ) t n = 1 L K n F ( n ln ( 2 ) t ) , (12)

where F ( ) is the Laplace transform of f ( t ) . The coefficients K n depends only on the (necessarily even) number of expansion terms, L , given by:

K n = ( 1 ) n + L / 2 k = n + 1 2 min ( n , L / 2 ) k L / 2 ( 2 k ) ! ( L / 2 k ) ! k ! ( k 1 ) ! ( n k ) ! ( 2 k n ) ! (13)

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 t = t 1 , t 2 , , t m . We define by f ˜ ( t ) = L 1 { C ( s ) } the numerical solution. Then E gives the root-mean-square deviation between the analytical and numerical solutions for the t values [19] :

E = ( i = 1 m ( f ( i / 2 ) f ˜ ( i / 2 ) ) 2 / m ) 1 / 2 , i = 1 , 2 , m (14)

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 precision 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.

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: N = 16 (double precision), N = 32 , N = 64 , N = 128 and N = 256 . 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:

C ( s ) = 1 s + 2 , with f ( t ) = e 2 t , (15)

where f ( t ) = e 2 t 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 5.66 × 10 6 .

Figure 1. Given C ( s ) = 1 s + 2 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = e 2 t .

calculations in double precision are at the order of 10 6 . In the plot, the inver-

sion f ˜ ( t ) = L 1 { 1 s + 2 } and the exact solution f ( t ) = e 2 t visibly overlapp-

ing even for N = 16 . The precision from N = 32 to N = 256 gives the sequence of improvements. High accuracy is borne out by the errors at the order 10 32 and 10 69 respectively to the precision N = 128 and N = 256 .

Example 2. The next test function inverted is

C ( s ) = 1 ( s + 0.2 ) 2 + 1 , with f ( t ) = e 0.2 t sin ( t ) , (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 10 2 . The precision level N = 32 and N = 64 show improvements in accuracy over double precision. In the plot, the inversion for N = 128 and the exact solution are visibly overlapping. High accuracy is borne out in Table 2 by the errors at the order 10 10 and 10 35 respectively to the precision level N = 128 and N = 256 .

Example 3. The Figure 3 provides the inverse results for the function

C ( s ) = 1 s 2 + 1 , with f ( t ) = sin ( t ) (17)

Steady improvement of the answer is observed through N = 128 . By N = 128 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 t .

Figure 2. Given C ( s ) = 1 ( s + 0.2 ) 2 + 1 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = e 0.2 t sin ( t ) .

Figure 3. Given C ( s ) = 1 s 2 + 1 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = sin ( t ) .

Example 4. Figure 4 shows the inverse results to the function

C ( s ) = s 2 1 ( s 2 + 1 ) 2 , with f ( t ) = t cos ( t ) (18)

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.

Figure 4. Given C ( s ) = s 2 1 ( s 2 + 1 ) 2 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = t cos ( t ) .

Example 5. Consider the inverse of the function

C ( s ) = tan 1 ( 1 s ) , with f ( t ) = sin ( t ) t (19)

In Figure 5 the precision level N = 32 and N = 64 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 N 64 and the exact solution are visibly overlapping.

Example 6. The inverse results to the function

C ( s ) = 1 1 + s 2 , with f ( t ) = J 0 ( t ) (20)

is plotted in Figure 6. Varying precision levels N = 32 , N = 64 , N = 64 and N = 128 show successive improvements in accuracy over double precision. The numerical inversion for N 128 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 N = 32 , increasing the number of terms to L = 64 . The result is slightly better, close to 2.26 × 10 2 , that for N = 64 and L = 64 . But using the precision level N = 32 and the number of terms in the approximation L = 128 (Figure 7), the error will be at the order of 10 39 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

C ( s ) = e 1 / s s 3 / 2 , with f ( t ) = sin ( 2 t ) π (21)

Figure 5. Given C ( s ) = tan 1 ( 1 s ) . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = sin ( t ) t .

Figure 6. Given C ( s ) = 1 1 + s 2 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = J 0 ( t ) , Bessel function of order 0.

In Figure 8 steady improvement is observed through N = 64 with the error at the order of 10 13 (13 decimal places accuracy, Table 2).

As can be seen from Table 2, all test functions have low level of accuracy for N = 16 , N = 32 and for N = 64 (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 C ( s ) = 1 1 + s 2 . The precision level is N = 32 and the number of expansion terms is L = 128 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = J 0 ( t ) . The algorithm was unable to provide even an order of magnitude estimate.

Figure 8. Given C ( s ) = e 1 / s s 3 / 2 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = sin ( 2 t ) π .

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, 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- Stehfest implementation for the function f ( t ) = e 0.2 t sin ( t ) having the La- place transform C ( s ) = ( ( s + 0.2 ) 2 + 1 ) 1 . The solutions are observed while increasing the number of expansion terms L . However, there is a limitation to adding additional terms.

Let N = 16 . Increasing the number of terms in the computation to L 64 it appears that the numerical inversion becomes unstable and functions are dominated by numerical errors. If N = 32 and N = 64 , the numerical inver-

Table 3. Numerical errors of the inverse Laplace transform f ˜ ( t ) = L 1 { ( ( s + 0.2 ) 2 + 1 ) 1 } as a function of the number of expansion terms L and precision N , Values 3.44e−2 3.44 × 10 2 .

Figure 9. Numerical errors E of the Laplace transform inversion f ˜ ( t ) = L 1 { ( ( s + 0.2 ) 2 + 1 ) 1 } as a function of the number of expansion terms L and precision N . Estimated data presented in terms of MAX ( Log 10 E , 0 ) , here 0 denotes the algorithm has failed.

sion becomes unstable for L 128 . If N = 128 , the numerical inversion becomes unstable for L 256 . On the other hand there is only slight improvement if the precision exceeds the number of terms, N > L .

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, L = 256 , N = 128 is a useless result while N = 256 gives high accuracy. The example, N = 16 , L = 16 gives the same accuracy as N = 256 , L = 16 .

Similar conclusions may be drawn from the results reported in Table 4 and in Figure 10.

Figure 10. Numerical errors E of the Laplace transform inversion f ˜ ( t ) = L 1 { ( e 1 / s / s 3 / 2 } as a function of the number of expansion terms L and precision N Estimated data presented in terms of MAX ( Log 10 E , 0 ) , here 0 denotes the algorithm has failed.

Table 4. Numerical errors of the inverse Laplace transform f ˜ ( t ) = L 1 { e 1 / s / s 3 / 2 } .

The error estimates are given for the example

f ( t ) = sin ( 2 t ) π , with C ( s ) = e 1 / s s 3 / 2 (22)

Let N = 16 . Increasing the number of terms in the computation to L 64 we can see that the numerical inversion becomes unstable. If N = 32 and N = 64 , the numerical inversion failed for L 128 . In case of N = 128 , the numerical inversion becomes unstable for L 256 . There is only a slight improvement 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.

7. Numerical Computation of the Direct Laplace Transform

The Laplace transform of a function f ( t ) is defined by (1) on the interval [ 0 , ] . The problem of an infinite upper limit of integration may be removed by the substitution t = ln ( u ) , d t = u 1 d u which replaces infinite by finite limits.

When t = 0 , u = 1 and when t , u 0 . Then

0 e s t f ( t ) d t = 0 e ln ( u s ) f ( ln ( u ) ) u 1 d u = 0 1 u s 1 f ( ln ( u ) ) d u (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 [ a , b ] , an even n is chosen such that the function is adequately smooth over each subinterval [ x j , x j + 1 ] where

x j = a + j h for all j { 0 , 1 , 2 , , n } with h = b a n . In particular, x 0 = a and x n = b . Then, the composite Simpson’s Rule is given by [25] :

a b f ( x ) d x h 3 [ f ( x 0 ) + 2 j = 1 n / 2 1 f ( x 2 j ) + 4 j = 1 n / 2 f ( x 2 j 1 ) + f ( x n ) ] (24)

Applying this to the transformed integrand from the Equation (23) we get u j = j h for all j { 0 , 1 , 2 , , n } with h = 1 n . Therefore,

C ( s ) 1 3 n [ 0 s 1 f ( ln ( 0 ) ) + 2 j = 1 n / 2 1 u 2 j s 1 f ( ln ( u 2 j ) ) + 4 j = 1 n / 2 u 2 j 1 s 1 f ( ln ( u 2 j 1 ) ) + 1 s 1 f ( ln ( 1 ) ) ] (25)

The basic Simpson’s rule formula divides the interval [ a , b ] of integration into two pieces. To apply the composite Simpson’s rule, the interval [ a , b ] must

be divided into an even number of subintervals n = 2 m . Then h = b a n = b a 2 m .

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 n 4 at worst. The error committed by the composite Simpson’s rule is bounded in absolute value by [25] .

E n = ( b a ) 5 180 n 4 max | f 4 ( ξ ) | , a < ξ < b (26)

For numerical Laplace transform in the obtaining approximate integral by this rule, some modifications must be made. First, the term 0 s 1 f ( ln ( 0 ) ) in (25) must be addressed. This term addresses the behaviour of the function at infinity. If the Laplace transform exists, the lim t e s t f ( t ) 0 , meaning that the exponential 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 0 s 1 f ( ln ( 0 ) ) vanishes and may be dropped from the formulation.

Next, we need to examine the last term, 1 s 1 f ( ln ( 1 ) ) = f ( 0 ) . Evaluating this term should require that the function be defined at t = 0 . This can however prove problematic since many applications of the Laplace transform result in t -domain functions that have singularities at t = 0 . As such, we change the domain of integration to be ( 0 , 1 ϵ ) as opposed to ( 0,1 ) , where ϵ is the machine epsilon depending on the precision level used. For example, in double precision ( N = 16 ) , ϵ 2.22 × 10 16 .

Therefore we have:

C ( s ) 1 3 n [ 2 j = 1 n / 2 1 u 2 j s 1 f ( ln ( u 2 j ) ) + 4 j = 1 n / 2 u 2 j 1 s 1 f ( ln ( u 2 j 1 ) ) + u n s 1 f ( ln ( u n ) ) ] , (27)

where n is the number of subintervals, h = 1 ϵ n , 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 u depend on s in the Equation (27). As such, the function evaluations, f ( ln ( u j ) ) , need not be evaluated every time a new s value is calculated. This is especially useful in the double transformation calculations because each evaluation of the function f ^ ( t ) is the numerical inversion of C ( s ) at some point t . 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 J 0 ( t ) 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 f ( t ) = n = 0 ( 1 ) n t 2 n 2 2 n ( n ! ) n . Compute C ˜ ( s ) = L { f ( t ) } . The determining function is C ( s ) = 1 1 + s 2 .

C ( s ) = L { J 0 ( t ) } and the exact solution C ( s ) = 1 1 + s 2 for Bessel function J 0 ( t ) is shown in Table 5.

Example 9. Consider the following Laplace transform example

f ( t ) = sin ( 2 t ) π , and C ( s ) = e 1 / s s 3 / 2 (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

f ( t ) = e 2 t , with C ( s ) = 1 s + 2 , (29)

where C ( s ) is the exact Laplace transform solution. Let the number of subintervals n = 2 m = 256 . From the Equation (26) we have

E n 1 180 × 256 4 max | f ( 4 ) ( u ) | (30)

The integral 0 1 u s 1 f ( ln ( u ) ) d u is used, where u = e t . Now f ( u ) = u 2 . This yields

0 1 u s 1 f ( ln ( u ) ) d u = 0 1 u s 1 u 2 d u = 0 1 u s + 1 d u (31)

The integrand f ( u ) = u s + 1 has four continuous derivatives and f ( 4 ) ( u ) = s ( s 1 ) ( s + 1 ) ( s 2 ) u s 3 .

Because 0 u 1 , then for example if s = 0.7 ,

E n 1 180 × 256 4 max | f ( 4 ) ( u ) | = 1 180 × 256 4 × 0.464 6 × 10 13

Table 5. Comparison of the Laplace transform C ˜ ( s ) = L { J 0 ( t ) } with the exact solution C ( s ) = 1 1 + s 2 for Bessel function J 0 ( t ) . Values 5.26e−2 5.26 × 10 2 .

Table 6. Comparison of the Laplace transform C ˜ ( s ) = L ( sin ( 2 t ) π ) with the exact solution C ( s ) = e 1 / s s 3 / 2 . Values 5.62E−6 5.62 × 10 6 .

Figure 12. Given f ( t ) = sin ( 2 t ) π . Compute C ˜ ( s ) = L ( sin ( 2 t ) π ) . The determining function is C ( s ) = e 1 / s s 3 / 2 .

Next we compare this theoretical error bound result with the numerical error obtaining by two step double transformation technique, on

C ˜ ˜ ( s ) = L { L 1 { 1 s + 2 } } (32)

Figure 13 displays the numerical error ε N ( s ) in varying precision levels for some range of s 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 t . Each time N is doubled, ε N ( 1 ) and ε N ( 2 ) are nearly halved. Examining lower values of s reveals that the difference between the various ε N ( s ) is not as pronounced as for s = 1 or s = 2 . 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 N = 64 since the double transformation errors are not showing a worthwhile improvement for the large increase in computational cost.

The error at s 0.7 ( log 2 s = 0.5 ) for N = 250 subintervals. For the precision level 64 the error is at the order of 10 17 , and for the precision level 128 the accuracy is much better, at the order 10 30 .

For the precision levels 16 and 32 the accuracy are at the order 10 7 and 10 12 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 ε N ( s ) = | C ( s ) C ˜ ˜ N ( s ) | for double transformation C ˜ ˜ N ( s ) = L { L 1 { 1 s + 2 } } in varying precision levels.

Example 11. This example illustrates two steps of the numerical double transformation calculation

C ˜ ˜ ( s ) = L { L 1 { e 1 / s s 3 / 2 } } (33)

Figure 14 depicts ε N ( s ) in varying precision levels for some range of s 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 ε N ( s ) = | C ( s ) C ˜ ˜ N ( s ) | for double transformation C ˜ ˜ N ( s ) = L { L 1 { e 1 / s s 3 / 2 } } in varying precision levels.

Table 7. Comparison of C ˜ ˜ ( s ) = L { L 1 { C ( s ) } } for Bessel function in Table 1 with the exact solution C ( s ) = 1 1 + s 2 . Values 7.3e−3 7.3 × 10 3 .

Figure 15. Given C ( s ) = 1 1 + s 2 . Compute C ˜ ˜ = L { L 1 { C ( s ) } } .

Example 13. Numerical double transformation results for the function C ( s ) = e 1 / s s 3 / 2 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 C ( s ) = e 1 / s s 3 / 2 .

Let s = 2 ( log 2 s = 1 ) and the precision level N = 32. The results illustrated in the four plots (Figure 17) correspond to different number of subintervals n in composite Simpson’s rule to approximate the Laplace transform. The error estimations are at the order 10 5 , 10 6 , 10 7 and 10 8 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 N . The number of subintervals used in the direct Laplace transform calculations is L = 512 . 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 N be the precision level. We approximate the relative CPU time t ( N ) for inverse and direct Laplace transform algorithms by the polynomial p ( N ) of degree 6. There are seven coefficients and the polynomial is:

p ( N ) = p 1 x 6 + p 2 ) x 5 + p 3 x 4 + p 4 x 3 + p 5 x 2 + p 6 x + p 7 (34)

For inverse Laplace transform algorithm, the coefficients in the polynomial p ( N ) of degree 6 are:

Figure 16. Given C ( s ) = e 1 / s s 3 / 2 . Compute C ˜ ˜ = L { L 1 { C ( s ) } } .

Figure 17. Given C ( s ) = e 1 / s s 3 / 2 . The precision level N = 32 . 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 C ˜ ( s ) = with original Laplace transform C ( s ) = e 1 / s s 3 / 2 . Values 4.3e−07 4.3 × 10 7 .

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

p ( 6 ) = ( p 1 , p 2 , p 3 , p 4 , p 5 , p 6 , p 7 ) = ( 4.673630 e 13 , 9.235484 e 10 , 4.371448 e 07 , 2.108954 e 04 , 1.350226 e 02 , 5.725313 e 01 , 6.348896 e + 00 ) .

For the direct Laplace transform, the coefficients in the approximation polynomial p ( N ) of degree 6 are:

p ( 6 ) = ( p 1 , p 2 , p 3 , p 4 , p 5 , p 6 , p 7 ) = ( 6.42651 e 17 , 2.43029 e 14 , 3.199151 e 10 , 4.079160 e 07 , 1.878244 e 04 , 9.7309 e 03 , 9.178160 e 01 ) .

In Figure 18 we plotted relative CPU time for inverse Laplace transform and for direct Laplace transform as a function of the precision level N .

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 N . 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 C ( s ) = s 2 1 ( s 2 + 1 ) 2 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = t cos ( t ) .

Figure 19. Given C ( s ) = s 2 1 ( s 2 + 1 ) 2 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = t cos ( t ) . 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 C ( s ) = s 2 1 ( s 2 + 1 ) 2 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = t cos ( t ) . Two screenshots are given for 16 and 256 precision levels.

Figure 21. Given C ( s ) = s 2 1 ( s 2 + 1 ) 2 . Compute f ˜ ( t ) = L 1 { C ( s ) } . The determining function is f ( t ) = t cos ( t ) . Two screenshots are given for 512 and 1000 precision levels.

Previously Figure 4 illustrated a good approximation for t values, up to 35, if the precision level N = 128 .

Let increase t up to 100. The results in Figure 19 leads to unstable solutions, for the number of expansion terms L = 512 , if the precision level much smaller than L , as illustrated for N = 128 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 N equals 16, 256, 512 and 1000. The same number of expansion terms L used as the precision level N . The accuracy is very poor for N is 16 and 256, at order at most roughly 10 2 . We see significant improvements in accuracy as precision level increased to 512 and 1000. They are at order at least roughly 10 25 and 10 139 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

Cite this paper
Krougly, Z. , Davison, M. and Aiyar, S. (2017) The Role of High Precision Arithmetic in Calculating Numerical Laplace and Inverse Laplace Transforms. Applied Mathematics, 8, 562-589. doi: 10.4236/am.2017.84045.
References
[1]   Cengel, Y. and Palm, W. (2013) Differential Equations for Engineers and Scientists. McGraw-Hill, New York.

[2]   Edwards, C. and Penney, C. (2015) Differential Equations and Boundary Value Problems. 5th Edition, Pearson, London.

[3]   Polyanin, A.D. (2016) Handbook of Linear Partial Differential Equations for Engineers and Scientists. 2nd Edition, Chapman and Hall/CRC, UK.
https://doi.org/10.1201/b19056

[4]   Cohen, A.M. (2007) Numerical Methods for Laplace Transform Inversion. Springer, Berlin.

[5]   Kreyszig, E. (2011) Advanced Engineering Mathematics. 10th Edition, Wiley, Hoboken.

[6]   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.

[7]   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.
https://doi.org/10.1080/03461230410000565

[8]   Kao, E. (1997) An Introduction to Stochastic Processes. Duxbury Press, Pacific Grove, CA.

[9]   Nadarajah, S. and Kotz, S. (2006) On the Laplace transform of the Pareto Distribution. Queueing Systems, 54, 243-244.
https://doi.org/10.1007/s11134-006-0299-1

[10]   Davison, M. and Doeschl, A. (2004) A Hyperbolic PDE with Parabolic Behavior. SIAM Reviev, 46, 115-127.
https://doi.org/10.1137/S0036144502409007

[11]   Cruz-Baez, D.I. and Gonazalex-Rodriguez, J.M. (2008) A Different Approach for Pricing Asian Options. Applied Mathematics Letters, 21, 303-306.

[12]   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.
https://doi.org/10.1155/2010/263451

[13]   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.

[14]   Kuhlman, K.L. (2013) Review of Inverse Laplace Transform Algorithms for Laplace-Space Numerical Approaches. Numerical Algorithms, 63, 339-355.
https://doi.org/10.1007/s11075-012-9625-3

[15]   Stehfest, H. (1970) Algorithm 368: Numerical Inversion of Laplace Transform. Communications of the ACM, 13, 47-49.
https://doi.org/10.1145/361953.361969

[16]   Kuznetsov, A. (2013) On the Convergence of the Gaver—Stehfest Algorithm. SIAM Journal of Numerical Analysis, 51, 2984-2998.
https://doi.org/10.1137/13091974X

[17]   Abate, A. and Valko, P. (2004) Multi-Precision Laplace Transform Inversion. International Journal for Numerical Methods in Engineering, 60, 979-993.
https://doi.org/10.1002/nme.995

[18]   Duffy, D.G. (2015) Green’s Functions with Applications. 2nd Edition, CRC Press, Boca Raton, FL.
https://doi.org/10.1201/b18159

[19]   Davies, B. and Martin, B. (1979) Numerical Inversion of the Laplace Transform, a Survey and Comparison of Methods. Journal of Computational Physics, 33, 1-32.

[20]   Murli, A. and Rizzardi, M. (1990) Algorithm 682: Talbot’s Method for the Laplace Inversion Problem. ACM Transactions on Mathematical Software, 16, 158-168.
https://doi.org/10.1145/78928.78932

[21]   Bailey, D.H., Hida, Y., Li, X.S and Thompson, B. (2002) ARPREC: An Arbitrary Precision Computation Package. Tech. Rep. LBNL-53651, Lawrence Berkley National Lab.

[22]   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.

[23]   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.

[24]   Abramowitz, M. and Stegun, I.A. (1964) Handbook of Mathematical Functions. National Bureau of Standards.

[25]   Atkinson, K.E. (1989) An Introduction to Numerical Analysis. 2nd Edition, John Wiley & Sons, Hoboken.

 
 
Top