Repository: commons-math Updated Branches: refs/heads/MATH_3_X 01673c2ce -> 777273e23
Fixed wrong event detection in case of close events pairs. JIRA: MATH-1226 Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/777273e2 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/777273e2 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/777273e2 Branch: refs/heads/MATH_3_X Commit: 777273e23f825d659695a4d69fca29eeba31cf37 Parents: 01673c2 Author: Luc Maisonobe <l...@apache.org> Authored: Tue May 19 13:26:33 2015 +0200 Committer: Luc Maisonobe <l...@apache.org> Committed: Tue May 19 13:26:33 2015 +0200 ---------------------------------------------------------------------- src/changes/changes.xml | 5 +- .../commons/math3/ode/events/EventState.java | 13 +++- .../math3/ode/events/EventStateTest.java | 74 ++++++++++++++++---- 3 files changed, 77 insertions(+), 15 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/777273e2/src/changes/changes.xml ---------------------------------------------------------------------- diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 3f1920b..d38f022 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -51,7 +51,10 @@ If the output is not quite correct, check for invisible trailing spaces! </properties> <body> <release version="3.6" date="XXXX-XX-XX" description=""> - <action dev="tn" type="fix"> <!-- backported to 3.6 --> + <action dev="luc" type="fix" issue="MATH-1226"> + Fixed wrong event detection in case of close events pairs. + </action> + <action dev="tn" type="fix"> Fix potential branching errors in "FastMath#pow(double, double)" when passing special values, i.e. infinity, due to erroneous JIT optimization. </action> http://git-wip-us.apache.org/repos/asf/commons-math/blob/777273e2/src/main/java/org/apache/commons/math3/ode/events/EventState.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/ode/events/EventState.java b/src/main/java/org/apache/commons/math3/ode/events/EventState.java index 76dae72..6d34266 100644 --- a/src/main/java/org/apache/commons/math3/ode/events/EventState.java +++ b/src/main/java/org/apache/commons/math3/ode/events/EventState.java @@ -296,7 +296,18 @@ public class EventState { ta = forward ? ta + convergence : ta - convergence; ga = f.value(ta); } while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb))); - --i; + + if (forward ^ (ta >= tb)) { + // we were able to skip this spurious root + --i; + } else { + // we can't avoid this root before the end of the step, + // we have to handle it despite it is close to the former one + // maybe we have two very close roots + pendingEventTime = root; + pendingEvent = true; + return true; + } } else if (Double.isNaN(previousEventTime) || (FastMath.abs(previousEventTime - root) > convergence)) { pendingEventTime = root; http://git-wip-us.apache.org/repos/asf/commons-math/blob/777273e2/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java b/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java index 6468fdc..4a8091f 100644 --- a/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java +++ b/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java @@ -27,6 +27,7 @@ import org.apache.commons.math3.ode.ExpandableStatefulODE; import org.apache.commons.math3.ode.FirstOrderDifferentialEquations; import org.apache.commons.math3.ode.SecondaryEquations; import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator; +import org.apache.commons.math3.ode.nonstiff.LutherIntegrator; import org.apache.commons.math3.ode.sampling.AbstractStepInterpolator; import org.apache.commons.math3.ode.sampling.DummyStepInterpolator; import org.junit.Assert; @@ -41,21 +42,9 @@ public class EventStateTest { final double r1 = 90.0; final double r2 = 135.0; final double gap = r2 - r1; - EventHandler closeEventsGenerator = new EventHandler() { - public void init(double t0, double[] y0, double t) { - } - public void resetState(double t, double[] y) { - } - public double g(double t, double[] y) { - return (t - r1) * (r2 - t); - } - public Action eventOccurred(double t, double[] y, boolean increasing) { - return Action.CONTINUE; - } - }; final double tolerance = 0.1; - EventState es = new EventState(closeEventsGenerator, 1.5 * gap, + EventState es = new EventState(new CloseEventsGenerator(r1, r2), 1.5 * gap, tolerance, 100, new BrentSolver(tolerance)); es.setExpandable(new ExpandableStatefulODE(new FirstOrderDifferentialEquations() { @@ -224,4 +213,63 @@ public class EventStateTest { } + @Test + public void testEventsCloserThanThreshold() + throws DimensionMismatchException, NumberIsTooSmallException, + MaxCountExceededException, NoBracketingException { + + FirstOrderDifferentialEquations equation = new FirstOrderDifferentialEquations() { + + public int getDimension() { + return 1; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + yDot[0] = 1.0; + } + }; + + LutherIntegrator integrator = new LutherIntegrator(20.0); + CloseEventsGenerator eventsGenerator = + new CloseEventsGenerator(9.0 - 1.0 / 128, 9.0 + 1.0 / 128); + integrator.addEventHandler(eventsGenerator, 1.0, 0.02, 1000); + double[] y = new double[1]; + double tEnd = integrator.integrate(equation, 0.0, y, 100.0, y); + Assert.assertEquals( 2, eventsGenerator.getCount()); + Assert.assertEquals( 9.0 + 1.0 / 128, tEnd, 1.0 / 32.0); + + } + + private class CloseEventsGenerator implements EventHandler { + + final double r1; + final double r2; + int count; + + public CloseEventsGenerator(final double r1, final double r2) { + this.r1 = r1; + this.r2 = r2; + this.count = 0; + } + + public void init(double t0, double[] y0, double t) { + } + + public void resetState(double t, double[] y) { + } + + public double g(double t, double[] y) { + return (t - r1) * (r2 - t); + } + + public Action eventOccurred(double t, double[] y, boolean increasing) { + return ++count < 2 ? Action.CONTINUE : Action.STOP; + } + + public int getCount() { + return count; + } + + } + }
