You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2015/12/09 17:16:51 UTC
[11/50] [abbrv] [math] One step towards immutable step interpolators.
One step towards immutable step interpolators.
The interpolators do not expect anymore the y and yDot arrays to
be shared with integrator and be updated by it.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/07e2b944
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/07e2b944
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/07e2b944
Branch: refs/heads/MATH_3_X
Commit: 07e2b9447b5950721f39b7d027a0ca95d0cd950e
Parents: 8180c9e
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue Dec 1 12:22:33 2015 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 1 12:22:33 2015 +0100
----------------------------------------------------------------------
.../ClassicalRungeKuttaFieldIntegrator.java | 7 +-
...lassicalRungeKuttaFieldStepInterpolator.java | 60 +-
.../DormandPrince54FieldIntegrator.java | 11 +-
.../DormandPrince54FieldStepInterpolator.java | 274 ++++----
.../DormandPrince853FieldIntegrator.java | 74 ++-
.../DormandPrince853FieldStepInterpolator.java | 630 +++++--------------
.../EmbeddedRungeKuttaFieldIntegrator.java | 31 +-
.../ode/nonstiff/EulerFieldIntegrator.java | 7 +-
.../nonstiff/EulerFieldStepInterpolator.java | 28 +-
.../math3/ode/nonstiff/GillFieldIntegrator.java | 7 +-
.../ode/nonstiff/GillFieldStepInterpolator.java | 62 +-
.../nonstiff/HighamHall54FieldIntegrator.java | 11 +-
.../HighamHall54FieldStepInterpolator.java | 61 +-
.../ode/nonstiff/LutherFieldIntegrator.java | 7 +-
.../nonstiff/LutherFieldStepInterpolator.java | 108 +---
.../ode/nonstiff/MidpointFieldIntegrator.java | 7 +-
.../nonstiff/MidpointFieldStepInterpolator.java | 40 +-
.../ode/nonstiff/RungeKuttaFieldIntegrator.java | 13 +-
.../RungeKuttaFieldStepInterpolator.java | 89 ++-
.../nonstiff/ThreeEighthesFieldIntegrator.java | 7 +-
.../ThreeEighthesFieldStepInterpolator.java | 48 +-
.../sampling/AbstractFieldStepInterpolator.java | 21 +-
.../commons/math3/ode/TestFieldProblem1.java | 5 +-
.../commons/math3/ode/TestFieldProblem5.java | 2 +-
.../math3/ode/TestFieldProblemAbstract.java | 18 +-
.../math3/ode/TestFieldProblemHandler.java | 2 +-
.../math3/ode/nonstiff/EulerIntegratorTest.java | 5 +-
27 files changed, 593 insertions(+), 1042 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
index eab1bbf..064b56d 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.util.MathArrays;
@@ -101,10 +100,8 @@ public class ClassicalRungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected ClassicalRungeKuttaFieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new ClassicalRungeKuttaFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new ClassicalRungeKuttaFieldStepInterpolator<T>(this, forward, mapper);
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
index 5bea4b9..fc56870 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaFieldStepInterpolator.java
@@ -21,7 +21,6 @@ import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
-import org.apache.commons.math3.util.MathArrays;
/**
* This class implements a step interpolator for the classical fourth
@@ -62,17 +61,13 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param rkIntegrator integrator being used
- * @param y reference to the integrator array holding the state at
- * the end of the step
- * @param yDotArray reference to the integrator array holding all the
- * intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
ClassicalRungeKuttaFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
- final T[] y, final T[][] yDotArray, final boolean forward,
+ final boolean forward,
final FieldEquationsMapper<T> mapper) {
- super(rkIntegrator, y, yDotArray, forward, mapper);
+ super(rkIntegrator, forward, mapper);
}
/** Copy constructor.
@@ -91,6 +86,7 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
}
/** {@inheritDoc} */
+ @SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
@@ -102,42 +98,28 @@ class ClassicalRungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
final T coeffDot1 = oneMinusTheta.multiply(oneMinus2Theta);
final T coeffDot23 = theta.multiply(oneMinusTheta).multiply(2);
final T coeffDot4 = theta.multiply(oneMinus2Theta).negate();
- final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
- final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
+ final T[] interpolatedState;
+ final T[] interpolatedDerivatives;
- if ((previousState != null) && (theta.getReal() <= 0.5)) {
- final T fourTheta2 = theta.multiply(theta).multiply(4);
- final T s = theta.multiply(h).divide(6.0);
- final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6));
- final T coeff23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
- final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));
- for (int i = 0; i < interpolatedState.length; ++i) {
- final T yDot1 = yDotK[0][i];
- final T yDot23 = yDotK[1][i].add(yDotK[2][i]);
- final T yDot4 = yDotK[3][i];
- interpolatedState[i] =
- previousState[i].add(coeff1.multiply(yDot1)).add(coeff23.multiply(yDot23)).add(coeff4.multiply(yDot4));
- interpolatedDerivatives[i] =
- coeffDot1.multiply(yDot1).add(coeffDot23.multiply(yDot23)).add(coeffDot4.multiply(yDot4));
- }
+ if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
+ final T fourTheta2 = theta.multiply(theta).multiply(4);
+ final T s = theta.multiply(h).divide(6.0);
+ final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(9)).add(6));
+ final T coeff23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
+ final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));
+ interpolatedState = previousStateLinearCombination(coeff1, coeff23, coeff23, coeff4);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot23, coeffDot23, coeffDot4);
} else {
- final T fourTheta = theta.multiply(4);
- final T s = oneMinusThetaH.divide(6);
- final T coeff1 = s.multiply(theta.multiply(fourTheta.negate().add(5)).subtract(1));
- final T coeff23 = s.multiply(theta.multiply(fourTheta.subtract(2)).subtract(2));
- final T coeff4 = s.multiply(theta.multiply(fourTheta.negate().subtract(1)).subtract(1));
- for (int i = 0; i < interpolatedState.length; ++i) {
- final T yDot1 = yDotK[0][i];
- final T yDot23 = yDotK[1][i].add(yDotK[2][i]);
- final T yDot4 = yDotK[3][i];
- interpolatedState[i] =
- currentState[i].add(coeff1.multiply(yDot1)).add(coeff23.multiply(yDot23)).add(coeff4.multiply(yDot4));
- interpolatedDerivatives[i] =
- coeffDot1.multiply(yDot1).add(coeffDot23.multiply(yDot23)).add(coeffDot4.multiply(yDot4));
- }
+ final T fourTheta = theta.multiply(4);
+ final T s = oneMinusThetaH.divide(6);
+ final T coeff1 = s.multiply(theta.multiply(fourTheta.negate().add(5)).subtract(1));
+ final T coeff23 = s.multiply(theta.multiply(fourTheta.subtract(2)).subtract(2));
+ final T coeff4 = s.multiply(theta.multiply(fourTheta.negate().subtract(1)).subtract(1));
+ interpolatedState = currentStateLinearCombination(coeff1, coeff23, coeff23, coeff4);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot23, coeffDot23, coeffDot4);
}
- return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
+ return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
index 7475c66..09bbc4a 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.util.MathArrays;
import org.apache.commons.math3.util.MathUtils;
@@ -93,7 +92,7 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double scalAbsoluteTolerance,
final double scalRelativeTolerance) {
- super(field, METHOD_NAME, true,
+ super(field, METHOD_NAME, 6,
minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
e1 = fraction( 71, 57600);
e3 = fraction( -71, 16695);
@@ -119,7 +118,7 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double[] vecAbsoluteTolerance,
final double[] vecRelativeTolerance) {
- super(field, METHOD_NAME, true,
+ super(field, METHOD_NAME, 6,
minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
e1 = fraction( 71, 57600);
e3 = fraction( -71, 16695);
@@ -190,10 +189,8 @@ public class DormandPrince54FieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected DormandPrince54FieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new DormandPrince54FieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new DormandPrince54FieldStepInterpolator<T>(this, forward, mapper);
}
/** {@inheritDoc} */
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
index d3fb208..23ba0e9 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54FieldStepInterpolator.java
@@ -21,7 +21,6 @@ import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
-import org.apache.commons.math3.util.MathArrays;
/**
* This class represents an interpolator over the last step during an
@@ -37,75 +36,63 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
extends RungeKuttaFieldStepInterpolator<T> {
/** Last row of the Butcher-array internal weights, element 0. */
- private static final double A70 = 35.0 / 384.0;
+ private final T a70;
// element 1 is zero, so it is neither stored nor used
/** Last row of the Butcher-array internal weights, element 2. */
- private static final double A72 = 500.0 / 1113.0;
+ private final T a72;
/** Last row of the Butcher-array internal weights, element 3. */
- private static final double A73 = 125.0 / 192.0;
+ private final T a73;
/** Last row of the Butcher-array internal weights, element 4. */
- private static final double A74 = -2187.0 / 6784.0;
+ private final T a74;
/** Last row of the Butcher-array internal weights, element 5. */
- private static final double A75 = 11.0 / 84.0;
+ private final T a75;
/** Shampine (1986) Dense output, element 0. */
- private static final double D0 = -12715105075.0 / 11282082432.0;
+ private final T d0;
// element 1 is zero, so it is neither stored nor used
/** Shampine (1986) Dense output, element 2. */
- private static final double D2 = 87487479700.0 / 32700410799.0;
+ private final T d2;
/** Shampine (1986) Dense output, element 3. */
- private static final double D3 = -10690763975.0 / 1880347072.0;
+ private final T d3;
/** Shampine (1986) Dense output, element 4. */
- private static final double D4 = 701980252875.0 / 199316789632.0;
+ private final T d4;
/** Shampine (1986) Dense output, element 5. */
- private static final double D5 = -1453857185.0 / 822651844.0;
+ private final T d5;
/** Shampine (1986) Dense output, element 6. */
- private static final double D6 = 69997945.0 / 29380423.0;
-
- /** First vector for interpolation. */
- private T[] v1;
-
- /** Second vector for interpolation. */
- private T[] v2;
-
- /** Third vector for interpolation. */
- private T[] v3;
-
- /** Fourth vector for interpolation. */
- private T[] v4;
-
- /** Initialization indicator for the interpolation vectors. */
- private boolean vectorsInitialized;
+ private final T d6;
/** Simple constructor.
* @param rkIntegrator integrator being used
- * @param y reference to the integrator array holding the state at
- * the end of the step
- * @param yDotArray reference to the integrator array holding all the
- * intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
DormandPrince54FieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
- final T[] y, final T[][] yDotArray, final boolean forward,
+ final boolean forward,
final FieldEquationsMapper<T> mapper) {
- super(rkIntegrator, y, yDotArray, forward, mapper);
- v1 = null;
- v2 = null;
- v3 = null;
- v4 = null;
- vectorsInitialized = false;
+ super(rkIntegrator, forward, mapper);
+ final T one = rkIntegrator.getField().getOne();
+ a70 = one.multiply( 35.0).divide( 384.0);
+ a72 = one.multiply( 500.0).divide(1113.0);
+ a73 = one.multiply( 125.0).divide( 192.0);
+ a74 = one.multiply(-2187.0).divide(6784.0);
+ a75 = one.multiply( 11.0).divide( 84.0);
+ d0 = one.multiply(-12715105075.0).divide( 11282082432.0);
+ d2 = one.multiply( 87487479700.0).divide( 32700410799.0);
+ d3 = one.multiply(-10690763975.0).divide( 1880347072.0);
+ d4 = one.multiply(701980252875.0).divide(199316789632.0);
+ d5 = one.multiply( -1453857185.0).divide( 822651844.0);
+ d6 = one.multiply( 69997945.0).divide( 29380423.0);
}
/** Copy constructor.
@@ -116,24 +103,17 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
DormandPrince54FieldStepInterpolator(final DormandPrince54FieldStepInterpolator<T> interpolator) {
super(interpolator);
-
- if (interpolator.v1 == null) {
-
- v1 = null;
- v2 = null;
- v3 = null;
- v4 = null;
- vectorsInitialized = false;
-
- } else {
-
- v1 = interpolator.v1.clone();
- v2 = interpolator.v2.clone();
- v3 = interpolator.v3.clone();
- v4 = interpolator.v4.clone();
- vectorsInitialized = interpolator.vectorsInitialized;
-
- }
+ a70 = interpolator.a70;
+ a72 = interpolator.a72;
+ a73 = interpolator.a73;
+ a74 = interpolator.a74;
+ a75 = interpolator.a75;
+ d0 = interpolator.d0;
+ d2 = interpolator.d2;
+ d3 = interpolator.d3;
+ d4 = interpolator.d4;
+ d5 = interpolator.d5;
+ d6 = interpolator.d6;
}
@@ -143,58 +123,13 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
return new DormandPrince54FieldStepInterpolator<T>(this);
}
-
- /** {@inheritDoc} */
- @Override
- public void storeState(final FieldODEStateAndDerivative<T> state) {
- super.storeState(state);
- vectorsInitialized = false;
- }
-
/** {@inheritDoc} */
+ @SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
- if (! vectorsInitialized) {
-
- if (v1 == null) {
- v1 = MathArrays.buildArray(time.getField(), previousState.length);
- v2 = MathArrays.buildArray(time.getField(), previousState.length);
- v3 = MathArrays.buildArray(time.getField(), previousState.length);
- v4 = MathArrays.buildArray(time.getField(), previousState.length);
- }
-
- // no step finalization is needed for this interpolator
-
- // we need to compute the interpolation vectors for this time step
- for (int i = 0; i < previousState.length; ++i) {
- final T yDot0 = yDotK[0][i];
- final T yDot2 = yDotK[2][i];
- final T yDot3 = yDotK[3][i];
- final T yDot4 = yDotK[4][i];
- final T yDot5 = yDotK[5][i];
- final T yDot6 = yDotK[6][i];
- v1[i] = yDot0.multiply(A70).
- add(yDot2.multiply(A72)).
- add(yDot3.multiply(A73)).
- add(yDot4.multiply(A74)).
- add(yDot5.multiply(A75));
- v2[i] = yDot0.subtract(v1[i]);
- v3[i] = v1[i].subtract(v2[i]).subtract(yDot6);
- v4[i] = yDot0.multiply(D0).
- add(yDot2.multiply(D2)).
- add(yDot3.multiply(D3)).
- add(yDot4.multiply(D4)).
- add(yDot5.multiply(D5)).
- add(yDot6.multiply(D6));
- }
-
- vectorsInitialized = true;
-
- }
-
// interpolate
final T one = theta.getField().getOne();
final T eta = one.subtract(theta);
@@ -202,35 +137,116 @@ class DormandPrince54FieldStepInterpolator<T extends RealFieldElement<T>>
final T dot2 = one.subtract(twoTheta);
final T dot3 = theta.multiply(theta.multiply(-3).add(2));
final T dot4 = twoTheta.multiply(theta.multiply(twoTheta.subtract(3)).add(1));
- final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
- final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
- if ((previousState != null) && (theta.getReal() <= 0.5)) {
- for (int i = 0; i < interpolatedState.length; ++i) {
- interpolatedState[i] = previousState[i].
- add(theta.multiply(h.multiply(v1[i].
- add(eta.multiply(v2[i].
- add(theta.multiply(v3[i].
- add(eta.multiply(v4[i])))))))));
- interpolatedDerivatives[i] = v1[i].
- add(dot2.multiply(v2[i])).
- add(dot3.multiply(v3[i])).
- add(dot4.multiply(v4[i]));
- }
+ final T[] interpolatedState;
+ final T[] interpolatedDerivatives;
+ if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
+ final T f1 = h.multiply(theta);
+ final T f2 = f1.multiply(eta);
+ final T f3 = f2.multiply(theta);
+ final T f4 = f3.multiply(eta);
+ final T coeff0 = f1.multiply(a70).
+ subtract(f2.multiply(a70.subtract(1))).
+ add(f3.multiply(a70.multiply(2).subtract(1))).
+ add(f4.multiply(d0));
+ final T coeff1 = theta.getField().getZero();
+ final T coeff2 = f1.multiply(a72).
+ subtract(f2.multiply(a72)).
+ add(f3.multiply(a72.multiply(2))).
+ add(f4.multiply(d2));
+ final T coeff3 = f1.multiply(a73).
+ subtract(f2.multiply(a73)).
+ add(f3.multiply(a73.multiply(2))).
+ add(f4.multiply(d3));
+ final T coeff4 = f1.multiply(a74).
+ subtract(f2.multiply(a74)).
+ add(f3.multiply(a74.multiply(2))).
+ add(f4.multiply(d4));
+ final T coeff5 = f1.multiply(a75).
+ subtract(f2.multiply(a75)).
+ add(f3.multiply(a75.multiply(2))).
+ add(f4.multiply(d5));
+ final T coeff6 = f4.multiply(d6).subtract(f3);
+ final T coeffDot0 = a70.
+ subtract(dot2.multiply(a70.subtract(1))).
+ add(dot3.multiply(a70.multiply(2).subtract(1))).
+ add(dot4.multiply(d0));
+ final T coeffDot1 = theta.getField().getZero();
+ final T coeffDot2 = a72.
+ subtract(dot2.multiply(a72)).
+ add(dot3.multiply(a72.multiply(2))).
+ add(dot4.multiply(d2));
+ final T coeffDot3 = a73.
+ subtract(dot2.multiply(a73)).
+ add(dot3.multiply(a73.multiply(2))).
+ add(dot4.multiply(d3));
+ final T coeffDot4 = a74.
+ subtract(dot2.multiply(a74)).
+ add(dot3.multiply(a74.multiply(2))).
+ add(dot4.multiply(d4));
+ final T coeffDot5 = a75.
+ subtract(dot2.multiply(a75)).
+ add(dot3.multiply(a75.multiply(2))).
+ add(dot4.multiply(d5));
+ final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
+ interpolatedState = previousStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
+ coeff4, coeff5, coeff6);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
+ coeffDot4, coeffDot5, coeffDot6);
} else {
- for (int i = 0; i < interpolatedState.length; ++i) {
- interpolatedState[i] = currentState[i].
- subtract(oneMinusThetaH.multiply(v1[i].
- subtract(theta.multiply(v2[i].
- add(theta.multiply(v3[i].
- add(eta.multiply(v4[i]))))))));
- interpolatedDerivatives[i] = v1[i].
- add(dot2.multiply(v2[i])).
- add(dot3.multiply(v3[i])).
- add(dot4.multiply(v4[i]));
- }
+ final T f1 = oneMinusThetaH.negate();
+ final T f2 = oneMinusThetaH.multiply(theta);
+ final T f3 = f2.multiply(theta);
+ final T f4 = f3.multiply(eta);
+ final T coeff0 = f1.multiply(a70).
+ subtract(f2.multiply(a70.subtract(1))).
+ add(f3.multiply(a70.multiply(2).subtract(1))).
+ add(f4.multiply(d0));
+ final T coeff1 = theta.getField().getZero();
+ final T coeff2 = f1.multiply(a72).
+ subtract(f2.multiply(a72)).
+ add(f3.multiply(a72.multiply(2))).
+ add(f4.multiply(d2));
+ final T coeff3 = f1.multiply(a73).
+ subtract(f2.multiply(a73)).
+ add(f3.multiply(a73.multiply(2))).
+ add(f4.multiply(d3));
+ final T coeff4 = f1.multiply(a74).
+ subtract(f2.multiply(a74)).
+ add(f3.multiply(a74.multiply(2))).
+ add(f4.multiply(d4));
+ final T coeff5 = f1.multiply(a75).
+ subtract(f2.multiply(a75)).
+ add(f3.multiply(a75.multiply(2))).
+ add(f4.multiply(d5));
+ final T coeff6 = f4.multiply(d6).subtract(f3);
+ final T coeffDot0 = a70.
+ subtract(dot2.multiply(a70.subtract(1))).
+ add(dot3.multiply(a70.multiply(2).subtract(1))).
+ add(dot4.multiply(d0));
+ final T coeffDot1 = theta.getField().getZero();
+ final T coeffDot2 = a72.
+ subtract(dot2.multiply(a72)).
+ add(dot3.multiply(a72.multiply(2))).
+ add(dot4.multiply(d2));
+ final T coeffDot3 = a73.
+ subtract(dot2.multiply(a73)).
+ add(dot3.multiply(a73.multiply(2))).
+ add(dot4.multiply(d3));
+ final T coeffDot4 = a74.
+ subtract(dot2.multiply(a74)).
+ add(dot3.multiply(a74.multiply(2))).
+ add(dot4.multiply(d4));
+ final T coeffDot5 = a75.
+ subtract(dot2.multiply(a75)).
+ add(dot3.multiply(a75.multiply(2))).
+ add(dot4.multiply(d5));
+ final T coeffDot6 = dot4.multiply(d6).subtract(dot3);
+ interpolatedState = currentStateLinearCombination(coeff0, coeff1, coeff2, coeff3,
+ coeff4, coeff5, coeff6);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot0, coeffDot1, coeffDot2, coeffDot3,
+ coeffDot4, coeffDot5, coeffDot6);
}
-
- return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
+ return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
index 6f6b436..569cba4 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.util.MathArrays;
import org.apache.commons.math3.util.MathUtils;
@@ -134,7 +133,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double scalAbsoluteTolerance,
final double scalRelativeTolerance) {
- super(field, METHOD_NAME, true,
+ super(field, METHOD_NAME, 12,
minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
e1_01 = fraction( 116092271.0, 8848465920.0);
e1_06 = fraction( -1871647.0, 1527680.0);
@@ -170,7 +169,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double[] vecAbsoluteTolerance,
final double[] vecRelativeTolerance) {
- super(field, METHOD_NAME, true,
+ super(field, METHOD_NAME, 12,
minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
e1_01 = fraction( 116092271.0, 8848465920.0);
e1_06 = fraction( -1871647.0, 1527680.0);
@@ -196,7 +195,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
final T sqrt6 = getField().getOne().multiply(6).sqrt();
- final T[] c = MathArrays.buildArray(getField(), 12);
+ final T[] c = MathArrays.buildArray(getField(), 15);
c[ 0] = sqrt6.add(-6).divide(-67.5);
c[ 1] = sqrt6.add(-6).divide(-45.0);
c[ 2] = sqrt6.add(-6).divide(-30.0);
@@ -209,6 +208,9 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
c[ 9] = fraction(6, 7);
c[10] = getField().getOne();
c[11] = getField().getOne();
+ c[12] = fraction(1.0, 10.0);
+ c[13] = fraction(1.0, 5.0);
+ c[14] = fraction(7.0, 9.0);
return c;
@@ -220,7 +222,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
final T sqrt6 = getField().getOne().multiply(6).sqrt();
- final T[][] a = MathArrays.buildArray(getField(), 12, -1);
+ final T[][] a = MathArrays.buildArray(getField(), 15, -1);
for (int i = 0; i < a.length; ++i) {
a[i] = MathArrays.buildArray(getField(), i + 1);
}
@@ -302,9 +304,8 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
a[10][ 9] = fraction(226716250.0, 18341897.0);
a[10][10] = fraction(1371316744.0, 2131383595.0);
- // this stage should be for interpolation only, but since it is the same
- // stage as the first evaluation of the next step, we perform it
- // here at no cost by specifying this is an fsal method
+ // the following stage is both for interpolation and the first stage in next step
+ // (the coefficients are identical to the B array)
a[11][ 0] = fraction(104257.0, 1920240.0);
a[11][ 1] = getField().getZero();
a[11][ 2] = getField().getZero();
@@ -318,6 +319,52 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
a[11][10] = fraction(171414593.0, 851261400.0);
a[11][11] = fraction(137909.0, 3084480.0);
+ // the following stages are for interpolation only
+ a[12][ 0] = fraction( 13481885573.0, 240030000000.0) .subtract(a[11][0]);
+ a[12][ 1] = getField().getZero();
+ a[12][ 2] = getField().getZero();
+ a[12][ 3] = getField().getZero();
+ a[12][ 4] = getField().getZero();
+ a[12][ 5] = getField().getZero() .subtract(a[11][5]);
+ a[12][ 6] = fraction( 139418837528.0, 549975234375.0) .subtract(a[11][6]);
+ a[12][ 7] = fraction( -11108320068443.0, 45111937500000.0) .subtract(a[11][7]);
+ a[12][ 8] = fraction(-1769651421925959.0, 14249385146080000.0).subtract(a[11][8]);
+ a[12][ 9] = fraction( 57799439.0, 377055000.0) .subtract(a[11][9]);
+ a[12][10] = fraction( 793322643029.0, 96734250000000.0) .subtract(a[11][10]);
+ a[12][11] = fraction( 1458939311.0, 192780000000.0) .subtract(a[11][11]);
+ a[12][12] = fraction( -4149.0, 500000.0);
+
+ a[13][ 0] = fraction( 1595561272731.0, 50120273500000.0) .subtract(a[11][0]);
+ a[13][ 1] = getField().getZero();
+ a[13][ 2] = getField().getZero();
+ a[13][ 3] = getField().getZero();
+ a[13][ 4] = getField().getZero();
+ a[13][ 5] = fraction( 975183916491.0, 34457688031250.0) .subtract(a[11][5]);
+ a[13][ 6] = fraction( 38492013932672.0, 718912673015625.0) .subtract(a[11][6]);
+ a[13][ 7] = fraction(-1114881286517557.0, 20298710767500000.0).subtract(a[11][7]);
+ a[13][ 8] = getField().getZero() .subtract(a[11][8]);
+ a[13][ 9] = getField().getZero() .subtract(a[11][9]);
+ a[13][10] = fraction( -2538710946863.0, 23431227861250000.0).subtract(a[11][10]);
+ a[13][11] = fraction( 8824659001.0, 23066716781250.0) .subtract(a[11][11]);
+ a[13][12] = fraction( -11518334563.0, 33831184612500.0);
+ a[13][13] = fraction( 1912306948.0, 13532473845.0);
+
+ a[14][ 0] = fraction( -13613986967.0, 31741908048.0) .subtract(a[11][0]);
+ a[14][ 1] = getField().getZero();
+ a[14][ 2] = getField().getZero();
+ a[14][ 3] = getField().getZero();
+ a[14][ 4] = getField().getZero();
+ a[14][ 5] = fraction( -4755612631.0, 1012344804.0) .subtract(a[11][5]);
+ a[14][ 6] = fraction( 42939257944576.0, 5588559685701.0) .subtract(a[11][6]);
+ a[14][ 7] = fraction( 77881972900277.0, 19140370552944.0) .subtract(a[11][7]);
+ a[14][ 8] = fraction( 22719829234375.0, 63689648654052.0) .subtract(a[11][8]);
+ a[14][ 9] = getField().getZero() .subtract(a[11][9]);
+ a[14][10] = getField().getZero() .subtract(a[11][10]);
+ a[14][11] = getField().getZero() .subtract(a[11][11]);
+ a[14][12] = fraction( -1199007803.0, 857031517296.0);
+ a[14][13] = fraction( 157882067000.0, 53564469831.0);
+ a[14][14] = fraction( -290468882375.0, 31741908048.0);
+
return a;
}
@@ -325,7 +372,7 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected T[] getB() {
- final T[] b = MathArrays.buildArray(getField(), 13);
+ final T[] b = MathArrays.buildArray(getField(), 16);
b[ 0] = fraction(104257, 1929240);
b[ 1] = getField().getZero();
b[ 2] = getField().getZero();
@@ -339,16 +386,17 @@ public class DormandPrince853FieldIntegrator<T extends RealFieldElement<T>>
b[10] = fraction( 171414593.0, 851261400.0);
b[11] = fraction( 137909.0, 3084480.0);
b[12] = getField().getZero();
+ b[13] = getField().getZero();
+ b[14] = getField().getZero();
+ b[15] = getField().getZero();
return b;
}
/** {@inheritDoc} */
@Override
protected DormandPrince853FieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new DormandPrince853FieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new DormandPrince853FieldStepInterpolator<T>(this, forward, mapper);
}
/** {@inheritDoc} */
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
index 3966121..7c1cdf5 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853FieldStepInterpolator.java
@@ -17,7 +17,6 @@
package org.apache.commons.math3.ode.nonstiff;
-import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.ode.AbstractFieldIntegrator;
@@ -38,268 +37,143 @@ import org.apache.commons.math3.util.MathArrays;
class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
extends RungeKuttaFieldStepInterpolator<T> {
- /** Propagation weights, element 1. */
- private final T b_01;
-
- // elements 2 to 5 are zero, so they are neither stored nor used
-
- /** Propagation weights, element 6. */
- private final T b_06;
-
- /** Propagation weights, element 7. */
- private final T b_07;
-
- /** Propagation weights, element 8. */
- private final T b_08;
-
- /** Propagation weights, element 9. */
- private final T b_09;
-
- /** Propagation weights, element 10. */
- private final T b_10;
-
- /** Propagation weights, element 11. */
- private final T b_11;
-
- /** Propagation weights, element 12. */
- private final T b_12;
-
- /** Time step for stage 14 (interpolation only). */
- private final T c14;
-
- /** Internal weights for stage 14, element 1. */
- private final T k14_01;
-
- // elements 2 to 5 are zero, so they are neither stored nor used
-
- /** Internal weights for stage 14, element 6. */
- private final T k14_06;
-
- /** Internal weights for stage 14, element 7. */
- private final T k14_07;
-
- /** Internal weights for stage 14, element 8. */
- private final T k14_08;
-
- /** Internal weights for stage 14, element 9. */
- private final T k14_09;
-
- /** Internal weights for stage 14, element 10. */
- private final T k14_10;
-
- /** Internal weights for stage 14, element 11. */
- private final T k14_11;
-
- /** Internal weights for stage 14, element 12. */
- private final T k14_12;
-
- /** Internal weights for stage 14, element 13. */
- private final T k14_13;
-
- /** Time step for stage 15 (interpolation only). */
- private final T c15;
-
-
- /** Internal weights for stage 15, element 1. */
- private final T k15_01;
-
- // elements 2 to 5 are zero, so they are neither stored nor used
-
- /** Internal weights for stage 15, element 6. */
- private final T k15_06;
-
- /** Internal weights for stage 15, element 7. */
- private final T k15_07;
-
- /** Internal weights for stage 15, element 8. */
- private final T k15_08;
-
- /** Internal weights for stage 15, element 9. */
- private final T k15_09;
-
- /** Internal weights for stage 15, element 10. */
- private final T k15_10;
-
- /** Internal weights for stage 15, element 11. */
- private final T k15_11;
-
- /** Internal weights for stage 15, element 12. */
- private final T k15_12;
-
- /** Internal weights for stage 15, element 13. */
- private final T k15_13;
-
- /** Internal weights for stage 15, element 14. */
- private final T k15_14;
-
- /** Time step for stage 16 (interpolation only). */
- private final T c16;
-
-
- /** Internal weights for stage 16, element 1. */
- private final T k16_01;
-
- // elements 2 to 5 are zero, so they are neither stored nor used
-
- /** Internal weights for stage 16, element 6. */
- private final T k16_06;
-
- /** Internal weights for stage 16, element 7. */
- private final T k16_07;
-
- /** Internal weights for stage 16, element 8. */
- private final T k16_08;
-
- /** Internal weights for stage 16, element 9. */
- private final T k16_09;
-
- /** Internal weights for stage 16, element 10. */
- private final T k16_10;
-
- /** Internal weights for stage 16, element 11. */
- private final T k16_11;
-
- /** Internal weights for stage 16, element 12. */
- private final T k16_12;
-
- /** Internal weights for stage 16, element 13. */
- private final T k16_13;
-
- /** Internal weights for stage 16, element 14. */
- private final T k16_14;
-
- /** Internal weights for stage 16, element 15. */
- private final T k16_15;
-
/** Interpolation weights.
* (beware that only the non-null values are in the table)
*/
private final T[][] d;
- /** Last evaluations. */
- private T[][] yDotKLast;
-
- /** Vectors for interpolation. */
- private T[][] v;
-
- /** Initialization indicator for the interpolation vectors. */
- private boolean vectorsInitialized;
-
/** Simple constructor.
* @param rkIntegrator integrator being used
- * @param y reference to the integrator array holding the state at
- * the end of the step
- * @param yDotArray reference to the integrator array holding all the
- * intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
DormandPrince853FieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
- final T[] y, final T[][] yDotArray, final boolean forward,
+ final boolean forward,
final FieldEquationsMapper<T> mapper) {
- super(rkIntegrator, y, yDotArray, forward, mapper);
- yDotKLast = null;
- v = null;
- vectorsInitialized = false;
-
- b_01 = fraction( 104257.0, 1920240.0);
- b_06 = fraction( 3399327.0, 763840.0);
- b_07 = fraction( 66578432.0, 35198415.0);
- b_08 = fraction( -1674902723.0, 288716400.0);
- b_09 = fraction(54980371265625.0, 176692375811392.0);
- b_10 = fraction( -734375.0, 4826304.0);
- b_11 = fraction( 171414593.0, 851261400.0);
- b_12 = fraction( 137909.0, 3084480.0);
- c14 = fraction(1.0, 10.0);
- k14_01 = fraction( 13481885573.0, 240030000000.0) .subtract(b_01);
- k14_06 = integrator.getField().getZero() .subtract(b_06);
- k14_07 = fraction( 139418837528.0, 549975234375.0) .subtract(b_07);
- k14_08 = fraction( -11108320068443.0, 45111937500000.0) .subtract(b_08);
- k14_09 = fraction(-1769651421925959.0, 14249385146080000.0).subtract(b_09);
- k14_10 = fraction( 57799439.0, 377055000.0) .subtract(b_10);
- k14_11 = fraction( 793322643029.0, 96734250000000.0) .subtract(b_11);
- k14_12 = fraction( 1458939311.0, 192780000000.0) .subtract(b_12);
- k14_13 = fraction( -4149.0, 500000.0);
- c15 = fraction(1.0, 5.0);
- k15_01 = fraction( 1595561272731.0, 50120273500000.0) .subtract(b_01);
- k15_06 = fraction( 975183916491.0, 34457688031250.0) .subtract(b_06);
- k15_07 = fraction( 38492013932672.0, 718912673015625.0) .subtract(b_07);
- k15_08 = fraction(-1114881286517557.0, 20298710767500000.0).subtract(b_08);
- k15_09 = integrator.getField().getZero() .subtract(b_09);
- k15_10 = integrator.getField().getZero() .subtract(b_10);
- k15_11 = fraction( -2538710946863.0, 23431227861250000.0).subtract(b_11);
- k15_12 = fraction( 8824659001.0, 23066716781250.0) .subtract(b_12);
- k15_13 = fraction( -11518334563.0, 33831184612500.0);
- k15_14 = fraction( 1912306948.0, 13532473845.0);
- c16 = fraction(7.0, 9.0);
- k16_01 = fraction( -13613986967.0, 31741908048.0) .subtract(b_01);
- k16_06 = fraction( -4755612631.0, 1012344804.0) .subtract(b_06);
- k16_07 = fraction( 42939257944576.0, 5588559685701.0) .subtract(b_07);
- k16_08 = fraction( 77881972900277.0, 19140370552944.0) .subtract(b_08);
- k16_09 = fraction( 22719829234375.0, 63689648654052.0) .subtract(b_09);
- k16_10 = integrator.getField().getZero() .subtract(b_10);
- k16_11 = integrator.getField().getZero() .subtract(b_11);
- k16_12 = integrator.getField().getZero() .subtract(b_12);
- k16_13 = fraction( -1199007803.0, 857031517296.0);
- k16_14 = fraction( 157882067000.0, 53564469831.0);
- k16_15 = fraction( -290468882375.0, 31741908048.0);
-
- /** Interpolation weights.
- * (beware that only the non-null values are in the table)
- */
- d = MathArrays.buildArray(integrator.getField(), 4, 12);
-
- d[0][ 0] = fraction( -17751989329.0, 2106076560.0);
- d[0][ 1] = fraction( 4272954039.0, 7539864640.0);
- d[0][ 2] = fraction( -118476319744.0, 38604839385.0);
- d[0][ 3] = fraction( 755123450731.0, 316657731600.0);
- d[0][ 4] = fraction( 3692384461234828125.0, 1744130441634250432.0);
- d[0][ 5] = fraction( -4612609375.0, 5293382976.0);
- d[0][ 6] = fraction( 2091772278379.0, 933644586600.0);
- d[0][ 7] = fraction( 2136624137.0, 3382989120.0);
- d[0][ 8] = fraction( -126493.0, 1421424.0);
- d[0][ 9] = fraction( 98350000.0, 5419179.0);
- d[0][10] = fraction( -18878125.0, 2053168.0);
- d[0][11] = fraction( -1944542619.0, 438351368.0);
-
- d[1][ 0] = fraction( 32941697297.0, 3159114840.0);
- d[1][ 1] = fraction( 456696183123.0, 1884966160.0);
- d[1][ 2] = fraction( 19132610714624.0, 115814518155.0);
- d[1][ 3] = fraction( -177904688592943.0, 474986597400.0);
- d[1][ 4] = fraction(-4821139941836765625.0, 218016305204281304.0);
- d[1][ 5] = fraction( 30702015625.0, 3970037232.0);
- d[1][ 6] = fraction( -85916079474274.0, 2800933759800.0);
- d[1][ 7] = fraction( -5919468007.0, 634310460.0);
- d[1][ 8] = fraction( 2479159.0, 157936.0);
- d[1][ 9] = fraction( -18750000.0, 602131.0);
- d[1][10] = fraction( -19203125.0, 2053168.0);
- d[1][11] = fraction( 15700361463.0, 438351368.0);
-
- d[2][ 0] = fraction( 12627015655.0, 631822968.0);
- d[2][ 1] = fraction( -72955222965.0, 188496616.0);
- d[2][ 2] = fraction( -13145744952320.0, 69488710893.0);
- d[2][ 3] = fraction( 30084216194513.0, 56998391688.0);
- d[2][ 4] = fraction( -296858761006640625.0, 25648977082856624.0);
- d[2][ 5] = fraction( 569140625.0, 82709109.0);
- d[2][ 6] = fraction( -18684190637.0, 18672891732.0);
- d[2][ 7] = fraction( 69644045.0, 89549712.0);
- d[2][ 8] = fraction( -11847025.0, 4264272.0);
- d[2][ 9] = fraction( -978650000.0, 16257537.0);
- d[2][10] = fraction( 519371875.0, 6159504.0);
- d[2][11] = fraction( 5256837225.0, 438351368.0);
-
- d[3][ 0] = fraction( -450944925.0, 17550638.0);
- d[3][ 1] = fraction( -14532122925.0, 94248308.0);
- d[3][ 2] = fraction( -595876966400.0, 2573655959.0);
- d[3][ 3] = fraction( 188748653015.0, 527762886.0);
- d[3][ 4] = fraction( 2545485458115234375.0, 27252038150535163.0);
- d[3][ 5] = fraction( -1376953125.0, 36759604.0);
- d[3][ 6] = fraction( 53995596795.0, 518691437.0);
- d[3][ 7] = fraction( 210311225.0, 7047894.0);
- d[3][ 8] = fraction( -1718875.0, 39484.0);
- d[3][ 9] = fraction( 58000000.0, 602131.0);
- d[3][10] = fraction( -1546875.0, 39484.0);
- d[3][11] = fraction( -1262172375.0, 8429834.0);
+ super(rkIntegrator, forward, mapper);
+
+ // interpolation weights
+ d = MathArrays.buildArray(integrator.getField(), 4, 16);
+
+ // this row is the same as the b array
+ d[0][ 0] = fraction(104257, 1929240);
+ d[0][ 1] = integrator.getField().getZero();
+ d[0][ 2] = integrator.getField().getZero();
+ d[0][ 3] = integrator.getField().getZero();
+ d[0][ 4] = integrator.getField().getZero();
+ d[0][ 5] = fraction( 3399327.0, 763840.0);
+ d[0][ 6] = fraction( 66578432.0, 35198415.0);
+ d[0][ 7] = fraction( -1674902723.0, 288716400.0);
+ d[0][ 8] = fraction( 54980371265625.0, 176692375811392.0);
+ d[0][ 9] = fraction( -734375.0, 4826304.0);
+ d[0][10] = fraction( 171414593.0, 851261400.0);
+ d[0][11] = fraction( 137909.0, 3084480.0);
+ d[0][12] = integrator.getField().getZero();
+ d[0][13] = integrator.getField().getZero();
+ d[0][14] = integrator.getField().getZero();
+ d[0][15] = integrator.getField().getZero();
+
+ d[1][ 0] = d[0][ 0].negate().add(1);
+ d[1][ 1] = d[0][ 1].negate();
+ d[1][ 2] = d[0][ 2].negate();
+ d[1][ 3] = d[0][ 3].negate();
+ d[1][ 4] = d[0][ 4].negate();
+ d[1][ 5] = d[0][ 5].negate();
+ d[1][ 6] = d[0][ 6].negate();
+ d[1][ 7] = d[0][ 7].negate();
+ d[1][ 8] = d[0][ 8].negate();
+ d[1][ 9] = d[0][ 9].negate();
+ d[1][10] = d[0][10].negate();
+ d[1][11] = d[0][11].negate();
+ d[1][12] = d[0][12].negate(); // really 0
+ d[1][13] = d[0][13].negate(); // really 0
+ d[1][14] = d[0][14].negate(); // really 0
+ d[1][15] = d[0][15].negate(); // really 0
+
+ d[2][ 0] = d[0][ 0].multiply(2).subtract(1);
+ d[2][ 1] = d[0][ 1].multiply(2);
+ d[2][ 2] = d[0][ 2].multiply(2);
+ d[2][ 3] = d[0][ 3].multiply(2);
+ d[2][ 4] = d[0][ 4].multiply(2);
+ d[2][ 5] = d[0][ 5].multiply(2);
+ d[2][ 6] = d[0][ 6].multiply(2);
+ d[2][ 7] = d[0][ 7].multiply(2);
+ d[2][ 8] = d[0][ 8].multiply(2);
+ d[2][ 9] = d[0][ 9].multiply(2);
+ d[2][10] = d[0][10].multiply(2);
+ d[2][11] = d[0][11].multiply(2);
+ d[2][12] = d[0][12].multiply(2); // really 0
+ d[2][13] = d[0][13].multiply(2); // really 0
+ d[2][14] = d[0][14].multiply(2); // really 0
+ d[2][15] = d[0][15].multiply(2); // really 0
+
+ d[3][ 0] = fraction( -17751989329.0, 2106076560.0);
+ d[3][ 1] = integrator.getField().getZero();
+ d[3][ 2] = integrator.getField().getZero();
+ d[3][ 3] = integrator.getField().getZero();
+ d[3][ 4] = integrator.getField().getZero();
+ d[3][ 5] = fraction( 4272954039.0, 7539864640.0);
+ d[3][ 6] = fraction( -118476319744.0, 38604839385.0);
+ d[3][ 7] = fraction( 755123450731.0, 316657731600.0);
+ d[3][ 8] = fraction( 3692384461234828125.0, 1744130441634250432.0);
+ d[3][ 9] = fraction( -4612609375.0, 5293382976.0);
+ d[3][10] = fraction( 2091772278379.0, 933644586600.0);
+ d[3][11] = fraction( 2136624137.0, 3382989120.0);
+ d[3][12] = fraction( -126493.0, 1421424.0);
+ d[3][13] = fraction( 98350000.0, 5419179.0);
+ d[3][14] = fraction( -18878125.0, 2053168.0);
+ d[3][15] = fraction( -1944542619.0, 438351368.0);
+
+ d[4][ 0] = fraction( 32941697297.0, 3159114840.0);
+ d[4][ 1] = integrator.getField().getZero();
+ d[4][ 2] = integrator.getField().getZero();
+ d[4][ 3] = integrator.getField().getZero();
+ d[4][ 4] = integrator.getField().getZero();
+ d[4][ 5] = fraction( 456696183123.0, 1884966160.0);
+ d[4][ 6] = fraction( 19132610714624.0, 115814518155.0);
+ d[4][ 7] = fraction( -177904688592943.0, 474986597400.0);
+ d[4][ 8] = fraction(-4821139941836765625.0, 218016305204281304.0);
+ d[4][ 9] = fraction( 30702015625.0, 3970037232.0);
+ d[4][10] = fraction( -85916079474274.0, 2800933759800.0);
+ d[4][11] = fraction( -5919468007.0, 634310460.0);
+ d[4][12] = fraction( 2479159.0, 157936.0);
+ d[4][13] = fraction( -18750000.0, 602131.0);
+ d[4][14] = fraction( -19203125.0, 2053168.0);
+ d[4][15] = fraction( 15700361463.0, 438351368.0);
+
+ d[5][ 0] = fraction( 12627015655.0, 631822968.0);
+ d[5][ 1] = integrator.getField().getZero();
+ d[5][ 2] = integrator.getField().getZero();
+ d[5][ 3] = integrator.getField().getZero();
+ d[5][ 4] = integrator.getField().getZero();
+ d[5][ 5] = fraction( -72955222965.0, 188496616.0);
+ d[5][ 6] = fraction( -13145744952320.0, 69488710893.0);
+ d[5][ 7] = fraction( 30084216194513.0, 56998391688.0);
+ d[5][ 8] = fraction( -296858761006640625.0, 25648977082856624.0);
+ d[5][ 9] = fraction( 569140625.0, 82709109.0);
+ d[5][10] = fraction( -18684190637.0, 18672891732.0);
+ d[5][11] = fraction( 69644045.0, 89549712.0);
+ d[5][12] = fraction( -11847025.0, 4264272.0);
+ d[5][13] = fraction( -978650000.0, 16257537.0);
+ d[5][14] = fraction( 519371875.0, 6159504.0);
+ d[5][15] = fraction( 5256837225.0, 438351368.0);
+
+ d[6][ 0] = fraction( -450944925.0, 17550638.0);
+ d[6][ 1] = integrator.getField().getZero();
+ d[6][ 2] = integrator.getField().getZero();
+ d[6][ 3] = integrator.getField().getZero();
+ d[6][ 4] = integrator.getField().getZero();
+ d[6][ 5] = fraction( -14532122925.0, 94248308.0);
+ d[6][ 6] = fraction( -595876966400.0, 2573655959.0);
+ d[6][ 7] = fraction( 188748653015.0, 527762886.0);
+ d[6][ 8] = fraction( 2545485458115234375.0, 27252038150535163.0);
+ d[6][ 9] = fraction( -1376953125.0, 36759604.0);
+ d[6][10] = fraction( 53995596795.0, 518691437.0);
+ d[6][11] = fraction( 210311225.0, 7047894.0);
+ d[6][12] = fraction( -1718875.0, 39484.0);
+ d[6][13] = fraction( 58000000.0, 602131.0);
+ d[6][14] = fraction( -1546875.0, 39484.0);
+ d[6][15] = fraction( -1262172375.0, 8429834.0);
}
@@ -312,74 +186,6 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
super(interpolator);
- if (interpolator.currentState == null) {
-
- yDotKLast = null;
- v = null;
- vectorsInitialized = false;
-
- } else {
-
- final int dimension = interpolator.currentState.length;
- final Field<T> field = interpolator.getGlobalPreviousState().getTime().getField();
-
- yDotKLast = MathArrays.buildArray(field, 3, dimension);
- for (int k = 0; k < yDotKLast.length; ++k) {
- System.arraycopy(interpolator.yDotKLast[k], 0, yDotKLast[k], 0,
- dimension);
- }
-
- v = MathArrays.buildArray(field, 7, dimension);
- for (int k = 0; k < v.length; ++k) {
- System.arraycopy(interpolator.v[k], 0, v[k], 0, dimension);
- }
-
- vectorsInitialized = interpolator.vectorsInitialized;
-
- }
-
- b_01 = interpolator.b_01;
- b_06 = interpolator.b_06;
- b_07 = interpolator.b_07;
- b_08 = interpolator.b_08;
- b_09 = interpolator.b_09;
- b_10 = interpolator.b_10;
- b_11 = interpolator.b_11;
- b_12 = interpolator.b_12;
- c14 = interpolator.c14;
- k14_01 = interpolator.k14_01;
- k14_06 = interpolator.k14_06;
- k14_07 = interpolator.k14_07;
- k14_08 = interpolator.k14_08;
- k14_09 = interpolator.k14_09;
- k14_10 = interpolator.k14_10;
- k14_11 = interpolator.k14_11;
- k14_12 = interpolator.k14_12;
- k14_13 = interpolator.k14_13;
- c15 = interpolator.c15;
- k15_01 = interpolator.k15_01;
- k15_06 = interpolator.k15_06;
- k15_07 = interpolator.k15_07;
- k15_08 = interpolator.k15_08;
- k15_09 = interpolator.k15_09;
- k15_10 = interpolator.k15_10;
- k15_11 = interpolator.k15_11;
- k15_12 = interpolator.k15_12;
- k15_13 = interpolator.k15_13;
- k15_14 = interpolator.k15_14;
- c16 = interpolator.c16;
- k16_01 = interpolator.k16_01;
- k16_06 = interpolator.k16_06;
- k16_07 = interpolator.k16_07;
- k16_08 = interpolator.k16_08;
- k16_09 = interpolator.k16_09;
- k16_10 = interpolator.k16_10;
- k16_11 = interpolator.k16_11;
- k16_12 = interpolator.k16_12;
- k16_13 = interpolator.k16_13;
- k16_14 = interpolator.k16_14;
- k16_15 = interpolator.k16_15;
-
d = MathArrays.buildArray(integrator.getField(), 4, -1);
for (int i = 0; i < d.length; ++i) {
d[i] = interpolator.d[i].clone();
@@ -403,72 +209,13 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
}
/** {@inheritDoc} */
- @Override
- public void storeState(final FieldODEStateAndDerivative<T> state) {
- super.storeState(state);
- vectorsInitialized = false;
- }
-
- /** {@inheritDoc} */
+ @SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH)
throws MaxCountExceededException {
- if (! vectorsInitialized) {
-
- if (v == null) {
- v = MathArrays.buildArray(time.getField(), 7, previousState.length);
- }
-
- // perform the last evaluations if they have not been done yet
- finalizeStep();
-
- // compute the interpolation vectors for this time step
- for (int i = 0; i < previousState.length; ++i) {
- final T yDot01 = yDotK[0][i];
- final T yDot06 = yDotK[5][i];
- final T yDot07 = yDotK[6][i];
- final T yDot08 = yDotK[7][i];
- final T yDot09 = yDotK[8][i];
- final T yDot10 = yDotK[9][i];
- final T yDot11 = yDotK[10][i];
- final T yDot12 = yDotK[11][i];
- final T yDot13 = yDotK[12][i];
- final T yDot14 = yDotKLast[0][i];
- final T yDot15 = yDotKLast[1][i];
- final T yDot16 = yDotKLast[2][i];
- v[0][i] = yDot01.multiply(b_01).
- add(yDot06.multiply(b_06)).
- add(yDot07.multiply(b_07)).
- add(yDot08.multiply(b_08)).
- add(yDot09.multiply(b_09)).
- add(yDot10.multiply(b_10)).
- add(yDot11.multiply(b_11)).
- add(yDot12.multiply(b_12));
- v[1][i] = yDot01.subtract(v[0][i]);
- v[2][i] = v[0][i].subtract(v[1][i]).subtract(yDotK[12][i]);
- for (int k = 0; k < d.length; ++k) {
- v[k+3][i] = yDot01.multiply(d[k][ 0]).
- add(yDot06.multiply(d[k][ 1])).
- add(yDot07.multiply(d[k][ 2])).
- add(yDot08.multiply(d[k][ 3])).
- add(yDot09.multiply(d[k][ 4])).
- add(yDot10.multiply(d[k][ 5])).
- add(yDot11.multiply(d[k][ 6])).
- add(yDot12.multiply(d[k][ 7])).
- add(yDot13.multiply(d[k][ 8])).
- add(yDot14.multiply(d[k][ 9])).
- add(yDot15.multiply(d[k][10])).
- add(yDot16.multiply(d[k][11]));
- }
- }
-
- vectorsInitialized = true;
-
- }
-
final T one = theta.getField().getOne();
final T eta = one.subtract(theta);
final T twoTheta = theta.multiply(2);
@@ -479,111 +226,32 @@ class DormandPrince853FieldStepInterpolator<T extends RealFieldElement<T>>
final T dot4 = theta2.multiply(theta.multiply(theta.multiply(5).subtract(8)).add(3));
final T dot5 = theta2.multiply(theta.multiply(theta.multiply(theta.multiply(-6).add(15)).subtract(12)).add(3));
final T dot6 = theta2.multiply(theta.multiply(theta.multiply(theta.multiply(theta.multiply(-7).add(18)).subtract(15)).add(4)));
- final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
- final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
-
- if ((previousState != null) && (theta.getReal() <= 0.5)) {
- for (int i = 0; i < interpolatedState.length; ++i) {
- interpolatedState[i] = previousState[i].
- add(theta.multiply(h.multiply(v[0][i].
- add(eta.multiply(v[1][i].
- add(theta.multiply(v[2][i].
- add(eta.multiply(v[3][i].
- add(theta.multiply(v[4][i].
- add(eta.multiply(v[5][i].
- add(theta.multiply(v[6][i])))))))))))))));
- interpolatedDerivatives[i] = v[0][i].
- add(dot1.multiply(v[1][i])).
- add(dot2.multiply(v[2][i])).
- add(dot3.multiply(v[3][i])).
- add(dot4.multiply(v[4][i])).
- add(dot5.multiply(v[5][i])).
- add(dot6.multiply(v[6][i]));
- }
+ final T[] interpolatedState;
+ final T[] interpolatedDerivatives;
+
+ if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
+ final T f0 = theta.multiply(h);
+ final T f1 = f0.multiply(eta);
+ final T f2 = f1.multiply(theta);
+ final T f3 = f2.multiply(eta);
+ final T f4 = f3.multiply(theta);
+ final T f5 = f4.multiply(eta);
+ final T f6 = f5.multiply(theta);
+ interpolatedState = previousStateLinearCombination(f0, f1, f2, f3, f4, f5, f6);
+ interpolatedDerivatives = derivativeLinearCombination(one, dot1, dot2, dot3, dot4, dot5, dot6);
} else {
- for (int i = 0; i < interpolatedState.length; ++i) {
- interpolatedState[i] = currentState[i].
- subtract(oneMinusThetaH.multiply(v[0][i].
- subtract(theta.multiply(v[1][i].
- add(theta.multiply(v[2][i].
- add(eta.multiply(v[3][i].
- add(theta.multiply(v[4][i].
- add(eta.multiply(v[5][i].
- add(theta.multiply(v[6][i]))))))))))))));
- interpolatedDerivatives[i] = v[0][i].
- add(dot1.multiply(v[1][i])).
- add(dot2.multiply(v[2][i])).
- add(dot3.multiply(v[3][i])).
- add(dot4.multiply(v[4][i])).
- add(dot5.multiply(v[5][i])).
- add(dot6.multiply(v[6][i]));
- }
+ final T f0 = oneMinusThetaH;
+ final T f1 = f0.multiply(theta).negate();
+ final T f2 = f1.multiply(theta);
+ final T f3 = f2.multiply(eta);
+ final T f4 = f3.multiply(theta);
+ final T f5 = f4.multiply(eta);
+ final T f6 = f5.multiply(theta);
+ interpolatedState = currentStateLinearCombination(f0, f1, f2, f3, f4, f5, f6);
+ interpolatedDerivatives = derivativeLinearCombination(one, dot1, dot2, dot3, dot4, dot5, dot6);
}
- return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
-
- }
-
- /** {@inheritDoc} */
- @Override
- protected void doFinalize() throws MaxCountExceededException {
-
- if (currentState == null) {
- // we are finalizing an uninitialized instance
- return;
- }
-
- T s;
- final T pT = getGlobalPreviousState().getTime();
- final T[] yTmp = MathArrays.buildArray(pT.getField(), currentState.length);
-
- // k14
- for (int j = 0; j < currentState.length; ++j) {
- s = yDotK[ 0][j].multiply(k14_01).
- add(yDotK[ 5][j].multiply(k14_06)).
- add(yDotK[ 6][j].multiply(k14_07)).
- add(yDotK[ 7][j].multiply(k14_08)).
- add(yDotK[ 8][j].multiply(k14_09)).
- add(yDotK[ 9][j].multiply(k14_10)).
- add(yDotK[10][j].multiply(k14_11)).
- add(yDotK[11][j].multiply(k14_12)).
- add(yDotK[12][j].multiply(k14_13));
- yTmp[j] = currentState[j].add(h.multiply(s));
- }
- yDotKLast[0] = integrator.computeDerivatives(pT.add(h.multiply(c14)), yTmp);
-
- // k15
- for (int j = 0; j < currentState.length; ++j) {
- s = yDotK[ 0][j].multiply(k15_01).
- add(yDotK[ 5][j].multiply(k15_06)).
- add(yDotK[ 6][j].multiply(k15_07)).
- add(yDotK[ 7][j].multiply(k15_08)).
- add(yDotK[ 8][j].multiply(k15_09)).
- add(yDotK[ 9][j].multiply(k15_10)).
- add(yDotK[10][j].multiply(k15_11)).
- add(yDotK[11][j].multiply(k15_12)).
- add(yDotK[12][j].multiply(k15_13)).
- add(yDotKLast[0][j].multiply(k15_14));
- yTmp[j] = currentState[j].add(h.multiply(s));
- }
- yDotKLast[1] = integrator.computeDerivatives(pT.add(h.multiply(c15)), yTmp);
-
- // k16
- for (int j = 0; j < currentState.length; ++j) {
- s = yDotK[ 0][j].multiply(k16_01).
- add(yDotK[ 5][j].multiply(k16_06)).
- add(yDotK[ 6][j].multiply(k16_07)).
- add(yDotK[ 7][j].multiply(k16_08)).
- add(yDotK[ 8][j].multiply(k16_09)).
- add(yDotK[ 9][j].multiply(k16_10)).
- add(yDotK[10][j].multiply(k16_11)).
- add(yDotK[11][j].multiply(k16_12)).
- add(yDotK[12][j].multiply(k16_13)).
- add(yDotKLast[0][j].multiply(k16_14)).
- add(yDotKLast[0][j].multiply(k16_15));
- yTmp[j] = currentState[j].add(h.multiply(s));
- }
- yDotKLast[2] = integrator.computeDerivatives(pT.add(h.multiply(c16)), yTmp);
+ return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
index 8da12a9..5f648f8 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/EmbeddedRungeKuttaFieldIntegrator.java
@@ -23,7 +23,6 @@ import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NoBracketingException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.ode.FieldExpandableODE;
import org.apache.commons.math3.ode.FieldODEState;
@@ -71,8 +70,8 @@ import org.apache.commons.math3.util.MathUtils;
public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
extends AdaptiveStepsizeFieldIntegrator<T> {
- /** Indicator for <i>fsal</i> methods. */
- private final boolean fsal;
+ /** Index of the pre-computed derivative for <i>fsal</i> methods. */
+ private final int fsal;
/** Time steps from Butcher array (without the first zero). */
private final T[] c;
@@ -98,7 +97,8 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
/** Build a Runge-Kutta integrator with the given Butcher array.
* @param field field to which the time and state vector elements belong
* @param name name of the method
- * @param fsal indicate that the method is an <i>fsal</i>
+ * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
+ * or -1 if method is not <i>fsal</i>
* @param minStep minimal step (sign is irrelevant, regardless of
* integration direction, forward or backward), the last step can
* be smaller than this
@@ -108,7 +108,7 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
* @param scalAbsoluteTolerance allowed absolute error
* @param scalRelativeTolerance allowed relative error
*/
- protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final boolean fsal,
+ protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
final double minStep, final double maxStep,
final double scalAbsoluteTolerance,
final double scalRelativeTolerance) {
@@ -132,7 +132,8 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
/** Build a Runge-Kutta integrator with the given Butcher array.
* @param field field to which the time and state vector elements belong
* @param name name of the method
- * @param fsal indicate that the method is an <i>fsal</i>
+ * @param fsal index of the pre-computed derivative for <i>fsal</i> methods
+ * or -1 if method is not <i>fsal</i>
* @param minStep minimal step (must be positive even for backward
* integration), the last step can be smaller than this
* @param maxStep maximal step (must be positive even for backward
@@ -140,7 +141,7 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
* @param vecAbsoluteTolerance allowed absolute error
* @param vecRelativeTolerance allowed relative error
*/
- protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final boolean fsal,
+ protected EmbeddedRungeKuttaFieldIntegrator(final Field<T> field, final String name, final int fsal,
final double minStep, final double maxStep,
final double[] vecAbsoluteTolerance,
final double[] vecRelativeTolerance) {
@@ -195,18 +196,11 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
protected abstract T[] getB();
/** Create an interpolator.
- * @param rkIntegrator integrator being used
- * @param y reference to the integrator array holding the state at
- * the end of the step
- * @param yDotArray reference to the integrator array holding all the
- * intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
* @return external weights for the high order method from Butcher array
*/
- protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(AbstractFieldIntegrator<T> rkIntegrator,
- T[] y, T[][] yDotArray,
- boolean forward,
+ protected abstract RungeKuttaFieldStepInterpolator<T> createInterpolator(boolean forward,
FieldEquationsMapper<T> mapper);
/** Get the order of the method.
* @return order of the method
@@ -248,7 +242,7 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
// set up an interpolator sharing the integrator arrays
final RungeKuttaFieldStepInterpolator<T> interpolator =
- createInterpolator(this, y0, yDotK, forward, equations.getMapper());
+ createInterpolator(forward, equations.getMapper());
interpolator.storeState(stepStart);
// set up integration control objects
@@ -328,11 +322,12 @@ public abstract class EmbeddedRungeKuttaFieldIntegrator<T extends RealFieldEleme
}
}
- final T stepEnd = stepStart.getTime().add(stepSize);
- final T[] yDotTmp = fsal ? yDotK[stages - 1] : computeDerivatives(stepEnd, yTmp);
+ final T stepEnd = stepStart.getTime().add(stepSize);
+ final T[] yDotTmp = (fsal >= 0) ? yDotK[fsal] : computeDerivatives(stepEnd, yTmp);
final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<T>(stepEnd, yTmp, yDotTmp);
// local error is small enough: accept the step, trigger events and step handlers
+ interpolator.setSlopes(yDotK);
interpolator.storeState(stateTmp);
System.arraycopy(yTmp, 0, y, 0, y0.length);
stepStart = acceptStep(interpolator, finalTime);
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java
index 7fb32a1..4c58e09 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.util.MathArrays;
@@ -86,10 +85,8 @@ public class EulerFieldIntegrator<T extends RealFieldElement<T>> extends RungeKu
/** {@inheritDoc} */
@Override
protected EulerFieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new EulerFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new EulerFieldStepInterpolator<T>(this, forward, mapper);
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java
index dba2312..9bff5f2 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/EulerFieldStepInterpolator.java
@@ -21,7 +21,6 @@ import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
-import org.apache.commons.math3.util.MathArrays;
/**
* This class implements a linear interpolator for step.
@@ -52,17 +51,13 @@ class EulerFieldStepInterpolator<T extends RealFieldElement<T>>
/** Simple constructor.
* @param rkIntegrator integrator being used
- * @param y reference to the integrator array holding the state at
- * the end of the step
- * @param yDotArray reference to the integrator array holding all the
- * intermediate slopes
* @param forward integration direction indicator
* @param mapper equations mapper for the all equations
*/
EulerFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
- final T[] y, final T[][] yDotArray, final boolean forward,
+ final boolean forward,
final FieldEquationsMapper<T> mapper) {
- super(rkIntegrator, y, yDotArray, forward, mapper);
+ super(rkIntegrator, forward, mapper);
}
/** Copy constructor.
@@ -81,22 +76,23 @@ class EulerFieldStepInterpolator<T extends RealFieldElement<T>>
}
/** {@inheritDoc} */
+ @SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
final T oneMinusThetaH) {
- final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
- if ((previousState != null) && (theta.getReal() <= 0.5)) {
- for (int i = 0; i < previousState.length; ++i) {
- interpolatedState[i] = previousState[i].add(theta.multiply(h).multiply(yDotK[0][i]));
- }
+ final T[] interpolatedState;
+ final T[] interpolatedDerivatives;
+ if ((getGlobalPreviousState() != null) && (theta.getReal() <= 0.5)) {
+ interpolatedState = previousStateLinearCombination(theta.multiply(h));
+ interpolatedDerivatives = derivativeLinearCombination(time.getField().getOne());
} else {
- for (int i = 0; i < previousState.length; ++i) {
- interpolatedState[i] = currentState[i].subtract(oneMinusThetaH.multiply(yDotK[0][i]));
- }
+ interpolatedState = currentStateLinearCombination(oneMinusThetaH.negate());
+ interpolatedDerivatives = derivativeLinearCombination(time.getField().getOne());
}
- return new FieldODEStateAndDerivative<T>(time, interpolatedState, yDotK[0]);
+ return new FieldODEStateAndDerivative<T>(time, interpolatedState, interpolatedDerivatives);
+
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java
index d8460c2..659eee8 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldIntegrator.java
@@ -19,7 +19,6 @@ package org.apache.commons.math3.ode.nonstiff;
import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
-import org.apache.commons.math3.ode.AbstractFieldIntegrator;
import org.apache.commons.math3.ode.FieldEquationsMapper;
import org.apache.commons.math3.util.MathArrays;
@@ -111,10 +110,8 @@ public class GillFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected GillFieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new GillFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new GillFieldStepInterpolator<T>(this, forward, mapper);
}
}