Prevent ODE integration to be stuck at initial point when it is restarted after

an event has stopped the previous integration
JIRA: MATH-421

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_X@1002827 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java

@@ -135,11 +135,8 @@ public class CombinedEventsManager {
             if (! initialized) {
 
                 // initialize the events states
-                final double t0 = interpolator.getPreviousTime();
-                interpolator.setInterpolatedTime(t0);
-                final double [] y = interpolator.getInterpolatedState();
                 for (EventState state : states) {
-                    state.reinitializeBegin(t0, y);
+                    state.reinitializeBegin(interpolator);
                 }
 
                 initialized = true;
@@ -170,6 +167,10 @@ public class CombinedEventsManager {
             return first != null;
 
         } catch (EventException se) {
+            final Throwable cause = se.getCause();
+            if ((cause != null) && (cause instanceof DerivativeException)) {
+                throw (DerivativeException) cause;
+            }
             throw new IntegratorException(se);
         } catch (ConvergenceException ce) {
             throw new IntegratorException(ce);

src/main/java/org/apache/commons/math/ode/events/EventState.java

@@ -140,18 +140,49 @@ public class EventState {
     }
 
     /** Reinitialize the beginning of the step.
-     * @param tStart value of the independent <i>time</i> variable at the
-     * beginning of the step
-     * @param yStart array containing the current value of the state vector
-     * at the beginning of the step
+     * @param interpolator valid for the current step
      * @exception EventException if the event handler
      * value cannot be evaluated at the beginning of the step
      */
-    public void reinitializeBegin(final double tStart, final double[] yStart)
+    public void reinitializeBegin(final StepInterpolator interpolator)
         throws EventException {
-        t0 = tStart;
-        g0 = handler.g(tStart, yStart);
-        g0Positive = g0 >= 0;
+        try {
+            // excerpt from MATH-421 issue:
+            // If an ODE solver is setup with an EventHandler that return STOP
+            // when the even is triggered, the integrator stops (which is exactly
+            // the expected behavior). If however the user want to restart the
+            // solver from the final state reached at the event with the same
+            // configuration (expecting the event to be triggered again at a
+            // later time), then the integrator may fail to start. It can get stuck
+            // at the previous event.
+
+            // The use case for the bug MATH-421 is fairly general, so events occurring
+            // less than epsilon after the solver start in the first step should be ignored,
+            // where epsilon is the convergence threshold of the event. The sign of the g
+            // function should be evaluated after this initial ignore zone, not exactly at
+            // beginning (if there are no event at the very beginning g(t0) and g(t0+epsilon)
+            // have the same sign, so this does not hurt ; if there is an event at the very
+            // beginning, g(t0) and g(t0+epsilon) have opposite signs and we want to start
+            // with the second one. Of course, the sign of epsilon depend on the integration
+            // direction (forward or backward). This explains what is done below.
+
+            final double ignoreZone = interpolator.isForward() ? getConvergence() : -getConvergence();
+            t0 = interpolator.getPreviousTime() + ignoreZone;
+            interpolator.setInterpolatedTime(t0);
+            g0 = handler.g(t0, interpolator.getInterpolatedState());
+            if (g0 == 0) {
+                // extremely rare case: there is a zero EXACTLY at end of ignore zone
+                // we will use the opposite of sign at step beginning to force ignoring this zero
+                final double tStart = interpolator.getPreviousTime();
+                interpolator.setInterpolatedTime(tStart);
+                g0Positive = handler.g(tStart, interpolator.getInterpolatedState()) <= 0;
+            } else {
+                g0Positive = g0 >= 0;
+            }
+
+        } catch (DerivativeException de) {
+            throw new EventException(de);
+        }
     }
 
     /** Evaluate the impact of the proposed step on the event handler.

src/site/xdoc/changes.xml

@@ -52,6 +52,9 @@ The <action> type attribute can be add,update,fix,remove.
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="2.2" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-421">
+        Fixed an error preventing ODE solvers to be restarted after they have been stopped by a discrete event
+      </action>
       <action dev="luc" type="add" issue="MATH-419">
         Added new random number generators from the Well Equidistributed Long-period Linear (WELL).
       </action>

src/site/xdoc/userguide/ode.xml

@@ -187,7 +187,12 @@ integrator.addStepHandler(stepHandler);
           event when its sign changes. The magnitude of the value is almost irrelevant,
           it should however be continuous (but not necessarily smooth) for the sake of
           root finding. The steps are shortened as needed to ensure the events occur
-          at step boundaries (even if the integrator is a fixed-step integrator).
+          at step boundaries (even if the integrator is a fixed-step integrator). Note that
+          g function signs changes at the very beginning of the integration (from t<sub>0</sub>
+          to t<sub>0</sub> + &#x3b5; where &#x3b5; is the events detection convergence threshold)
+          are explicitely ignored. This prevents having the integration stuck at its
+          initial point when a new integration is restarted just at the same point a
+          previous one had been stopped by an event.
         </p>
         <p>
           When an event is triggered, the event time, current state and an indicator

src/test/java/org/apache/commons/math/ode/events/EventStateTest.java

@@ -49,11 +49,16 @@ public class EventStateTest {
         final double tolerance = 0.1;
         EventState es = new EventState(closeEventsGenerator, 1.5 * gap, tolerance, 10);
 
-        double t0 = r1 - 0.5 * gap;
-        es.reinitializeBegin(t0, new double[0]);
         AbstractStepInterpolator interpolator =
             new DummyStepInterpolator(new double[0], new double[0], true);
-        interpolator.storeTime(t0);
+        interpolator.storeTime(r1 - 2.5 * gap);
+        interpolator.shift();
+        interpolator.storeTime(r1 - 1.5 * gap);
+        es.reinitializeBegin(interpolator);
+
+        interpolator.shift();
+        interpolator.storeTime(r1 - 0.5 * gap);
+        Assert.assertFalse(es.evaluateStep(interpolator));
 
         interpolator.shift();
         interpolator.storeTime(0.5 * (r1 + r2));