MATLAB Answers

Thresholds in ODE solvers

10 views (last 30 days)
Daniel Schwartz
Daniel Schwartz on 28 May 2018
I am trying to implement the extinction threshold solution given in:
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?
Daniel Schwartz
Daniel Schwartz on 31 May 2018
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.

Sign in to comment.

Answers (2)

Steven Lord
Steven Lord on 28 May 2018
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.
  1 Comment
Daniel Schwartz
Daniel Schwartz on 28 May 2018
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...

Sign in to comment.

Francisco de Castro
Francisco de Castro on 16 Aug 2018
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!