1. INTRODUCTION
Birth–death processes with immigration and emigration are major
models for the study of the population process of biological and
ecological systems, and their transient analysis is important in the
understanding of the structural behavior of such systems (see, e.g.,
Keyfitz [14]). However, analytic
transient solutions for such systems have been derived only for some
special cases (see, e.g., Kendall [13], Karlin and McGregor [11], Bailey [1], Bartlett [2], Iosifescu and Tautu [7], and Ismail, Letessier, and Valent
[8], among others). In this
article, we are concerned with the transient behavior of the general
linear birth–death processes with both immigration and
emigration.
Transient analysis of birth and death processes has also attracted
much interest because of its connection with queuing models. Morse
[20] studied the
M/M/1 queue and obtained the time-dependent
probabilities for the number of customers at any time t (see
also Kleinrock [15]). Using the
spectral measure method, Karlin and McGregor [10] analyzed the transient solutions of
multiple-server queues. Saaty [26]
considered the M/M/s queue and
derived the Laplace transform for the distribution of the number of
customers in a system at any time t (see also Kelton and Law
[12]). Rothkopf and Oren
[25] studied some generalized
M/M/s queues with time-dependent
arrival or service rates, and they obtained approximation results for
transient solutions. The spectral function of the continuous time
transient solution for M/M/s queue
was studied in van Doorn [27,
Chap.6], and further papers relevant to transient solution for
M/M/s queues include Whitt
[29], Halfin and Whitt [6], and Pegden and Rosenshine [21]. All these queuing systems are modeled as
and analyzed using simple birth–death processes.
The spectral method has been widely used in the literature in
analyzing birth–death processes (see Karlin and McGregor
[11]). In this article, we employ
an alternative approach, the method of characteristics, to
study birth and death processes. As will be seen in this article, this
method not only provides a unified approach to the study of
birth–death processes with immigration and emigration, but it
also yields explicit solutions that have not been explicitly reported
in the literature. The method can be considered as the dual method and
it is complementary to the spectral method of Karlin and McGregor
[11], further refined in Ismail et
al. [8]. Karlin and McGregor's
method is essentially a separation of variables which converts a
time-dependent partial differential equation to a time-independent
(stationary) problem; that is, the time variable is removed first. In
contrast, our method converts a time-dependent partial differential
equation to a set of ordinary differential equations with respect to
time, in which the spatial variable x is a benign parameter.
Put simply, our method is to remove x and Karlin and
McGregor's is to remove t. Karlin and McGregor's
eigenvalues and eigenfunctions approach is powerful. Nonetheless our
method of characteristics gives a comprehensive integral representation
of solutions for all the linear birth and death processes (see Theorem
1 in Section 2) and yields many explicit solutions (Section 4).
Furthermore, the spectral method requires the computation of measures
that depend heavily on Laguerre and Meixner polynomials, whereas our
method is more elementary.
In this article, we present a transient analysis for the general
linear birth–death processes with both immigration and emigration
by the method of characteristics. Our analysis involves solving partial
differential equations for the moment generating functions of the
transient distribution. For the general case, we present a
Volterra-type integral equation for the chance of extinction
P0(t) (i.e., the probability that the
system is empty at time t) and an explicit formula for the
z-transform of the transient distribution in terms of
P0(t). Furthermore, an explicit-form
infinite-series solution can be found for
P0(t). These results allow us to obtain
explicit solutions for the transient solutions of several special cases
that have not been previously presented in the literature.
This article is organized as follows. In the following section, we
present the transient analysis for the general linear birth–death
processes with immigration and emigration. In Section 3, we discuss the
chance of extinction P0(t) and its
asymptotic behavior. In Section 4, we present the closed-form solutions
for several special cases. We conclude the article with a discussion in
Section 5.
2. MODEL AND ANALYSIS
Consider a linear birth–death process with immigration and
emigration. The state of the system is the population size. When the
state of the system is n,n = 0,1,…, the
immigration rate is ν, the emigration rate is θ, the birth
rate is nλ, and the death rate is nμ.
Let X(t) be the population size at time t.
Then, {X(t); t ≥ 0} is a continuous-time
Markov process with transition rates
We are concerned with the probability distribution of the population
size at any time t ≥ 0 [i.e., of
X(t)] . Assume that the initial population
distribution is
and let
Note that if the system is initially empty, then
h(z) ≡ 1.
Let
then the Kolmogorov forward differential equation for
Pn(t) is (see, e.g., Ross
[24] and Bailey [1])
The initial condition for these differential equations is
Pn(0) = hn,
n = 0,1,….
Let
be the z-transform of
Pn(t),n ≥ 0, then it
follows from (2) that P(t,z) satisfies the
partial differential equation (PDE)
with initial-boundary conditions
Thus, the problem boils down to solving the PDE (3) with initial
condition (4). We remark that the appearance of the boundary term
P(t,0) in the equation makes the problem unusual from
the view point of standard PDEs. Thus, a key step is to determine
P(t,0) ≡ P0(t),
which is the probability that the population size at time t is
zero.
To solve the differential equation, we first study its
characteristics. The characteristic equation is (see, e.g., John
[9])
We note that this characteristic equation is the same as that studied
in Chao and Zheng [3]. Solving this
simple ordinary differential equation gives the characteristics of PDE
(3):
See Chao and Zheng [3, Figs. 1–5]
for the different scenarios of the characteristic curve. Some characteristic
curves intersect the t-axis if and only if μ > 0. We will
assume μ > 0 in this article.
Along any characteristic curve, the PDE becomes an ordinary
differential equation (ODE):
For convenience, we let
We now introduce the integration factor
where τ > 0 and z0 ∈ [0,1]
are two independent variables and the function z(ζ)
depends on z0. Then, the solution of the ODE can be
expressed implicitly as
where the parameter z0 is an arbitrary number in
[0,1] and z(0) = z0.
To obtain P(t,z) at a given point
(t,z), we let z(t) = z in
(5) and solve for z0 to obtain
for λ ≠ μ and
for λ = μ. Thus, the z-transform
P(t,z) is given by
where z0 is given by (8) for λ ≠ μ and
by (9) for λ = μ, and z(τ) is given by (5) (where
t is replaced by τ), with z0 as just
defined. The remaining key issue is to obtain
a(t).
If we let z0 > 0 be less than both
μ/λ and 1, then the characteristic curve starting at
z0 must intersect the t axis at a finite
time. See Figure 1. Denote this finite time
by t*. By setting z(t) = 0 in the
characteristic solution (5), we find that
as t → t*. Thus, the integration factor (6)
satisfies
as t → t* (assuming θ > 0). Since
P(t,z(t)) is expected to approach
P(t*,0) as t → t*, we conclude
that the second factor in (10) must approach zero for
z0 ∈ [0,min{1,μ/λ}); that
is,
This is the equation that determines a(t), which is
a key observation in our derivation.
Characteristics connecting the two axes.
The following is the main result of this section. The reader is
referred to [5], [22], and [23] for the concept of Volterra integral
equations (IEs).
Theorem 1: Suppose μ > 0 and θ >
0. Then, the following hold:
(a) The solution P(t,z) of PDE
(3) is given by (10) with z(t) given by (5),
Λ(t; x) given by (6), z0 given
by (8) and (9), and P(t,0) ≡ a(t)
given by integral equation (11).
(b) The integral equation (11) can be simplified to a singular
Volterra integral equation of the first kind
for 0 ≤ y < 1/(λ −
μ)+, where (λ − μ)+ =
max{λ − μ,0},
for y ≥ 0 and H(y) = 0 for
y < 0,
(c) Once A(s) is obtained from Volterra
integral equation (12), the solution P(t,z) of
PDE (3) given in (10) can be simplified to
where s = s(t) is given in
(14).
Proof:
Our earlier analysis has shown that any continuous solution
P(t,z) with a finite P(t,0)
yields a solution of the integral equation (11). Conversely, a finite solution
of the integral equation means that the second factor of (10) vanishes as
t → t*. Since the integral factor (6) goes to zero as
t → t*, we can use L'Hôpital's rule on
the ratio of the second factor over the integral factor, which yields that
the right-hand side of (10) approaches the value P(t*,0),
which shows that P(t,z) is continuous up to
the boundary z = 0. The problem is trivial when z
> 0. The proof of part (a) is thus completed.
(b) We first simplify (11). We use the substitution s =
s(t) defined in (14). The inverse formula is
and t = s if λ = μ. The differential
relation is ds = (1 − (λ − μ)s)
dt. From the formula for t*, we find the corresponding
s*:
We can express the characteristic curves in the new variable
s as
where, of course, s(τ) is as defined in (14). With
the above simplification, the integral factor becomes
and the integral equation (11) becomes
Substituting τ into s of (14) we
obtain
Finally, let us apply the change of variables
Notice that the range 0 ≤ z0 <
min{1,μ/λ} corresponds to 0 ≤ y < 1/(λ
− μ) for λ > μ or all y ≥ 0 for λ ≤
μ. For this range of y, the IE transforms to the pivotal
equation (12).
(c) Assume that A(s) has been obtained from
the Volterra integral equation (12). For any given (t,z),
we let z(t) = z in (5), from which we solve for
z0 to obtain (8) and (9). Using this
z0, we find (10) in which a(t) =
A(s(t)). Using (11) in (10) to eliminate
h(z0), we find
From the expression of the integral factor (18) and (16), we do
a change of variable to find
It follows from (8) and (9) that
Let the integral in (20) be denoted by L,
Transforming the variable w →
y,
This completes the proof of Theorem 1. █
It follows from Theorem 1 that the solution P(t,
z) is obtained once we solve the singular Volterra integral
equation (12). From the theory of integral equations, we know that the
singular integral equation (12) can be reduced to a regular equation
and this first-kind equation can be brought to the second kind, for
which the iteration method can be used to construct an infinite-series
solution. The details will not be given here and the reader is referred
to [5], [22], or [23].
In the following sections, we consider the asymptotic behavior of the
solution to this equation and we also consider several cases for which
we are able to find closed-form solutions to this equation.
Remark 1: The probability P0(t)
= a(t) is determined by (11). Once
P0(t) is obtained, recursive formulas can
be developed for Pn(t),
n = 1,2,…, using (2), which can then be used for
calculating explicit numerical solutions. However, for many other
purposes, we prefer to obtain a closed-form explicit formula for
Pn(t), which, among other things,
shows more explicitly its dependencies on the systems parameters.
Indeed, very few stochastic systems possess closed-form solutions for
transient probabilities, and in Section 4, we focus on finding
closed-form solutions for the process.
3. CHANCE OF EXTINCTION
The probability that the population size is zero at time t,
i.e., P0(t), is referred to as the chance
of extinction (Kendall [13]). In
this section we apply Theorem 1 to analyze the limiting behavior of
P0(t) as time tends to infinity.
Theorem 2: As time tends to infinity, we have the
following:
Proof: We notice that (12) can be written as
For the case λ > μ, we let y → 1/(λ
− μ) in (21) to find
If ν/λ ≥ θ/μ, the term (1 −
s)−ν/λ+θ/μ−1
is not integrable on [0,1] . Thus,
A(s/(λ − μ)) must vanish as
s → 1; hence, A(1/(λ − μ)) =
0. In this case, we have
This proves the first half of part (a).
If λ = μ, we let y → ∞ in (21) and use
h(1) = 1 to find
This implies A(∞) = 0 if ν ≥ θ and
A(∞) = (θ − ν)/θ if ν <
θ. This proves part (b) and completes part (a).
If λ < μ, we similarly let y → ∞ in
(21) and use h(1) = 1 to find
Thus,
If ν = 0, we find A(∞) = 1. Otherwise, it is a
number between zero and one since the factor (1 −
λs/μ)−ν/λ is greater
than one. This proves parts (c) and (d). The proof of Theorem 2 is
complete. █
Remark 2: The limit of P0(t)
for the case λ > μ and ν/λ < θ/μ
is generally dependent on the initial condition h(z).
For example, a(t) → 1 if h(z)
≡ 1 for the case λ > μ, ν = 0, any θ > 0,
and μ > 0. However, a(t) → 1 −
((λ − μ)/λ)2 if h(z)
= z for the case λ > μ, ν = 0, and θ =
μ. These results can be obtained from the explicit formulas of
Corollaries 3 and 4 of the next section.
Remark 3: It follows from parts (b) and (d) of Theorem 2
that if λ ≤ μ and ν = 0, then
limt→∞
P0(t) = 1; that is, the population is
“almost certain” to die out for linear birth–death
process with emigration when λ ≤ μ. This result was proved
in Kendall [13] for the special
case ν = θ = 0.
4. CLOSED-FORM SOLUTIONS
Birth and death processes have been extensively studied in the
literature and some cases have been solved explicitly a long time ago.
For example, when λ = μ = 0, we obtain an
M/M/1 queue and its transient solution can be
found, for example, in Kleinrock [15]. The case λ = 0 and θ = 0 is
equivalent to an M/M/∞ queue and its
transient solution can also be found in many texts on stochastic
processes (e.g., Example 4.5 of Cox and Miller [4] and Massey and Whitt [19]). The case ν = 0 and θ = 0 is the
typical linear birth–death process and the case θ = 0 is a
linear birth–death process with immigration, also referred to as
the Kendall process (Kendall [13]).
The transient solutions to all of these cases have been solved in
closed forms (see, e.g., Bailey [1]).
Based on the analysis in Section 2, we present in this section the
closed-form transient solutions for four cases. Our first case is
covered by Ismail et al. [8], but
we choose to include it here since our presentation is simple and more
elementary. The transient solutions to the other cases appear to be
new.
We first consider the case λ = μ. In this case, the integral
equation (12) becomes
We recognize that the right-hand side of this equation is a
convolution. In what follows, we use f * g to
represent the convolution of f and g. Let
then,
By using the Laplace transform ℒ, we have
This analysis yields the following result.
Corollary 1: The integral equation (12) has closed-form
solutions if λ = μ. In this case, a(t) =
A(t) is given by (27) via Laplace transform, and the
solution is
We next consider the case ν = 0, resulting in a linear
birth–death process with emigration. In this case, the integral
equation (12) becomes
where 0 ≤ y < 1/(λ − μ)+.
We note that since ν = 0, we can assume that (28) holds for all
y ≥ 0 (in this case, the term (1 + μy −
λs)−ν/λ is not present).
Note that the right-hand side of the IE (28) is also a convolution. Let
Then
Using the Laplace transform ℒ, we have
More specifically, we find
where Γ is the gamma function. Using αΓ(α) =
Γ(α + 1), we have
If θ/μ = m is a positive integer, then
since H(n)(0) = 0 for all n =
0,1,2,…,m − 1. If θ/μ = m
+ α, where m is a nonnegative integer and α ∈
(0,1), then
since H(n)(0) = 0 for all n =
0,1,2,…,m. Using the convolution formula, we find
We note that (31) is valid even for α = 0; thus, (31) includes
(30) as a special case.
We comment in passing that the equation for A(s)
for the case θ/μ = α ∈ (0,1) is an Abel's
equation; see [5], [22], [23],
or [28]. In particular, if
θ/μ = ½, then we have
Corollary 2: If ν = 0, μ = λ, and
θ = mμ for some positive integer m, we
have
where
and
with
We remark that the function G(y) is very similar to
H(y). We prefer to use G(y) here.
Also, the notation G(k)(λt)
denotes the kth derivative of G(y) with
respect to the generic variable y at λt. One can
see that λ and t appear as a product in the formula, which
is consistent with the equation.
Proof: The function a(t) follows from (30)
and a simple change of variable. From (15), we obtain that
Through Taylor expansion and the beta function, we find
That simplifies to
Using the expansion
we find for all n ≥ 1 that
where the sum is over all k ≥ 1 such that l +
k = n, and it simplifies to yield the result stated
in Corollary 2. The expressions in (33) can be verified easily through
either integration by parts on (34) or inspection of their Taylor
expansion. █
Corollary 3: For the case ν = 0, μ ≠ λ,
and θ = mμ for some positive integer m, we
have
where G(y) = h(y/(1 +
y))ym is the same as in (32),
and
where s = s(t) is defined in (14),
with
Proof: Once we make use of the variable
s(t), the proof will be parallel to that of Corollary
2. We omit the details. █
Corollary 4: Suppose ν = 0, θ/μ =
m + α > 0 for some nonnegative integer m and 0
≤ α < 1, and μ > 0, λ > 0, we
have
where s = s(t) is given in (14). The
corresponding Pn(t) for all
n ≥ 1 are
Proof: It is a simple consequence of Theorem 1 and (31).
█
We next present an alternative method. As used in (21) in the proof
of Theorem 2, we notice that (12) can always be written as
where we have defined
When ν = 0, we have
It follows from (39) that the nth derivative of F
at zero is
for all nonnegative integers n. Thus, for all n, we
obtain A(n)(0) as
where β is the beta function. Then, we find
Since a(τ(s)) = A(s), we have
This method can be used to deal with the general case. We can use
Taylor expansion on F(y), A(sy),
and (1 + (μ −
λs)y)−ν/λ at
y = 0 in (37) to derive a recursive formula for
A(n)(0):
where α = −ν/λ. In this way we can find an
expression for a(t); however, since it is very
complicated, we choose not to present its formula.
5. DISCUSSION
We presented the use of the method of characteristics for a transient
analysis of the linear birth–death processes with immigration and
emigration. For the general case, we presented a Volterra-type integral
equation for extinction probability P0(t)
and an explicit formula for P(t,z) in terms
of P0(t), which is ready for expansion to
yield Pn(t). For several special
cases that have not been previously solved, we presented closed-form
solutions for P(t,z) and
Pn(t) for all n ≥ 0.
We also give semiexplicit closed-form solutions P(t,
z) and P0(t) for the case λ =
μ.
Recently, there has been several studies that incorporate total
catastrophes in the stochastic systems (see, e.g., Kumar and
Arivudainambi [16], and Kyriakidis
[17,18],
Chao and Zheng [3], among others).
The catastrophes arrive according to a Poisson process and its arriving
effect is to reduce the total population size to zero. Let the
catastrophes arrival rate be γ; then PDE (3) becomes
It is not hard to see from our analysis that the procedure used in
solving PDE (3) equally applies to solving PDE (41). Thus, adding total
catastrophes does not add difficulty in solving the problem. The
details will not be given here.
