Help with differential equation Kolmogorov in queuing theory

Hi everybody,

I have a problem with differential equation Kolmogorov in queuing theory. Now I need to write a code to solve the equation system (1) with the condition (2) and (3) according to Euler method. I wrote the code but the output didn't seem right (please see images). Any one can see my code and give me some advices. Thanks so much

clear all
tspan = [0:0.01:5];
n=6;
ic=zeros(1,n+1);
for ii=1:n+1
    ic(1,1)=1;
end
[t, p] = ode45(@odeFun, tspan, ic);
figure ('name','xac suat theo time')
plot(t, p)
legend({'p7'})

function dpdt = odeFun(t, p)
lambda = 4.8;
mu = 2;
A=[-1 0 0 0 0 0 0;
  1 -1 0 0 0 0 0;
  0 1 -1 0 0 0 0;
  0 0 1 -1 0 0 0;
  0 0 0 1 -1 0 0;
  0 0 0 0 1 -1 0;
  0 0 0 0 0 1 0];

B=[0 1 0 0 0 0 0;
   0 -1 0 0 0 0 0;
   0 0 -2 2 0 0 0;
   0 0 0 -2 2 0 0;
   0 0 0 0 -2 2 0;
   0 0 0 0 0 -2 2;
   0 0 0 0 0 0 -2];

dpdt = (lambda.*A + mu.*B)*p;

end

.

Answer 1

[Alan Stevens]

Your code doesn't maintain condition (2) for all times. You should eliminate p6 from equations (1) using condition (2), then use the ode solver to solve for the others. You can subsequently calculate p6 for each timestep.

[Le Duc Long]

Thank you very much. But when replacing p6 at (1) by condition (2), the system (1) will no longer be the first order differential system? I am a new member so can you explain more clearly?

[Le Duc Long]

I mean, can you show me new code when replacing p6 by condition (2). Thanks so much. Best regards!

[Alan Stevens]

You need to replace p6 in the equation for dp5/dt by 1-(p0+p1+...+p5). Eliminate dp6/dt and modify A and B accordingly (so they are 6x6 rather than 7x7 matrices).

[Alan Stevens]

Try it yourself first.

[Le Duc Long]

Thanks so much. I will try to do it first. After that if i have problem, please help me. Best regards!

[Le Duc Long]

Hi Alan Stevens

If replace p6=1-(p0+p1+p2+p3+p4+p5) for dp5/dt I will get the non-homogeneous system of differential equations?

[Alan Stevens]

Try this:

tspan = 0:0.01:5;
n=6;
ic=zeros(1,n);
for ii=1:n
    ic(1,1)=1;
end
[t, p] = ode45(@odeFun, tspan, ic);
p6 = (1 - sum(p,2));
figure ('name','xac suat theo time')
plot(t, p)
legend('p0','p1','p2','p3','p4','p5')

function dpdt = odeFun(~, p)
lambda = 4.8;
mu = 2;
A=[-1 0 0 0 0 0 ;
  1 -1 0 0 0 0 ;
  0 1 -1 0 0 0 ;
  0 0 1 -1 0 0 ;
  0 0 0 1 -1 0 ;
  0 0 0 0 1 -1 ];
  

B=[0 1 0 0 0 0 ;
   0 -1 2 0 0 0 ;
   0 0 -2 2 0 0 ;
   0 0 0 -2 2 0 ;
   0 0 0 0 -2 2 ;
   -2 -2 -2 -2 -2 -4] ;
v = [0 0 0 0 0 2]';
   

dpdt = (lambda.*A + mu.*B + mu*v)*p;

end

Notice that the book plots only p0 to p5. It's still not identical to the book! I don't know if this is because the book used a simple Euler method (as you mentioned in your original post) and MATLAB is using a more accurate Runge-Kutta method, or not. I'll leave you to investigate!

[Le Duc Long]

Thanks for your help. I do not know why the condition p0(t)+p1(t)+...+pn(t) is diffrent 1. The condition (2) is true with any time value. Can you explain more? Best regards!