Proper way to write ode45

I am having some problems with my results using ode 45 and i dont know if i am using the right way to run the function.I send to you a few lines of the functions and the scripts so that you can get the general idea.I am using so many inputs to the function to avoid using global and if you know any other way to make it run faster i would be very interested in hearing it.

so the function is

• function [ dydt ] = adm1_no_w5 (t,y,Ssu_in,Saa_in,Sfa_in,Sva_in,Sbu_in,Spro_in,Sac_in,Sh2_in,Sch4_in,Sic_in,Sin_in,Si_in,Xc_in,Xch_in,Xpr_in,Xli_in,Xsu_in,Xaa_in,Xfa_in,Xc4_in,Xpro_in,Xac_in,Xh2_in,Xi_in,Scation_in,Sanion_in,kLa,T_base,T_op,Ka_va,Ka_bu,Ka_pro,Ka_ac,Ka_co2,Ka_in,Ka_Bva,Ka_Bbu,Ka_Bpro,Ka_Bac,Ka_Bco2,Ka_Bin,Patm,kp,KH_co2,KH_ch4,KH_h2,Kh_h20_base,Vgas,R,Kw,Vliq,qin,Kdis,Khyd_ch,Khyd_pr,Khyd_li,Kdec,pHul_aa,pHll_aa,Km_su,Ks_su,Ysu,Km_aa,Ks_aa,Yaa,Km_fa,Ks_fa,Yfa,Klh2_fa,Km_c4,Ks_c4,Yc4,Klh2_c4,Km_pro,Ks_pro,Ypro,Klh2_pro,Km_ac,Ks_ac,Yac,pHul_ac,pHll_ac,Klnh3,Km_h2,Ks_h2,Yh2,pHul_h2,pHll_h2,Ks_in,fsi_xc,fxi_xc,fch_xc,fpr_xc,fli_xc,Nxc,ffa_li,fh2_su,fbu_su,fpro_su,fac_su,fh2_aa,Naa,fva_aa,fbu_aa,fpro_aa,fac_aa,Ni,Nbac,C_xc,C_si,C_ch,C_pr,C_li,C_xi,C_su,C_aa,C_fa,C_bu,C_pro,C_ac,C_bac,C_va,C_ch4)

dydt(1)=(((qin/Vliq)*(Ssu_in-y(1))+r2+(1-ffa_li)*r4)-r5);

dydt(2)=((qin/Vliq)*(Saa_in-y(2)))+r3-r6;

end

the general type is

𝑑𝑆/dt= 𝑞𝑖𝑛/Vliq ∗(𝑆in −𝑆 )+Σr

And the way i am writting the ode solver is

• options = odeset('RelTol',1e-2,'AbsTol',1e-2);

tspan=0.000001:0.1:5;

[t,y]=ode45(@(t,y) adm1_no_w5 (t,y,Ssu_in,Saa_in,Sfa_in,Sva_in,Sbu_in,Spro_in,Sac_in,Sh2_in,Sch4_in,Sic_in,Sin_in,Si_in,Xc_in,Xch_in,Xpr_in,Xli_in,Xsu_in,Xaa_in,Xfa_in,Xc4_in,Xpro_in,Xac_in,Xh2_in,Xi_in,Scation_in,Sanion_in,kLa,T_base,T_op,Ka_va,Ka_bu,Ka_pro,Ka_ac,Ka_co2,Ka_in,Ka_Bva,Ka_Bbu,Ka_Bpro,Ka_Bac,Ka_Bco2,Ka_Bin,Patm,kp,KH_co2,KH_ch4,KH_h2,Kh_h20_base,Vgas,R,Kw,Vliq,qin,Kdis,Khyd_ch,Khyd_pr,Khyd_li,Kdec,pHul_aa,pHll_aa,Km_su,Ks_su,Ysu,Km_aa,Ks_aa,Yaa,Km_fa,Ks_fa,Yfa,Klh2_fa,Km_c4,Ks_c4,Yc4,Klh2_c4,Km_pro,Ks_pro,Ypro,Klh2_pro,Km_ac,Ks_ac,Yac,pHul_ac,pHll_ac,Klnh3,Km_h2,Ks_h2,Yh2,pHul_h2,pHll_h2,Ks_in,fsi_xc,fxi_xc,fch_xc,fpr_xc,fli_xc,Nxc,ffa_li,fh2_su,fbu_su,fpro_su,fac_su,fh2_aa,Naa,fva_aa,fbu_aa,fpro_aa,fac_aa,Ni,Nbac,C_xc,C_si,C_ch,C_pr,C_li,C_xi,C_su,C_aa,C_fa,C_bu,C_pro,C_ac,C_bac,C_va,C_ch4),tspan,y0,options);

