Trouble solving ODE equations
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
I am trying to solve these differential equations in order to plot graphs, however, when running the code I am receiving multiple errors but I'm unsure why as I can't see why the code is not running smoothly? I'm looking for some insight on why this might be happening and how I can solve the problem?
Equations and Parameter values:
function dydt = Dunster(t,y)
%Parameters and their values%
K1a = 20;
Y1a = 1;
K1b = 0.25;
K2a = 14;
K2am = 72;
K2b = 2.6;
K2bm = 10;
K2c = 24;
K2cm = 20;
K3a = 10;
K3b = 0.05;
K3c = 24;
K3cm = 20;
K4a = 2.3;
K4am = 58;
K4b = 2000;
K4bm = 210;
K4c = 1.3;
K5a = 0.0014;
K5b = 0.35;
K6 = 2000;
K6m = 2000;
Kx = 50;
Kb = 0.01;
Lb = 1;
Lx = 1;
Bxva = 0;
H = K1a*Y1a*exp(-Y1a*t);
bl = 0.5*((Kb+Lb+Bxva)-sqrt((Kb+Lb+Bxva)^2 - 4*Lb*Bxva));
%Equation 1%
%Generation and inactivation of factor Xa%
F1= y(1);
%y(1)= -K6*F1;
% Factor I
%Equation 2%
F5=y(2);
%y(2)=(-K2a*F2a*F5)/(F5+K2am(1+F1/K6m) - ((K2b*F10a*F5)/F5+K2bm(1+F2/K4am)));
%Equation 3%
F2=y(3);
%y(3)= ((-K4a*F1*F2)/F2+K4am(1+F5/K2bm))-((K4b*b1*F2)/(F2+K4bm));
%Equation 4%
C=y(4);
%y(4)= (-K5a*C);
%Equation 5%
F10a=y(5);
%y(5)= H +(K3c*Ca*Bxva)/(Bxva+K3cm) - K1b*F10a - K3a*F10a*F5a;
%Equation 6%
F5a=y(6);
%y(6)= (K2a*F2a*F5)/(F5+K2am(1+F1/K6m)) + (K2b*F10a*F5a)/(F5+K2bm(1+F1/K4am))+ K3b*Bxva - (K2c*Ca*F5a)/(F5a+K2am) - K3a*F10a*F5a;
%Equation 7%
Bxva=y(7);
%y(7)= K3a*F10a*F5a - K3b*Bxva - (K3c*Ca*Bxva)/(Bxva+K3cm);
%Equation 8%
F2a=y(8);
%y(8)= (K4a*F1*F2)/(F2+K4am(1+F5/K2bm)) + (K4b*BF2)/(F2+K4bm) - K4c*F2a;
%Equation 9%
Ca=y(9);
%y(9)= K5a*C - K5b*Ca;
%Equation 10%
F10i=y(10);
%y(10)= K1b*F10a + K3b*Bxva;
%Equation 11%
F5i=y(11);
%y(11) = (K2c*Ca*F5a)/(F5a+K2cm) + (K3c*Ca*Bxva/(Bxva+K3cm));
%Equation 12%
F2i=y(12);
%y(12) = K4c*F2a;
%Equation 13%
Ci=y(13);
%y(13) = K5b*Ca;
%Equation 14%
F1a=y(14);
%y(14) = K6*F1;
dydt=[
-K6*F1; % d F1 / dt
((-K2a*F2a*F5)/(F5+K2am*(1+F1/K6m)) - ((K2b*F10a*F5)/F5+K2bm*(1+F2/K4am))); % d F5 / dt
((-K4a*F1*F2)/(F2+K4am*(1+F5/K2bm))-((K4b*bl*F2)/(F2+K4bm))); % d F2 / dt
-K5a*C; % d C / dt
H+(K3c*Ca*Bxva)/(Bxva+K3cm) - K1b*F10a - K3a*F10a*F5a ;
% d F10a / dt
(K2a*F2a*F5)/(F5+K2am(1+F1/K6m)) + (K2b*F10a*F5a)/(F5+K2bm(1+F1/K4am))
+ K3b*Bxva - (K2c*Ca*F5a)/(F5a+K2am) - K3a*F10a*F5a; % d F5a / dt
K3a*F10a*F5a - K3b*Bxva - (K3c*Ca*Bxva)/(Bxva+K3cm); % d Bxva / dt
(K4a*F1*F2)/(F2+K4am(1+F5/K2bm)) + (K4b*BF2)/(F2+K4bm) - K4c*F2a ;
% d F2a / dt
K5a*C - K5b*Ca ; % d Ca / dt
K1b*F10a + K3b*Bxva ; % d F10i / dt
((K2c*Ca*F5a)/(F5a+K2cm) + (K3c*Ca*Bxva/(Bxva+K3cm))) ; % d F5i / dt
K4c*F2a ; % d F2i / dt
K5b*Ca ; % d Ci / dt
K6*F1 ; % d F1a / dt
];
end
Running Part of the Code:
tspan=[0 60];
F1=10000;
F5=30;
F2=1000;
C=100;
F10a=0;
F5a=0;
Bxva=0;
F2a=0;
Ca=0;
F10i=0;
F5i=0;
F2i=0;
Ci=0;
F1a=0;
y0=[F1;F5;F2;C;F10a;F5a;Bxva;F2a;Ca;F10i;F5i;F2i;Ci;F1a];
[t,y]=ode45(@Dunster, tspan, y0);%,[], pars);
plot(t/60, y(:,14)*1000, "r")
xlabel("Time [Min]")
ylabel("Concentration [nM]")
Risposte (1)
Star Strider
il 11 Apr 2021
There were several missing multiplication operators, and a reference to ‘BF2’ that I corrected to ‘F2’, since there is no ‘B’ that I can find, so no missing multiplication operator there. With those edits, and concatenating an additional 0 to ‘y0’ to make its length equal to the number of differential equations, this runs without error:
function dydt = Dunster(t,y)
%Parameters and their values%
K1a = 20;
Y1a = 1;
K1b = 0.25;
K2a = 14;
K2am = 72;
K2b = 2.6;
K2bm = 10;
K2c = 24;
K2cm = 20;
K3a = 10;
K3b = 0.05;
K3c = 24;
K3cm = 20;
K4a = 2.3;
K4am = 58;
K4b = 2000;
K4bm = 210;
K4c = 1.3;
K5a = 0.0014;
K5b = 0.35;
K6 = 2000;
K6m = 2000;
Kx = 50;
Kb = 0.01;
Lb = 1;
Lx = 1;
Bxva = 0;
H = K1a*Y1a*exp(-Y1a*t);
bl = 0.5*((Kb+Lb+Bxva)-sqrt((Kb+Lb+Bxva)^2 - 4*Lb*Bxva));
%Equation 1%
%Generation and inactivation of factor Xa%
F1= y(1);
%y(1)= -K6*F1;
% Factor I
%Equation 2%
F5=y(2);
%y(2)=(-K2a*F2a*F5)/(F5+K2am(1+F1/K6m) - ((K2b*F10a*F5)/F5+K2bm(1+F2/K4am)));
%Equation 3%
F2=y(3);
%y(3)= ((-K4a*F1*F2)/F2+K4am(1+F5/K2bm))-((K4b*b1*F2)/(F2+K4bm));
%Equation 4%
C=y(4);
%y(4)= (-K5a*C);
%Equation 5%
F10a=y(5);
%y(5)= H +(K3c*Ca*Bxva)/(Bxva+K3cm) - K1b*F10a - K3a*F10a*F5a;
%Equation 6%
F5a=y(6);
%y(6)= (K2a*F2a*F5)/(F5+K2am(1+F1/K6m)) + (K2b*F10a*F5a)/(F5+K2bm(1+F1/K4am))+ K3b*Bxva - (K2c*Ca*F5a)/(F5a+K2am) - K3a*F10a*F5a;
%Equation 7%
Bxva=y(7);
%y(7)= K3a*F10a*F5a - K3b*Bxva - (K3c*Ca*Bxva)/(Bxva+K3cm);
%Equation 8%
F2a=y(8);
%y(8)= (K4a*F1*F2)/(F2+K4am(1+F5/K2bm)) + (K4b*BF2)/(F2+K4bm) - K4c*F2a;
%Equation 9%
Ca=y(9);
%y(9)= K5a*C - K5b*Ca;
%Equation 10%
F10i=y(10);
%y(10)= K1b*F10a + K3b*Bxva;
%Equation 11%
F5i=y(11);
%y(11) = (K2c*Ca*F5a)/(F5a+K2cm) + (K3c*Ca*Bxva/(Bxva+K3cm));
%Equation 12%
F2i=y(12);
%y(12) = K4c*F2a;
%Equation 13%
Ci=y(13);
%y(13) = K5b*Ca;
%Equation 14%
F1a=y(14);
%y(14) = K6*F1;
dydt=[-K6*F1; % d F1 / dt
((-K2a*F2a*F5)/(F5+K2am*(1+F1/K6m)) - ((K2b*F10a*F5)/F5+K2bm*(1+F2/K4am))); % d F5 / dt
((-K4a*F1*F2)/(F2+K4am*(1+F5/K2bm))-((K4b*bl*F2)/(F2+K4bm))); % d F2 / dt
-K5a*C; % d C / dt
H+(K3c*Ca*Bxva)/(Bxva+K3cm) - K1b*F10a - K3a*F10a*F5a ;
% d F10a / dt
(K2a*F2a*F5)/(F5+K2am*(1+F1/K6m)) + (K2b*F10a*F5a)/(F5+K2bm*(1+F1/K4am))
+ K3b*Bxva - (K2c*Ca*F5a)/(F5a+K2am) - K3a*F10a*F5a; % d F5a / dt
K3a*F10a*F5a - K3b*Bxva - (K3c*Ca*Bxva)/(Bxva+K3cm); % d Bxva / dt
(K4a*F1*F2)/(F2+K4am*(1+F5/K2bm)) + (K4b*F2)/(F2+K4bm) - K4c*F2a ;
% d F2a / dt
K5a*C - K5b*Ca ; % d Ca / dt
K1b*F10a + K3b*Bxva ; % d F10i / dt
((K2c*Ca*F5a)/(F5a+K2cm) + (K3c*Ca*Bxva/(Bxva+K3cm))) ; % d F5i / dt
K4c*F2a ; % d F2i / dt
K5b*Ca ; % d Ci / dt
K6*F1 ; % d F1a / dt
];
end
% % Running Part of the Code:
tspan=[0 60];
F1=10000;
F5=30;
F2=1000;
C=100;
F10a=0;
F5a=0;
Bxva=0;
F2a=0;
Ca=0;
F10i=0;
F5i=0;
F2i=0;
Ci=0;
F1a=0;
y0=[F1;F5;F2;C;F10a;F5a;Bxva;F2a;Ca;F10i;F5i;F2i;Ci;F1a;0];
[t,y]=ode45(@Dunster, tspan, y0);%,[], pars);
plot(t/60, y(:,14)*1000, "r")
xlabel("Time [Min]")
ylabel("Concentration [nM]")
.
6 Commenti
Star Strider
il 11 Apr 2021
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!