Asymptotic Analysis of a Loss Model with Trunk Reservation I : Trunks Reserved for Fast Traffic

We consider a model for a single link in a circuit-switched network. The link has C circuits, and the input consists of offered calls of two types, that we call primary and secondary traffic. Of the C links, R are reserved for primary traffic. We assume that both traffic types arrive as Poisson arrival streams. Assuming that C is large and R O 1 , the arrival rate of primary traffic is O C , while that of secondary traffic is smaller, of the order O √ C . The holding times of the primary calls are assumed to be exponentially distributed with unit mean. Those of the secondary calls are exponentially distributed with a large mean, that is, O √ C . Thus, the primary calls have fast arrivals and fast service, compared to the secondary calls. The loads for both traffic types are comparable O C , and we assume that the system is “critically loaded”; that is, the system’s capacity is approximately equal to the total load. We analyze asymptotically the steady state probability that n1 resp., n2 circuits are occupied by primary resp., secondary calls. In particular, we obtain two-term asymptotic approximations to the blocking probabilities for both traffic types.


Introduction
A classic model in teletraffic is the Erlang loss model.Here, we have C servers or circuits , and customers telephone calls arrive as a Poisson process with rate parameter λ.The arriving customer takes one of the circuits if one is available, and if they are all occupied then the call is blocked and lost.When occupying a circuit, the customer has an exponentially distributed holding time whose mean we take as the unit of time.It is well known that the steady state probability that n circuits are occupied is the truncated Poisson distribution, that is, Kλ n e −λ /n!, 0 ≤ n ≤ C, with 1/K C n 0 λ n e −λ /n!.This model dates back to circa 1918 1 .When n C, we obtain the steady state blocking probability.The transient probability distribution is much more complicated, but it can be computed in terms of special functions see 2 .
Over the years, many generalizations of the basic model have been analyzed, including networks of such loss models see 3, 4 .One important extension is that of trunk reservation, which is fundamental in the analysis of circuit-switched communication networks.Here, we consider a model with C circuits that are used by the two types of customers or offered calls .We refer to these as primary or high-priority calls and secondary or low-priority calls.They arrive as Poisson arrival streams with respective rates λ and ν.Of the C circuits, R are reserved for primary calls.Thus, if a high-priority call arrives, it is blocked if all C circuits are busy, while a low-priority call is blocked if at least C − R circuits are busy.All calls are assumed to have independent and exponentially distributed holding times, with respective means 1 and 1/κ.The total load on the system is λ ν/κ.If this exceeds C, then typically all the circuits are busy an overloaded link , while if λ ν/κ < C, typically some circuits are free an underloaded link .An interesting situation is when C is large and λ ν/κ ≈ C; we refer to this as "critical loading." Previous work on this and related models includes Mitra and Gibbens 5 who considered the asymptotic regime λ, C → ∞ with C λ − O √ λ , R O √ λ , and ν O √ λ thus, secondary calls are less frequent than primary calls .We thus have C/λ 1 O 1/ √ λ ; so this is an example of critical loading.They analyzed a single link and used their results to obtain approximations for more complicated loss networks with a distributed, state-dependent, dynamic routing strategy.Related work appears in 6, 7 , and optimization and control policies for such problems were analyzed by Hunt and Laws 8 .
Of fundamental importance in this model is the probability B 1 resp., B 2 that a primary resp., secondary call is blocked and lost in the steady state.Roberts 9, 10 obtained approximations to these blocking probabilities, which are based on a certain recursion which is exact for special cases of the model parameters, but not for all cases.Morrison 11 investigated this model for R O √ λ and R O 1 , and obtained the blocking probabilities as asymptotic series in powers of 1/ √ λ.This led to a better understanding of the asymptotic validity of Roberts' approximation s .However, the coefficients in the asymptotic series in 11 were not explicit, as their calculation still involves recursively solving an infinite system of differential equations.But, if it is further assumed that γ ν/ √ λ is small, the blocking probabilities were obtained more explicitly in terms of parabolic cylinder functions.Also, if R O 1 rather than R O √ λ , explicit results are obtained without the small γ assumption.
In 12 , we analyzed the case R O 1 , with κ O 1 , and with the arrival rates λ and ν both O C .Expressions for the blocking probabilities were obtained for the overloaded and underloaded cases.In the first case, both blocking probabilities remain O 1 as λ → ∞, while in the second case they are exponentially small.In this paper, we investigate the case of critical loading, where again C, λ → ∞ but now with λ ν/κ ∼ C. We will also assume that ν O √ λ secondary calls are less frequent than primary ones but now with κ O 1/ √ λ ; that is, secondary calls have large holding times.Thus, primary calls have faster arrivals and faster service.Note that the loads due to the primary and secondary calls remain asymptotically comparable with this scaling.In this asymptotic regime, we are able to obtain explicit analytic expressions for the first two terms in the expansions in powers of 1/ √ λ for the blocking probabilities, which involve readily evaluated definite integrals.
We are currently investigating situations where the secondary calls are the ones with fast arrivals and service 13 .Here, the asymptotic structure of the problem turns out to be quite different.
We comment that the basic problem to be solved is a two-dimensional difference equation cf.2.1 , with discontinuities of the coefficient functions at various boundaries and an interface.Such a problem appears to be very difficult or impossible to solve exactly.From a numerical point of view, the problem corresponds to solving roughly N C 2 /2 linear equations.A good general method such as Gaussian elimination has computational complexity O N 3 or O C 6 .Some methods that use the sparseness of the system and some iteration procedures may improve this to O N 2 or O C 4 .The purpose of our asymptotic analysis is to obtain reasonable approximations whose numerical evaluation has computational complexity that is independent of C, and also to obtain explicit formulas that show the dependence of the stationary distribution and blocking probabilities on the model parameters i.e., C, R, and the arrival and service rates .
The paper is organized as follows.In Section 2, the problem is stated more precisely and the basic equations are obtained.Here, we summarize our main results, which are derived in detail in Sections 3 and 5.In Section 3, we obtain the leading terms for the blocking probabilities.In Section 4, we relate the present results to the ones in 11, 12 using asymptotic matching.The first-order correction terms to the blocking probabilities are derived in Section 5, while in Section 6 we present some numerical studies to assess the accuracy of the asymptotics.

Statement of the problem and summary of results
We denote by N 1 t the number of servers serving high-priority customers, and by N 2 t the number of servers serving low-priority ones.The total number of servers circuits is C of which R are reserved for the high-priority customers.Thus, if N 1 N 2 C, a newly arriving high-priority customer call is lost; if N 1 N 2 ≥ C − R, then a newly arriving low-priority call is lost.The high-and low-priority customers arrive as independent Poisson processes, with respective rates λ and ν.The service times are exponentially distributed with respective means 1 and 1/κ.Thus, the unit of time is taken as the service rate of the high-priority customers.
We denote the steady state joint distribution of the numbers of servers used by the two priority classes by p n 1 , n 2 lim t→∞ Pr N 1 t n 1 , N 2 t n 2 .We let I{A} be the indicator function on the event A. Then, from the description of the model, we obtain the following balance equation:

2.1
This applies over the domain Thus, we may view the problem as solving a second-order difference equation in two variables, over the triangle with the two subdomains separated by the "interface" { n 1 , n 2 : n 1 n 2 C − R}.There are also boundary conditions inherent in 2.1 , along n 2 0, n 1 0, Of particular interest are the blocking probabilities for the high-priority customers, defined by and for the low-priority ones, defined by

2.5
Note that we clearly have 0 We analyze the problem in the asymptotic limit where We furthermore assume that the arrival rate λ of high-priority customers is large, of the same magnitude as C, and then scale the other rate parameters as Thus, the arrival rate of low-priority customers is large, but only of the order O √ λ .The service times of these customers however are also large, and the total load due to low-priority customers is ν/κ γλ/μ.Also, we have λ ν/κ ∼ C so that the total load due to all customers is roughly equal to the capacity of the system.Hence, this asymptotic limit may certainly be considered as "heavy traffic" or "critical loading." Once we input the model parameters λ, ν, κ, C, and R, we can compute γ, μ, and ω from 2.7 .For some of our numerical studies, it is desirable to fix C, R, γ, μ, and ω and then vary the original rate parameters, which are computed by inverting 2.7 using and ν and κ are obtained from 2.7 once λ is known.We next scale n 1 and n 2 as and x and are taken as O 1 .We consider the scaled state space with X, Y C −1 n 1 , n 2 .Then, since C is large and R is O 1 , the domain in the X, Y plane is the triangle 0 ≤ X Y ≤ 1; X, Y ≥ 0. The scaling 2.9 corresponds to a small neighborhood of the point However, this is where most of the probability mass accumulates in this asymptotic limit, and the analysis of this range is sufficient to obtain the blocking probabilities B j .We will obtain these as asymptotic series in powers of 1/ √ λ.Using 2.9 and 2.10 in 2.1 , we obtain

