[jira] [Comment Edited] (MATH-937) NoBracketingException after event was found

["Luc Maisonobe (JIRA)" <ji...@apache.org>]

    [ https://issues.apache.org/jira/browse/MATH-937?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=13590613#comment-13590613 ] 

Luc Maisonobe edited comment on MATH-937 at 3/1/13 3:12 PM:
------------------------------------------------------------

I have analyzed the problem.

As you identified yourself, the function is not continuous near the zero. In fact, it is both discontinuous as you wrote in the comment above, and it reverts its slope. I'm not sure we can identify this in all cases. This is exactly the reason why we have written in the javadoc that continuity is required in the *neighborhood* of the root, i.e. both just before and just after the event. I will improve the documentation in this case.

What happens under the hood is that we use a solver to locate the instant at which the function crosses 0. For this, we need continuity before the event is triggered. Then, if the event did not trigger a stop, the integrator continues after the event. In order to detect the next event properly, it needs to have a state corresponding to what happens just after the event. If computations were perfect, the value g0 = g(t0) should be exactly 0, but due to convergence criteria, it may be slightly non-zero. So we store t0, g0 and a boolean g0Positive that indicate that indicates the sign just after t0. This helps keeping everything consistent.

In this issue, the event occurs as the g function decreases from positive to negative values. We find a very good estimate of the root at t0 = 1.4278431229270647, which corresponds to g0 = g(t0) = -2.7755575615628914E-16. It is very close to zero and the solver properly selected a value after the real root, i.e. a negative one. As additional protection the g0Positive was properly set to false, i.e. we are really sure after the event the function should (and in fact is) negative.

However, the user code explicitly changes the slope at this time and the function starts to increase again. This is of course what the user wants because we want the function to represent something bouncing on the floor, i.e. when the altitude decreases down to zero, it should increase again. So the function does not really *crosses* the zero line, it reaches it and then we change it. As numerical inaccuracies arise, we introduce a discontinuity near the zero and in fact create another which really occurs before the one we have already detected. Then the algorithm gets lost as it did not expect this discontinuity.

We have already identified this problem when we developed TestProblem4 and we introduced the sign attribute in the Bounce event class exactly for this reason. This sign attribute helps preserving consistency of the g function on both sides of the event.

So I would say we cannot improve the code here, but we should probably explain better the problem and show the sign trick to users. We could write something like:

{panel:title=improved javadoc}
The discrete events are generated when the sign of this
switching function changes. The integrator will take care to change
the stepsize in such a way these events occur exactly at step boundaries.

The switching function must be continuous in its roots neighborhood
(but not necessarily smooth), as the integrator will need to find its
roots to locate precisely the events.

Also note that the integrator expect that once an event has occurred, the
sign of the switching function at the start of the next step (i.e. just
after the event) is the opposite of the sign just before the event. This
consistency between the steps *must* be preserved, otherwise exceptions
related to root not being bracketed will occur.

This need for consistency is sometimes tricky to achieve. A typical example
is using an event to model a ball bouncing on the floor. The first idea to
represent this would be to have g(t) = h(t) where h is the height above the floor
at time t. When g(t) reaches 0, the ball is on the floor, so it should bounce
and the typical way to do this is to reverse its vertical velocity. However,
this would mean that before the event g(t) was decreasing from positive values
to 0, and after the event g(t) would be increasing from 0 to positive values
again. Consistency is broken here! The solution here is to have g(t) = sign * h(t),
where sign is a variable that starts at +1. Then, each time eventOccurred is called,
sign should be reset to -sign. This allows the g(t) function to remain continuous
(and even smooth) even across events, despite h(t) is not. 
{panel}

Would this be OK for you?

                
      was (Author: luc):
    I have analyzed the problem.

As you identified yourself, the function is not continuous near the zero. In fact, it is both discontinuous as you wrote in the comment above, and it reverts its slope. I'm not sure we can identify this in all cases. This is exactly the reason while we have written in the javadoc that continuity is required in the *neighborhood* of the root, i.e. both just before and just after the event. I will improve the documentation in this case.

What happens under the hood is that we use a solver to locate the instant at which the function crosses 0. For this, we need continuity before the event is triggered. Then, if the event did not trigger a stop, the integrator continues after the event. In order to detect the next event properly, it needs to have a state corresponding to what happens just after the event. If computations were perfect, the value g0 = g(t0) should be exactly 0, but due to convergence criteria, it may be slightly non-zero. So we store t0, g0 and a boolean g0Positive that indicate that indicates the sign just after t0. This helps keeping everything consistent.

In this issue, the event occurs as the g function decreases from positive to negative values. We find a very good estimate of the root at t0 = 1.4278431229270647, which corresponds to g0 = g(t0) = -2.7755575615628914E-16. It is very close to zero and the solver properly selected a value after the real root, i.e. a negative one. As additional protection the g0Positive was properly set to false, i.e. we are really sure after the event the function should (and in fact is) negative.

However, the user code explicitly changes the slope at this time and the function starts to increase again. This is of course what the user wants because we want the function to represent something bouncing on the floor, i.e. when the altitude decreases down to zero, it should increase again. So the function does not really *crosses* the zero line, it reaches it and then we change it. As numerical inaccuracies arise, we introduce a discontinuity near the zero and in fact create another which really occurs before the one we have already detected. Then the algorithm gets lost as it did not expect this discontinuity.

We have already identified this problem when we developed TestProblem4 and we introduced the sign attribute in the Bounce event class exactly for this reason. This sign attribute helps preserving consistency of the g function on both sides of the event.

So I would say we cannot improve the code here, but we should probably explain better the problem and show the sign trick to users. We could write something like:

{panel:title=improved javadoc}
The discrete events are generated when the sign of this
switching function changes. The integrator will take care to change
the stepsize in such a way these events occur exactly at step boundaries.

The switching function must be continuous in its roots neighborhood
(but not necessarily smooth), as the integrator will need to find its
roots to locate precisely the events.

Also note that the integrator expect that once an event has occurred, the
sign of the switching function at the start of the next step (i.e. just
after the event) is the opposite of the sign just before the event. This
consistency between the steps *must* be preserved, otherwise exceptions
related to root not being bracketed will occur.

This need for consistency is sometimes tricky to achieve. A typical example
is using an event to model a ball bouncing on the floor. The first idea to
represent this would be to have g(t) = h(t) where h is the height above the floor
at time t. When g(t) reaches 0, the ball is on the floor, so it should bounce
and the typical way to do this is to reverse its vertical velocity. However,
this would mean that before the event g(t) was decreasing from positive values
to 0, and after the event g(t) would be increasing from 0 to positive values
again. Consistency is broken here! The solution here is to have g(t) = sign * h(t),
where sign is a variable that starts at +1. Then, each time eventOccurred is called,
sign should be reset to -sign. This allows the g(t) function to remain continuous
(and even smooth) even across events, despite h(t) is not. 
{panel}

Would this be OK for you?

                  
> NoBracketingException after event was found
> -------------------------------------------
>
>                 Key: MATH-937
>                 URL: https://issues.apache.org/jira/browse/MATH-937
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 3.1.1
>            Reporter: Christoph Höger
>         Attachments: ApacheCommonsBouncingBall.java, ApacheCommonsBouncingBallTest.java
>
>
> The BracketingNthOrderBrentSolver used by the EmbeddedRungeKuttaIntegrator fails, if an event is detected twice with a NoBracketingException.
> The problem lies in line EventState.java line 262 (version 3.1.1). Here the event detection function f is applied to an arbitrary choosen value of time. If the event detector crosses zero before this time, the solver throws the mentioned exception.

--
This message is automatically generated by JIRA.
If you think it was sent incorrectly, please contact your JIRA administrators
For more information on JIRA, see: http://www.atlassian.com/software/jira