Thresholds in ODE solvers

I am trying to implement the extinction threshold solution given in:

https://www.mathworks.com/matlabcentral/answers/7710-thresholds-in-ode-solvers.

However, I cannot get it to work.

Following the example I wrote this function:

function dn= expdecay(t,x) dn= -0.1*x; dead = x < 10; assignin('caller','dead',dead); evalin('caller','y(dead) = 0;');

when I call it using this code line:

[t,x]= ode45('expdecay',[1 100],[100]);

I get the following error messages:

Error using assignin Attempt to add "dead" to a static workspace. See Variables in Nested and Anonymous Functions.

Error in expdecay (line 6) assignin('caller','dead',dead);

Error in odefcncleanup>@(t,y)oldFcnFun(t,y) (line 17) newFcn = @(t,y) oldFcnFun(t,y);

Error in ode45 (line 299) f2 = odeFcn_main(t2, y2);

Error in expdecay_ODE (line 5) [t,x]= ode45('expdecay',[1 100],[100]);

What am I doing wrong?

[Francisco de Castro]

I'm the author of the original question on threshold in solvers. I run your example and I get no errors whatsoever. Have you tried to call the solver with anonymous function syntax? I don't know if that would help.

[Daniel Schwartz]

It seems to be a version issue. I have ver9.4.0.813654 (R2018a) installed and the code I posted doesn't work there. I tried the same code on ver9.2.0.556344 (R2017a) and it does work there.

Answer 1

[Steven Lord]

I recommend using an events function to stop the ODE solver when the population drops below 10. See the ballode example for a demonstration of this technique. While ballode stops when the ball's height crosses 0, you'd stop when the population crosses 10.

[Daniel Schwartz]

I already implemented that in my own code (not the example that I asked about). That works for the time point at which the threshold is crossed:

• I use the event to stop the solver
• I update the y values (in the example changing 10 to 0)
• I start again from the time point the solver was terminated by the event.

However, my model also includes another source for x, and I want to keep it extinct unless it rises above the extinction threshold.

I think the solution given to the question I referred to in my own question is what I need. I just can't get it to work...

Answer 2

[Francisco de Castro]

I found a partial (and not very elegant) solution for this, albeit only for some of the ODE solvers.

1. Write the example ode function as:

    function dn = expdecay(t,x)   
    thresh= 1;    
    largenegrate= -1E2;

    dn= -0.1*x;
    dn(x < thresh)= largenegrate;

2. Set option non-negative for all components of the solution (only 1 in this example):

opt= odeset('NonNegative',1);

3. Test:

[t,x]= ode45(@expdecay,[1 50],10,opt); plot(t,x), grid on

In this case, the extinction threshold is 1, and -1E2 is just a negative number much larger (in abs. value) than the typical rates in your system. However, only works with the solvers that admit the NonNegative option. At least some solvers of stiff systems do not take this option and you'll get negative values in the solution.