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

Improved Dormand-Prince 8(5,3) step interpolator accuracy at step start. 

The previous step is preserved and if interpolation time is in the first
half ot the step, then interpolation is based on step start, otherwise
is it based on step end. Previously, interpolation was always performed
with respect to step end.

If this trick proves useful, it will be extended to other Runge-Kutta
type step interpolators.

JIRA: MATH-705

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1204270 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2011-11-20 21:41:46 +00:00
parent 6d877485ed
commit 645d642b8f
19 changed files with 166 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java
View File

@ -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

2
src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
View File

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

42
src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
View File

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

3
src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
View File

@ -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) {

2
src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java
View File

@ -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

2
src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java
View File

@ -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

2
src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java
View File

@ -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

2
src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java
View File

@ -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

2
src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
View File

@ -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

38
src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java
View File

@ -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) {

2
src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java
View File

@ -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

82
src/test/java/org/apache/commons/math/ode/events/ReappearingEventTest.java Normal file
View File

@ -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.
}
}
}

4
src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java
View File

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

4
src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolatorTest.java
View File

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

4
src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolatorTest.java
View File

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

4
src/test/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolatorTest.java
View File

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

4
src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolatorTest.java
View File

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

4
src/test/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolatorTest.java
View File

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

4
src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java
View File

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