diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java index 42796c755..baf236a52 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java @@ -48,7 +48,7 @@ class ClassicalRungeKuttaStepInterpolator extends RungeKuttaStepInterpolator { /** Serializable version identifier. */ - private static final long serialVersionUID = 20110928L; + private static final long serialVersionUID = 20111120L; /** Simple constructor. * This constructor builds an instance that is not usable yet, the diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java index af32c4c01..1dac8ab30 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java @@ -72,7 +72,7 @@ class DormandPrince54StepInterpolator private static final double D6 = 69997945.0 / 29380423.0; /** Serializable version identifier. */ - private static final long serialVersionUID = 20110928L; + private static final long serialVersionUID = 20111120L; /** First vector for interpolation. */ private double[] v1; diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java index 91442bbf2..64dd7d77b 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java @@ -39,7 +39,7 @@ class DormandPrince853StepInterpolator extends RungeKuttaStepInterpolator { /** Serializable version identifier. */ - private static final long serialVersionUID = 20110928L; + private static final long serialVersionUID = 20111120L; /** Propagation weights, element 1. */ private static final double B_01 = 104257.0 / 1920240.0; @@ -368,18 +368,34 @@ class DormandPrince853StepInterpolator final double dot5 = theta2 * (3 + theta * (-12 + theta * (15 - 6 * theta))); final double dot6 = theta2 * theta * (4 + theta * (-15 + theta * (18 - 7 * theta))); - for (int i = 0; i < interpolatedState.length; ++i) { - interpolatedState[i] = currentState[i] - - oneMinusThetaH * (v[0][i] - - theta * (v[1][i] + - theta * (v[2][i] + - eta * (v[3][i] + - theta * (v[4][i] + - eta * (v[5][i] + - theta * (v[6][i]))))))); - interpolatedDerivatives[i] = v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] + - dot3 * v[3][i] + dot4 * v[4][i] + - dot5 * v[5][i] + dot6 * v[6][i]; + if ((previousState != null) && (theta <= 0.5)) { + for (int i = 0; i < interpolatedState.length; ++i) { + interpolatedState[i] = previousState[i] + + theta * h * (v[0][i] + + eta * (v[1][i] + + theta * (v[2][i] + + eta * (v[3][i] + + theta * (v[4][i] + + eta * (v[5][i] + + theta * (v[6][i]))))))); + interpolatedDerivatives[i] = v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] + + dot3 * v[3][i] + dot4 * v[4][i] + + dot5 * v[5][i] + dot6 * v[6][i]; + } + } else { + for (int i = 0; i < interpolatedState.length; ++i) { + interpolatedState[i] = currentState[i] - + oneMinusThetaH * (v[0][i] - + theta * (v[1][i] + + theta * (v[2][i] + + eta * (v[3][i] + + theta * (v[4][i] + + eta * (v[5][i] + + theta * (v[6][i]))))))); + interpolatedDerivatives[i] = v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] + + dot3 * v[3][i] + dot4 * v[4][i] + + dot5 * v[5][i] + dot6 * v[6][i]; + } } } diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java index 87fd71602..34d2c0073 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java @@ -202,7 +202,7 @@ public abstract class EmbeddedRungeKuttaIntegrator final double[] y = y0.clone(); final int stages = c.length + 1; final double[][] yDotK = new double[stages][y.length]; - final double[] yTmp = new double[y.length]; + final double[] yTmp = y0.clone(); final double[] yDotTmp = new double[y.length]; // set up an interpolator sharing the integrator arrays @@ -294,6 +294,7 @@ public abstract class EmbeddedRungeKuttaIntegrator System.arraycopy(yTmp, 0, y, 0, y0.length); System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length); stepStart = acceptStep(interpolator, y, yDotTmp, t); + System.arraycopy(y, 0, yTmp, 0, y.length); if (!isLastStep) { diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java index 99a604f3b..d2807a465 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java @@ -42,7 +42,7 @@ class EulerStepInterpolator extends RungeKuttaStepInterpolator { /** Serializable version identifier. */ - private static final long serialVersionUID = 20110928L; + private static final long serialVersionUID = 20111120L; /** Simple constructor. * This constructor builds an instance that is not usable yet, the diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java index e956b200e..537775556 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java @@ -54,7 +54,7 @@ class GillStepInterpolator private static final double TWO_PLUS_SQRT_2 = 2 + FastMath.sqrt(2.0); /** Serializable version identifier. */ - private static final long serialVersionUID = 20110928L; + private static final long serialVersionUID = 20111120L; /** Simple constructor. * This constructor builds an instance that is not usable yet, the diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java index ba112d449..b928b53ab 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java @@ -33,7 +33,7 @@ class HighamHall54StepInterpolator extends RungeKuttaStepInterpolator { /** Serializable version identifier */ - private static final long serialVersionUID = 20110928L; + private static final long serialVersionUID = 20111120L; /** Simple constructor. * This constructor builds an instance that is not usable yet, the diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java index 1c7648365..25f77c543 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java @@ -44,7 +44,7 @@ class MidpointStepInterpolator extends RungeKuttaStepInterpolator { /** Serializable version identifier */ - private static final long serialVersionUID = 20110928L; + private static final long serialVersionUID = 20111120L; /** Simple constructor. * This constructor builds an instance that is not usable yet, the diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java index fb82d0c9e..c51620c17 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java @@ -107,7 +107,7 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator { for (int i = 0; i < stages; ++i) { yDotK [i] = new double[y0.length]; } - final double[] yTmp = new double[y0.length]; + final double[] yTmp = y0.clone(); final double[] yDotTmp = new double[y0.length]; // set up an interpolator sharing the integrator arrays diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java index 987dfb1f3..55146c32c 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java @@ -38,6 +38,9 @@ import org.apache.commons.math.ode.sampling.AbstractStepInterpolator; abstract class RungeKuttaStepInterpolator extends AbstractStepInterpolator { + /** Previous state. */ + protected double[] previousState; + /** Slopes at the intermediate points */ protected double[][] yDotK; @@ -55,9 +58,9 @@ abstract class RungeKuttaStepInterpolator * uninitialized model and latter initializing the copy. */ protected RungeKuttaStepInterpolator() { - super(); - yDotK = null; - integrator = null; + previousState = null; + yDotK = null; + integrator = null; } /** Copy constructor. @@ -82,16 +85,16 @@ abstract class RungeKuttaStepInterpolator super(interpolator); if (interpolator.currentState != null) { - final int dimension = currentState.length; + + previousState = interpolator.previousState.clone(); yDotK = new double[interpolator.yDotK.length][]; for (int k = 0; k < interpolator.yDotK.length; ++k) { - yDotK[k] = new double[dimension]; - System.arraycopy(interpolator.yDotK[k], 0, - yDotK[k], 0, dimension); + yDotK[k] = interpolator.yDotK[k].clone(); } } else { + previousState = null; yDotK = null; } @@ -129,10 +132,18 @@ abstract class RungeKuttaStepInterpolator final EquationsMapper primaryMapper, final EquationsMapper[] secondaryMappers) { reinitialize(y, forward, primaryMapper, secondaryMappers); + this.previousState = null; this.yDotK = yDotArray; this.integrator = rkIntegrator; } + /** {@inheritDoc} */ + @Override + public void shift() { + previousState = currentState.clone(); + super.shift(); + } + /** {@inheritDoc} */ @Override public void writeExternal(final ObjectOutput out) @@ -143,6 +154,10 @@ abstract class RungeKuttaStepInterpolator // save the local attributes final int n = (currentState == null) ? -1 : currentState.length; + for (int i = 0; i < n; ++i) { + out.writeDouble(previousState[i]); + } + final int kMax = (yDotK == null) ? -1 : yDotK.length; out.writeInt(kMax); for (int k = 0; k < kMax; ++k) { @@ -165,6 +180,15 @@ abstract class RungeKuttaStepInterpolator // read the local attributes final int n = (currentState == null) ? -1 : currentState.length; + if (n < 0) { + previousState = null; + } else { + previousState = new double[n]; + for (int i = 0; i < n; ++i) { + previousState[i] = in.readDouble(); + } + } + final int kMax = in.readInt(); yDotK = (kMax < 0) ? null : new double[kMax][]; for (int k = 0; k < kMax; ++k) { diff --git a/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java b/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java index 731ec444c..fdc9d757c 100644 --- a/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java +++ b/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java @@ -49,7 +49,7 @@ class ThreeEighthesStepInterpolator extends RungeKuttaStepInterpolator { /** Serializable version identifier */ - private static final long serialVersionUID = 20110928L; + private static final long serialVersionUID = 20111120L; /** Simple constructor. * This constructor builds an instance that is not usable yet, the diff --git a/src/test/java/org/apache/commons/math/ode/events/ReappearingEventTest.java b/src/test/java/org/apache/commons/math/ode/events/ReappearingEventTest.java new file mode 100644 index 000000000..239854a92 --- /dev/null +++ b/src/test/java/org/apache/commons/math/ode/events/ReappearingEventTest.java @@ -0,0 +1,82 @@ +/* + * 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 static org.junit.Assert.assertEquals; + +import java.util.Arrays; + +import org.apache.commons.math.analysis.solvers.PegasusSolver; +import org.apache.commons.math.ode.FirstOrderDifferentialEquations; +import org.apache.commons.math.ode.FirstOrderIntegrator; +import org.apache.commons.math.ode.events.EventHandler; +import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator; +import org.apache.commons.math.ode.nonstiff.GraggBulirschStoerIntegrator; +import org.junit.Test; + +public class ReappearingEventTest { + @Test + public void testDormandPrince() { + double tEnd = test(1); + assertEquals(10.0, tEnd, 1e-7); + } + + @Test + public void testGragg() { + double tEnd = test(2); + assertEquals(10.0, tEnd, 1e-7); + } + + public double test(int integratorType) { + double e = 1e-15; + FirstOrderIntegrator integrator; + integrator = (integratorType == 1) + ? new DormandPrince853Integrator(e, 100.0, 1e-7, 1e-7) + : new GraggBulirschStoerIntegrator(e, 100.0, 1e-7, 1e-7); + PegasusSolver rootSolver = new PegasusSolver(e, e); + integrator.addEventHandler(new Event(), 0.1, e, 1000, rootSolver); + double t0 = 6.0; + double tEnd = 10.0; + double[] y = {2.0, 2.0, 2.0, 4.0, 2.0, 7.0, 15.0}; + return integrator.integrate(new Ode(), t0, y, tEnd, y); + } + + private static class Ode implements FirstOrderDifferentialEquations { + public int getDimension() { + return 7; + } + + public void computeDerivatives(double t, double[] y, double[] yDot) { + Arrays.fill(yDot, 1.0); + } + } + + /** State events for this unit test. */ + protected static class Event implements EventHandler { + public double g(double t, double[] y) { + return y[6] - 15.0; + } + + public Action eventOccurred(double t, double[] y, boolean increasing) { + return Action.STOP; + } + + public void resetState(double t, double[] y) { + // Never called. + } + } +} diff --git a/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java b/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java index 57bbf7a9f..13b4eb28f 100644 --- a/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java @@ -61,8 +61,8 @@ public class ClassicalRungeKuttaStepInterpolatorTest { oos.writeObject(handler); } - Assert.assertTrue(bos.size () > 750000); - Assert.assertTrue(bos.size () < 800000); + Assert.assertTrue(bos.size () > 880000); + Assert.assertTrue(bos.size () < 900000); ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray()); ObjectInputStream ois = new ObjectInputStream(bis); diff --git a/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolatorTest.java b/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolatorTest.java index 51cba08e7..22ec5f301 100644 --- a/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolatorTest.java @@ -73,8 +73,8 @@ public class DormandPrince54StepInterpolatorTest { oos.writeObject(handler); } - Assert.assertTrue(bos.size () > 125000); - Assert.assertTrue(bos.size () < 130000); + Assert.assertTrue(bos.size () > 135000); + Assert.assertTrue(bos.size () < 145000); ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray()); ObjectInputStream ois = new ObjectInputStream(bis); diff --git a/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolatorTest.java b/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolatorTest.java index 0a8ea97c8..e2c0f45ba 100644 --- a/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolatorTest.java @@ -73,8 +73,8 @@ public class DormandPrince853StepInterpolatorTest { oos.writeObject(handler); } - Assert.assertTrue(bos.size () > 85000); - Assert.assertTrue(bos.size () < 95000); + Assert.assertTrue(bos.size () > 90000); + Assert.assertTrue(bos.size () < 100000); ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray()); ObjectInputStream ois = new ObjectInputStream(bis); diff --git a/src/test/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolatorTest.java b/src/test/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolatorTest.java index 7cbb992f1..fdf27dfd4 100644 --- a/src/test/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolatorTest.java @@ -61,8 +61,8 @@ public class GillStepInterpolatorTest { oos.writeObject(handler); } - Assert.assertTrue(bos.size () > 750000); - Assert.assertTrue(bos.size () < 800000); + Assert.assertTrue(bos.size () > 880000); + Assert.assertTrue(bos.size () < 900000); ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray()); ObjectInputStream ois = new ObjectInputStream(bis); diff --git a/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolatorTest.java b/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolatorTest.java index 691ce4a68..42787d2e3 100644 --- a/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolatorTest.java @@ -73,8 +73,8 @@ public class HighamHall54StepInterpolatorTest { oos.writeObject(handler); } - Assert.assertTrue(bos.size () > 170000); - Assert.assertTrue(bos.size () < 175000); + Assert.assertTrue(bos.size () > 185000); + Assert.assertTrue(bos.size () < 195000); ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray()); ObjectInputStream ois = new ObjectInputStream(bis); diff --git a/src/test/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolatorTest.java b/src/test/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolatorTest.java index 3b71c961c..5efd99033 100644 --- a/src/test/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolatorTest.java @@ -62,8 +62,8 @@ public class MidpointStepInterpolatorTest { oos.writeObject(handler); } - Assert.assertTrue(bos.size () > 120000); - Assert.assertTrue(bos.size () < 125000); + Assert.assertTrue(bos.size () > 135000); + Assert.assertTrue(bos.size () < 145000); ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray()); ObjectInputStream ois = new ObjectInputStream(bis); diff --git a/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java b/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java index 6ffad4fdd..e60e0656e 100644 --- a/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java +++ b/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java @@ -61,8 +61,8 @@ public class ThreeEighthesStepInterpolatorTest { oos.writeObject(handler); } - Assert.assertTrue(bos.size () > 750000); - Assert.assertTrue(bos.size () < 800000); + Assert.assertTrue(bos.size () > 880000); + Assert.assertTrue(bos.size () < 900000); ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray()); ObjectInputStream ois = new ObjectInputStream(bis);