Thresholds in ODE solvers

14 Ansichten (letzte 30 Tage)
Daniel Schwartz
Daniel Schwartz am 28 Mai 2018
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?
  2 Kommentare
Francisco de Castro
Francisco de Castro am 30 Mai 2018
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
Daniel Schwartz am 31 Mai 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.

Melden Sie sich an, um zu kommentieren.

Antworten (2)

Steven Lord
Steven Lord am 28 Mai 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 Kommentar
Daniel Schwartz
Daniel Schwartz am 28 Mai 2018
Bearbeitet: Daniel Schwartz am 28 Mai 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...

Melden Sie sich an, um zu kommentieren.


Francisco de Castro
Francisco de Castro am 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!

Translated by