2.12
Note that in this asymptotic scaling, the indicator functions I{n 1 ≥ 1}, I{n 2 ≥ 1}, and I{n 2 1 ≤ C −R} may be replaced by one, since these correspond to boundaries that are far from the point in 2.11 .However, the interface n 1 n 2 C − R and the boundary n 1 n 2 C are evident in 2.12 and will play a large part in the analysis.In Section 3, we will analyze 2.12 , and then also consider a second scale where is large and negative.This second scale will lead to a diffusion equation in two variables.
The main results are as follows. For where dx.

2.15
The correction term in 2.13 takes the form p For the latter function, we have φ 1 x, y ψ 1 x, x y , where ψ 1 is given by 5.5 in terms of g 0 and g 1 .Then, g 1 is given by 5.66 -5.68 with Λ ξ and E x defined in 5.6 and 5.56 , and the constant A 1 in 5.66 can be obtained by using 5.78 in 5.83 .
On the x, y scale, where where φ 0 x, y e −x 2 /2 g 0 x y and φ 1 is as above.The analysis of these two scales leads to two-term approximations to the blocking probabilities in 2.4 and 2.5 .To leading order, these are where the integral is evaluated using 2.14 and 2.15 see also 3.49 .The O λ −1 correction terms follow from 5.80 and 5.81 , using also 5.74 .We note that the numerical evaluation of the leading order asymptotic results involves only the integrals in 2.15 and 2.17 .The correction terms involve numerically evaluating some double or triple integrals cf.5.78 , but the computational complexity of evaluating the asymptotic results is independent of C or λ .

Asymptotic analysis: leading terms
We consider 2.12 and assume that for λ → ∞ the probabilities have the expansion

3.1
In this section, we focus on the leading term, but its calculation will necessitate that we also analyze the problem for p 1 x .This correction term is calculated completely in Section 5. We will also need to couple the analysis of the scale O 1 to that where is large and negative, with − O √ λ .Using 3.1 in 2.12 and equating coefficients of powers of 1/ √ λ, we obtain to leading order and this applies for all ≤ R.This is a simple difference equation with a boundary condition at R. Its most general solution is where F 0 x is to be determined.For 1 ≤ ≤ R, we obtain from 3.1 , 3.3 , and 2.12 the problem whose solution is Here, x and hence Here, we imposed continuity between 3.5 and 3.7 along 0, and J 1 x is another function not yet determined.By setting 0 in 2.12 and comparing terms of order O 1/ √ λ , we obtain

3.8
Using 3.7 to compute p 1 −1 and 3.5 for p 1 1 and p 1 0 , we obtain 3.9 so that 3.7 becomes We next consider the problem 2.12 for → −∞, with the scaling

3.11
Note that this still corresponds to a local approximation near the point in 2.11 .In terms of x, y , we let 12 and 2.12 , upon multiplying by λ, becomes

3.13
We note that ± 1 corresponds to y ∓ 1/ √ λ, and that for y > 0 the indicator functions in 2.12 can all be replaced by one.
We assume that φ has an expansion in the form

3.14
Using 3.14 in 3.13 , we obtain to leading order the PDE This is a parabolic PDE whose solution is facilitated by the change of variables x ξ, x y η, φ 0 x, y ψ 0 ξ, η .

3.19
We will determine g 0 shortly.
We observe that on the scale, with y − / √ λ, expansion 3.14 becomes

3.20
Comparing this to 3.1 , for < 0 we conclude from 3.3 that and, from 3.10 that

3.22
It follows that which is a boundary condition for the PDE 3.15 along y 0. In terms of ξ, η , this becomes ψ 0 ξ ξ, ξ ξψ 0 ξ, ξ 0, but this holds automatically for any g 0 in view of 3.18 .To determine g 0 , we must analyze the correction term in 3.14 .
From 3.13 and 3.14 , we find that the first correction term φ 1 satisfies the PDE

3.24
Switching to the ξ, η variables and then using 3.18 , we get

