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.