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 2009/12/07 00:04:55 UTC
svn commit: r887794 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/ode/events/EventState.java
site/xdoc/changes.xml test/java/org/apache/commons/math/ode/events/
test/java/org/apache/commons/math/ode/events/EventStateTest.java
Author: luc
Date: Sun Dec 6 23:04:55 2009
New Revision: 887794
URL: http://svn.apache.org/viewvc?rev=887794&view=rev
Log:
Fixed an error in handling of very close events during ODE integration
JIRA: MATH-322
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/
commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java (with props)
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
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=887794&r1=887793&r2=887794&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 Sun Dec 6 23:04:55 2009
@@ -19,6 +19,7 @@
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
+import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.solvers.BrentSolver;
import org.apache.commons.math.ode.DerivativeException;
@@ -187,6 +188,26 @@
if (g0Positive ^ (gb >= 0)) {
// there is a sign change: an event is expected during this step
+ if (ga * gb > 0) {
+ // this is a corner case:
+ // - there was an event near ta,
+ // - there is another event between ta and tb
+ // - when ta was computed, convergence was reached on the "wrong side" of the interval
+ // this implies that the real sign of ga is the same as gb, so we need to slightly
+ // shift ta to make sure ga and gb get opposite signs and the solver won't complain
+ // about bracketing
+ final double epsilon = (forward ? 0.25 : -0.25) * convergence;
+ for (int k = 0; (k < 4) && (ga * gb > 0); ++k) {
+ ta += epsilon;
+ interpolator.setInterpolatedTime(ta);
+ ga = handler.g(ta, interpolator.getInterpolatedState());
+ }
+ if (ga * gb > 0) {
+ // this should never happen
+ throw MathRuntimeException.createInternalError(null);
+ }
+ }
+
// variation direction, with respect to the integration direction
increasing = gb >= ga;
@@ -205,16 +226,9 @@
final BrentSolver solver = new BrentSolver();
solver.setAbsoluteAccuracy(convergence);
solver.setMaximalIterationCount(maxIterationCount);
- double root;
- try {
- root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta);
- } catch (IllegalArgumentException iae) {
- // the interval did not really bracket a root
- root = Double.NaN;
- }
- if (Double.isNaN(root) ||
- ((Math.abs(root - ta) <= convergence) &&
- (Math.abs(root - previousEventTime) <= convergence))) {
+ final double root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta);
+ if ((Math.abs(root - ta) <= convergence) &&
+ (Math.abs(root - previousEventTime) <= convergence)) {
// we have either found nothing or found (again ?) a past event, we simply ignore it
ta = tb;
ga = gb;
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=887794&r1=887793&r2=887794&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Sun Dec 6 23:04:55 2009
@@ -39,6 +39,9 @@
</properties>
<body>
<release version="2.1" date="TBD" description="TBD">
+ <action dev="luc" type="fix" issue="MATH-322" >
+ Fixed an error in handling very close events in ODE integration.
+ </action>
<action dev="psteitz" type="fix" issue="MATH-305" due-to="Erik van Ingen">
Fixed an overflow error in MathUtils.distance that was causing KMeansPlusPlusClusterer
to fail with a NullPointerException when component distances between points
Added: 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=887794&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java Sun Dec 6 23:04:55 2009
@@ -0,0 +1,71 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.events;
+
+import junit.framework.Assert;
+
+import org.apache.commons.math.ConvergenceException;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
+import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
+import org.junit.Test;
+
+public class EventStateTest {
+
+ // JIRA: MATH-322
+ @Test
+ public void closeEvents()
+ throws EventException, ConvergenceException, DerivativeException {
+
+ final double r1 = 90.0;
+ final double r2 = 135.0;
+ final double gap = r2 - r1;
+ EventHandler closeEventsGenerator = new EventHandler() {
+ public void resetState(double t, double[] y) {
+ }
+ public double g(double t, double[] y) {
+ return (t - r1) * (r2 - t);
+ }
+ public int eventOccurred(double t, double[] y, boolean increasing) {
+ return CONTINUE;
+ }
+ };
+
+ 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], true);
+ interpolator.storeTime(t0);
+
+ interpolator.shift();
+ interpolator.storeTime(0.5 * (r1 + r2));
+ Assert.assertTrue(es.evaluateStep(interpolator));
+ Assert.assertEquals(r1, es.getEventTime(), tolerance);
+ es.stepAccepted(es.getEventTime(), new double[0]);
+
+ interpolator.shift();
+ interpolator.storeTime(r2 + 0.4 * gap);
+ Assert.assertTrue(es.evaluateStep(interpolator));
+ Assert.assertEquals(r2, es.getEventTime(), tolerance);
+
+ }
+
+}
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/events/EventStateTest.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision