Events not consistently triggering with ODE45

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?

Answer 1

[Mischa Kim]

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;

[Conrad]

Thanks for responding.The ode45 solver doesn't return the values of te and ye. The function is called with the following code:

options = odeset('Events',@(t,Y) cdevent(t,Y,lakeBottom));

[t, H, te, ye, ie] = ode45(@(t, Y) drainageODE(t, Y, K, L, Qin, areaVec, heightVec, lakeBottom), [0,3600], [Hl, Hs],options)

and then I receive the error.

The 0.1 is the threshold so that the ode terminates 10cm before the lake bottom.

[Mischa Kim]

Of course, I need to read more carefully. Could you attach the entire code (see paper clip)?

[Conrad]

I've attached 4 files. cdevent.m is the event function, drainageODE contains the system of ODEs to be solved, and odeTest.m, loads the data in Param.mat and sets the ODE solver in motion.

[Mischa Kim]

Hello Conrad, reduce the error tolerances (in the main file). That should do it.

options = odeset('Events',@(t,Y) cdevent(t,Y,lakeBottom),'RelTol',1e-10,'AbsTol',1e-10);

You might want to experiment with different tolerance settings if the simulation time gets inhibitively long.

[Conrad]

Hello Mischa, Thanks, that does the trick. However, I'm perplexed why the ode solver doesn't trigger the event with a lower tolerance, since the documentation says matlab looks for the value crossing zero. Since I'm investigating the cumulative impact of many lake drainages rather then the behavior of any one lake, I'm looking for the fastest solution. Is there a way to get the event caught with the lower tolerance? Thanks!

[Mischa Kim]

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.

[Conrad]

The idea of a trade-off between step size and accuracy makes sense. But that doesn't address the issue that if you accept less accurate results, why you can't catch the event. The magnitude or accuracy of the change of a value, shouldn't prevent matlab from being able to detect when it goes from positive to negative. I'm struggling to understand why when in the event function, value is positive at one time step, and negative at the next, matlab can't trigger the event. Is there a way to put a conditional test in the event function, and use that to terminate ode45?

Answer 2

[Darren Wethington]

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.