Matlab understanding ode solver

Question

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

How are you calling ode23? That does affect the results that are returned. If you specify just an initial time and final time, you will get all the intermediate time points in the results, but if you specify a list of time points, Matlab will still calculate the same intermediate time points as if you had just specified initial/final times, then it will interpolate the results onto the time points you give it. Is this what you were asking about? — David
"[With] the debugging function [, I] noticed that it starts to jump back in 'time' around this point": you might be observing the adaptive time-step algorithm. Since your function's derivative is, I'm guessing, discontinuous at a point in time, the solver may have to backtrack to ensure it's convergence criteria are satisfied and only the converged state value is stored. — TroyHaskin
Bad idea to put a discontinuity into the integration function itself – that's not what they're for. It will make your equations stiff, slower to simulate, and likely produce inaccurate or completely wrong results. If you change the value of a parameter, you should stop the integration and restart it with the new parameter specified. You'll need to use event detection in your case to do this properly – see my answer here. — horchler
The function is not discontinuous, max(0,y(2)) is continuous in the function value. However, the derivative is discontinuous and that brings an order reduction in the local discretization error. — Lutz Lehmann
Horchler, thanks for pointing this out to me. I will give it a try! Also, I clarified my problem in the original question (see "//Edit") and in a comment to solution below, if you want to have a look and tell me whether event detection will help. — Dahlai

Lutz Lehmann Lutz Lehmann · Accepted Answer · 2015-04-15T06:47:30

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.

Matlab understanding ode solver

1 Answers