Ayuda ODE45 con condición If

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.

Answer 1

[Star Strider]

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));

Answer 2

[Cesar García]

Thank you for your response, i am traying to pass from euler to ode45 solution technique. Attached is the eulers method.

I appreciate your help.

[Star Strider]

I cannot make sense out of that. I would appreciate your uploading a description of your problem as a PDF file or other description of your problem.

[Cesar García]

Hi,

In that script i am solving a set of non-linear differential equations. I have the fed-batch model with the states (X,S,N,P,O2L, CO2 and V), these states are function of time and i have a polymerization model where we have the states (M,ESHM,P,P1 and P2) this polymerization model is solved by a discretization technique called fixed pivot technique. This is function of time "i" and the discretization "j" as well. I have already the matrix A(j,k) and B(j,k) and the script is working. I have solved this set of equations by the euler´s method, but, the polymerization model arise some problems with the numerical method, this is the reason why i want to use ode.

Attached is the paper that i am trying to reproduce. (A combined metabolic/polymerization kinetic model on the microbial production of poly(3-hydroxybutyrate)-Giannis Penloglou et al-2010)

I´ll appreciate all your help. And if is possible, i can explain you in a better way by skype.

I really need help with this.

Best regards

[Star Strider]

It would take me weeks to understand that, even considering that I have a relatively strong (although not current) background in chemistry and biochemistry, so I will leave the details to you. My background is not sufficient to understand that paper.

What parts do you need help on? It is probably best for me to help you with the MATLAB coding. The paper is necessary because it is easier for me to read the relevant parts and ask you about them. It will take a bit of time regardless.

[Cesar García]

Hi,

Thank you for your response.

I understand that it takes many time to understand, so far for me.

In the code, i need to write the states of fed-batch model (X,S,N,P,O2L,CO2 ans V) and the states of polymerization model (M,ESHM,P,P1 and D) into an ODE, this is because i can use an stiff solution for this method. They actually are solved by the euler´s method, but it takes many computational time and a big step size in order to reach a good solution.

At the end, i want to reproduce de Molecular weight distribution during time of fermentation but the real problem is to apply the ODE for stiff problem like this.

If you want, i can send you the .m at your e-mail.

Best Regards.

[Star Strider]

What equations in the paper do you want to integrate? (Please be aware that biochemical engineering is very far from my areas of expertise.)

I will look through the code you posted for the constants.

[Cesar García]

Hi,

In the paper, equations 8 for [M], 9 for [ESHM], and equations A2, A3 and A4 for [P],[P1] and [D] in the code.

[Star Strider]

Eq (8) is probably straightforward if we can assume the summation of Pn is a constant. What is Jm(t)? Do you have a function expression for it as defined in (12)?

I have to admit that I do not understand (A2)-(A4) at all. They apparently have to be iterated over j, and the summation done at each step would be a constant. You would have to put these in a loop with a fixed vector of integration times, and then store each vector of results, most easily in a cell array although a matrix would work as well. It will obviously not be possible to iterate (A3) an infinite number of times.

The ode15s function would seem to me to be the obvious choice for these. but they will not be straightforward to code.

[Cesar García]

Actually, i already have working all the states as they should. In the code i have defined Jm (production rate term) as:

Jf(i)=(Csx*u(i)*X(i))+(Rcsx*X(i))+Csp*((k1*u(i)*X(i))+(k2*X(i)));

Jm(i)=(Jf(i)*Yms);

Yms=Calculated constant.

A2-A4 are the discret solution of the equations 5,6 and 7. This is using the fixed pivot technique (this is working ok)

The idea is try to apply an ode15s. Do you have any idea how to do it?

[Star Strider]

It should be relatively easy to apply it to (8) and (9) if you already have the code (as functions of time) for the rest of it.

I have been working on this intermittently and cannot find some necessary information. I cannot figure out how to get an expression in time for ‘Jf’. The ODE functions return time to the ODE you are integrating, so writing a function to return a time-dependent value for ‘Jf’ (and therefore ‘Jm’) would be straightforward if I already had an expression in time for the functions you use to calculate it.

The infinite sum for ‘Pn’ is another problem, but lacking a function for ‘Jf’ renders it irrelevant.

So while I have an idea of how to do it, without the requisite background information, I cannot write a function to use in ode15s to integrate it.

Also, this would likely be easier to implement in SimBiology if you have access to it.

[Cesar García]

Hi,

well, Jf(t) is a function of the states (X,S, N and so on), we get a value based in the solution of the dynamics states.

i will check the simbiology.

Thank you very much