3.25
Integrating 3.25 with respect to ξ yields

3.27
We will show that 3.27 leads to a differential equation for g 0 x .But, we must first consider φ 1 y x, 0 .In view of 3.20 , this term arises as a part of the O 1/λ term in the expansion on the scale.We therefore return to 3.1 and 2.12 .
For 1 ≤ ≤ R, we obtain from 2.12 and 3.1 , at order O 1/λ , the equation

10
Journal of Applied Mathematics and Stochastic Analysis Using 3.3 and 3.5 , we find that 3.28 has a solution in the form

3.31
It follows that

3.32
We note that p

3.34
Using 3.3 and 3.10 , we find after some calculation that 3.34 simplifies to and hence Setting 0 in 2.12 , we obtain at order O 1/λ the interface relation

3.37
Using 3.36 , with 0 and −1, and 3.33 , we obtain from 3.37 after some calculation

3.38
Now, comparing the − /λ terms in 3.20 and 3.1 with 3.36 , we conclude that .27 involves F 0 x only, and since F 0 and g 0 are related via 3.21 , we obtain from 3.27 the following ODE for g 0 x : e −u 2 /2 du γg 0 x μ ω x g 0 x μg 0 x 0.

3.41
Integrating once and requiring g 0 to vanish as x → ∞ yield and hence

3.43
where A 0 is a constant, which will be fixed by normalization.We use 2.10 and 3.12 in the normalization sum 2.3 , and then use the Euler-MacLaurin formula to approximate sums by integrals.To leading order, this yields 3.44 We use 3.19 and 3.43 and evaluate one of the two integrals in 3.44 using integration by parts, with the result

3.45
To summarize the calculation, we have shown that on the x, y scale in 3.11 we have which applies for −∞ < ≤ R and is independent of .To calculate correction terms to 3.46 and 3.47 , we would need to find F 1 x in 3.5 and 3.10 , and solve 3.24 for φ 1 x, y .This ultimately involves calculating the O 1/λ and O 1/λ 3/2 terms in 3.1 and the O 1/λ term in 3.14 ; this is done in Section 5.
We next calculate the blocking probabilities in 2.4 and 2.5 .Evaluating these sums requires the expansion on the x, scale.Again, approximating sums by integrals and using the scaling 2.9 , we obtain as p 0 x F 0 x .From 3.21 , 3.43 , and 3.45 , we then obtain The numerical accuracy of 3.48 is investigated in Section 6.This completes the analysis of the leading terms.

Consistency with previous results
In where We show that 4.1 matches asymptotically to 3.48 , in the limit where μ → ∞.When μ → ∞, we can simplify the integrals in both the numerator and denominator in 3.49 , as the factor exp −μ ω x 2 /2γ has the effect of freezing the remaining integrands at x −ω.Thus, by the Laplace method, we obtain We consider first the asymptotic matching of 3.48 -3.49 to the underloaded case in 12 .We note that the parameters ρ and σ are related to the current ones by For the matching, we must thus let ω → −∞ with |ω| √ λ, and ρ/σ ↑ 1.The results in 12 for B 1 and B 2 were as follows: where In view of 4.4 , we have, in the matching region, ρ ∼ σ, a ∼ 1, κ/σ ∼ μ/ μ γ , and Hence, and 4.5 becomes

Journal of Applied Mathematics and Stochastic Analysis
We show that 4.9 agrees with 3.48 -3.49 when the latter is expanded for ω → −∞.In 3.49 , the main contribution to the integral in the denominator comes from x −ω |ω|, and in the numerator from x μ|ω|/ γ μ , where d/dx x 2 /2 μ/2γ ω x 2 0. Thus, we can approximate ∞ −x e −u 2 /2 du by √ 2π everywhere, and obtain
Next, we consider the overloaded case in 12 , for which we obtained where and ζ is the solution of We note that by using 4.12 in 4.13 , ζ is a particular root of a polynomial.For the matching, we will take ρ/σ ↓ 1 in 4.11 and ω → ∞ in 3.48 -3.49 .As ρ/σ → 1, we will have ζ → 1 and this will allow us to simplify 4.11 .Let us set

4.16
Using 4.14 -4.16 in 4.13 , we find after some calculation that θ ∼ μω R 1 γ μ 4.17 and then 4.11 , for ρ/σ ↓ 1, simplify to We expand 3.49 for ω → ∞.Scaling x ωt, we obtain For t > 0, we use the asymptotic approximation and conclude that both integrals in 4.19 have their major contribution from where But then by the Laplace method, we have and 3.48 agrees with 4.18 .This completes the matching verifications.

