Matlab, Fel i koden?

Jag ska använda Gillespies algoritm för att simulera en matematisk modell.

Följande kod har jag gjort:



%Begynnelsevärden
DA0 = 1;
DR0 = 1;
DpA0 = 0;
DpR0 = 0;
MA0 = 0;
MR0 = 0;
A0 = 0;
R0 = 0;
C0 = 0;

%Begynnelsevärden sätts in i en vektor
x = [A0 C0 DA0 DpA0 DR0 DpR0 MA0 MR0 R0];
t = 0;
slut = 200;
tid =t:0.5:slut;
Aresult = zeros(1,length(tid));
Rresult = zeros(1,length(tid));

i = 1;
Aresult(i) = x(1);
Rresult(i) = x(9);
tid(i) = t;

nr = nr_vilar;



while t<slut
w = prop_vilar(x); Här är det något fel
a0 = sum(w);
u1 = rand;
tao = -(log(1-u1))/a0;
P = w./a0;
F = cumsum(P);
u2 =rand;
r = find(F>=u2,1,'first');
x = x + nr(r,;
t = t+tao;
if t>tid(i)
i = i +1;
tid(i) = t;
Aresult(i) = x(1);
Rresult(i) = x(9);
end
end

plot(tid,Aresult,tid,Rresult)
title('Cirkadiska rytmer och stokastiska modeller')
legend('A','R')
xlabel('tid/h')
ylabel('mängd/antal')



Men jag får felmeddelandet "Error in (line 32) w = prop_vilar(x);"

nr_vilar och prop_vilar (stökiometrimatris och propensiteter) var angivna som m-filer från början.


De ser ut på följande sätt:

Se fortsättning på nästa post...

...Fortsättning från föregående post..

nr_vilar filen

function nr = nr_vilar( )% nr = nr_vilar()

% Stoichiometry matrix, nr, for the vilar oscillator.
The variables (corresponding to the columns in nr) are ordered as:
A C D_A D_A' D_R D_R' M_A M_R R



nr = zeros(18,9);

nr(1,: ) = [ 0 0 1 -1 0 0 0 0 0];
nr(2,: ) = [-1 0 -1 1 0 0 0 0 0];
nr(3,: ) = [1 0 0 0 1 -1 0 0 0];
nr(4,: ) = [-1 0 0 0 -1 1 0 0 0];
nr(5,: ) = [0 0 0 0 0 0 0 1 0];
nr(6,: ) = [0 0 0 0 0 0 0 1 0];
nr(7,: ) = [0 0 0 0 0 0 0 -1 0];
nr(8,: ) = [0 0 0 0 0 0 1 0 0];

osv



prop_vilar filen

function w = prop_vilar(u, p)
%
% w = prop_vilar(u, p)
% Propensities, w, for the Vilar oscillator.
%
% Input: u - the current state.
% p - list of parameters
%
% The current state variables (u) are ordered as:
% A C D_A D_A' D_R D_R' M_A M_R R
% The parameters (in p) are ordered as:
% alfa_A alfa'_A alfa_R alfa'_R beta_A beta_R teta_A teta_R ...
% gamma_A gamma_R gamma_C delta_M_R delta_M_A delta_A delta_R
%

alfaA = p(1); alfapA = p(2);
alfaR = p(3); alfapR = p(4);
betaA = p(5); betaR = p(6);
tetaA = p(7); tetaR = p(8);
gammaA = p(9); gammaR = p(10); gammaC = p(11);
deltaMR = p(12); deltaMA = p(13); deltaA = p(14);
deltaR = p(15);

w = zeros(18,1);
w(1 ) = tetaA*u(4);
w(2 ) = gammaA*u(1)*u(3);
w(3 ) = tetaR*u(6);
w(4 ) = gammaR*u(5)*u(1);
w(5 ) = alfapR*u(6);
w(6 ) = alfaR*u(5);
w(7 ) = deltaMR*u(8);
w(8 ) = alfapA*u(4);
w(9 ) = alfaA*u(3);
w(10) = deltaMA*u(7);
w(11) = betaR*u(8);
w(12) = deltaR*u(9);
w(13) = deltaA*u(2);
w(14) = betaA*u(7);
w(15) = tetaA*u(4);
w(16) = tetaR*u(6);
w(17) = deltaA*u(1);
w(18) = gammaC*u(1)*u(9);





Tack på förhand

Enligt koden för prop_vilar så tar den två argument, du anropar bara med ett argument:

% w = prop_vilar(u, p)
% Propensities, w, for the Vilar oscillator.
%
% Input: u - the current state.
% p - list of parameters
%

alltså förväntar sig prop_vilar att få u och p medan du bara ger x som ser ut att vara parametrarna. Det är som om en matematisk funktion vara:

f(x,y) = x+y och så försöker du anropa f(3) till exempel men det kan programmet/funktionen inte använda eftersom den behöver två variabler.

Men vad ska jag stoppa in som u och p egentligen?

Det står i filen:

function w = prop_vilar(u, p)
%
% w = prop_vilar(u, p)
% Propensities, w, for the Vilar oscillator.
%
% Input: u - the current state.
% p - list of parameters
%
% The current state variables (u) are ordered as:
% A C D_A D_A' D_R D_R' M_A M_R R
% The parameters (in p) are ordered as:
% alfa_A alfa'_A alfa_R alfa'_R beta_A beta_R teta_A teta_R ...
% gamma_A gamma_R gamma_C delta_M_R delta_M_A delta_A delta_R
%

Så först current state u och sen parametrarna p.

Men jag förstår fortfarande inte vad jag ska stoppa in som u och p?