CC BY
10
Accuracy and Precision Requirements in Probability Models
Zinovi Krougly
� Department of Applied Mathematics, Western University, London, Ontario, Canada N6A5B7 e-mail: [email protected]
Abstract
Numerical Laplace transform and inverse Laplace transform is a challenging task in queueing theory and others probability models. A double transformation approach is used to find stable, accurate, and computationally efficient methods for performing the numerical Laplace and inverse Laplace transform. To validate and improve the inversion solution obtained, direct Laplace transforms are taken of the numerically inverted transforms to compare with the original function. Algorithms provide increasing accuracy as precision level increases. The most promising methods were applied to computational probability models, when there are no closed-form solutions of the Laplace transform inversion.The computational efficiency compared to precision levels is demonstrated for different service models in ../../1 queuing systems.
Keywords: Numerical Laplace Transform, Numerical Laplace Transform Inversion, High Precision Computation, Applications in Probability Models
1 Introduction
Numerical inverting of Laplace transform to get various performance measures is an important techniques in queueing theory and related stochastic models [1], [6], [16]. Laplace transform techniques may simplify the task of solving systems of differential equations [5], and can be considered in terms of typical applications [4], [8].
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 [1], [6], [16].
Several algorithms have been proposed for numerical Laplace transforms inversion, see for instance the surveys in [4] and [13]. The Gaver-Stehfest algorithm [18] is one of the most powerful algorithms for this purpose. The convergence of this algorithm has been examined in [14]. Unfortunately despite its theoretical advantages, in many practical applications, numerical approximation often encounters numerical accuracy problems [1], [9], [11], [12], [13], [15]. As such, small rounding errors in computation in standard double arithmetic may significantly corrupt the results, rendering these algorithms impractical to apply. With extended precision, we are able to add additional significant figures, and produce results that converge to the solution.
We used C++ and Matlab numerical class library [10], [12] and applied ARPREC, Arbitrary precision computation package [3], for numerical implementation of Laplace transform and inversions.
In [9] has been introduced double transformation approach which perform computationally efficient methods for inverse Laplace transform. Challenging numerical examples involving periodic and oscillatory functions, are investigated. It was found that the number of
expansion terms and precision level selected must be in a harmonious balance in order for correct and stable results to be obtained. In this paper we examine the stability and accuracy of the Laplace transform inversion by the Gaver-Stehfest algorithm [18]. The numerical results were first compared to known analytical solutions. The most interesting methods were then applied to computational probability models, where the solution requires numerical Laplace transform inversion. For computation efficiency a composite Simpson�s rule is implemented for the numerical direct Laplace transform [9]. The numerical examples illustrate the computational accuracy and efficiency of the direct Laplace transform and its inverse due to increasing precision level (N) and the number of terms (L) included in the expansion. The remainder paper is organized as follows. In general, we use lowercase letters 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, we compare the original ..(..) and numerical solution ....(..) after double transformation. The results are illustrated in the plots and error estimations.
In Section 2, a brief description and notation of the underlying theory is given. In section 3 introduced numerical computation of the direct Laplace transform by composite Simpson�s rule. In section 4 defined numerical Laplace double transformation technique. Section 5, 6 and 7 illustrated challenges numerical examples and the role of high precision arithmetic in application for probability models. Numerical Laplace transform and their inverses, particular for applications in M/D/1 and M/G/1 queue, are given in sections 8, 9 and 10. We examine the stability and accuracy of the Laplace transform inversion and the role that the number of expansion terms and precision level play in the numerical approximation. We discuss the numerical double transformation technique to confirm agreement of the numerical inversion results. Section 11 demonstrates double transformation technique and precision requirement for approximation of waiting time distribution in M/D/1 queue.
2 Numerical Laplace transforms and their inverses
Let ..(..) be a function defined for ...0. Then the integral
.{..(..)}=..0.........(..)..... (2.1)
is said to be the Laplace transform of ..(..), provided the integral converges. The symbol . is the Laplace transformation operator, which act on the function ..(..) and generates a new function, ..(..)=.{..(..)}.
If ..(..) represents the Laplace transform of a function ..(..), that is, .{..(..)}=..(..), then ..(..) is the inverse Laplace transform of ..(..) and ..(..)=..1{..(..)}. The inverse Laplace transform ..1{..(..)} is uniquely determined in the sense that if ..(..)=..(..) and ..(..) and ..(..) are continuous functions, then ..(..)=..(..).
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 ..(..).
If .. is the random variable with the probability density function .. and the cumulative distribution function .., this gives
..(0)=..0...........(..)=..0.........(..)....=1 (2.2)
3 Numerical computation of the direct Laplace transform
To validate and improve the inversion solution obtained by Gaver-Stehfest algorithm, the numerical direct Laplace transform are used for this inversion, to compare it with the original Laplace transform. To insure high accuracy of approximation, numerical direct Laplace transform are implemented [9] by the composite Simpson�s Rule [2]. We used a composite Simpson�s Rule calculation with large number of subintervals to ensure high accuracy.
The Laplace transform of a function ..(..) is defined by (2.1) on the interval [0,.]. The problem of an infinite upper limit of integration may be removed by the substitution ..=.ln(..),....=...1.... which replaces infinite by finite limits.
When .. = 0, .. = 1 and when ..>., ..>0. Then
..0.........(..)....=.10..ln(....)..(.ln(..))...1....=.10.....1..(.ln(..)).... (3.1)
The behaviour of the function to be transformed, needs to be considered at the new limits, and the exponential function inside the integral requires special examination in terms of high accuracy.
3.1 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 [....,....+1] where ....=..+... for all ...{0,1,2,...,..} with .=........ In particular, ..0=.. and ....=... Then, the composite Simpson�s Rule is given by [2]:
.......(..)......3[..(..0)+2.../2.1..=1..(..2..)+4.../2..=1..(..2...1)+..(....)] (3.2)
Applying this to the transformed integrand from the equation (3.1) we get ....=... for all ...{0,1,2,...,..} with .=1... Therefore,
..(..).13..[0...1..(.ln(0))+2.../2.1..=1..2.....1..(.ln(..2..))+4.../2..=1..2...1...1..(.ln(..2...1))+1...1..(.ln(1))] (3.3)
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 ..=2... Then .=.......=.....2...
4 Numerical Laplace double transformation technique
We define the following double transformation technique for the Laplace transform of the inversion [9]:
....(..)=.{..1{..(..)}} (4.1)
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 precision level and the number of terms in the approximation. Setting precision level to ..1, we get
.....1(..)=...1.1{..(..)} (4.2)
3. Take the Numerical Laplace Transform of .....1(..), resulting in
.{.....1(..)}=......1(..) (4.3)
4. Compare the functions ..(..) and ......1(..), and define the error-function as:
....1(..)=|..(..).......1(..)| (4.4)
5. Repeat the process with some other precision level ..2.
6. Compare ....1(..) and ....2(..). 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 precision level.
To validate and improve the inversion solution obtained by Gaver-Stehfest algorithm, the numerical direct Laplace transform are used for this inversion, to compare it with the original Laplace transform. To insure high accuracy of approximation, numerical direct Laplace transform are implemented [9] by the composite Simpson�s Rule [2]. We used a composite Simpson�s Rule calculation with large number of subintervals to ensure high accuracy.
5 Testing numerical inversion algorithms in arbitrary precision
This demonstration applies the inverse Laplace transforms of the test functions (Table 1) to various type numerical accuracy.
Given ..(..) we want to find ..(..) such that the following must hold:
..(..)=..0.........(..)..... (5.1)
Table 1: Laplace and inverse transforms for test functions used in the numerical calculations
Example 1. Find ...01=..1{..01(..)}, where
..01(..)=1(1+../..).. (..>0and..>0) (5.2)
then we know that
..01(..)=.....(..).....1....... (5.3)
where .(..) is the gamma function.
We are reached higher accuracy for the numerical inverting of the function ..01(..)=1(1+../..).. by multiple precision calculations. The exact inversion is ..01(..)=.....(..).....1........ The results shown in Fig.1 correspond to the parameters .. = 1 and .. = 20. Fig. 1 illustrate a poor approximation for double precision level (N = 16). The numerical inversion evaluated also with precision level and the number of terms (..,..) = (32, 32). Two screenshots are presented in Fig. 2 for the same function, as in Fig. 1. We see significant improvements in accuracy as precision level increased to 256, at order at least roughly 10.73.
Figure 1: Inverse Laplace transform of the function ..01(..)=1/(1+../..).. evaluated in double and multiple precision. The exact and numerical solution with precision level N = 32 are indistinguishable to the eye
Figure 2: Two screenshots for inverse Laplace transform of the function ..01(..)=1/(1+../..).. evaluated in double and multiple precision
Numerous examples exists where there are no closed-form solution of the Laplace transform inversion. For such problems we compare numerical solution .... (s) after double transformation technique (4.1) with the original Laplace transform ..(..).
First we illustrate that double transformation technique in Fig. 3 for the gamma distribution with .. = 1 (exponential distribution) and .. = 2.5. Both Laplace transform and inverting worked extremely well. The errors .. are the following: 2.5.10.5 and 7.45.10.4 corresponding to the plots on the left and the right sides.
Figure 3: Given the original Laplace transform ..01(..)=1/(1+../..)... Compute its numerical approximation ....01(..)=.{..1{1/(1+../..)..}} evaluated for .. = 1.0 and .. = 1.0 (left plot) and .. = 2.5 (right plot). The original and numerical solution are indistinguishable to the eye
6 Numerical inverse Laplace transform of the incomplete Gamma function
The folloving example is quite different from the previous as we cannot express the inverse Laplace transform analytically. The lower incomplete gamma function .. and the upper incomplete gamma function .. are defined by
..(..,..)=1.(..)...0.....1......... (6.1)
..(..,..)=1.(..)..0.....1......... (6.2)
We used the normalized definition of the incomplete gamma function, where ..(..,..)+..(..,..)=1
Example 2. We determine ...02(..)=..1{..02(..)} and ...02(..)=.{..1{..02(..)}}, where
..02(..)=..(..,..)=1.(..)...0.....1........., (6.3)
We obtained the approximation (Fig. 4) for the inverting of the function (6.3) with the parameters .. = 1.0. As easy to prove, the exact solution of inverse Laplace transform is ...(...1), where ..(..) is the Dirac delta function (7.1).
Improvements can be effected with increasing the number of precision .. from double precision (left plot) to precision level 32 and 64 (right plot). We use the number of terms in approximation equals to precision level, .. = ... To obtain a more accurate estimate precision levels 500 and 1000 used as displayed in (Fig. 5).
Figure 4: Inverse Laplace transform of the function ..02(..)=1.(..)...0.....1......... in double precision (left plot) and precision level 32 and 64 (right plot). Note that two plots use different scale
Figure 5: Inverse Laplace transform of the function ..02(..)=1.(..)...0.....1......... in precision level 512 (left plot) and 1000 (right plot). Note that two plots use different scale
The original ..02(..) compares with numerical solution ....02(..)=.{..1{..02(..)}} evaluated after double transformation. So, ...02(..) computed as numerical inversion of ..02(..). Then Laplace transform ....02(..) of ...02(..) compares with the original function ..02(..). The original Laplace transform ..02(..) (Exact) and numerical approximation (Numerical) of this double transformation seen in Fig. 6. The following parameters used: .. = 0.5, 1, 3 and 5.
Figure 6: Given the incomplete gamma function ..02(..)=1.(..)...0.....1.......... Compute its numerical approximation ....02(..)=.{..1{1.(..)...0.....1.........}} for values .. = 0.5, 1, 3 and 5.
7 Approximation of the Dirac delta function
The Dirac delta function [5] can be loosely thought of as a function on the real line which is zero everywhere except at the origin, where it is infinite,
..(..)=(+. (..=0)0 (...0) (7.1)
and which is also constrained to satisfy the identity
......(..)....=1 (7.2)
This is merely a heuristic characterization. The Dirac delta is not a function in the traditional sense as no function defined on the real numbers has these properties. This function can be rigorously defined either as a distribution or as a measure.
Noting that the Dirac delta function can be defined as the limit (in the sense of distributions) of the sequence of zero-centered normal distributions,
....(..)=1..v.......2/..2, (7.3)
as ..> 0.
By analytic continuation of the Fourier transform, the Laplace transform of the delta function is found to be [5],
..0.........(.....)....=......., (7.4)
consistent with the definition of the Laplace transform of ..(.....) as ........
Example 3. Find ...03(..)=..1{..03(..)} and ...03(..)=.{..1{..03(..)}}, where
..03(..)=exp(.......) (..>0and...(0,1]) (7.5)
The expression of the inverse Laplace transform, in terms of standard mathematical functions, is unknown. We can handle the inverse Laplace transform and double transformation technique, including Dirac delta function and its shifted form.
So, if ..=1, then ..(..)=..(.....), where ..(..) is the Dirac delta function [5].
Mathematica gives for ..=0.5 and ..=0.5
..1{exp(.......)}=0.14104739588693907exp(.0.0625/..)..3/2. (7.6)
Numerical inversion of Laplace transform ..03(..)=....... is known to be equivalent to the approximation of the Dirac delta function. Fig. 7 illustrates the approximation of the Dirac delta function with parameter .. = 1 evaluated in double precision. On the left plot used equal number of terms and precision level ..=... This approximation take negative values while the delta function is strictly positive. On the right plot .. = 16 and .. = 64. When looking at numerical inversion, it is important to note the accuracy with a varying number of expasion terms and precision level. We compare the inverses using Gaver-Stehfest implementation and observe the accuracy of the inversions as we increase the number of the expansion terms and precision level. However, there exists a limit to adding additional terms [12]. As we increase the number of expansion terms using .. = 64 we quickly discover that the numerical inversion becomes unstable and our function is dominated by numerical error (right plot).
Figure 7: Approximation of the Dirac delta function, .. = 1, in double precision. Precision level .. and the number of expansion terms .. are (16, 16) (left plot) and (16, 64) (right plot)
Extended precisions (Fig. 8) allow to combat the numerical limitation that we experience when dealing with double precision. Thus, we can use a larger number of terms. These examples use equal number of terms and precision level ..=... To improve the accuracy of the approximation we extended precision level to 32, 64, 128, 512 and 1000.
Figure 8: Approximation of the Dirac delta function, .. = 1. Precision level .. and number of terms .. are equal, ..=..: 32, 64, 128, 256, 512 and 1000. Note the difference in scales seen in four plots
In Fig. 9 are shown Inverse Laplace transform (left plot) and Numerical Laplace double transformation (right plot) of ..03(..)=exp(.......) evaluated for .. = 0.5 and .. = 0.5. The analytical solution of the inverse Laplace transform are given by (7.6). The approximation are given in double precision, and the errors are: ..=6.6.10.4 and 1.3.10.2 corresponding to the left and to the right plot. In Fig. 10 presented the screenshot for Laplace transform inversion of ..03(..)=exp(.......) evaluated for the same .. and .. with precision level 64. Note the accuracy of approximation improved at the order of at least 10.18.
Figure 9: Given the original Laplace transform ..03(..)=exp(.......). Compute its inverse Laplace transform (left plot) and numerical approximation ....03(..)=.{..1{exp(.......)}} (right plot)
Figure 10: The screenshot for inverse Laplace transform of the function ..03=exp(.......), evaluated for .. = 0.5 and .. = 0.5 with precision level 64
Tab. 2 displays numerical error ..03(..)=|..03(..).....03(..)| for the numerical approximation ....03(..)=.{..1{exp(.......)}, evaluated for .. = 0.5 and .. = 0.5, in varying precision levels of 16 and 64.
Table 2: Numerical error ..03(..)=|..03(..).....03(..)| for the numerical approximation ....03(..)=.{..1{exp(.......)}, evaluated for .. = 0.5 and .. = 0.5, in varying precision levels of 16 and 64. Values 9.23e-6 .9.23.10.6
8 Waiting times distribution in ../../.. queue
The ../../1 model assumes Poisson arrival at rate .. to a single server with generally distributed service time ... The coefficient of variation .... of .. defined by ....=../.., where ..=..[..] and .. are the expectation and the standard deviation of ...
If ....=1, we have M/M/1 model with ...(..)=../(..+..).
The PDF and CDF of waiting time simplified to [17]
....(..)=..(1...).....(1...)..,...0, ......, ../../1 (8.1)
....(..)=1........(1...)..,...0, ......, ../../1 (8.2)
If the coefficient ....>1 an approximation can be effective provided by hyperexponential distribution using for it the definition by parallel stages.
Let the service time .. follows a hyperexponential ..2 distribution whose PDF is given by (9.2), and the CDF is defined by (9.3).
The Laplace transform of the PDF is given by (9.1). This gives the Laplace transform of the CDF
..(..)=..(..)/..=1...2..=1............+.. (8.3)
First and second derivatives of ..(..) are:
....(..)....=..2..=1........(....+..)2 (8.4)
..2..(..)....2=2.2..=1........(....+..)3 (8.5)
The expectation ..[..]=.....(..)....|..=0 and the variance ......[..]=..[..2].(..[..])2 of the random variable .. are the following:
..[..]=..1..1+..2..2 (8.6)
......[..]=..1(2...1)..12+..2(2...2)..22.2..1..2..1..2 (8.7)
..[..2]=..2..(..)....2|..=0=2..1..12+2..2..22 (8.8)
To satisfy the condition (8.6), let
..1=2..1..[..], ..2=2..2..[..] (8.9)
Substituting (8.6), (8.8) and (8.9) in ....2=(..[..])2(..[..])2=..[..2].(..[..])2(..[..])2, this gives the unique ..2 PDF [6]:
..1=12(1+v....2.1....2+1),..2=1...1,..1=2..1..[..],..2=2..2..[..] (8.10)
9 Solving performance measures in ../..../.. queue
Example 4. Two cases consider for this ../../1 queue example.
Case 1. Laplace transform are given for PDF of the service time distribution.
Find ...04(..)=..1{..04(..)} and ...04(..)=.{..1{..04(..)}}, where
..04(..)=.2..=1............+.., (0...1.1,0...2.1,..1+..2=1), (9.1)
..04(..)=.2..=1................. (..>0). (9.2)
Case 2. Identical to case 1, but now for CDF.
Find ...04(..)=..1{..04(..)} and ...04(..)=.{..1{..04(..)}}, where
..04(..)=..04(..)/..=1...2..=1............+.. (9.3)
..04(..)=1..2..=1............., (..>0). (9.4)
Several cases were considered for ../../1 queue example. First is for Laplace transform of the PDF for service time distribution. Second is identical but for CDF. The ../../1 queue describes with ..=0.8, the expectation ..[..]=1, and the coefficient of variation .... = 1.5, 2.5 and 4.5. In Fig. 11 are shown inverse Laplace transform of the function ..04(..), evaluated for the PDF (left plot), and inverse Laplace transform ..04(..)/.. evaluated for the CDF (right plot). The errors are: ..=3.52.10.5 (left plot) and ..=1.2.10.5 (right plot).
Figure 11: Inverse Laplace transform of the functions ..04(..), evaluated for the PDF (left plot), and for ..04(..)/.. evaluated for the CDF (right plot). The ../../1 queue has ..=0.8, the expectation ..[..]=1.0, and the coefficients of variation .... = 1.5, 2.5 and 4.5.
10 Waiting time distribution in M/D/1 queues
Let service times has .... density with mean 1/.. and PDF
..(..)=(....).......1.........(...1)! (0<..<.) (10.1)
and LST
...(..)=(......+....).. (10.2)
The model M/D/1 can be considered as a limited case of ../..../1. As ..>. and ..>. in such a way that .....1>.. (0<..<.) then .... service times are deterministic with the constant ... The traffic intensity ....<1. Now ...>....... as ..>.. PDF of waiting time has the Laplace transform [17].
....(..)=(1...).......[1........] (10.3)
Example 5. Numerically estimate waiting time distribution ....(..) for different service models in M/G/1 queue. The Laplace transform of ....(..) is given by Pollaczek-Khinchine (P-K) transform equation [17]
....(..)=....(..)..=(1...).....[1....(..)],where (10.4)
...(..)=..0....(..)=..0.........(..)...., (10.5)
...(..) is the Laplace-Stieltjes transform (LST) of ..(..), ..(..) is the CDF of the service time distribution, .. and .. are the averages of arrival rate and service tine distribution, and ..=.... is the traffic intensity.
As with ../../1 queue the following service models considered:
../../1 ...(..)=....+.. (10.6)
../..../1 ...(..)=(....+..).. (10.7)
../../1 ...(..)=....... (10.8)
../..2/1 ...(..)=.2..=1............+.. (10.9)
Consider ../../1 model for calculating waiting time distributions. For the system ../..2/1 the arrival rate ..=5.0. The service time distribution is ..2 evaluated for ..=6, ..=1/.. and ....=1.5. For this ../..2/1 system the traffic intensity ..=../.., and the CDF at the time 0 is ..(0)=1....
Gaver-Sthefest algorithm used for inverting Laplace transform for variety ...(..) in (10.4). The CDF distribution of the waiting time in queue distribution CDF for ../../1,../..2/1 and ../../1 are shown in Fig. 12.
Figure 12: PDF and CDF waiting time distributions for ../../1, ../../1 and ../..2/1
For ../../1 queue service times are deterministic and equal to the value ... For deterministic service ...(..) are given by (10.8). We compare ....(..) calculated by inverting of (10.4) with ....(..) also known and analytically given by [17]:
....(..)=(1...).[../..]..=0.....(.......)(.......)....!...., (10.10)
where [..] is the greatest integer less than or equal to ...
Figure 13 shows results for this queue with .. = 5.0, .. = 6.0 (..=1/.. and ..=../.. = 0.8). The figure shows errors, by analytical estimate, in logarithmic scale Log....(..) (left plot). The estimate for ....(..) is then obtained by Gaver-Stehfest algorithm. We failed to get the numerical analytical solution for precision level: .. = 16 if ..>10; .. = 64 if ..>26; .. = 256 if ..>90. Only for ..=512 and ..=1000 we get the correct result for all range 0<...100. Gaver-Stehfest inversion ....(..) are shown in the right plot, even for double precision we have accurate solution, and can not distinguish by eyes the curves with different precision .. = 16, 64, 256 and 512.
Figure 13: ../../1 outputs: Log of ....(..) error (left plot) and ....(..)
by Gaver-Stehfest inversion (right plot)
Fig. 14 shows numerical results in double precision (N = 16) for ../../1 waiting time distribution ....(..) by analytical solution (10.10) and Gaver-Stehfest inversion. The analytical solution is dominated by noise after ..>9.
Figure 14: ....(..) for ../../1,CDF waiting time distribution
by analytical solution and Gaver-Stehfest inversion
The most common service time distribution is the exponential, in which case the waiting time distribution are available in closed-form. For more general cases, closed-form solutions of (P-K) ../../1 transform equation are mathematically intractable. The following specific example used to compare analytical solution and Gaver-Stehfest inversion. Consider the system ../..2/1 with pdf for service time distribution [7]
..(..)=14.........+34(2..)...2...., (10.11)
where .. is the arrival rate, ..=5/(8..), and .....=....=5/8. For numerical solutions we used ..=5. The corresponding laplace transform ...(..) is
...(..)=(14)....+..+(34)2..2..+.. (10.12)
Using ...(..) and (P-K) transform equation (10.4) ....(..) and ....(..) for the waiting time density [7]:
.....(..)=(1...)[1+../4(3/2)..+..+3../4(1/2)..+..] (10.13)
....(..)=38..0(..)+3..32...(3/2)....+9..32...(1/2).... ...0, (10.14)
where ..0(..) is unit impulse function.
The analytical solution for CDF of waiting time distribution can be easily found as:
.....(..)=(1...)[1..+../4((3/2)..+..)..+3../4((1/2)..+..)..] (10.15)
....(..)=(1...)[83.16...(3/2).....32...(1/2)....] ...0 (10.16)
Fig. 15 shows numerical results in double precision (N = 16) for ../..2/1 waiting time distribution for PDF (on the left) and CDF (on the right). To recognize by eyes the distinguish an analytical and Gaver-Stehfest inversion outputs is almost impossible.
Figure 15: Waiting time distribution for ../..2/1 by analytical solution and Gaver-Stehfest inversion for PDF (left plot) and CDF (right plot)
11 Accuracy and precision requirements of the ../../.. large .. analysis
The results of double transformation approximation for waiting time distribution ....(..) in ../../1 are shown in Figs. 16 and 17.
The following double transformation technique used for the Laplace transform of the inversion:
......(..)=.{..1{....(..)}} (11.1)
The exact solution is original ....(..) which compares with ......(..) after double transformation technique. Laplace transform inversion implemented by the Gaver-Stehfest algorithm, and a composite Simpson�s rule is performed for the numerical direct Laplace transform. Approximation for the waiting time distribution in ../../1 queue is convenient for light traffic intensity and small .., providing suitable approximation in double precision. The Figures fit the wide-range small Laplace transform parameter .., corresponding to large number ...
Figure 16: Double transformation approximation for waiting time distribution ....(..) of the ../../1: double precision (.. = 16) with different number of subinterval .. = 500, 5000, 50000, 150000 and ranges [0.1, 10] (left plot), and [0.1, 5] (right plot)
Figure 17: Double transformation approximation for waiting time distribution ....(..) of the ../../1 with multiple precision calculation with the number of subintervals .. = 500
On the Fig. 16 all curves are given for double precision level (.. = 16) with different number of subinterval .. = 500, 5000, 50000 and 150000 in numerical direct Laplace calculation. It looks like the is no effect of .. increased for .. on the interval [0.1, 10] (left plot), but curve are different for .. on the smaller interval [0.1, 0.5] (right plot).
The Fig.17 demonstrate the impact of precision level ... The number of subintervals ..=500. In double precision the method does not work well, and significant improvements in accuracy illustrated as precision level increased up to 256.
12 Conclusions
Accuracy and stability of numerical Laplace transform and inversion are crucial for many applications in computation probability models. In this work, we proposed and evaluated different numerical implementations of the Laplace transform and inversion in multiple precision arithmetic systems. We are concerned with transformation that can occur in two ways.
If the examples cover functions with known inverses, accuracy of multiple precision models can be asserted by comparison to the exact solution.
The most realistic and challenging problems cover functions with analytically unknown inverses. So double transformation approach proposed to find computationally efficient methods for performing the numerical Laplace transform and inversion. In this approach numerical Laplace transform inversion used directly as input into numerical Laplace transform. The accuracy can be asserted by comparison to the original Laplace planform. We observe the accuracy of the inversions as we increase the number of the expansion terms and precision, which lead to stable solutions. The computational efficiency compared to precision levels is demonstrated for waiting times distribution in ../../1 queuing systems.
References
[1] J. Abate, P. Valko, Multi-precision Laplace transform inversion, International Journal for Numerical Methods in Engineering 60 (2004) 979-993.
[2] K. E. Atkinson, An Introduction to Numerical Analysis, second ed., John Wiley & Sons, 1989.
[3] D. H. Bailey, Y. Hida, X. S. Li, B. Thompson, ARPREC: An Arbitrary Precision Computation Package, Tech. Rep. LBNL-53651, Lawrence Berkley National Lab., 2002.
[4] A. M. Cohen, Numerical Methods for Laplace Transform Inversion, Springer, 2007.
[5] C. Edwards, C. Penney, Differential Equations and Boundary Value Problems: Computing and Modeling, 5th ed., 2015.
[6] E. Kao, An Introduction to Stochastic Processes, Duxbury Press, 1997.
[7] L. Kleinrock, Queueing Systems, Volume I, Hardcover, 1975
[8] E. Kreyszig, Advanced Engineering Mathematics, 10th ed., 2011.
[9] Z. Krougly, M. Davison, S. Aiyar, The role of high precision arithmetic in calculating numerical Laplace and inverse Laplace transforms, Applied Mathematics 8 (2017) 562-589.
[10] Z. L. Krougly, D. J. Jeffrey, Implementation and application of extended precision in Matlab, in: N. Mastorakis et al (Eds.), Proc. of the Applied Computing Conference ACC�09, WSEAS Press (2009) 103-108.
[11] Z. L. Krougly, D. A. Stanford, Iterative algorithms for performance evaluation of closed network models, Performance Evaluation 61 (2005), 41-64.
[12] Z. L. Krougly, D. J. Jeffrey, D. Tsarapkina, Software implementation of numerical algorithms in arbitrary precision, in: 15th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC 2013), N. Bjorner et al. (Eds.), IEEE Computer Society (2014) 132-138.
[13] K. L. Kuhlman, Review of inverse Laplace transform algorithms for Laplace-space numerical approaches, Numerical Algorithms 63 (2013) 339-355.
[14] A. Kuznetsov, On the convergence of the Gaver�Stehfest algorithm, SIAM J. Numer. Anal. 51 (2013) 2984-2998.
[15] A. Murli, M. Rizzardi, Algorithm 682 Talbot�s method for the Laplace inversion problem, ACM Transactions on Mathematical Software 16 (1990) 158-168.
[16] S. Nadarajah, S. Kotz, On the Laplace transform of the Pareto distribution, Queueing System 54 (2006) 243-244.
[17] J.F. Shortle, J.F. Thompson, D. Gross, C. M. Harris. Fundamentals of Queueing Theory, 5th ed., Wiley, New York, 2018.
[18] H. Stehfest, Algorithm 368: Numerical Inversion of Laplace Transform, Communications of the ACM, 13(1), 1970: 47-49.