Correction terms
We will compute the O 1/ √ λ terms in expansions 3.1 and 3.14 , and then obtain O 1/λ corrections to the blocking probabilities.
We first consider 3.13 .Using the relations φ x, y we obtain from 3.13 where L i are the operators

5.4
Before analyzing 5.3 , we obtain a more complete description of ψ 1 .We recall that ψ 0 is known completely, in view of 3.18 , 3.43 , and 3.45 .But, we computed ψ 1 only partially as the combination ψ 1 ξ ξψ 1 in 3.26 .We integrate 3.26 to get where 1 − e ξt t dt.

5.8
The function g 1 will be determined shortly actually, not very shortly, but only after a lengthy calculation .
We evaluate the right side of 5.3 more explicitly.Some terms are expressible as derivatives with respect to ξ, while the ones involving derivatives in η may be evaluated using 3.18 and 5.5 .Then, 5.3 becomes μ ω η g 0 η μg 0 η .

5.9
From 3.18 , we have and by direct calculation

5.11
With the above, we integrate 5.9 to get where

5.13
We will show that the calculation of g 1 will require only that we evaluate 5.12 along η ξ.However, we must first reconsider scale O 1 and analyze at least partly the term p 3 x in 3.1 i.e., the coefficient of λ −3/2 in the series .
Returning to 2.12 with the expansion 3.1 , we find that for 1 ≤ ≤ R, p 3 satisfies

5.14
We recall that for 0 ≤ ≤ R, p 0 is given by 3.3 , p 1 by 3.5 , and p 2 by 3.29 .Let us write

Journal of Applied Mathematics and Stochastic Analysis
We then rewrite 5.14 as

5.17
After some rearrangement, we obtain

5.18
where all derivatives are with respect to x. Setting 0 in 2.12 and using again expansion 3.1 , we obtain at O λ −3/2 the relation

5.19
Next, we consider 2.12 for < 0, and recall that p 0 is given by 3.3 for all , p 1 is given by 3.7 or 3.10 , and where J 2 is in 3.27 , and 3.36 yields

5.22
From 2.12 , for ≤ − 1, we obtain the following problem for p 3 :

5.23
Here, we used I{ ≤ R − 1} 1 in this range of , and all functions in 5.23 are evaluated at x.After rearranging 5.23 and using 5.20 , 3.3 , and 3.7 to evaluate p j for j 0, 1, 2, we are led to

5.24
We sum 5.24 from −m to −1 with m ≥ 1 to obtain

5.25
Note that 5.25 remains true if m 0. Using 5.22 in 5.25 gives us

5.26
Using our previous results for p

5.27
Adding 5.26 and 5.27 , we obtain an explicit expression for p and K 2 which is quadratic in m.By summing from m 0 to m n − 1, we get

5.28
This holds for all n ≥ 0. We write 5.28 as where F 3 x p 3 0 x , and J 3 , K 3 , and L 3 may be identified from 5.28 .We recall that p −n x φ x, n/ √ λ for n ≥ 0, and the coefficient of λ −3/2 in the expansion 3.20 is Comparing 5.29 to 5.30 yields

5.31
In view of 3.19 , we have e −x 2 /2 g 0 x , 5.32 while 5.28 shows that

5.33
But, from 5.21 and 3.21 , we get K 2 x 1/2 e −x 2 /2 g 0 x , and from 3.9 , J 1 x −e −x 2 /2 g 0 x .Then, we can easily verify that 5.32 is consistent with 5.33 .Also, from 5.28 , we find that

5.34
We show that this is the same as −φ 1 yy x, 0 .We recall that φ 1 x, y ψ 1 x, x y is given by 5.5 and J 2 is expressed in terms of F 1 in 3.38 .From 3.38 and 3.21 , it follows that

5.35
where D x is the operator in 5.13 , with η replaced by x.The second equality in 5.35 follows from 3.40 and 5.6 .Using 3.9 , we obtain which when combined with 5.35 gives

5.37
We use 5.37 in 5.34 , also noting that which follows from 5.5 .Then, 5.34 becomes

