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

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!

Answer 1

[Alan Stevens]

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

[Jack95]

Thank you so much it looks perfect! but what if I wanted the hydrogen to be 300 again at different times (for example at day 15 and 30)? because in my experiments I supplied hydrogen periodically. Thank you so much!

[Alan Stevens]

Just use a more complicated if statement.

[Jack95]

Hi Alan, I'm still trying but putting in "if" t=20 C(2)=C(0) doesn't work. Even if I write C(2)=300 it doesn't work. Can you give me one more suggestion? or maybe if you can show me an example,

Thanks again

[Alan Stevens]

Don't test for exact equality with floating point numbers. Try:

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;

    if abs(t-20)<10^-6
        C(2) = 300;
    end

    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

[Jack95]

Hi Alan, I tried this code, but the new hydrogen peak at 300 on day 20 does not appear. How can I do it?

[Alan Stevens]

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

[Jack95]

wonderful! You've been a great help, thank you so much!