diff --git a/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java b/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java index d58fa3f45..96893e06a 100644 --- a/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java +++ b/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java @@ -22,6 +22,7 @@ import org.apache.commons.math.ode.AbstractIntegrator; import org.apache.commons.math.ode.DerivativeException; import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.IntegratorException; +import org.apache.commons.math.ode.events.CombinedEventsManager; import org.apache.commons.math.ode.sampling.AbstractStepInterpolator; import org.apache.commons.math.ode.sampling.DummyStepInterpolator; import org.apache.commons.math.ode.sampling.StepHandler; @@ -106,19 +107,20 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator { } interpolator.storeTime(t0); - // recompute the step - long nbStep = Math.max(1l, Math.abs(Math.round((t - t0) / step))); - boolean lastStep = false; + // set up integration control objects stepStart = t0; - stepSize = (t - t0) / nbStep; + stepSize = step; for (StepHandler handler : stepHandlers) { handler.reset(); } - for (long i = 0; ! lastStep; ++i) { + CombinedEventsManager manager = addEndTimeChecker(t, eventsHandlersManager); + boolean lastStep = false; + + // main integration loop + while (!lastStep) { interpolator.shift(); - boolean needUpdate = false; for (boolean loop = true; loop;) { // first stage @@ -148,11 +150,10 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator { yTmp[j] = y[j] + stepSize * sum; } - // Discrete events handling + // discrete events handling interpolator.storeTime(stepStart + stepSize); - if (eventsHandlersManager.evaluateStep(interpolator)) { - needUpdate = true; - stepSize = eventsHandlersManager.getEventTime() - stepStart; + if (manager.evaluateStep(interpolator)) { + stepSize = manager.getEventTime() - stepStart; } else { loop = false; } @@ -162,12 +163,8 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator { // the step has been accepted final double nextStep = stepStart + stepSize; System.arraycopy(yTmp, 0, y, 0, y0.length); - eventsHandlersManager.stepAccepted(nextStep, y); - if (eventsHandlersManager.stop()) { - lastStep = true; - } else { - lastStep = (i == (nbStep - 1)); - } + manager.stepAccepted(nextStep, y); + lastStep = manager.stop(); // provide the step data to the step handler interpolator.storeTime(nextStep); @@ -176,19 +173,14 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator { } stepStart = nextStep; - if (eventsHandlersManager.reset(stepStart, y) && ! lastStep) { + if (manager.reset(stepStart, y) && ! lastStep) { // some events handler has triggered changes that // invalidate the derivatives, we need to recompute them equations.computeDerivatives(stepStart, y, yDotK[0]); } - if (needUpdate) { - // an event handler has changed the step - // we need to recompute stepsize - nbStep = Math.max(1l, Math.abs(Math.round((t - stepStart) / step))); - stepSize = (t - stepStart) / nbStep; - i = -1; - } + // make sure step size is set to default before next step + stepSize = step; } diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index 4a7af0f2a..13394b7c7 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -39,6 +39,10 @@ The type attribute can be add,update,fix,remove. + + Replaced size adjustment of all steps of fixed steps Runge-Kutta integrators by + a truncation of the last step only. + The ODE integrators now support several step handlers at once, instead of just one. This is more consistent with event handlers management. diff --git a/src/test/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java b/src/test/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java index 1bf1ecd0a..036c7b1c0 100644 --- a/src/test/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java +++ b/src/test/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java @@ -20,6 +20,7 @@ package org.apache.commons.math.ode.nonstiff; import junit.framework.*; import org.apache.commons.math.ode.DerivativeException; +import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderIntegrator; import org.apache.commons.math.ode.IntegratorException; import org.apache.commons.math.ode.events.EventHandler; @@ -191,6 +192,36 @@ public class ClassicalRungeKuttaIntegratorTest private TestProblem3 pb; } + public void testStepSize() + throws DerivativeException, IntegratorException { + final double step = 1.23456; + FirstOrderIntegrator integ = new ClassicalRungeKuttaIntegrator(step); + integ.addStepHandler(new StepHandler() { + private static final long serialVersionUID = 0L; + public void handleStep(StepInterpolator interpolator, boolean isLast) { + if (! isLast) { + assertEquals(step, + interpolator.getCurrentTime() - interpolator.getPreviousTime(), + 1.0e-12); + } + } + public boolean requiresDenseOutput() { + return false; + } + public void reset() { + } + }); + integ.integrate(new FirstOrderDifferentialEquations() { + private static final long serialVersionUID = 0L; + public void computeDerivatives(double t, double[] y, double[] dot) { + dot[0] = 1.0; + } + public int getDimension() { + return 1; + } + }, 0.0, new double[] { 0.0 }, 5.0, new double[1]); + } + public static Test suite() { return new TestSuite(ClassicalRungeKuttaIntegratorTest.class); } diff --git a/src/test/org/apache/commons/math/ode/nonstiff/EulerIntegratorTest.java b/src/test/org/apache/commons/math/ode/nonstiff/EulerIntegratorTest.java index a8809415d..14ddebe28 100644 --- a/src/test/org/apache/commons/math/ode/nonstiff/EulerIntegratorTest.java +++ b/src/test/org/apache/commons/math/ode/nonstiff/EulerIntegratorTest.java @@ -20,10 +20,13 @@ package org.apache.commons.math.ode.nonstiff; import junit.framework.*; import org.apache.commons.math.ode.DerivativeException; +import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderIntegrator; import org.apache.commons.math.ode.IntegratorException; import org.apache.commons.math.ode.events.EventHandler; import org.apache.commons.math.ode.nonstiff.EulerIntegrator; +import org.apache.commons.math.ode.sampling.StepHandler; +import org.apache.commons.math.ode.sampling.StepInterpolator; public class EulerIntegratorTest extends TestCase { @@ -123,7 +126,37 @@ public class EulerIntegratorTest assertEquals(0, handler.getMaximalTimeError(), 1.0e-12); } - + + public void testStepSize() + throws DerivativeException, IntegratorException { + final double step = 1.23456; + FirstOrderIntegrator integ = new EulerIntegrator(step); + integ.addStepHandler(new StepHandler() { + private static final long serialVersionUID = 0L; + public void handleStep(StepInterpolator interpolator, boolean isLast) { + if (! isLast) { + assertEquals(step, + interpolator.getCurrentTime() - interpolator.getPreviousTime(), + 1.0e-12); + } + } + public boolean requiresDenseOutput() { + return false; + } + public void reset() { + } + }); + integ.integrate(new FirstOrderDifferentialEquations() { + private static final long serialVersionUID = 0L; + public void computeDerivatives(double t, double[] y, double[] dot) { + dot[0] = 1.0; + } + public int getDimension() { + return 1; + } + }, 0.0, new double[] { 0.0 }, 5.0, new double[1]); + } + public static Test suite() { return new TestSuite(EulerIntegratorTest.class); } diff --git a/src/test/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java b/src/test/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java index 251344533..e5e6344de 100644 --- a/src/test/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java +++ b/src/test/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java @@ -20,6 +20,7 @@ package org.apache.commons.math.ode.nonstiff; import junit.framework.*; import org.apache.commons.math.ode.DerivativeException; +import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderIntegrator; import org.apache.commons.math.ode.IntegratorException; import org.apache.commons.math.ode.events.EventHandler; @@ -181,6 +182,36 @@ public class GillIntegratorTest private TestProblem3 pb; } + public void testStepSize() + throws DerivativeException, IntegratorException { + final double step = 1.23456; + FirstOrderIntegrator integ = new GillIntegrator(step); + integ.addStepHandler(new StepHandler() { + private static final long serialVersionUID = 0L; + public void handleStep(StepInterpolator interpolator, boolean isLast) { + if (! isLast) { + assertEquals(step, + interpolator.getCurrentTime() - interpolator.getPreviousTime(), + 1.0e-12); + } + } + public boolean requiresDenseOutput() { + return false; + } + public void reset() { + } + }); + integ.integrate(new FirstOrderDifferentialEquations() { + private static final long serialVersionUID = 0L; + public void computeDerivatives(double t, double[] y, double[] dot) { + dot[0] = 1.0; + } + public int getDimension() { + return 1; + } + }, 0.0, new double[] { 0.0 }, 5.0, new double[1]); + } + public static Test suite() { return new TestSuite(GillIntegratorTest.class); } diff --git a/src/test/org/apache/commons/math/ode/nonstiff/MidpointIntegratorTest.java b/src/test/org/apache/commons/math/ode/nonstiff/MidpointIntegratorTest.java index f1e2b2942..26fbaffe6 100644 --- a/src/test/org/apache/commons/math/ode/nonstiff/MidpointIntegratorTest.java +++ b/src/test/org/apache/commons/math/ode/nonstiff/MidpointIntegratorTest.java @@ -20,10 +20,13 @@ package org.apache.commons.math.ode.nonstiff; import junit.framework.*; import org.apache.commons.math.ode.DerivativeException; +import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderIntegrator; import org.apache.commons.math.ode.IntegratorException; import org.apache.commons.math.ode.events.EventHandler; import org.apache.commons.math.ode.nonstiff.MidpointIntegrator; +import org.apache.commons.math.ode.sampling.StepHandler; +import org.apache.commons.math.ode.sampling.StepInterpolator; public class MidpointIntegratorTest extends TestCase { @@ -124,6 +127,36 @@ public class MidpointIntegratorTest } + public void testStepSize() + throws DerivativeException, IntegratorException { + final double step = 1.23456; + FirstOrderIntegrator integ = new MidpointIntegrator(step); + integ.addStepHandler(new StepHandler() { + private static final long serialVersionUID = 0L; + public void handleStep(StepInterpolator interpolator, boolean isLast) { + if (! isLast) { + assertEquals(step, + interpolator.getCurrentTime() - interpolator.getPreviousTime(), + 1.0e-12); + } + } + public boolean requiresDenseOutput() { + return false; + } + public void reset() { + } + }); + integ.integrate(new FirstOrderDifferentialEquations() { + private static final long serialVersionUID = 0L; + public void computeDerivatives(double t, double[] y, double[] dot) { + dot[0] = 1.0; + } + public int getDimension() { + return 1; + } + }, 0.0, new double[] { 0.0 }, 5.0, new double[1]); + } + public static Test suite() { return new TestSuite(MidpointIntegratorTest.class); } diff --git a/src/test/org/apache/commons/math/ode/nonstiff/StepInterpolatorAbstractTest.java b/src/test/org/apache/commons/math/ode/nonstiff/StepInterpolatorAbstractTest.java index 11cfccc13..8a61c0191 100644 --- a/src/test/org/apache/commons/math/ode/nonstiff/StepInterpolatorAbstractTest.java +++ b/src/test/org/apache/commons/math/ode/nonstiff/StepInterpolatorAbstractTest.java @@ -48,6 +48,10 @@ public class StepInterpolatorAbstractTest extends TestCase { final double h = 0.001 * (interpolator.getCurrentTime() - interpolator.getPreviousTime()); final double t = interpolator.getCurrentTime() - 300 * h; + if (Math.abs(h) < 10 * Math.ulp(t)) { + return; + } + interpolator.setInterpolatedTime(t - 4 * h); final double[] yM4h = interpolator.getInterpolatedState().clone(); interpolator.setInterpolatedTime(t - 3 * h); @@ -73,6 +77,9 @@ public class StepInterpolatorAbstractTest extends TestCase { 32 * (yP3h[i] - yM3h[i]) + -168 * (yP2h[i] - yM2h[i]) + 672 * (yP1h[i] - yM1h[i])) / (840 * h); + if (Math.abs(approYDot - yDot[i]) >= threshold) { + System.out.println("gotcha!"); + } assertEquals(approYDot, yDot[i], threshold); } diff --git a/src/test/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java b/src/test/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java index d865efd3b..dd4c85c62 100644 --- a/src/test/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java +++ b/src/test/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java @@ -20,6 +20,7 @@ package org.apache.commons.math.ode.nonstiff; import junit.framework.*; import org.apache.commons.math.ode.DerivativeException; +import org.apache.commons.math.ode.FirstOrderDifferentialEquations; import org.apache.commons.math.ode.FirstOrderIntegrator; import org.apache.commons.math.ode.IntegratorException; import org.apache.commons.math.ode.events.EventHandler; @@ -178,6 +179,36 @@ public class ThreeEighthesIntegratorTest } + public void testStepSize() + throws DerivativeException, IntegratorException { + final double step = 1.23456; + FirstOrderIntegrator integ = new ThreeEighthesIntegrator(step); + integ.addStepHandler(new StepHandler() { + private static final long serialVersionUID = 0L; + public void handleStep(StepInterpolator interpolator, boolean isLast) { + if (! isLast) { + assertEquals(step, + interpolator.getCurrentTime() - interpolator.getPreviousTime(), + 1.0e-12); + } + } + public boolean requiresDenseOutput() { + return false; + } + public void reset() { + } + }); + integ.integrate(new FirstOrderDifferentialEquations() { + private static final long serialVersionUID = 0L; + public void computeDerivatives(double t, double[] y, double[] dot) { + dot[0] = 1.0; + } + public int getDimension() { + return 1; + } + }, 0.0, new double[] { 0.0 }, 5.0, new double[1]); + } + public static Test suite() { return new TestSuite(ThreeEighthesIntegratorTest.class); }