Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fixed an error in handling of very close events during ODE integration 

JIRA: MATH-322

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@887794 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-12-06 23:04:55 +00:00
parent 335572b99e
commit f7c0b403d2
3 changed files with 98 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -19,6 +19,7 @@ package org.apache.commons.math.ode.events;
import org.apache.commons.math.ConvergenceException; import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException; 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.UnivariateRealFunction;
import org.apache.commons.math.analysis.solvers.BrentSolver; import org.apache.commons.math.analysis.solvers.BrentSolver;
import org.apache.commons.math.ode.DerivativeException; import org.apache.commons.math.ode.DerivativeException;
@ -187,6 +188,26 @@ public class EventState {
if (g0Positive ^ (gb >= 0)) { if (g0Positive ^ (gb >= 0)) {
// there is a sign change: an event is expected during this step // 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 // variation direction, with respect to the integration direction
increasing = gb >= ga; increasing = gb >= ga;
@ -205,16 +226,9 @@ public class EventState {
final BrentSolver solver = new BrentSolver(); final BrentSolver solver = new BrentSolver();
solver.setAbsoluteAccuracy(convergence); solver.setAbsoluteAccuracy(convergence);
solver.setMaximalIterationCount(maxIterationCount); solver.setMaximalIterationCount(maxIterationCount);
double root; final double root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta);
try { if ((Math.abs(root - ta) <= convergence) &&
root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta); (Math.abs(root - previousEventTime) <= convergence)) {
} 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))) {
// we have either found nothing or found (again ?) a past event, we simply ignore it // we have either found nothing or found (again ?) a past event, we simply ignore it
ta = tb; ta = tb;
ga = gb; ga = gb;

3
src/site/xdoc/changes.xml
View File

@ -39,6 +39,9 @@ The <action> type attribute can be add,update,fix,remove.
</properties> </properties>
<body> <body>
<release version="2.1" date="TBD" description="TBD"> <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"> <action dev="psteitz" type="fix" issue="MATH-305" due-to="Erik van Ingen">
Fixed an overflow error in MathUtils.distance that was causing KMeansPlusPlusClusterer Fixed an overflow error in MathUtils.distance that was causing KMeansPlusPlusClusterer
to fail with a NullPointerException when component distances between points to fail with a NullPointerException when component distances between points

71
src/test/java/org/apache/commons/math/ode/events/EventStateTest.java Normal file
View File

@ -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);
}
}