We introduce and investigate a framework for constructing algorithms to invert Laplace transforms numerically. Given a Laplace transform \hat{f} of a complex-valued function of a nonnegative real-variable, f, the function f is approximated by a finite linear combination of the transform values; i.e., we use the inversion formula f(t) \approx f_n (t) \equiv \frac{1}{t} \sum_{k = 0}^{n}\omega_{k}\hat{f}\biggl(\frac{\alpha_{k}}{t}\biggr),\quad 0 < t < \infty, where the weights ω k and nodes α k are complex numbers, which depend on n, but do not depend on the transform \hat{f} or the time argument t. Many different algorithms can be put into this framework, because it remains to specify the weights and nodes. We examine three one-dimensional inversion routines in this framework: the Gaver-Stehfest algorithm, a version of the Fourier-series method with Euler summation, and a version of the Talbot algorithm, which is based on deforming the contour in the Bromwich inversion integral. We show that these three building blocks can be combined to produce different algorithms for numerically inverting two-dimensional Laplace transforms, again all depending on the single parameter n. We show that it can be advantageous to use different one-dimensional algorithms in the inner and outer loops.

It is often possible to effectively calculate probability density functions (pdf's) and cumulative distribution functions (cdf's) by numerically inverting Laplace transforms. However, to do so it is necessary to compute the Laplace transform values. Unfortunately, convenient explicit expressions for required transforms are often unavailable for component pdf's in a probability model. In that event, we show that it is sometimes possible to find continued-fraction representations for required Laplace transforms that can serve as a basis for computing the transform values needed in the inversion algorithm. This property is very likely to prevail for completely monotone pdf's, because their Laplace transforms have special continued fractions called S fractions, which have desirable convergence properties. We illustrate the approach by considering applications to compute first-passage-time cdf's in birth-and-death processes and various cdf's with non-exponential tails, which can be used to model service-time cdf's in queueing models. Included among these cdf's is the Pareto cdf.

The Laguerre method for numerically inverting Laplace transforms is an old established method based on the 1935 Tricomi–Widder theorem, which shows (under suitable regularity conditions) that the desired function can be represented as a weighted sum of Laguerre functions, where the weights are coefficients of a generating function constructed from the Laplace transform using a bilinear transformation. We present a new variant of the Laguerre method based on: (i) using our previously developed variant of the Fourier-series method to calculate the coefficients of the Laguerre generating function, (ii) developing systematic methods for scaling, and (iii) using Wynn's ϵ-algorithm to accelerate convergence of the Laguerre series when the Laguerre coefficients do not converge to zero geometrically fast. These contributions significantly expand the class of transforms that can be effectively inverted by the Laguerre method. We provide insight into the slow convergence of the Laguerre coefficients as well as propose a remedy. Before acceleration, the rate of convergence can often be determined from the Laplace transform by applying Darboux's theorem. Even when the Laguerre coefficients converge to zero geometrically fast, it can be difficult to calculate the desired functions for large arguments because of roundoff errors. We solve this problem by calculating very small Laguerre coefficients with low relative error through appropriate scaling. We also develop another acceleration technique for the case in which the Laguerre coefficients converge to zero geometrically fast. We illustrate the effectiveness of our algorithm through numerical examples.

We present a simple algorithm for numerically inverting Laplace transforms. The algorithm is designed especially for probability cumulative distribution functions, but it applies to other functions as well. Since it does not seem possible to provide effective methods with simple general error bounds, we simultaneously use two different methods to confirm the accuracy. Both methods are variants of the Fourier-series method. The first, building on Dubner and Abate (Dubner, H., J. Abate. 1968. Numerical inversion of Laplace transforms by relating them to the finite Fourier cosine transform. JACM 15 115–123.) and Simon, Stroot, and Weiss (Simon, R. M., M. T. Stroot, G. H. Weiss. 1972. Numerical inversion of Laplace transforms with application to percentage labeled experiments. Comput. Biomed. Res. 6 596–607.), uses the Bromwich integral, the Poisson summation formula and Euler summation; the second, building on Jagerman (Jagerman, D. L. 1978. An inversion technique for the Laplace transform with applications. Bell System Tech. J. 57 669–710 and Jagerman, D. L. 1982. An inversion technique for the Laplace transform. Bell System Tech. J. 61 1995–2002.), uses the Post-Widder formula, the Poisson summation formula, and the Stehfest (Stehfest, H. 1970. Algorithm 368. Numerical inversion of Laplace transforms. Comm. ACM 13 479–490 (erratum: 13 624).) enhancement. The resulting program is short and the computational experience is encouraging. INFORMS Journal on Computing, ISSN 1091-9856, was published as ORSA Journal on Computing from 1989 to 1995 under ISSN 0899-1499.

Methods are discussed for computing transient performance measures for the M/M/1 queue. These performance measures are often expressed in terms of modified Bessel functions without any discussion about computation. In fact, a common expression for the probability transition function of the M/M/1 queue length process has an infinite sum of modified Bessel functions. For actually generating numbers, however, it is convenient to use numerical integration with associated integral representations, as was first pointed out by P.M. Morse (Oper. Res., vol.3, p.255-61, 1955).>

This paper analyzes the performance, for a particular strategy, of a head-per-track type auxiliary storage system in a real-time environment. The prototype system of either fixed-head disks or drums incorporates a hardware queuer. This device attempts to always select first that request for a data transfer which will incur the shortest rotational latency relative to all other possible waiting file access requests. The queueing analysis is concerned with the tradeoff which is experienced in practice between throughput of a stochastic service device and the response time for each service request. An exact analysis of the system is shown to be totally unmanageable. Therefore, the approach is to estimate the results via physical arguments. We obtain an approximate stationary distribution for the waiting time of a file access request in the system.

In this paper the use of the techniques of queueing theory in analyzing the performance of a mass storage device in a real-time environment is demonstrated; concern is with the tradeoff experienced in practice between throughput of a stochastic service device and the response time for each service request. For concreteness, the analysis is applied to the IBM 2314 disk storage facility. The results are presented in a series of graphs showing the file system response time versus the throughput for several distributions of record length and arm movement. The queueing model and the theoretical tools used are described in sufficient detail to permit the reader to apply the techniques to other systems. In particular, any disk whose seek time characteristic can be approximated by a piecewise linear continuous function may be analyzed by the methods presented.

In this paper the problem of readily determining the inverse Laplace transform numerically by a method which meets the efficiency requirements of automatic digital computation is discussed. Because the result inverse function is given as a Fourier cosine series, the procedure requires only about ten FORTRAN statements. Furthermore, it does not require the use of involved algorithms for the generation of any special functions, but uses only cosines and exponentials. The basis of the method hinges on the fact that in evaluating the inverse Laplace transform integral there exists a freedom in choosing the contour of integration. Given certain restrictions, the contour may be any vertical line in the right-half plane. Specifying a line, the integral to be evaluated is essentially a Fourier integral. However, the method is concerned with determining the proper line, so that when the integral (along this line) is approximated, the error is as small as desired by virtue of having chosen the particular contour.