Answer 1

[Torsten]

The usual way is to use arrays in which you can write and from which you can read variables at certain positions.

Or you can use structures

var.p1 = p1;

var.p2 = p2

etc

and you only pass the structure "var".

[Michael Skyvalos]

you mean istead of reading the variables from the excell sheet?

[Torsten]

I mean that it is not practical to have inputs to a function which go over several lines.

Therefore, you should not pass every variable, but save them in an array or in a structure and pass this array or this structure between the functions.

But still, I'm not sure what your question is. If you have problems with the solver, you should supply values of the relevant variables for testing.

[Michael Skyvalos]

my main question is for the diferential equation that i have which is 𝑑𝑆/dt= 𝑞𝑖𝑛/Vliq ∗(𝑆in −𝑆 )+Σr

so i wrote it in the function like this

dydt(1)=(((qin/Vliq)*(Ssu_in-y(1))+r2+(1-ffa_li)*r4)-r5);

dydt(2)=((qin/Vliq)*(Saa_in-y(2)))+r3-r6;

which is the right way to write the ode45 in the script to solve this function

I am asking because i dont know how exactly ode45 works and i am getting wrong results so i think maybe i am writting and using it wrong.

As for the variables it was a side question and thank you very much for the answer it is very helpful!

[Torsten]

If the r_i are reaction rates, they usually depend on y(1) and y(2) and have to be recalculated in each call to your function.

At the moment, you use constant values that are transfered to the function "adm1_no_w5" from the calling program.

[Michael Skyvalos]

you are so right!how do i recalculate them?i mean how do i have to write them so each time it gets the new value

[Torsten]

function [ dydt ] = adm1_no_w5(....)
  r2 = something depending on y(1) and y(2);
  r4 = something depending on y(1) and y(2);  
  %Same for the other model constants
  
  dydt = zeros(2,1);
  dydt(1) = ...;
  dydt(2) = ...;
end

[Michael Skyvalos]

i have already done that and this also solved the error for the function not returning a column vector.

Do you know if maybe the problem is in the way i call the ode

[t,dydt]=ode45(@(t,y) adm1_no_w5(t,y,....),tspan,y0,options);

am i writting this the right way.?

