Thresholds in ODE solvers

7 次查看（过去 30 天）
Daniel Schwartz
Daniel Schwartz 2018-5-28
回答： Francisco de Castro 2018-8-16
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 个评论
Francisco de Castro
Francisco de Castro 2018-5-30
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 2018-5-31
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Steven Lord
Steven Lord 2018-5-28
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 个评论
Daniel Schwartz
Daniel Schwartz 2018-5-28
编辑：Daniel Schwartz 2018-5-28
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...

请先登录，再进行评论。

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.

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by