MATLAB: Ayuda ODE45 con condición If

if statement

Saludos,
Estoy tratando de resolver en set de ODEs usando ODE45, el problema es que algunos estados que dependen del tiempo tienes restricciones las cuales debo definir con if.
condiciones iniciales:
x0=[0.43 5 0.58 0.14 0.0024 0.01 4];
Condición inicial para el estado que no depende del tiempo
u(1)=um*((x0(3)/x0(2))/((x0(3)/x0(2))+Ksr))*(1-((x0(3)/x0(2))/Sm)^nk)*((x0(5)/(Kox*x0(1)+x0(5))));
Time definition
tsim=50;
dt=0.0001;
t(1)=0;
i=1;
Ni=round(tsim/dt);
for i=1:Ni
if t(i)<24
            F1=0;
            F2=0;
            elseif t(i)>24 && t(i)<30
            F1=80*(1/1000);
            F2=0;%70*(1/1000);
             else
            F1=0;
            F2=0;
end
dx=@(t,x)    [u(i)*x(1)-((F1+F2)/x(7))*x(1);
             -((Csx*u(i)*x(1))+(Rcsx*x(1))+Csp*((k1*u(i)*x(1))+(k2*x(1))))+(F1/x(7))*Sin-((F1+F2)/x(7))*x(2);
             -((Cnx*u(i)*x(1))+(Rcnx*x(1)))+(F2/x(7))*Nin-((F1+F2)/x(7))*x(3);
             ((k1*u(i)*x(1))+(k2*x(1)))-((F1+F2)/x(7))*x(4);
              ((KL*(O2Leq-x(5)))-(k3*u(i)*x(1))-((k4*k1*u(i)*x(1))+(k4*k2*x(1))))+(F3/x(7))*O2in-((F1+F2)/x(7))*x(5);
              ((alfa1*u(i)+alfa2)*x(1))+(alfa3)-((F1+F2)/x(7))*x(6);
              (F1+F2)];
[t,x]=ode45(dx,tspan,x0);     
u(i+1)=um*((x(3)/x(2))/((x(3)/x(2))+Ksr))*(1-((x(3)/x(2))/Sm)^nk)*((x(5)/(Kox*x(1)+x(4))));
        if u(i+1)<0
            u(i+1)=0;
        end
         t(i+1)=t(i)+dt;
          i=i+1;
  end
condiciones para los estados que dependen del tiempo
if x(3)<0.15
            x(3)=0.15;
        end
  if x (5)<0.002432
            x(5)=0.002432;
        end
aprecio su ayuda.
Gracias
Best Regards.

Best Answer

I did not run your code, and I do not follow what you are doing. However, looking at it I have two observations:
1. I do not know what tspan is, but the ODE solvers do not do well with discontinuities with respect to time. For example, I would break your ODE integration time into:
tspan1 = [ 0 24];
tspan2 = [24 30];
tspan3 = [30 60];
or whatever you want the third interval to be. Use the last values of ‘x’ of the previous interval as the initial values for the next interval. That way, ode45 does not have to integrate over the discontinuities.
2. If you want to limit ‘x(3)’ and ‘x(5)’, it is easier to use the max function than if blocks:
x(3) = max(x(3), 0.015);
x(5) = max(0.002432, x(5));
Related Question