also in the function there is nowhere the `time factor is this maybe the problem?

[Michael Skyvalos]

Here is the whole function if it helps

function [ dydt ] = adm1_no_w5 (t,y,Ssu_in,Saa_in,Sfa_in,Sva_in,Sbu_in,Spro_in,Sac_in,Sh2_in,Sch4_in,Sic_in,Sin_in,Si_in,Xc_in,Xch_in,Xpr_in,Xli_in,Xsu_in,Xaa_in,Xfa_in,Xc4_in,Xpro_in,Xac_in,Xh2_in,Xi_in,Scation_in,Sanion_in,kLa,T_base,T_op,Ka_va,Ka_bu,Ka_pro,Ka_ac,Ka_co2,Ka_in,Ka_Bva,Ka_Bbu,Ka_Bpro,Ka_Bac,Ka_Bco2,Ka_Bin,Patm,kp,KH_co2,KH_ch4,KH_h2,Kh_h20_base,Vgas,R,Kw,Vliq,qin,Kdis,Khyd_ch,Khyd_pr,Khyd_li,Kdec,pHul_aa,pHll_aa,Km_su,Ks_su,Ysu,Km_aa,Ks_aa,Yaa,Km_fa,Ks_fa,Yfa,Klh2_fa,Km_c4,Ks_c4,Yc4,Klh2_c4,Km_pro,Ks_pro,Ypro,Klh2_pro,Km_ac,Ks_ac,Yac,pHul_ac,pHll_ac,Klnh3,Km_h2,Ks_h2,Yh2,pHul_h2,pHll_h2,Ks_in,fsi_xc,fxi_xc,fch_xc,fpr_xc,fli_xc,Nxc,ffa_li,fh2_su,fbu_su,fpro_su,fac_su,fh2_aa,Naa,fva_aa,fbu_aa,fpro_aa,fac_aa,Ni,Nbac,C_xc,C_si,C_ch,C_pr,C_li,C_xi,C_su,C_aa,C_fa,C_bu,C_pro,C_ac,C_bac,C_va,C_ch4)

dydt=zeros(length(y),1);

fprintf('t=%f\n',t)

theta=y(25)+y(11)-y(32) -y(31) -(y(30) / 64)-(y(29)/112)-(y(28)/160)-(y(27) /208)-y(26) ; %ola ta S excell

Shcation=(-theta*0.5)+(0.5*(sqrt((theta^2)+(4*Kw))));

pHLim_aa=10^(-(pHul_aa+pHll_aa)/2);

pHLim_ac=10^(-(pHul_ac + pHll_ac)/2);

pHLim_h2=10^(-(pHul_h2 + pHll_h2)/2);

k_aa=3/(pHul_aa - pHll_aa);

k_ac=3/(pHul_ac - pHll_ac);

k_h2=3/(pHul_h2 - pHll_h2);

IpH_aa=(pHLim_aa^k_aa)/((Shcation^k_aa)+(pHLim_aa^k_aa));

IpH_ac=(pHLim_ac^k_ac)/((Shcation^k_ac)+(pHLim_ac^k_ac));

IpH_h2=(pHLim_h2^k_h2)/((Shcation^k_h2)+(pHLim_h2^k_h2));

param1=-C_xc+(fsi_xc*C_si)+(fch_xc*C_ch)+(fpr_xc*C_pr)+(fli_xc*C_li)+(fxi_xc*C_xi);

param2=-C_ch+C_su;

param3=-C_pr+C_aa;

param4=-C_li+((1-ffa_li)*C_su)+(ffa_li*C_fa);

param5=-C_su+((1-Ysu)*((fbu_su*C_bu)+(fpro_su*C_pro)+(fac_su*C_ac)))+(Ysu*C_bac);

param6=-C_aa+((1-Yaa)*((fva_aa*C_va)+(fbu_aa*C_bu)+(fpro_aa*C_pro)+(fac_aa*C_ac)))+(Yaa*C_bac);

param7=-C_fa+((1-Yfa)*0.7*C_ac)+(Yfa*C_bac);

param8=-C_va+((1-Yc4)*0.54*C_pro)+((1-Yc4)*0.31*C_ac)+(Yc4*C_bac);

param9=-C_bu+((1-Yc4)*0.8*C_ac)+(Yc4*C_bac);

param10=-C_pro+((1-Ypro)*0.57*C_ac)+(Ypro*C_bac);

param11=-C_ac+((1-Yac)*C_ch4)+(Yac*C_bac);

param12=((1-Yh2)*C_ch4)+(Yh2*C_bac);

param13=-C_bac+C_xc;

Iin_lim=1/(1+(Ks_in/y(11))) ;

Ih2_fa=1/(1+(y(8)/Klh2_fa)) ;

Ih2_pro=1/(1+(y(8)/Klh2_pro));

Ih2_c4=1/(1+(y(8)/Klh2_c4));

Inh3=1/(1+(y(32)/Klnh3));

I1=IpH_aa*Iin_lim;

I2=I1*Ih2_fa;

I3=I1*Ih2_c4;

I4=I1*Ih2_pro;

I5=IpH_ac*Iin_lim*Inh3;

I6=IpH_h2*Iin_lim;

r1=Kdis*y(13);

r2=Khyd_ch*y(14);

r3=Khyd_pr*y(15);

r4=Khyd_li*y(16);

r5=Km_su*(y(1)/(Ks_su+y(1)))*y(17)*I1;

r6=Km_aa*(y(2)/(Ks_aa+y(2)))*y(18)*I1;

r7=Km_fa*(y(3)/(Ks_fa+y(3)))*y(19)*I2;

r8=Km_c4 *(y(4)/(Ks_c4+y(4)))*y(20)*(y(4)/(y(5)+y(4)+0.000001))*I3;

r9=Km_c4*(y(5)/(Ks_c4+y(5)))*y(20)*(y(5)/(y(4)+y(5)+0.000001))*I3; excell

r10=Km_pro*(y(6)/(Ks_pro+y(6)))*y(21)*I4;

r11=Km_ac*(y(7)/(Ks_ac+y(7)))*y(22)*I5 ;

r12=Km_h2*(y(8)/(Ks_h2+y(8)))*y(23)*I6;

r13=Kdec*y(17);

r14=Kdec*y(18);

r15=Kdec*y(19);

r16=Kdec*y(20);

r17=Kdec*y(21);

r18=Kdec*y(22);

r19=Kdec*y(23);

rt8=kLa*(y(8)-(16*KH_h2*(y(33)*R*T_op/16)));

rt9=kLa*(y(9)-(64*KH_ch4*(y(34)*R*T_op/64)));

rt10=kLa*(y(10)-y(31)-(KH_co2*(y(35)*R*T_op)));

%%

ra4=Ka_Bva*(y(27)*(Ka_va+Shcation)-(Ka_va*y(4)));

ra5=Ka_Bbu*(y(28)*(Ka_bu+Shcation)-(Ka_bu*y(5)));

ra6=Ka_Bpro*(y(29)*(Ka_pro+Shcation)-(Ka_pro*y(6)));

ra7=Ka_Bac*(y(30)*(Ka_ac+Shcation)-(Ka_ac*y(7)));

ra10=Ka_Bco2*(y(31)*(Ka_co2+Shcation)-(Ka_co2*y(10)));

ra11=Ka_Bin*(y(32)*(Ka_in+Shcation)-(Ka_in*y(11)));

Pgas_h2=y(33)*R*T_op/16;

Pgas_ch4=y(34)*R*T_op/64;

Pgas_co2=y(35)*R*T_op;

Pgas_h20=Kh_h20_base*exp(5290*(1/T_base-1/T_op));

Pgas=Pgas_h2+Pgas_ch4+Pgas_co2+Pgas_h20;

qgas=kp*(Pgas-Patm)*Pgas/Patm;

if qgas<0;

qgas=0;

end

dydt(1)=(((qin/Vliq)*(Ssu_in-y(1))+r2+(1-ffa_li)*r4)-r5);

dydt(2)=((qin/Vliq)*(Saa_in-y(2)))+r3-r6;

dydt(3)=((qin/Vliq)*(Sfa_in-y(3)))+(ffa_li*r4)-r7;

dydt(4)=((qin/Vliq)*(Sva_in-y(4)))+((1-Yaa)*fva_aa*r6)-r8;

dydt(5)=((qin/Vliq)*(Sbu_in-y(5)))+((1-Ysu)*fbu_su*r5)+((1-Yaa)*fbu_aa*r6)-r9;

dydt(6)=((qin/Vliq)*(Spro_in-y(6))+((1-Ysu)*fpro_su*r5)+((1-Yaa)*fpro_aa*r6)+((1-Yc4)*0.54*r8)-r10);

dydt(7)=((qin/Vliq)*(Sac_in-y(7))+((1-Ysu)*fac_su*r5)+((1-Yaa)*fac_aa*r6)+((1-Yfa)*0.7*r7)+((1-Yc4)*0.31*r8)+((1-Yc4)*0.8*r9)+((1-Ypro)*0.57*r10)-r11);

dydt(8)=((qin/ Vliq)*(Sh2_in-y(8))+((1-Ysu)*fh2_su*r5)+((1-Yaa)*fh2_aa*r6)+((1-Yfa)*0.3*r7)+((1-Yc4)*0.15*r8)+((1-Yc4)*0.2*r9)+((1-Ypro)*0.43*r10)-r12-rt8);

dydt(9)=((qin/Vliq)*(Sch4_in-y(9))+((1-Yac)*r11)+((1-Yh2)*r12))-rt9;

dydt(10)=((qin/Vliq)*(Sic_in-y(10))-((param1*r1)+(param2*r2)+(param3*r3)+(param4*r4)+(param5*r5)+(param6*r6)+(param7*r7)+(param8*r8)+(param9*r9)+(param10*r10)+(param11*r11)+(param12*r12)+(param13*r13)+(param13*r14)+(param13*r15)+(param13*r16)+(param13*r17)+(param13*r18)+(param13*r19))-(rt10));

dydt(11)=((qin/Vliq)*(Sin_in-y(11))-(Ysu*Nbac*r5)+((Naa-(Yaa*Nbac))*r6)-(Yfa*Nbac*r7)-(Yc4*Nbac*r8)-(Yc4*Nbac*r9)-(Ypro*Nbac*r10)-(Yac*Nbac*r11)-(Yh2*Nbac*r12)+((Nbac-Nxc)*(r13+r14+r15+r16+r17+r18+r19))+(Nxc-(fxi_xc*Ni)-(fsi_xc*Ni)-(fpr_xc*Naa))*r1);

dydt(12)=((qin/Vliq)*(Si_in-y(12))+(fsi_xc*r1));

dydt(13)=((qin/Vliq)*(Xc_in-y(13))-r1+r13+r14+r15+r16+r17+r18+19);

dydt(14)=((qin/Vliq)*(Xch_in-y(14))+(fch_xc*r1)-r2);

dydt(15)=((qin/Vliq)*(Xpr_in-y(15))+(fpr_xc*r1)-r3);

dydt(16)=((qin/Vliq)*(Xli_in-y(16))+(fli_xc*r1)-r4);

dydt(17)=((qin/Vliq)*(Xsu_in-y(17)))+(Ysu*r5)-r13;

dydt(18)=((qin/Vliq)*(Xaa_in-y(18))+(Yaa*r6)-r14);

dydt(19)=((qin/Vliq)*(Xfa_in-y(19))+(Yfa*r7)-r15);

dydt(20)=((qin/Vliq)*(Xc4_in-y(20))+(Yc4*r8)+(Yc4*r9)-r16);

dydt(21)=((qin/Vliq)*(Xpro_in-y(21)))+(Ypro*r10)-r17;

dydt(22)=((qin/Vliq)*(Xac_in-y(22)))+(Yac*r11)-r18;

dydt(23)=((qin/Vliq)*(Xh2_in-y(23))+(Yh2*r12)-r19);

dydt(24)=((qin/Vliq)*(Xi_in-y(24))+(fxi_xc*r1));

dydt(26)=((qin/Vliq)*(Sanion_in-y(26)));

dydt(27)=-ra4;

dydt(28)=-ra5;

dydt(29)=-ra6;

dydt(30)=-ra7;

dydt(31)=-ra10;

dydt(32)=-ra11;

dydt(33)=(-y(33)*(qgas/Vgas))+(rt8*(Vliq/Vgas));

dydt(34)=(-y(34)*(qgas/Vgas))+(rt9*(Vliq/Vgas));

dydt(35)=(-y(35)*(qgas/Vgas))+(rt10*(Vliq/Vgas));

end

[Torsten]

Help for what ?

[Michael Skyvalos]

no so that you get a better idea of the function