Repository: commons-math Updated Branches: refs/heads/master cf571ba2a -> c44bfe000
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/c44bfe00 Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/c44bfe00 Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/c44bfe00 Branch: refs/heads/master Commit: c44bfe000cd84d0b48401f8c7107bc37c09ec622 Parents: cf571ba Author: Luc Maisonobe <l...@apache.org> Authored: Tue May 19 13:18:32 2015 +0200 Committer: Luc Maisonobe <l...@apache.org> Committed: Tue May 19 13:18:32 2015 +0200 ---------------------------------------------------------------------- src/changes/changes.xml | 3 + .../commons/math4/ode/events/EventState.java | 13 +++- .../math4/ode/events/EventStateTest.java | 74 ++++++++++++++++---- 3 files changed, 76 insertions(+), 14 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-math/blob/c44bfe00/src/changes/changes.xml ---------------------------------------------------------------------- diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 66f0175..228a7f3 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -54,6 +54,9 @@ If the output is not quite correct, check for invisible trailing spaces! </release> <release version="4.0" date="XXXX-XX-XX" description=""> + <action dev="luc" type="fix" issue="MATH-1226"> <!-- backported to 3.6 --> + Fixed wrong event detection in case of close events pairs. + </action> <action dev="luc" type="add" > Reimplemented pow(double, double) in FastMath, for better accuracy in integral power cases and trying to fix erroneous JIT optimization again. http://git-wip-us.apache.org/repos/asf/commons-math/blob/c44bfe00/src/main/java/org/apache/commons/math4/ode/events/EventState.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math4/ode/events/EventState.java b/src/main/java/org/apache/commons/math4/ode/events/EventState.java index fe3039a..1908440 100644 --- a/src/main/java/org/apache/commons/math4/ode/events/EventState.java +++ b/src/main/java/org/apache/commons/math4/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/c44bfe00/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java ---------------------------------------------------------------------- diff --git a/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java b/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java index b2041dc..bbb2bc5 100644 --- a/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java +++ b/src/test/java/org/apache/commons/math4/ode/events/EventStateTest.java @@ -29,6 +29,7 @@ import org.apache.commons.math4.ode.SecondaryEquations; import org.apache.commons.math4.ode.events.EventHandler; import org.apache.commons.math4.ode.events.EventState; import org.apache.commons.math4.ode.nonstiff.DormandPrince853Integrator; +import org.apache.commons.math4.ode.nonstiff.LutherIntegrator; import org.apache.commons.math4.ode.sampling.AbstractStepInterpolator; import org.apache.commons.math4.ode.sampling.DummyStepInterpolator; import org.junit.Assert; @@ -43,21 +44,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() { @@ -226,4 +215,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; + } + + } + }
