CODEWITHSUNDEEP
matlabode
Dahlai
Dahlai

Reputation: 705

Matlab understanding ode solver

I have a system of linked differential equations that I am solving with the ode23 solver. When a certain threshold is reached one of the parameters changes which reverses the slope of my function.

I followed the behavior of the ode with the debugging function and noticed that it starts to jump back in "time" around this point. Basically it generates more data points.However, these are not all represented in the final solution vector.

Can somebody explain this behavior, especially why not all calculated values find their way into the solution vector?

//Edit: To clarify, the behavior starts when v changes from 0 to any other value. (When I write every value of v to a vector it has more than a 1000 components while the ode solver solution only has ~300).

Find the code of my equations below:

%chemostat model, based on:
%DCc=-v0*Cc/V + umax*Cs*Cc/(Ks+Cs)-rd
%Dcs=(v0/V)*(Cs0-Cs) - Cc*(Ys*umax*Cs/(Ks+Cs)-m)
function dydt=systemEquationsRibose(t,y,funV0Ribose,V,umax,Ks,rd,Cs0,Ys,m)
     v=funV0Ribose(t,y); %funV0Ribose determines v dependent on y(1)

if y(2)<0
    y(2)=0
end
      dydt=[-(v/V)*y(1)+(umax*y(1)*y(2))/(Ks+y(2))-rd; 
           (v/V)*(Cs0-y(2))-((1/Ys)*(umax*y(2)*y(1))/(Ks+y(2)))];

Thanks in advance!

Cheers, dahlai

Upvotes: 0

Views: 357

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25972

The first conditional can also be expressed as

y(2) = max(0, y(2)).

As one can see, this is still a continuous function, but with a kink, i.e., a discontinuity in the first derivative. One can this also interpret as a point with curvature radius 0, i.e., infinite curvature.

ode23 uses an order 2 method to integrate, an order 3 method to estimate the error and probably the order 1 Euler step to estimate stiffness.

An integration step over the kink renders all discretization errors to be order 1 (or 2, depending on the convention), confounding the logic of the step size control. This forces a rather radical step-size reduction, but since that small step then falls, most probably, short of the kink, the correct orders are found again, resulting in a step-size increase in the next step which could again go over the kink etc.

The return array only contains successful integration steps, not the failed attempts of the step-size control.

Upvotes: 1

Related Questions