Abstract

We develop an algorithm for numerically inverting multi-dimensional transforms. Our algorithm applies to any number of continuous variables (Laplace transforms) and discrete variables (generating functions). We use the Fourier-series method; i.e., the inversion formula is the Fourier series of a periodic function constructed by aliasing. This amounts to an application of the Poisson summation formula. By appropriately exponentially damping the given function, we control the aliasing error. We choose the periods of the multi-dimensional periodic function so that each infinite series is a finite sum of nearly alternating infinite series; then we apply the Euler transformation to compute the infinite series from finitely many terms. The multidimensional inversion algorithm enables us, evidently for the first time, to quickly and accurately calculate probability distributions from several classical transforms in queueing theory. For example, we apply our algorithm to invert the two-dimensional transforms of the joint distribution of the duration of a busy period and the number served in that busy period, and the time-dependent of the transient queue-length and workload distributions, in the M/G/1 queue. In other related work, we have applied the inversion algorithms here to calculate time-dependent distributions in the transient BMAP/G/1 queue (with a batch Markovian arrival process) and the piecewise-stationary M t / G t /1 queue.

1. Introduction

In this paper we present an algorithm for numerically inverting multi-dimensional transforms. We are motivated by the desire to compute probability distributions of interest in queues and related stochastic models, but of course there are many other applications. We even allow the inverse transform to be complex-valued. However, in our error analysis we exploit the fact that the modulus of our functions have known bounds, so that the algorithm is particularly appropriate for probability transforms (where the bound is 1). This algorithm is evidently the first multidimensional inversion algorithm in the queueing literature. However, we have learned that a different multi-dimensional inversion algorithm intended for queueing models has recently been developed in Russia by Frolov and Kitaev (1992). Their algorithm evidently is similar to a multi-dimensional version of the POST-WIDDER algorithm in Abate and Whitt (1992a). Of course, there is substantial literature on numerical transform inversion, as reviewed in Abate and Whitt (1992a). However, relatively little attention has been given to inversion of multidimensional transforms; for some instances, see Singhal, Vlach and Vlach (1975), Huntley and Zinober (1979) and Shephard (1991).

We consider both continuous variables (Laplace transforms) and discrete variables (generating functions). We thus consider three types of two-dimensional transforms: (i) continuouscontinuous, (ii) continuous-discrete and (iii) discrete-discrete. We also show how the formulas can be generalized to more than two dimensions with any number of continuous and discrete variables.

The multi-dimensional inversions obviously allows us to compute multivariate probability distributions, as we illustrate here. However, the multi-dimensional inversions also allow us to calculate time-dependent probability distributions in queueing models that are not in steady state. As an example, in this paper we invert the classical double transform expressions for the transient workload and queue-length distributions in the M/G/1 queue; see Takacs (1962). The special case of the M/M/1 transient queue length has been widely studied in the literature; e.g., see Abate and Whitt (1989). We show that our algorithm in this case is comparable in speed and accuracy to the numerical integration of the integral representation in Abate and Whitt (1989). Moreover, our algorithm applies equally well to the case of non-exponential service times, with no loss of speed or accuracy. In fact, we have applied the algorithm here to calculate time-dependent distributions in the transient BMAP/G/1 queue (with a batch Markovian arrival process) in Lucantoni, Choudhury and Whitt (1994) and the piecewise-stationary M t / G t / 1 queue in Choudhury, Lucantoni and Whitt (1993). We plan to report on applications of the multidimensional inversion algorithm to other important queueing problems in the future.

Our algorithms here is a multivariate generalization of the algorithms EULER and LATTICE-POISSON in Abate and Whitt (1992a). We also introduce an enhancement of those algorithms to be able to simultaneously control the aliasing and roundoff errors. As in Abate and Whitt (1992a), we exploit the Fourier-series method. The general approach goes back at least to Fettis (1955); see Abate and Whitt (1992a) for a review. For the multi-dimensional transforms, this means that we apply the multivariate version of the Poisson summation formula, as given for the two-dimensional continuous-continuous case in (5.47) of Abate and Whitt (1992a); also see Good (1962). The approach is closely related to the fast Fourier transform (FFT). The idea is relatively simple: Just as in the one-dimensional case, in the two-dimensional continuouscontinuous case we damp the given function by multiplying by a two-dimensional decaying exponential function and then approximate the damped function by a periodic function constructed by aliasing. We use the two exponential parameters to control the aliasing error in this approximation by the periodic function. The inversion formula is then the two-dimensional Fourier series of the periodic function. This yields what we want, because the transform values are the two-dimensional Fourier coefficients. Moreover, the two periods of the periodic function can be chosen so that the two-dimensional Fourier series is a series nested within a second series, each of them being nearly an alternating series. Hence, we can efficiently calculate each infinite series from finitely many terms by exploiting the Euler transformation (or summation). In practice, this usually means that it suffices to compute 100 or fewer terms of each infinite series to achieve a truncation error of the order 10^-13 or less; see Abate and Whitt (1992a), Johnsonbaugh (1979) and Wimp (1981). When the inverse transform is real, as with probabilities, the overall computation can be reduced by a factor of two.

. . .Continue to read rest of article (PDF).