5.39
Here, we also used Λ x xΛ x − 1.Now, from 5.21 and 3.9 , we obtain x γ g 0 x .

5.40
Using 5.38 to compute J 2 xJ 2 and 5.40 , we get In view of 5.5 , the above is the same as 2K 3 x φ 1 yy x, 0 ψ 1 ηη ξ, ξ .We next examine the relation J 3 x −φ 2 y x, 0 −ψ 2 η ξ, ξ .We will use this to ultimately obtain g 1 η in 5.5 , which will complete the determination of the O 1/ √ λ correction terms in 3.1 and 3.14 .We first note from 5.28 and 5.29 that

5.42
We solve 5.42 for the combination J 3 F 2 xF 2 , that we rewrite as where

5.44
Next we note that Z and U may be integrated explicitly and we write where

5.47
It follows that x, x xψ 2 x, x .

5.48
The right side of 5.48 was computed in 5.12 .
Using F 1 x ψ 1 x, x , 5.38 , and the identities F 0 x e −x 2 /2 g 0 x , J 1 x − e −x 2 /2 g 0 x , and K 2 x 1/2 e −x 2 /2 g 0 x , we evaluate W x in terms of g 0 x , and then use 5.12 .After some simplification, this leads to

5.49
Since this must be equal to − Z RU , we try to write the right side of 5.49 as a perfect derivative.To this end, we note that

5.50
Adding Z RU to 5.49 , we rewrite that equation as

5.51
We next evaluate Z in terms of g 0 , and U in terms of g 0 and g 1 , and then we integrate 5.51 and thus explicitly obtain g 1 .Since J 1 −xF 0 − F 0 and K 2 is in 5.21 , we have

5.52
We thus note that Z x is canceled by the bracketed term that precedes it in 5.51 .Using 5.46 , we explicitly calculate U x , recalling that G 1 is given by 5.15 , and 2 can be computed from 3.30 and 3.32 .After some cancellation of terms, we obtain

5.53
Using 5.52 and 5.53 , we integrate 5.51 , subject to the condition that the solution decays exponentially as x → ∞.Hence,

5.54
What remains is a linear first-order ordinary differential equation for g 1 , which is readily solved by multiplying by the integrating factor: .

5.55
We introduce the notation

5.57
Then, we have

5.58
From 3.43 , we have and thus

5.60
With 5.60 , we have

5.65
Here, A 1 is a constant that will be fixed by normalization.We thus write g 1 as where

5.68
where D v g 0 v γg 0 v μ ω v g 0 v μg 0 v is as in 5.13 .We next determine A 1 by normalization and then obtain correction terms to the blocking probabilities B 1 and B 2 .This requires that we evaluate the integrals , η dη dξ.From 5.5 , we have 5.69 and from 3.43 , we calculate D x g 0 and obtain

5.70
Consider the contribution to ∞ −∞ φ 1 x, 0 dx that comes from the term Λ x D x g 0 x in 5.69 and the part of g 1 that is proportional to N x cf.5.66 .We obtain

5.71
Here, we wrote and integrated by parts.We also have

5.73
Using 5.71 and 5.73 , we integrate 5.69 and get

5.74
Here, we also used 5.70 to eliminate D x g 0 from the expression.

5.83
In view of 5.78 , 5.83 may be viewed as a linear equation for A 1 , and thus all the correction terms are now known fully.
To summarize the calculations in this section, we have determined g 1 x in 5.66 -5.68 , with A 1 computed from 5.78 and 5.83 .In terms of g 1 , the second term in the expansion on the x, y or ξ, η x, x y scale is given by 5.5 .Then, F 1 x φ 1 x, 0 ψ 1 ξ, ξ and the second terms on the scale are given by 3.5 for 0 ≤ ≤ R, and by 3.10 for < 0 with F 0 x in 3.21 and 3.43 .

