commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject 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
Date Sun, 06 Dec 2009 23:04:55 GMT
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



Mime
View raw message