Events not consistently triggering with ODE45
10 ビュー (過去 30 日間)
古いコメントを表示
I am calculating lake drainage through a stream using ode45, and set a threshold -- 10cm above the lake bottom -- using the 'Events' property when the solver should terminate. However, the event function doesn't seem to be triggered consistently by the lake height dropping below the threshold, resulting in the following error:
------------------------------
Error using odezero (line 50)
odezero: an event disappeared (internal error)
Error in ode45 (line 353)
[te,ye,ie,valt,stop] = ...
------------------------------
This is my Events function:
-----------------------------
function [value,isterminal,direction] = cdevent(t,Y, lakeBottom)
value = Y(1) - lakeBottom - .1;
disp(value)
isterminal = 1;
direction = -1;
------------------------------
As the program runs, disp(value) gives the following output:
3.2317
3.0090
2.7688
2.4971
2.1645
1.7984
1.3130
0.6979
-0.1054
NaN
The NaN is the result when the lake level drops below the assigned bottom, and properties such as area are no longer defined. I am trying to figure out why on the step that value goes from 0.69 to -0.1 the ODE doesn't terminate?
0 件のコメント
採用された回答
Mischa Kim
2014 年 1 月 30 日
編集済み: Mischa Kim
2014 年 1 月 30 日
Please verify that the solver returns the correct event values, ye, te (simply remove the semi-colon at the end of the ode45 command). Also, what is the 0.1 you are subtracting in
value = Y(1) - lakeBottom - .1;
7 件のコメント
Mischa Kim
2014 年 1 月 30 日
The solver works with a default error tolerance (1e-3, or something like that) that typically needs to be adjusted for more complex, numerically more demanding problems. The error tolerance determines the step size the solver can and will take to deliver the result. Greater step size means shorter sim time and less accurate results, which may result in issues similar to the ones you observed.
その他の回答 (1 件)
Darren Wethington
2019 年 3 月 7 日
I've been digging around Matlab's ODE Events and searching documentation this week when I came across this question. I was able to find the solution for your problem, if you still are curious several years later.
When you execute disp(value), what you're seeing is something along the lines of:
"Value's positive"
"Value's positive"
...
"Value is negative! Let's go back and figure out where it's exactly 0."
"Value is NAN. We don't know what to do."
I've tracked it down to some specific values: if you call drainageODE(2880,[40.4271,35.4175],...), then A=max(...,0) becomes 0. As a result, hldot=...+Qin/A becomes NAN, and you get your error. Conditional breakpoints helped me to find this result, if you'd like to try to find it yourself as well.
0 件のコメント
参考
カテゴリ
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!