event location questions in MATLAB

Question

Suppose in matlab the following:

[t, x, te, xe, ie] = ode15s(@myfunc, [tStart tFinal], x0, odeset('Events', @events));

Question 1

1a) The function events is called only after a successful step of the solver. Is this true?

1b) Just after the solver has made a successful step, is it possible that the last call of myfunc not be the call that lead to the successful step?

1c) If the events function contains multiple terminal events and upon a successful step it is detected that two of them (not just one) have occured, what would be the behavior of the solver?

Question 2

Suppose that myfunc contains the following code

if (check(x) > 2)
    dx(3) = x(1)*x(2);
else
    dx(3) = x(2)^2;
end

where check is some function of x.

One way to solve this problem is not use an events function. In my experience the ode solver can solve such problems that way.

The other way to solve this problem is use an events function to locate check(x) - 2 == 0, one terminal event for direction = 1 and another one for direction = -1. After the solver stops on either of the events a global variable eg myvar is set appropriately to distinguish between the two events, and then the simulation continues from where it stopped. In that case the code in myfunc will be

if (myvar == 1)
    dx(3) = x(1)*x(2);
else
    dx(3) = x(2)^2;
end

Both way yield correct results in simple cases. However I am trying to solve a very complex problem (additional events than the above and discontinous right-hand parts of the differential equations that are proven to be solvable at some cases) and I am trying to find out if the first way would yield different results than the second one.

One might tell that the ode would either fail to return a solution before tFinal or return a correct solution, but due to the discontinuity of the right-hand part the solver might not return a solution while a solution exists.

So in some sense, the question is: what is the practical-theoretical difference between using the first way and the second way?

Answer 1

Since I've spent some effort on these questions, I'm posting back some feedback.

Question 1

1a) Yes this is true. For reference see for example the ode15s.m matlab solver. However note that before the solver continues solving, the events function may be called several times for the sake of a more accurate te value.

1b) Yes this is also true.

1c) In this case the solver would terminate returning an ie vector containing the two (or even more) indexes of the events that stopped the solver. In that case, the te vector would contain equal elements (te(1) == te(2) will always return true). This is the only way to distiguish "double events" (meaning events that simultaneously stopped the solver after the same successful step) from "fake" events that are recorded when the ode solver continues solving after a terminal event (to understand better what I'm saying also read ode solver event location index in MATLAB ).

Tracing the odezero function will make 1c answer very clear.

Question 2

Now this is a tricky answer. **In general* both ways return correct results. However (and most naturally) they are not bound to return the exact solution points at the exact time points with the exact number of steps.

The notable difference between the two ways is that in the second one, we have a branch change only when a check(x)-2 sign change occurs using only the currently active branch. For example, suppose that the currently active branch is the first one. When the solver notices a sign change in check(x) - 2 after a successful step that was produced using only that branch, only then changes to the second branch. In simple words, both successful and unsuccessful steps are calculated using the very same branch before a use of the other branch can occur. However if we use the first way we may notice a use of the non-active branch during (for example) an unsucessful step.

With these in mind comes the verdict; the most general and strictly correct way to choose is the second one (using events). The first way should also return correct results. HOWEVER, due to the difference between the two ways, the first one could fail in very specific/extreme problems. I'm very tempted to provide information about my case, in which ONLY the second way could be safely used, but it truly is a long way to go.