# The following statements compute the sojourn time trivariate # cumulative distribution function associated with the # five cases for the first three customers in an M/M/1 # queue that is initially empty and idle. For each case # the partial derivatives are taken resulting in the # trivariate probability distribution functions. These # PDFs are integrated over their respective regions to # calculate their contribution to the trivariate # distribution. The contributions are summed for all # five cases providing the total of 1 as required. restart; assume(lambda > 0); assume(mu > 0); # Case A FA := simplify((1 - exp(-(lambda + mu) * t1)) * (1 - exp(-(lambda + mu) * t2)) * (1 - exp(-mu * t3))): # Case B # # FBtop is the joint cdf of t1, t2, t3 when t2 < t3 # FBbot is the joint cdf of t1, t2, t3 when t2 > t3 FBtop := simplify(int( (1 - exp(-(lambda + mu) * t1)) * (1 - exp(-(lambda + mu) * (t2 - y))) * (1 - exp(-mu * (t3 - y))) * mu * exp(-mu * y), y = 0 .. t2)): FBbot := simplify(int( (1 - exp(-(lambda + mu) * t2)) * (1 - exp(-(lambda + mu) * (t1 - y))) * (1 - exp(-mu * (t3 - y))) * mu * exp(-mu * y), y = 0 .. t3)): # Case C # # FCtop is the joint cdf of t1, t2, t3 when t1 < t2 # FCbot is the joint cdf of t1, t2, t3 when t1 > t2 FCtop := simplify(int( (1 - exp(-(lambda + mu) * (t1 - y))) * (1 - exp(-(lambda + mu) * (t2 - y))) * (1 - exp(-mu * t3)) * (lambda + mu) * exp(-(lambda + mu) * y), y = 0 .. t1)): FCbot := simplify(int( (1 - exp(-(lambda + mu) * (t1 - y))) * (1 - exp(-(lambda + mu) * (t2 - y))) * (1 - exp(-mu * t3)) * (lambda + mu) * exp(-(lambda + mu) * y), y = 0 .. t2)): # Case D # # FD1 is the joint cdf of t1, t2, t3 in Case 1 # FD2 is the joint cdf of t1, t2, t3 in Case 2 # FD3 is the joint cdf of t1, t2, t3 in Case 3 # FD4 is the joint cdf of t1, t2, t3 in Case 4 # FD5 is the joint cdf of t1, t2, t3 in Case 5 integrand := (1 - exp(-(lambda + mu) * (t1 - y))) * (1 - exp(-(lambda + mu) * (t2 - y - z))) * (1 - exp(-mu * (t3 - z))) * (lambda + mu) * exp(-(lambda + mu) * y) * mu * exp(-mu * z): FD1 := simplify(int(int( integrand, z = 0 .. (t2 - y)), y = 0 .. t2)): FD2 := simplify(int(int( integrand, z = 0 .. (t2 - y)), y = 0 .. t1)): FD3 := simplify(int(int( integrand, y = 0 .. (t2 - z)), z = 0 .. t3)): FD4 := simplify(int(int( integrand, z = 0 .. t3), y = 0 .. (t2 - t3)) + int(int( integrand, z = 0 .. (t2 - y)), y = (t2 - t3) .. t1)): FD5 := simplify(int(int( integrand, z = 0 .. t3), y = 0 .. t1)): # Case E # # FE1 is the joint cdf of t1, t2, t3 in Case 1 # FE2 is the joint cdf of t1, t2, t3 in Case 2 # FE3 is the joint cdf of t1, t2, t3 in Case 3 # FE4 is the joint cdf of t1, t2, t3 in Case 4 # FE5 is the joint cdf of t1, t2, t3 in Case 5 integrand := (1 - exp(-(lambda + mu) * (t1 - y - z))) * (1 - exp(-mu * (t3 - z - w))) * (lambda + mu) * exp(-(lambda + mu) * y) * mu * exp(-mu * z) * mu * exp(-mu * w): FE1 := simplify(int(int(int( integrand, w = 0 .. (t2 - y - z)), z = 0 .. (t2 - y)), y = 0 .. t2)): FE2 := simplify(int(int(int( integrand, w = 0 .. (t2 - y - z)), z = 0 .. (t1 - y)), y = 0 .. t1)): FE3 := simplify(int(int(int( integrand, y = 0 .. (t2 - w - z)), w = 0 .. (t3 - z)), z = 0 .. t3)): FE4a := simplify(int(int(int( integrand, w = 0 .. (t3 - z)), z = 0 .. (t1 - y)), y = 0 .. (t2 - t3)) + int(int(int( integrand, w = 0 .. (t2 - y - z)), z = 0 .. (t1 - y)), y = (t2 - t3) .. t1)): FE4b := simplify(int(int(int( integrand, w = 0 .. (t3 - z)), z = 0 .. (t1 - y)), y = 0 .. t1)): FE5a := simplify(int(int(int( integrand, y = 0 .. (t1 - z)), z = 0 .. (t3 - w)), w = 0 .. (t2 - t1)) + int(int(int( integrand, y = 0 .. (t2 - w - z)), z = 0 .. (t3 - w)), w = (t2 - t1) .. t3)): FE5b := simplify(int(int(int( integrand, y = 0 .. (t1 - z)), z = 0 .. (t3 - w)), w = 0 .. t3)): > # Check the PDF in the A case fA := simplify(diff(diff(diff(FA, t1), t2), t3)): simplify(int(int(int(fA, t1 = 0 .. infinity), t2 = 0 .. infinity), t3 = 0 .. infinity)): # Check the PDF in the B cases fBtop := simplify(diff(diff(diff(FBtop, t1), t2), t3)): fBbot := simplify(diff(diff(diff(FBbot, t1), t2), t3)): pBtop := simplify(int(int(int(fBtop, t2 = 0 .. t3), t3 = 0 .. infinity), t1 = 0 .. infinity)): pBbot := simplify(int(int(int(fBbot, t3 = 0 .. t2), t2 = 0 .. infinity), t1 = 0 .. infinity)): simplify(pBtop + pBbot): # Check the PDF in the C cases fCtop := simplify(diff(diff(diff(FCtop, t1), t2), t3)): fCbot := simplify(diff(diff(diff(FCbot, t1), t2), t3)): pCtop := simplify(int(int(int(fCtop, t1 = 0 .. t2), t2 = 0 .. infinity), t3 = 0 .. infinity)): pCbot := simplify(int(int(int(fCbot, t2 = 0 .. t1), t1 = 0 .. infinity), t3 = 0 .. infinity)): simplify(pCtop + pCbot): # Check the PDF in the D cases fD1 := simplify(diff(diff(diff(FD1, t1), t2), t3)): fD2 := simplify(diff(diff(diff(FD2, t1), t2), t3)): fD3 := simplify(diff(diff(diff(FD3, t1), t2), t3)): fD4 := simplify(diff(diff(diff(FD4, t1), t2), t3)): fD5 := simplify(diff(diff(diff(FD5, t1), t2), t3)): pD1a := simplify(int(int(int(fD1, t2 = 0 .. t1), t1 = 0 .. t3), t3 = 0 .. infinity)): pD1b := simplify(int(int(int(fD1, t2 = 0 .. t3), t3 = 0 .. t1), t1 = 0 .. infinity)): pD1 := simplify(pD1a + pD1b): pD2 := simplify(int(int(int(fD2, t1 = 0 .. t2), t2 = 0 .. t3), t3 = 0 .. infinity)): pD3 := simplify(int(int(int(fD3, t3 = 0 .. t2), t2 = 0 .. t1), t1 = 0 .. infinity)): pD4a := simplify(int(int(int(fD4, t2 = t3 .. (t1 + t3)), t1 = 0 .. t3), t3 = 0 .. infinity)): pD4b := simplify(int(int(int(fD4, t2 = t1 .. (t1 + t3)), t3 = 0 .. t1), t1 = 0 .. infinity)): pD4 := simplify(pD4a + pD4b): pD5 := simplify(int(int(int(fD5, t2 = (t1 + t3) .. infinity), t1 = 0 .. infinity), t3 = 0 .. infinity)): simplify(pD1 + pD2 + pD3 + pD4 + pD5): # Check the PDF in the E cases fE1 := simplify(diff(diff(diff(FE1, t1), t2), t3)): fE2 := simplify(diff(diff(diff(FE2, t1), t2), t3)): fE3 := simplify(diff(diff(diff(FE3, t1), t2), t3)): fE4a:= simplify(diff(diff(diff(FE4a, t1), t2), t3)): fE4b:= simplify(diff(diff(diff(FE4b, t1), t2), t3)): fE5a := simplify(diff(diff(diff(FE5a, t1), t2), t3)): fE5b := simplify(diff(diff(diff(FE5b, t1), t2), t3)): pE1a := simplify(int(int(int(fE1, t2 = 0 .. t1), t1 = 0 .. t3), t3 = 0 .. infinity)): pE1b := simplify(int(int(int(fE1, t2 = 0 .. t3), t3 = 0 .. t1), t1 = 0 .. infinity)): pE1 := simplify(pE1a + pE1b): pE2 := simplify(int(int(int(fE2, t1 = 0 .. t2), t2 = 0 .. t3), t3 = 0 .. infinity)): pE3 := simplify(int(int(int(fE3, t3 = 0 .. t2), t2 = 0 .. t1), t1 = 0 .. infinity)): pE4a:= simplify(int(int(int(fE4a, t2 = t3 .. (t1 + t3)), t1 = 0 .. t3), t3 = 0 .. infinity)): pE4b:= simplify(int(int(int(fE4b, t2 = (t1 + t3) .. infinity), t1 = 0 .. t3), t3 = 0 .. infinity)): pE4 := simplify(pE4a + pE4b): pE5a := simplify(int(int(int(fE5a, t2 = t1 .. (t3 + t1)), t3 = 0 .. t1), t1 = 0 .. infinity)): pE5b := simplify(int(int(int(fE5b, t2 = t3 .. (t1 + t3)), t1 = 0 .. t3), t3 = 0 .. infinity)): pE5 := simplify(pE5a + pE5b): simplify(pE1 + pE2 + pE3 + pE4 + pE5):