How can I add two differential equations to the system in a given time interval?

1 次查看（过去 30 天）
Jack95
Jack95 2023-2-1
评论： Jack95 2023-2-3
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
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
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!

请先登录，再回答此问题。

采纳的回答

Alan Stevens
Alan Stevens 2023-2-1
Ran in:
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 个评论
显示 5更早的评论
Alan Stevens
Alan Stevens 2023-2-3
Ran in:
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Lighting, Transparency, and Shading 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by