# copyright Barry L. Nelson and Michael R. Taaffe 1999
with(linalg):
with(plots):
PhPhInfMoments := proc(A, lambda, B, mu, lower, upper)
#
# this procedure returns a plot of the mean and variance
# of number in queue for the Ph_t/Ph_t/Infinity queue as given
# in Section 4 of Nelson and Taaffe (1999)
#
# initial conditions for the queue are empty and idle
#
# this version is designed to work with Maple V release 5
#
# last update 8/1/99
#
# Inputs:
# A = Markov chain transition matrix for the arrival process
# B = Markov chain transition matrix for the service process
# mdeqs = numerical procedure for moment differential equations
# lower = lower time limit for plot
# upper = upper time limit for plot
#
#
# Output:
# odeplot of mean and variance of number in queue for
# lower < t < upper
local mA, mB, N, NN, k, i, j, g1, a, b, c, d, mdeqs;
#
# create index mapping function
g1 := proc(i,j)
RETURN(N.i.N.j);
end;
#
# obtain procedure to evaluate moment differential equations
mdeqs := PhPhInfNumericEmpty(A, lambda, B, mu):
#
# determine number of phases in phase representations
mA := rowdim(A) -1;
mB := rowdim(B) -1;
# create the N vector
N := matrix(mB, 1, [N.(1..mB)]);
N := N(t);
# create the NN matrix
NN := vector(mA);
for k from 1 to mA do
NN[k] := map(cat, matrix(mB,mB,g1), A.k)(t);
od;
# create the mean and variance plot
odeplot(mdeqs, [[t, sum(N[a,1], a=1..mB)],
[t, sum( sum( sum(NN[b][c,d], b=1..mA) - N[c,1]*N[d,1], c=1..mB), d=1..mB)]], lower..upper, color=black);
end;