Modèle Gillespie - Problème de plot

Fermé
Thom_Lart - 25 sept. 2016 à 18:21
Bonjour,

J'ai un projet à rendre qui s'inspire du modèle gillespie appliqué aux horloges circadiennes.
Je dois développer sous MatLab est mon souci est que j'ai dû utiliser cet outil deux fois dans ma vie.
Dans mon cas je souhaite juste imprimer a en fonction du temps mais cela ne me retourne rien.
Si quelqu'un pouvait m'expliquer pourquoi ce serait sympa. Et désolé si la question paraît vraiment bête :
Voici le code :
clear all;

% Set initial state
a = 1;
b = 1;
c = 0;
d = 0;
e = 0;
t = 0;
tend = 1000;

%Mesoscopic :
r(1) = 4; %A->ARNA
r(2) = 1; %ARNA->ARNA + A
r(3) = 1; %ARNA->
r(4) = 0.1; %A->
r(5) = 100; %A + R->AR (AR->A + R)
r(6) = 0.01; %R->
r(7) = 0.1; %RRNA->R
r(8) = 0.02; %RRNA->
r(9) = 0.001; %A->RRNA

while t<tend
u1 = r(1)*a;
u2 = r(2)*c;
u3 = r(3)*c;
u4 = r(4)*a;
u5 = r(5)*a*b;
u6 = r(6)*b;
u7 = r(7)*d;
u8 = r(8)*d;
u9 = r(9)*a;

lambda = u1+u2+u3+u4+u5+u6+u7+u8+u9;
Y = rand;
tau = -log(Y)/lambda;
t=t+tau;

x1 = u1/lambda;
x2 = (u1+u2)/lambda;
x3 = (u1+u2+u3)/lambda;
x4 = (u1+u2+u3+u4)/lambda;
x5 = (u1+u2+u3+u4+u5)/lambda;
x6 = (u1+u2+u3+u4+u5+u6)/lambda;
x7 = (u1+u2+u3+u4+u5+u6+u7)/lambda;
x8 = (u1+u2+u3+u4+u5+u6+u7+u8)/lambda;

z = rand;

if z<x1
a = a - 1;
c = c +1;

elseif z<x2
a = a + 1;

elseif z<x3
c = c - 1;

elseif z<x4
a = a - 1;

elseif z<x5
a = a - 1;
b = b - 1;
c = c + 1;

elseif z<x6
b = b-1;

elseif z<x7
b = b + 1;
d = d - 1;

elseif z<x8
d = d - 1;

else
a = a - 1;
d = d-1;

end
end

plot(t,a);
xlabel('Time');
ylabel('Number of particles A');
plot(b,t);