How can I add two differential equations to the system in a given time interval?
1 次查看(过去 30 天)
显示 更早的评论
Hello everyone! I am i'm modeling a process in which we initially have 3 equations: acetogen biomass production C(1), hydrogen substrate consumption C(2), acetate production C(3).
The code is as follows:
Tspan = [0:30];
C = [50 300 100];
Miumax=0.07;
Ks=10;
Yxs=0.1;
alfa=1-Yxs;
betha=0;
[t,C]=ode45(@(t,C)ParameterJack(t,C,Miumax,Ks,Yxs,alfa,betha),Tspan,C);
figure
hold on
yyaxis left
plot(t,C(:,1),'-g',t,C(:,2),'-b')
ylabel('acetogens (mgCOD/L), hidrogen (mgCOD/L)')
ylim([min(C(:,1)) max(ylim)])
ylim([min(C(:,2)) max(ylim)])
yyaxis right
plot(t,C(:,3));
ylabel('acetate (mgCOD/L)')
ylim([min(C(:,3)) max(ylim)])
xlabel('tempo (giorni)')
grid
legend('acetogeni','idrogeno','acetato', 'Location','best')
function dCndt=ParameterJack(t,C,Miumax,Ks,Yxs,alfa,betha)
mu = (Miumax*C(2))/(Ks+C(2));
alfa=1-Yxs;
Miumax=0.07;
Ks=10;
Yxs=0.175
dCndt(1,:) = mu*C(1);
dCndt(2,:) = -C(1)*(mu/Yxs);
dCndt(3,:) = (C(1)*((alfa*mu)+betha))
end
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1280420/image.jpeg)
after 20 days I have that competing organisms (methanogens) C(4) are activated and produce methane (5), and the code becomes:
Tspan = [0:30];
C0 = [50 300 100 10 0];
Miumax=0.07;
Ks=10;
Yxs=0.1;
alfa=9;
betha=0;
y=0.2
km=100;
ka=3
a2=0.8
a1=0.559;
u=0.1
[t,C]=ode45(@(t,C)ParameterJack(t,C,Miumax,Ks,km,ka,a2,a1,Yxs,u,y,alfa,betha),Tspan,C0);
figure
hold on
yyaxis left
plot(t,C(:,1),'-g',t,C(:,2),'-b',t,C(:,4),'-m');
ylabel('acetogens (mgCOD/L), hidrogen (mgCOD/L),metanogens')
ylim([min(C(:,1)) max(ylim)])
ylim([min(C(:,2)) max(ylim)])
yyaxis right
plot(t,C(:,3),'r', t,C(:,5),'-o');
ylabel('acetate (mgCOD/L), methane (mgCOD/L)')
ylim([min(C(:,3)) max(ylim)])
hold off
xlabel('tempo (giorni)')
grid
legend('acetogens','Hidrogen','metanogens','acetate','methane', 'Location','best')
function dCdt=ParameterJack(t,C,Miumax,Ks,km,u,ka,y,a2,a1,Yxs,alfa,betha)
mum = (0.1)*(C(2))/(km+C(2));
mu = (Miumax*C(2))/(Ks+C(2));
alfa=9;
Miumax=0.07;
u=0.1;
ka=50;
Ks=10;
Yxs=0.175
y=0.2;
km=100;
a2=0.0121
a1=0.159
dCdt(1,:) = mu*C(1);
dCdt(2,:) = -C(1)*(mu/Yxs)-(C(4)*(mum/y));
dCdt(3,:) = (C(1)*((alfa*mu)+betha))-C(3)*a2*C(4);
dCdt(4,:) = mum*C(4)+(u)*((C(3)*(1-a2)/((C(3)*(1-a2))+ka))*C(4));
dCdt(5,:) = (C(3)*a2)+(C(2)*a1);
end
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1280425/image.jpeg)
I would like to understand how I can represent this thing in the same graph, where the growth of methanogens and the production of methane starts on day 20 and not on day 0.
Thanks everyone for the help!
0 个评论
采纳的回答
Alan Stevens
2023-2-1
Like this? (though I'm not sure I've interpreted which constants apply during which interval correctly!)
Tspan = 0:30;
C0 = [50 300 100 10 0];
Miumax=0.07;
Ks=10;
Yxs=0.1;
alfa=9;
betha=0;
y=0.2;
km=100;
ka=3;
a2=0.8;
a1=0.559;
u=0.1;
[t,C]=ode45(@(t,C)ParameterJack(t,C,betha),Tspan,C0);
figure
hold on
yyaxis left
plot(t,C(:,1),'-g',t,C(:,2),'-b',t,C(:,4),'-m');
ylabel('acetogens (mgCOD/L), hidrogen (mgCOD/L),metanogens')
ylim([min(C(:,1)) max(ylim)])
ylim([min(C(:,2)) max(ylim)])
yyaxis right
plot(t,C(:,3),'r', t,C(:,5),'-o');
ylabel('acetate (mgCOD/L), methane (mgCOD/L)')
ylim([0 max(ylim)])
hold off
xlabel('tempo (giorni)')
grid
legend('acetogens','Hidrogen','metanogens','acetate','methane', 'Location','best')
function dCdt=ParameterJack(t,C,betha)
Miumax=0.07;
alfa=9;
u=0.1;
Ks = 10;
y=0.2;
km=100;
Yxs = 0.175;
mum = (0.1)*(C(2))/(km+C(2));
mu = (Miumax*C(2))/(Ks+C(2));
if t<20
term2 = 0;
term3 = 0;
term4 = 0;
term5 = 0;
else
ka = 50;
a2=0.0121;
a1=0.159;
term2 = (C(4)*(mum/y));
term3 = C(3)*a2*C(4);
term4 = mum*C(4)+(u)*((C(3)*(1-a2)/((C(3)*(1-a2))+ka))*C(4));
term5 = (C(3)*a2)+(C(2)*a1);
end
dCdt(1,:) = mu*C(1);
dCdt(2,:) = -C(1)*(mu/Yxs)-term2;
dCdt(3,:) = (C(1)*((alfa*mu)+betha))-term3;
dCdt(4,:) = term4;
dCdt(5,:) = term5;
end
7 个评论
Alan Stevens
2023-2-3
Sorry, i did it too quickly, without thinking! Since you want to change things at times, then the following should work:
Tspan1 = 0:20;
Tspan2 = 20:30;
C01 = [50 300 100 10 0];
Miumax=0.07;
Ks=10;
Yxs=0.1;
alfa=9;
betha=0;
y=0.2;
km=100;
ka=3;
a2=0.8;
a1=0.559;
u=0.1;
[t1,C1] = ode45(@(t,C)ParameterJack(t,C,betha),Tspan1,C01);
C02 = C1(end,:);
C02(1,2) = 300;
[t2,C2] = ode45(@(t,C)ParameterJack(t,C,betha),Tspan2,C02);
t = [t1; t2];
C = [C1; C2];
figure
hold on
yyaxis left
plot(t,C(:,1),'-g',t,C(:,2),'-b',t,C(:,4),'-m');
ylabel('acetogens (mgCOD/L), hidrogen (mgCOD/L),metanogens')
ylim([min(C(:,1)) max(ylim)])
ylim([min(C(:,2)) max(ylim)])
yyaxis right
plot(t,C(:,3),'r', t,C(:,5),'-o');
ylabel('acetate (mgCOD/L), methane (mgCOD/L)')
ylim([0 max(ylim)])
hold off
xlabel('tempo (giorni)')
grid
legend('acetogens','Hidrogen','metanogens','acetate','methane', 'Location','best')
function dCdt=ParameterJack(t,C,betha)
Miumax=0.07;
alfa=9;
u=0.1;
Ks = 10;
y=0.2;
km=100;
Yxs = 0.175;
mum = (0.1)*(C(2))/(km+C(2));
mu = (Miumax*C(2))/(Ks+C(2));
if t<20
term2 = 0;
term3 = 0;
term4 = 0;
term5 = 0;
else
ka = 50;
a2=0.0121;
a1=0.159;
term2 = (C(4)*(mum/y));
term3 = C(3)*a2*C(4);
term4 = mum*C(4)+(u)*((C(3)*(1-a2)/((C(3)*(1-a2))+ka))*C(4));
term5 = (C(3)*a2)+(C(2)*a1);
end
dCdt(1,:) = mu*C(1);
dCdt(2,:) = -C(1)*(mu/Yxs)-term2;
dCdt(3,:) = (C(1)*((alfa*mu)+betha))-term3;
dCdt(4,:) = term4;
dCdt(5,:) = term5;
end
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!