You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2010/09/29 21:50:24 UTC

svn commit: r1002829 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ode/events/ site/xdoc/ site/xdoc/userguide/ test/java/org/apache/commons/math/ode/events/

Author: luc
Date: Wed Sep 29 19:50:24 2010
New Revision: 1002829

URL: http://svn.apache.org/viewvc?rev=1002829&view=rev
Log:
Prevent ODE integration to be stuck at initial point when it is restarted after
an event has stopped the previous integration
JIRA: MATH-421

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java?rev=1002829&r1=1002828&r2=1002829&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java Wed Sep 29 19:50:24 2010
@@ -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;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java?rev=1002829&r1=1002828&r2=1002829&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java Wed Sep 29 19:50:24 2010
@@ -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.

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1002829&r1=1002828&r2=1002829&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Wed Sep 29 19:50:24 2010
@@ -78,6 +78,9 @@ The <action> type attribute can be add,u
       </action>
     </release>
     <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>

Modified: commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml?rev=1002829&r1=1002828&r2=1002829&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml Wed Sep 29 19:50:24 2010
@@ -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

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java?rev=1002829&r1=1002828&r2=1002829&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java Wed Sep 29 19:50:24 2010
@@ -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));