From 8ad7c0a62e3fbe8a882ac3e4dac030939f000a7e Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Wed, 29 Sep 2010 19:49:34 +0000 Subject: [PATCH] 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 --- .../ode/events/CombinedEventsManager.java | 9 ++-- .../commons/math/ode/events/EventState.java | 47 +++++++++++++++---- src/site/xdoc/changes.xml | 3 ++ src/site/xdoc/userguide/ode.xml | 7 ++- .../math/ode/events/EventStateTest.java | 11 +++-- 5 files changed, 61 insertions(+), 16 deletions(-) diff --git a/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java b/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java index 0831e4ab8..08d873862 100644 --- a/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java +++ b/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); diff --git a/src/main/java/org/apache/commons/math/ode/events/EventState.java b/src/main/java/org/apache/commons/math/ode/events/EventState.java index 173265084..e8c10f285 100644 --- a/src/main/java/org/apache/commons/math/ode/events/EventState.java +++ b/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 time 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. diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index b7a560856..909b85091 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -52,6 +52,9 @@ The type attribute can be add,update,fix,remove. If the output is not quite correct, check for invisible trailing spaces! --> + + Fixed an error preventing ODE solvers to be restarted after they have been stopped by a discrete event + Added new random number generators from the Well Equidistributed Long-period Linear (WELL). diff --git a/src/site/xdoc/userguide/ode.xml b/src/site/xdoc/userguide/ode.xml index 20db5dd09..318484cb9 100644 --- a/src/site/xdoc/userguide/ode.xml +++ b/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 t0 + to t0 + ε where ε 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.

When an event is triggered, the event time, current state and an indicator diff --git a/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java b/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java index e4baf3bd8..dff6d6b5b 100644 --- a/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java +++ b/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));