1 Case One ()
Set , and
A sequence of functions is defined iteratively as follows:
where the empty sum convention holds for and it is understood that for . Since is a degree polynomial in with coefficients of the form , where , are rational numbers, the integrals , , all possess closedform expressions. Therefore the differential equation for can be solved exactly. Stitching the fragments together gives rise to the density function
pictured in Figure 1, where denotes the greatest integer .
Let us illustrate in greater detail:
thus
implying
Continuing:
thus
implying
The pattern is maintained:
until (the equation length becomes fixed and replaces the rightmost ):
We count terms in each upon expansion, at least for .
When numerically evaluating large symbolic expressions, it is important to ensure that the working precision of floating point quantities is suitably high. To evaluate may require several hundred decimal digits because, for instance, the first two of the 348 numerator terms comprising are
upon setting . The subtraction of nearly equal numbers, such as these, will lead to a horrific loss of precision unless appropriate care is taken.
2 Case Two ()
Set , and
A sequence of functions is defined iteratively as follows:
where is the Kronecker delta. Stitching the fragments together gives rise to
pictured in Figure 2.
Let us illustrate in greater detail:
thus
implying
Beyond this point, the derivatives at match:
thus
implying
3 Case Three ()
Set and
A sequence of functions is defined iteratively as follows:
The role of
here is more pronounced than in Section 2: a jump discontinuity occurs in the density at
(as opposed to merely a sharp corner) . Stitching the fragments together gives rise topictured in Figure 3.
4 Personal Notes
I taught a semesterlong class in statistical programming at Harvard University (as a preceptor) for nearly ten years. A favorite set of problems began with an M/M/1 example [10] involving a hospital emergency room (ER).
Patients (clients) arrive according to a Poisson process with rate . One doctor (server) is available to treat them. The doctor, when busy, treats patients with rate
. More precisely: interarrival times are exponentially distributed with mean
and treatment lengths are exponentially distributed with mean. The ER is open 24 hours a day, 7 days a week. Patients must wait until the doctor is free and are treated in the order via which they arrive. Simulate the performance of the ER over many weeks. What can be said about waiting time in queue per patient (excluding service time)? Determine the mean, variance, mode and median of
to as high accuracy as possible. Assume for this purpose that the expected number of arriving patients per hour is and that the expected number of treatment completions per hour (for a continuously busy ER) is .
Putting aside experiment in favor of pristine theory, the Laplace transform of service density is . Formulas for the mean and variance of follow immediately from
consequently
giving and respectively. An expression for the density of :
shows that the mode is ; integrating and solving the equation
implies
giving the median to be .
The aforementioned problem set continued with an M/G/1 example, involving the same ER parameters, but with Uniform[] treatment lengths. Since was required, I arbitrarily chose and .
Perform exactly the same simulation as before, except assume treatment lengths (in minutes) are uniformly distributed on the interval
. Less is known about this scenario than the preceding (with exponential service times).
The thought of choosing or
did not occur to me until later. I had imagined that numerical inversion of Laplace transforms was the only avenue available to reliably estimate the mode and median.
Just as
is the first moment of treatment lengths,
are the corresponding second and third moments. Again, favoring pristine theory over experiment,
giving and respectively. The mode (location of the density maximum, excluding ) occurs when
i.e.,
solving the equation
gives the median to be .
My classroom example deviates from the direction of research [13], which emphasizes heavytailed service time distributions. Ramsay [14] discovered a remarkably compact formula – a single integral of a nonoscillating function over the real line – in connection with Pareto service times. My formulas for uniform service times are sprawling by comparison.
The recursive solution of delaydifferential equations is certainly not new [15] but application of such to queueing theory does not seem to be widespread. Counterexamples include [16, 17]; surely there are more that I’ve missed. Symbolics are mentioned in [13].
The waiting time probability for M/D/1 seems to decay exponentially (as do the other two cases). The CramérLundberg approximation [18] is applicable; alternatively, we have asymptotics [11, 19]
where is the unique root of .
Returning finally to M/U/1, let denote the number of patients in the system (both queue and service). Define
Under equilibrium, with , and , , we have
variance  
giving and respectively. The final term for the variance is missing in [21]. The probability generating function seems to decay geometrically, but details surrounding the exact limit of successive ratios have not been verified [20].
I am grateful to innumerable software developers, as my “effective” formulas are too lengthy to be studied in any traditional sense. Mathematica routines NDSolve for DDEs and InverseLaplaceTransform (for Mma version ) plus ILTCME [22] assisted in numerically confirming many results. R steadfastly remains my favorite statistical programming language. A student asked in 2006 for my help in writing a relevant R simulation, leading to the computational exercises described here and to my abiding interest in queues.
References
 [1] J. F. Shortle, J. M. Thompson, D. Gross and C. M. Harris, Fundamentals of Queueing Theory, 5 ed., Wiley, 2018, pp. 261–268, 271–275, 443–444; MR3791493.
 [2] M. F. Neuts, Generalizations of the PollaczekKhinchin integral equation in the theory of queues, Adv. in Appl. Probab. 18 (1986) 952–990; MR0867095.
 [3] K. H. Lundberg, H. R. Miller and D. L. Trumper, Initial conditions, generalized functions, and the Laplace transform troubles at the origin, IEEE Control Systems Magazine 27 (2007) 22–35.
 [4] J. L. Schiff, The Laplace Transform: Theory and Applications, SpringerVerlag, 1999, pp. 88–89; MR1716143.
 [5] A. Henderson, Extension of the initial value theorem, Internat. J. Electrical Engineering & Education 15 (1978) 237–241.
 [6] A. K. Erlang, The theory of probabilities and telephone conversations (in Danish), Nyt Tidsskrift for Matematik 20B (1909) 33–39; Engl. transl. in E. Brockmeyer, H. L. Halstrøm and A. Jensen, eds., The Life and Works of A. K. Erlang, Trans. Danish Acad. Tech. Sci., 1948, pp. 131137; MR0027976.
 [7] V. B. Iversen and L. Staalhagen, Waiting time distribution in M/D/1 queueing systems, Electronics Letters 35 (1999) 2184–2185.
 [8] N. U. Prabhu, Stochastic Storage Processes. Queues, Insurance Risk, Dams, and Data Communication, 2 ed., SpringerVerlag, 1998, pp. 36–37, 75–78; MR1492990.
 [9] H. Tijms, New and old results for the M/D/c queue, AEU  Internat. J. Electronics and Communications 60 (2006) 125–130.
 [10] F. S. Hillier and G. J. Lieberman, Introduction to Operations Research, 3 ed., HoldenDay, 1980, pp. 400–457 (chapter entitled “Queueing Theory”); MR0569591.
 [11] H. C. Tijms, A First Course in Stochastic Models, Wiley, 2003, pp. 327–331, 348–353, 381–384; MR2190630.
 [12] X. Zhao, J. Hou and K. Gilbert, Measuring the variance of customer waiting time in service operations, Management Decision 52 (2014) 296–312.
 [13] J. F. Shortle, P. H. Brill, M. J. Fischer, D. Gross and D. M. B. Masi, An algorithm to compute the waiting time distribution for the M/G/1 queue, INFORMS J. Comput. 16 (2004) 152–161; MR2065995.
 [14] C. M. Ramsay, Exact waiting time and queue size distributions for equilibrium M/G/1 queues with Pareto service, Queueing Systems 57 (2007) 147–155; MR2369214.
 [15] S. R. Finch, Components and cycles of random mappings, arXiv:2205.05579.
 [16] M. Vlasiou, LindleyType Recursions, Ph.D. thesis, Eindhoven University of Technology, 2006, http://research.tue.nl/en/publications/lindleytyperecursions.
 [17] O. J. Boxma and M. Vlasiou, On queues with service and interarrival times depending on waiting times, Queueing Systems 56 (2007) 121–132; arXiv:1404.5549; MR2336100.
 [18] T. Sakurai, Approximating M/G/1 waiting time tail probabilities, Stochastic Models 20 (2004) 173–191; MR2048251.
 [19] H. Tijms and K. Staats, Negative probabilities at work in the M/D/1 queue, Probab. Engrg. Inform. Sci. 21 (2007) 67–76; MR2289220.
 [20] M. L. Chaudhry, U. C. Gupta and M. Agarwal, Exact and approximate numerical solutions to steadystate singleserver queues: M/G/1 – a unified approach, Queueing Systems 10 (1992) 351–379; MR1150007.
 [21] J. D. Griffiths, The coefficient of variation of queue size for heavy traffic, J. Operational Research Society 47 (1996) 10711076.

[22]
G. Horváth, I. Horváth, M. Telek, S. AlDeen Almousa
and Z. Talyigás, Inverse Laplace Transform with Concentrated
MatrixExponential Functions, http://inverselaplace.org.
Steven Finch MIT Sloan School of Management Cambridge, MA, USA steven_finch@harvard.edu