Fixed wrong event detection in case of close events pairs.

JIRA: MATH-1226

src/changes/changes.xml

@@ -54,6 +54,9 @@ If the output is not quite correct, check for invisible trailing spaces!
     </release>
 
     <release version="4.0" date="XXXX-XX-XX" description="">
+      <action dev="luc" type="fix" issue="MATH-1226"> <!-- backported to 3.6 -->
+        Fixed wrong event detection in case of close events pairs.
+      </action>
       <action dev="luc" type="add" >
         Reimplemented pow(double, double) in FastMath, for better accuracy in
         integral power cases and trying to fix erroneous JIT optimization again.

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

@@ -296,7 +296,18 @@ public class EventState {
                             ta = forward ? ta + convergence : ta - convergence;
                             ga = f.value(ta);
                         } while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb)));
-                        --i;
+
+                        if (forward ^ (ta >= tb)) {
+                            // we were able to skip this spurious root
+                            --i;
+                        } else {
+                            // we can't avoid this root before the end of the step,
+                            // we have to handle it despite it is close to the former one
+                            // maybe we have two very close roots
+                            pendingEventTime = root;
+                            pendingEvent = true;
+                            return true;
+                        }
                     } else if (Double.isNaN(previousEventTime) ||
                                (FastMath.abs(previousEventTime - root) > convergence)) {
                         pendingEventTime = root;

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

@@ -29,6 +29,7 @@ import org.apache.commons.math4.ode.SecondaryEquations;
 import org.apache.commons.math4.ode.events.EventHandler;
 import org.apache.commons.math4.ode.events.EventState;
 import org.apache.commons.math4.ode.nonstiff.DormandPrince853Integrator;
+import org.apache.commons.math4.ode.nonstiff.LutherIntegrator;
 import org.apache.commons.math4.ode.sampling.AbstractStepInterpolator;
 import org.apache.commons.math4.ode.sampling.DummyStepInterpolator;
 import org.junit.Assert;
@@ -43,21 +44,9 @@ public class EventStateTest {
         final double r1  = 90.0;
         final double r2  = 135.0;
         final double gap = r2 - r1;
-        EventHandler closeEventsGenerator = new EventHandler() {
-            public void init(double t0, double[] y0, double t) {
-            }
-            public void resetState(double t, double[] y) {
-            }
-            public double g(double t, double[] y) {
-                return (t - r1) * (r2 - t);
-            }
-            public Action eventOccurred(double t, double[] y, boolean increasing) {
-                return Action.CONTINUE;
-            }
-        };
 
         final double tolerance = 0.1;
-        EventState es = new EventState(closeEventsGenerator, 1.5 * gap,
+        EventState es = new EventState(new CloseEventsGenerator(r1, r2), 1.5 * gap,
                                        tolerance, 100,
                                        new BrentSolver(tolerance));
         es.setExpandable(new ExpandableStatefulODE(new FirstOrderDifferentialEquations() {
@@ -226,4 +215,63 @@ public class EventStateTest {
 
     }
 
+    @Test
+    public void testEventsCloserThanThreshold()
+        throws DimensionMismatchException, NumberIsTooSmallException,
+               MaxCountExceededException, NoBracketingException {
+
+        FirstOrderDifferentialEquations equation = new FirstOrderDifferentialEquations() {
+            
+            public int getDimension() {
+                return 1;
+            }
+            
+            public void computeDerivatives(double t, double[] y, double[] yDot) {
+                yDot[0] = 1.0;
+            }
+        };
+
+        LutherIntegrator integrator = new LutherIntegrator(20.0);
+        CloseEventsGenerator eventsGenerator =
+                        new CloseEventsGenerator(9.0 - 1.0 / 128, 9.0 + 1.0 / 128);
+        integrator.addEventHandler(eventsGenerator, 1.0, 0.02, 1000);
+        double[] y = new double[1];
+        double tEnd = integrator.integrate(equation, 0.0, y, 100.0, y);
+        Assert.assertEquals( 2, eventsGenerator.getCount());
+        Assert.assertEquals( 9.0 + 1.0 / 128, tEnd, 1.0 / 32.0);
+
+    }
+
+    private class CloseEventsGenerator implements EventHandler {
+
+        final double r1;
+        final double r2;
+        int count;
+
+        public CloseEventsGenerator(final double r1, final double r2) {
+            this.r1    = r1;
+            this.r2    = r2;
+            this.count = 0;
+        }
+
+        public void init(double t0, double[] y0, double t) {
+        }
+
+        public void resetState(double t, double[] y) {
+        }
+
+        public double g(double t, double[] y) {
+            return (t - r1) * (r2 - t);
+        }
+
+        public Action eventOccurred(double t, double[] y, boolean increasing) {
+            return ++count < 2 ? Action.CONTINUE : Action.STOP;
+        }
+
+        public int getCount() {
+            return count;
+        }
+
+    }
+
 }