0

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

Dahlai
  • 695
  • 1
  • 4
  • 17
  • 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 Apr 14 '15 at 23:36
  • "[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 Apr 14 '15 at 23:54
  • 1
    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](http://www.mathworks.com/help/matlab/ref/odeset.html?refresh=true#f92-1017470) in your case to do this properly – [see my answer here](http://stackoverflow.com/a/28436049/2278029). – horchler Apr 15 '15 at 00:19
  • 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 Apr 15 '15 at 05:00
  • 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 Apr 15 '15 at 14:00
  • okay the way I understand it, the event detection can only help with zero crossings of the function, which is still something I will have a look into. In my final simulation, however, I will never have zero crossing. Mainly I am looking for a way to keep track of my v over time compared to Cc and Cs. – Dahlai Apr 15 '15 at 15:21
  • @Dahlai: You define what zero means in the zeros crossings. It doesn't mean that the state itself need to cross zero. When you write your Events function you define exactly the arbitrary conditions that will be met. When those conditions are zero, an event is triggered. Also, please use the "@" prefix when referring to usernames or they won't be notified. – horchler Apr 15 '15 at 22:40
  • @Horcher: Making sure I understand this correctly, the change in v leads to a discontinuity in my derivatives. Therefore, I have to restart the solver at that point, using event detection. Yes/No? – Dahlai Apr 16 '15 at 11:03

1 Answers1

1

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Thanks for the clarification Lutzl! I clarified my original question. I forgot to mention that this behavior can be observed around the switch of v (from 0 to some other positive value), rather then when y(2) becomes 0. I think your answer still applies, though! I will give event detection a try to see whether this behavior can be circumvented. I actually want to track and overlay v with Cc and Cs (generated from the differential equations) but at the moment I get about three times more v values than data points from the solver. – Dahlai Apr 15 '15 at 13:56
  • 1
    This is equally bad as using an `if` statement (it is a nice trick for speeding up code in some cases though) and not a substitute for solving the problem properly with event detection. Just because integrations steps are "successful" does not mean that the solution is correct. It may work here in some cases, but it's a bad idea to use in general I'd never trust code with this in it without a lot more detailed experimentation. Lastly, one should probably use a stiff solver like `ode15s` in this case. – horchler Apr 15 '15 at 22:35