Numerical studies
We test the numerical accuracy of our asymptotic expansions, focusing on the blocking probabilities B 1 and B 2 .The numerical results are obtained by solving the linear system 2.1 with the normalization 2.3 .We simply omitted the equation with n 1 n 2 0 in 2.1 and replaced it by 2.3 , thus obtaining an inhomogeneous problem with a unique solution.
We solved 2.1 by two different methods.First, we simply used the program MAPLE to solve the linear system numerically.We also tried an iteration method of the form p n 1 , n 2 ; M 1 p n 1 , n 2 ; M T/MAX • Lp n 1 , n 2 ; M , where Lp 0 is the basic equation 2.1 .Starting from some initial guess p n 1 , n 2 ; 0 and iterating up to M MAX − 1 correspond to solving approximately for the transient solution for this model, from time t 0 to t T .We verified that choosing T sufficiently large leads to the same results as the MAPLE solution of 2.1 .
Since the asymptotic results are expansions in powers of 1/ √ λ with coefficients expressed in terms of γ, μ, ω , we input the five parameters C, R, γ, μ, ω , calculate λ from 2.8 , then ν and κ from 2.7 , and solve 2.1 numerically.In Table 1, we have R 2, γ 1, μ 1, and ω 1, and we compare the exact numerical results for B 1 and B 2 with the one-and two-term asymptotic approximations.We give various values of C and also tabulate the corresponding We next consider a different purely numerical approach to estimating the coefficients in the expansions of the B j .We choose some C0, and for C C0 − 1, C0, and C0 1, we equate 6.1 Note that B i B i C and λ λ C , for fixed values of R, γ, μ, ω .Thus, 6.1 may be viewed as 3 × 3 systems of linear equations for the T i and S i , respectively.This allows us to numerically estimate the first three coefficients in the asymptotic series.In Table 4, we consider C0 in the range of 5 to 70, and give the T i and S i , fixing R, γ, μ, ω 2, 1, 1, 1 .We see that the sequence of T 1 and S 1 does appear to converge as C0 → ∞; the convergence of T 2 and S 2 is slower, and that of T 3 and S 3 is even slower.The asymptotic results in Sections 3 and 5 show that for these parameter values This is in good agreement with Table 4.The data in Table 4 also give a rough estimate of the third O λ −3/2 terms in the expansions of the blocking probabilities.These can be computed analytically by continuing our expansions further, but the calculations are too foreboding.
In Tables 5 and 6, we again give the T i and S i for C0 between 5 and 70, but now with ω 0 and ω −1, respectively.For these values, our asymptotic analysis predicts that Again this is in good agreement with the apparent limiting values of T 1 , T 2 , S 1 , and S 2 as C0 → ∞.The data in Tables 4-6 show that the expansions do indeed appear to be in powers of 1/ √ λ, and that we correctly computed the leading two terms.Note that in each case, the second coefficient T 2 and S 2 is negative, while the first and third ones are positive.This is consistent with the fact that in Tables 1-3 the leading terms always overestimate the exact answer, while the two-term approximations underestimate it.As we decrease ω, the ratio |T 2 /T 1 | increase, as do | S 2 /S 1 | though these are always larger .Hence, we expect that decreasing ω leads to further cancellation between the first and second terms in the asymptotic series, and this again is in agreement with the data in Tables 1-3.It also explains why the two-term asymptotic approximations to B 1 sometimes lead to negative answers, for moderate C values.
To summarize, we have shown that the asymptotic approximations are reasonably accurate, though certainly not excellent, and that there is merit to computing the O 1/λ correction terms unless C is quite small.For small C, however, the one-term approximations may be superior, as the two-term approximations may lead to negative answers.The accuracy of the asymptotic approximations presumably increases as C increases further, and the two-term approximations are presumably better than the one-term approximations.However, limitations of the available computing facilities have so far prevented the evaluation of the exact numerical results for larger values of C.

2 − 1
from both sides of 5.19 and substitute 5.18 into the resulting equation.We then use 3.3 , 3.10 , 5.20 , 5.15 , and 3.29 , and after some calculation, obtain

Table 1 :
R 2, γ 1, μ 1, ω 1. R, leads to the third term in 5.82 .For ≤ 0, the expansion on the y scale contains that on the scale.The leading term in 5.82 regains 3.45 and determines A 0 .The O 1/

Table 3 :
R 2, γ 1, μ 1, ω −1. as computed from 2.8 .We see that the one-term approximations always overestimate the true values, while the two-term approximations underestimate them.The two-term approximations are more accurate especially for the second blocking probability B 2 and for larger values of C. In Table2, we have ω 0, and in Table3, ω −1, with the other parameter values unchanged.With decreasing ω which corresponds to increasing the total load cf.2.7 , we get similar results, but the overall asymptotics both one-and two-term are getting somewhat worse.Also note that the two-term approximations may sometimes lead to negative answers, and this is explained in what follows.