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> + ε 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.
</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));
