Added a way to specify a custom root solver to find events in ODE integration.

JIRA: MATH-586

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

src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java

@@ -28,6 +28,8 @@ import java.util.TreeSet;
 
 import org.apache.commons.math.ConvergenceException;
 import org.apache.commons.math.MaxEvaluationsExceededException;
+import org.apache.commons.math.analysis.solvers.BrentSolver;
+import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
 import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.exception.util.LocalizedFormats;
 import org.apache.commons.math.ode.events.EventException;
@@ -123,7 +125,18 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
                                 final double maxCheckInterval,
                                 final double convergence,
                                 final int maxIterationCount) {
-        eventsStates.add(new EventState(handler, maxCheckInterval, convergence, maxIterationCount));
+        addEventHandler(handler, maxCheckInterval, convergence,
+                        maxIterationCount, new BrentSolver(convergence));
+    }
+
+    /** {@inheritDoc} */
+    public void addEventHandler(final EventHandler handler,
+                                final double maxCheckInterval,
+                                final double convergence,
+                                final int maxIterationCount,
+                                final UnivariateRealSolver solver) {
+        eventsStates.add(new EventState(handler, maxCheckInterval, convergence,
+                                        maxIterationCount, solver));
     }
 
     /** {@inheritDoc} */

src/main/java/org/apache/commons/math/ode/ODEIntegrator.java

@@ -19,6 +19,8 @@ package org.apache.commons.math.ode;
 
 import java.util.Collection;
 
+import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
+import org.apache.commons.math.analysis.solvers.BrentSolver;
 import org.apache.commons.math.ode.events.EventHandler;
 import org.apache.commons.math.ode.sampling.StepHandler;
 
@@ -62,7 +64,9 @@ public interface ODEIntegrator  {
      */
     void clearStepHandlers();
 
-    /** Add an event handler to the integrator.
+    /** Add an event handler to the integrator. Uses a {@link BrentSolver}
+     * with an absolute accuracy equal to the given convergence threshold,
+     * as root-finding algorithm to detect the state events.
      * @param handler event handler
      * @param maxCheckInterval maximal time interval between switching
      * function checks (this interval prevents missing sign changes in
@@ -76,6 +80,23 @@ public interface ODEIntegrator  {
     void addEventHandler(EventHandler handler, double maxCheckInterval,
                          double convergence, int maxIterationCount);
 
+    /** Add an event handler to the integrator.
+     * @param handler event handler
+     * @param maxCheckInterval maximal time interval between switching
+     * function checks (this interval prevents missing sign changes in
+     * case the integration steps becomes very large)
+     * @param convergence convergence threshold in the event time search
+     * @param maxIterationCount upper limit of the iteration count in
+     * the event time search
+     * @param solver The root-finding algorithm to use to detect the state
+     * events.
+     * @see #getEventHandlers()
+     * @see #clearEventHandlers()
+     */
+    void addEventHandler(EventHandler handler, double maxCheckInterval,
+                         double convergence, int maxIterationCount,
+                         UnivariateRealSolver solver);
+
     /** Get all the event handlers that have been added to the integrator.
      * @return an unmodifiable collection of the added events handlers
      * @see #addEventHandler(EventHandler, double, double, int)

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

@@ -19,7 +19,7 @@ package org.apache.commons.math.ode.events;
 
 import org.apache.commons.math.ConvergenceException;
 import org.apache.commons.math.analysis.UnivariateRealFunction;
-import org.apache.commons.math.analysis.solvers.BrentSolver;
+import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
 import org.apache.commons.math.exception.MathInternalError;
 import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.ode.sampling.StepInterpolator;
@@ -81,6 +81,9 @@ public class EventState {
     /** Next action indicator. */
     private int nextAction;
 
+    /** Root-finding algorithm to use to detect state events. */
+    private final UnivariateRealSolver solver;
+
     /** Simple constructor.
      * @param handler event handler
      * @param maxCheckInterval maximal time interval between switching
@@ -89,13 +92,16 @@ public class EventState {
      * @param convergence convergence threshold in the event time search
      * @param maxIterationCount upper limit of the iteration count in
      * the event time search
+     * @param solver Root-finding algorithm to use to detect state events
      */
     public EventState(final EventHandler handler, final double maxCheckInterval,
-                      final double convergence, final int maxIterationCount) {
+                      final double convergence, final int maxIterationCount,
+                      final UnivariateRealSolver solver) {
         this.handler           = handler;
         this.maxCheckInterval  = maxCheckInterval;
         this.convergence       = FastMath.abs(convergence);
         this.maxIterationCount = maxIterationCount;
+        this.solver            = solver;
 
         // some dummy values ...
         t0                = Double.NaN;
@@ -235,7 +241,6 @@ public class EventState {
                             }
                         }
                     };
-                    final BrentSolver solver = new BrentSolver(convergence);
 
                     if (ga * gb >= 0) {
                         // this is a corner case:

src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java

@@ -17,6 +17,7 @@
 
 package org.apache.commons.math.ode.nonstiff;
 
+import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
 import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.IntegratorException;
@@ -347,8 +348,10 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
   public void addEventHandler(final EventHandler function,
                               final double maxCheckInterval,
                               final double convergence,
-                              final int maxIterationCount) {
+                              final int maxIterationCount,
-    super.addEventHandler(function, maxCheckInterval, convergence, maxIterationCount);
+                              final UnivariateRealSolver solver) {
+    super.addEventHandler(function, maxCheckInterval, convergence,
+                          maxIterationCount, solver);
 
     // reinitialize the arrays
     initializeArrays();

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="3.0" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-586" due-to="Dennis Hendriks">
+        Added a way to specify a custom root solver to find events in ODE integration.
+      </action>
       <action dev="psteitz" type="fix" issue="MATH-540">
         Fixed error in javadoc describing Integer distribution inverse cumulative
         probability API.

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

@@ -19,6 +19,7 @@ package org.apache.commons.math.ode.events;
 
 
 import org.apache.commons.math.ConvergenceException;
+import org.apache.commons.math.analysis.solvers.BrentSolver;
 import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
 import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
@@ -47,7 +48,9 @@ public class EventStateTest {
         };
 
         final double tolerance = 0.1;
-        EventState es = new EventState(closeEventsGenerator, 1.5 * gap, tolerance, 10);
+        EventState es = new EventState(closeEventsGenerator, 1.5 * gap,
+                                       tolerance, 10,
+                                       new BrentSolver(tolerance));
 
         AbstractStepInterpolator interpolator =
             new DummyStepInterpolator(new double[0], new double[0], true);