SteadyStates1

182 days ago by sma5jlw

var('a,a1,a2,y,i,n,N,E,L1,L2,p,f,g') a = 1.79111 a1 = 0.81787 a2 = 0.97323 n=2 N = sum(((a^i)/factorial(i)),i,0,n) E = (a^n/factorial(n))/N L1 = 0.5*(1+(a1/n) + (a2/n)*(1-p)-((1+(a1/n) + (a2/n)*(1-p))^2-(4*a1/n))^(1/2)) L2 = 0.5*(1+(a1/n) + (a2/n)*(1-p)+((1+(a1/n) + (a2/n)*(1-p))^2-(4*a1/n))^(1/2)) 
       
f(p,y) = ((n-a)*(E)/((n-a)+(a*E)))*(((1-p)*(1-(L2*p))*(L1^y))) show(f(p,y)) # This is the generating function 
        
\newcommand{\Bold}[1]{\mathbf{#1}}0.0883823990144534 \, {\left(p - 1\right)} {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{y}
\newcommand{\Bold}[1]{\mathbf{#1}}0.0883823990144534 \, {\left(p - 1\right)} {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{y}
# To find P(0) Cat A's in the system, we need to find coefficient of p^0 before differentated (with any y power)...is this correct methodoogy? But remember this function is only relevant when there is a queue!! # To find P(1) would need to find the coefficient of p^1 - can find this by differentiating, and finding the coefficient of p(0) # To find P(2) would need to find the coefficient of p^2 - can find this by differentiating twice, and finding the coefficient of p(0)/2! # To find P(3) would need to find the coefficient of p^3 - can find this by differentiating three times, and finding the coefficient of p(0)/3! # And to find P(3 cat A's, 2 Cat B for example, we need to find function of coefficient of p^3, y^2???) # There may even be a simpler way of finding the coefficient...see next box below! g(p,y) = diff(f(p,y),p,y) # This is the function differentiated once - the coefficient of p(0) here should give the probability of 1 Cat A customer in the system if the formula still applied here?! (Long way?!) i(p,y) = diff(f(p,y),p,y,2) # This is the function differentiated twice?! The coefficent of p(0)/2! here would be the prob of 1 cat A in the system? show(g(p,y)) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}0.0883823990144534 \, {\left(-\frac{0.250000000000000 \, {\left(0.473588316450000 \, p - 1.84480612650000\right)}}{\sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000}} - 0.243307500000000\right)} {\left(p - 1\right)} {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} y {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{{\left(y - 1\right)}} \log\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right) + 0.0883823990144534 \, {\left(-\frac{0.250000000000000 \, {\left(0.473588316450000 \, p - 1.84480612650000\right)}}{\sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000}} - 0.243307500000000\right)} {\left(p - 1\right)} {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{{\left(y - 1\right)}} + 0.0883823990144534 \, {\left(p - 1\right)} {\left({\left(\frac{0.250000000000000 \, {\left(0.473588316450000 \, p - 1.84480612650000\right)}}{\sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000}} - 0.243307500000000\right)} p - 0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{y} \log\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right) + 0.0883823990144534 \, {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{y} \log\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)
\newcommand{\Bold}[1]{\mathbf{#1}}0.0883823990144534 \, {\left(-\frac{0.250000000000000 \, {\left(0.473588316450000 \, p - 1.84480612650000\right)}}{\sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000}} - 0.243307500000000\right)} {\left(p - 1\right)} {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} y {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{{\left(y - 1\right)}} \log\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right) + 0.0883823990144534 \, {\left(-\frac{0.250000000000000 \, {\left(0.473588316450000 \, p - 1.84480612650000\right)}}{\sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000}} - 0.243307500000000\right)} {\left(p - 1\right)} {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{{\left(y - 1\right)}} + 0.0883823990144534 \, {\left(p - 1\right)} {\left({\left(\frac{0.250000000000000 \, {\left(0.473588316450000 \, p - 1.84480612650000\right)}}{\sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000}} - 0.243307500000000\right)} p - 0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{y} \log\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right) + 0.0883823990144534 \, {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{y} \log\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)
h(p,y) = simplify(f(p,y)) show(h(p,y)) # An initial attempt to find the coefficient of p(0) - probability of 0 Cat A's in system (if gen fn worked for this number...on;y really valid if there is a queue!) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}0.0883823990145 \, {\left(p - 1\right)} {\left({\left(-0.2433075 \, p + 0.5 \, \sqrt{{\left(-0.486615 \, p + 1.89555\right)}^{2} - 1.63574} + 0.947775\right)} p - 1\right)} {\left(-0.2433075 \, p - 0.5 \, \sqrt{{\left(-0.486615 \, p + 1.89555\right)}^{2} - 1.63574} + 0.947775\right)}^{y}
\newcommand{\Bold}[1]{\mathbf{#1}}0.0883823990145 \, {\left(p - 1\right)} {\left({\left(-0.2433075 \, p + 0.5 \, \sqrt{{\left(-0.486615 \, p + 1.89555\right)}^{2} - 1.63574} + 0.947775\right)} p - 1\right)} {\left(-0.2433075 \, p - 0.5 \, \sqrt{{\left(-0.486615 \, p + 1.89555\right)}^{2} - 1.63574} + 0.947775\right)}^{y}
simplify(g(p,y)) # An initial attempt to find the coefficient of p(0)/1! - probability of 1 Cat A's in system (if gen fn worked for this number...on;y really valid if there is a queue!) 
        
0.0883823990145*((-0.118397079113*p + 0.461201531625)/sqrt((-0.486615*p
+ 1.89555)^2 - 1.63574) - 0.2433075)*(p - 1)*((-0.2433075*p +
0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)*p -
1)*y*(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) +
0.947775)^(y - 1)*log(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2
- 1.63574) + 0.947775) + (0.0883823990145*(-0.118397079113*p +
0.461201531625)/sqrt((-0.486615*p + 1.89555)^2 - 1.63574) -
0.0215041005482)*(p - 1)*((-0.2433075*p + 0.5*sqrt((-0.486615*p +
1.89555)^2 - 1.63574) + 0.947775)*p - 1)*(-0.2433075*p -
0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)^(y - 1) +
0.0883823990145*(p - 1)*(((0.118397079113*p -
0.461201531625)/sqrt((-0.486615*p + 1.89555)^2 - 1.63574) - 0.2433075)*p
- 0.2433075*p + 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) +
0.947775)*(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574)
+ 0.947775)^y*log(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 -
1.63574) + 0.947775) + 0.0883823990145*((-0.2433075*p +
0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)*p -
1)*(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) +
0.947775)^y*log(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 -
1.63574) + 0.947775)
0.0883823990145*((-0.118397079113*p + 0.461201531625)/sqrt((-0.486615*p + 1.89555)^2 - 1.63574) - 0.2433075)*(p - 1)*((-0.2433075*p + 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)*p - 1)*y*(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)^(y - 1)*log(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775) + (0.0883823990145*(-0.118397079113*p + 0.461201531625)/sqrt((-0.486615*p + 1.89555)^2 - 1.63574) - 0.0215041005482)*(p - 1)*((-0.2433075*p + 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)*p - 1)*(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)^(y - 1) + 0.0883823990145*(p - 1)*(((0.118397079113*p - 0.461201531625)/sqrt((-0.486615*p + 1.89555)^2 - 1.63574) - 0.2433075)*p - 0.2433075*p + 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)*(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)^y*log(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775) + 0.0883823990145*((-0.2433075*p + 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)*p - 1)*(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)^y*log(-0.2433075*p - 0.5*sqrt((-0.486615*p + 1.89555)^2 - 1.63574) + 0.947775)
#By using the series function, we can see the function up to a certain order. e.g. below shows series up to power of p=0 j(p,y) = (h(p,y)).series(p,3) show(i(p,y)) #For example, if we only wanted to see the coefficient of p3 (without being confused by lower powers of p) i(p,y) = (h(p,y)).series(p,3)-(h(p,y)).series(p,2) show(i(p,y)) #But we would also need to think about the powers of y 
        
\newcommand{\Bold}[1]{\mathbf{#1}}{\left(0.00534606075709 \, 0.2482448358^{y} y^{2} - 0.0741369798657 \, 0.2482448358^{y} y + 0.196232222473 \, 0.2482448358^{y}\right)} p^{2}
\newcommand{\Bold}[1]{\mathbf{#1}}{\left(0.00534606075709 \, 0.2482448358^{y} y^{2} - 0.0741369798657 \, 0.2482448358^{y} y + 0.196232222473 \, 0.2482448358^{y}\right)} p^{2}
simplify(i(p,y)) # To see if this is any clearer - but it's not! Can we also find specific powers of y??! 
        
(0.00534606075709*0.2482448358^y*y^2 - 0.0741369798657*0.2482448358^y*y
+ 0.196232222473*0.2482448358^y)*p^2
(0.00534606075709*0.2482448358^y*y^2 - 0.0741369798657*0.2482448358^y*y + 0.196232222473*0.2482448358^y)*p^2
show(f(p,y)) # Reminder! 
        
\newcommand{\Bold}[1]{\mathbf{#1}}0.0883823990144534 \, {\left(p - 1\right)} {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{y}
\newcommand{\Bold}[1]{\mathbf{#1}}0.0883823990144534 \, {\left(p - 1\right)} {\left({\left(-0.243307500000000 \, p + 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)} p - 1\right)} {\left(-0.243307500000000 \, p - 0.500000000000000 \, \sqrt{{\left(-0.486615000000000 \, p + 1.89555000000000\right)}^{2} - 1.63574000000000} + 0.947775000000000\right)}^{y}
f(p,y).series(p,2) # Reminder of generating function expanded up to point p(0) 
        
(0.0883823990144534*0.248244835800485^y) +
(0.0307407766651731*0.248244835800485^y*y -
0.233975181335305*0.248244835800485^y)*p + Order(p^2)
(0.0883823990144534*0.248244835800485^y) + (0.0307407766651731*0.248244835800485^y*y - 0.233975181335305*0.248244835800485^y)*p + Order(p^2)
f(0,1) #This gives the coefficent of p^0, y^1 (i.e. 0 cat A's, 1 Cat B in system?? (Check system!) If we multiply out above this IS correct! 
        
0.0219404741309960
0.0219404741309960
f(8,7) # This also seems to work for higher powers of y and p - haven't checked here! # Is there some way of carrying this function out like a macro, and exporting results?! 
        
94.4669253808455
94.4669253808455