Mailing-List: contact issues-help@commons.apache.org; run by ezmlm
Precedence: bulk
Reply-To: issues@commons.apache.org
Date: Fri, 1 Mar 2013 15:13:12 +0000 (UTC)
From: "Luc Maisonobe (JIRA)" <jira@apache.org>
To: issues@commons.apache.org
Message-ID: <JIRA.12634597.1362052108101.364970.1362150792911@arcas>
In-Reply-To: <JIRA.12634597.1362052108101@arcas>
References: <JIRA.12634597.1362052108101@arcas>
Subject: [jira] [Comment Edited] (MATH-937) NoBracketingException after
 event was found
MIME-Version: 1.0
Content-Type: text/plain; charset=utf-8
Content-Transfer-Encoding: quoted-printable


    [ https://issues.apache.org/jira/browse/MATH-937?page=3Dcom.atlassian.j=
ira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=3D135906=
13#comment-13590613 ]=20

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. I=
n 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 re=
quired in the *neighborhood* of the root, i.e. both just before and just af=
ter 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 a=
t which the function crosses 0. For this, we need continuity before the eve=
nt 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 n=
eeds to have a state corresponding to what happens just after the event. If=
 computations were perfect, the value g0 =3D 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 afte=
r t0. This helps keeping everything consistent.

In this issue, the event occurs as the g function decreases from positive t=
o negative values. We find a very good estimate of the root at t0 =3D 1.427=
8431229270647, which corresponds to g0 =3D g(t0) =3D -2.7755575615628914E-1=
6. 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 fun=
ction should (and in fact is) negative.

However, the user code explicitly changes the slope at this time and the fu=
nction starts to increase again. This is of course what the user wants beca=
use 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 w=
e 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 expec=
t 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 expl=
ain better the problem and show the sign trick to users. We could write som=
ething like:

{panel:title=3Dimproved 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) =3D h(t) where h is the height above t=
he floor
at time t. When g(t) reaches 0, the ball is on the floor, so it should boun=
ce
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 val=
ues
to 0, and after the event g(t) would be increasing from 0 to positive value=
s
again. Consistency is broken here! The solution here is to have g(t) =3D si=
gn * h(t),
where sign is a variable that starts at +1. Then, each time eventOccurred i=
s called,
sign should be reset to -sign. This allows the g(t) function to remain cont=
inuous
(and even smooth) even across events, despite h(t) is not.=20
{panel}

Would this be OK for you?

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

As you identified yourself, the function is not continuous near the zero. I=
n 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 a=
t which the function crosses 0. For this, we need continuity before the eve=
nt 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 n=
eeds to have a state corresponding to what happens just after the event. If=
 computations were perfect, the value g0 =3D 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 afte=
r t0. This helps keeping everything consistent.

In this issue, the event occurs as the g function decreases from positive t=
o negative values. We find a very good estimate of the root at t0 =3D 1.427=
8431229270647, which corresponds to g0 =3D g(t0) =3D -2.7755575615628914E-1=
6. 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 fun=
ction should (and in fact is) negative.

However, the user code explicitly changes the slope at this time and the fu=
nction starts to increase again. This is of course what the user wants beca=
use 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 w=
e 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 expec=
t 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 expl=
ain better the problem and show the sign trick to users. We could write som=
ething like:

{panel:title=3Dimproved 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) =3D h(t) where h is the height above t=
he floor
at time t. When g(t) reaches 0, the ball is on the floor, so it should boun=
ce
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 val=
ues
to 0, and after the event g(t) would be increasing from 0 to positive value=
s
again. Consistency is broken here! The solution here is to have g(t) =3D si=
gn * h(t),
where sign is a variable that starts at +1. Then, each time eventOccurred i=
s called,
sign should be reset to -sign. This allows the g(t) function to remain cont=
inuous
(and even smooth) even across events, despite h(t) is not.=20
{panel}

Would this be OK for you?

                 =20
> 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=C3=B6ger
>         Attachments: ApacheCommonsBouncingBall.java, ApacheCommonsBouncin=
gBallTest.java
>
>
> The BracketingNthOrderBrentSolver used by the EmbeddedRungeKuttaIntegrato=
r fails, if an event is detected twice with a NoBracketingException.
> The problem lies in line EventState.java line 262 (version 3.1.1). Here t=
he event detection function f is applied to an arbitrary choosen value of t=
ime. 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 administrato=
rs
For more information on JIRA, see: http://www.atlassian.com/software/jira