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/01 12:23:45 UTC
[1/4] [math] Fixed wrong state reset in field ode.
Repository: commons-math
Updated Branches:
refs/heads/field-ode a40d822c0 -> 87664da98
Fixed wrong state reset in field ode.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/8180c9e5
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/8180c9e5
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/8180c9e5
Branch: refs/heads/field-ode
Commit: 8180c9e5e444c661ec99a41444590f95bafcb6b9
Parents: a40d822
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue Dec 1 12:19:24 2015 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 1 12:19:24 2015 +0100
----------------------------------------------------------------------
.../math3/ode/AbstractFieldIntegrator.java | 20 ++++++++++----------
.../math3/ode/events/FieldEventState.java | 18 ++++++++++++------
2 files changed, 22 insertions(+), 16 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/8180c9e5/src/main/java/org/apache/commons/math3/ode/AbstractFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/AbstractFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/AbstractFieldIntegrator.java
index c15b827..b39ed3b 100644
--- a/src/main/java/org/apache/commons/math3/ode/AbstractFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/AbstractFieldIntegrator.java
@@ -340,17 +340,17 @@ public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> imp
return eventState;
}
- boolean needReset = false;
+ FieldODEState<T> newState = null;
for (final FieldEventState<T> state : eventsStates) {
- needReset = needReset || state.reset(eventState);
- }
- if (needReset) {
- // some event handler has triggered changes that
- // invalidate the derivatives, we need to recompute them
- final T[] y = equations.getMapper().mapState(eventState);
- final T[] yDot = computeDerivatives(eventState.getTime(), y);
- resetOccurred = true;
- return equations.getMapper().mapStateAndDerivative(eventState.getTime(), y, yDot);
+ newState = state.reset(eventState);
+ if (newState != null) {
+ // some event handler has triggered changes that
+ // invalidate the derivatives, we need to recompute them
+ final T[] y = equations.getMapper().mapState(newState);
+ final T[] yDot = computeDerivatives(newState.getTime(), y);
+ resetOccurred = true;
+ return equations.getMapper().mapStateAndDerivative(newState.getTime(), y, yDot);
+ }
}
// prepare handling of the remaining part of the step
http://git-wip-us.apache.org/repos/asf/commons-math/blob/8180c9e5/src/main/java/org/apache/commons/math3/ode/events/FieldEventState.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/events/FieldEventState.java b/src/main/java/org/apache/commons/math3/ode/events/FieldEventState.java
index 2210b4c..044b889 100644
--- a/src/main/java/org/apache/commons/math3/ode/events/FieldEventState.java
+++ b/src/main/java/org/apache/commons/math3/ode/events/FieldEventState.java
@@ -23,6 +23,7 @@ import org.apache.commons.math3.analysis.solvers.AllowedSolution;
import org.apache.commons.math3.analysis.solvers.BracketedRealFieldUnivariateSolver;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.ode.FieldODEState;
import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
import org.apache.commons.math3.ode.sampling.FieldStepInterpolator;
import org.apache.commons.math3.util.FastMath;
@@ -325,22 +326,27 @@ public class FieldEventState<T extends RealFieldElement<T>> {
/** Let the event handler reset the state if it wants.
* @param state state at the beginning of the next step
- * @return true if the integrator should reset the derivatives too
+ * @return reset state (may by the same as initial state if only
+ * derivatives should be reset), or null if nothing is reset
*/
- public boolean reset(final FieldODEStateAndDerivative<T> state) {
+ public FieldODEState<T> reset(final FieldODEStateAndDerivative<T> state) {
if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) {
- return false;
+ return null;
}
+ final FieldODEState<T> newState;
if (nextAction == Action.RESET_STATE) {
- handler.resetState(state);
+ newState = handler.resetState(state);
+ } else if (nextAction == Action.RESET_DERIVATIVES) {
+ newState = state;
+ } else {
+ newState = null;
}
pendingEvent = false;
pendingEventTime = null;
- return (nextAction == Action.RESET_STATE) ||
- (nextAction == Action.RESET_DERIVATIVES);
+ return newState;
}
[3/4] [math] One step towards immutable step interpolators.
Posted by lu...@apache.org.
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/field-ode
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);
}
}
[4/4] [math] Starting tests on field ODE!
Posted by lu...@apache.org.
Starting tests on field ODE!
At least we can now run some code ...
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/87664da9
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/87664da9
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/87664da9
Branch: refs/heads/field-ode
Commit: 87664da98fb8586058973f4f9a50a5d231944f14
Parents: 07e2b94
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue Dec 1 12:23:20 2015 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 1 12:23:20 2015 +0100
----------------------------------------------------------------------
.../ode/nonstiff/EulerFieldIntegratorTest.java | 247 +++++++++++++++++++
1 file changed, 247 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/87664da9/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegratorTest.java
new file mode 100644
index 0000000..2e41788
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerFieldIntegratorTest.java
@@ -0,0 +1,247 @@
+/*
+ * 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.math3.ode.nonstiff;
+
+
+import java.lang.reflect.Array;
+
+import org.apache.commons.math3.Field;
+import org.apache.commons.math3.RealFieldElement;
+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.FieldExpandableODE;
+import org.apache.commons.math3.ode.FieldFirstOrderDifferentialEquations;
+import org.apache.commons.math3.ode.FieldFirstOrderIntegrator;
+import org.apache.commons.math3.ode.FieldODEState;
+import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
+import org.apache.commons.math3.ode.TestFieldProblem1;
+import org.apache.commons.math3.ode.TestFieldProblem2;
+import org.apache.commons.math3.ode.TestFieldProblem3;
+import org.apache.commons.math3.ode.TestFieldProblem4;
+import org.apache.commons.math3.ode.TestFieldProblem5;
+import org.apache.commons.math3.ode.TestFieldProblem6;
+import org.apache.commons.math3.ode.TestFieldProblemAbstract;
+import org.apache.commons.math3.ode.TestFieldProblemHandler;
+import org.apache.commons.math3.ode.events.FieldEventHandler;
+import org.apache.commons.math3.ode.sampling.FieldStepHandler;
+import org.apache.commons.math3.ode.sampling.FieldStepInterpolator;
+import org.apache.commons.math3.util.Decimal64Field;
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.MathArrays;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class EulerFieldIntegratorTest {
+
+ @Test(expected=DimensionMismatchException.class)
+ public void testDimensionCheck()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ doTestDimensionCheck(Decimal64Field.getInstance());
+ }
+
+ private <T extends RealFieldElement<T>> void doTestDimensionCheck(Field<T> field)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ TestFieldProblemAbstract<T> pb = new TestFieldProblem1<T>(field);
+ FieldExpandableODE<T> ode = new FieldExpandableODE<T>(pb);
+ FieldFirstOrderIntegrator<T> integ =
+ new EulerFieldIntegrator<T>(field, field.getZero().add(0.01));
+ FieldODEState<T> start =
+ new FieldODEState<T>(field.getZero().add(0.0),
+ MathArrays.buildArray(field, pb.getDimension() + 10));
+ integ.integrate(ode, start, field.getZero().add(1.0));
+ }
+
+ @Test
+ public void testDecreasingSteps()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ dotTestDecreasingSteps(Decimal64Field.getInstance());
+ }
+
+ private <T extends RealFieldElement<T>> void dotTestDecreasingSteps(Field<T> field)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ @SuppressWarnings("unchecked")
+ TestFieldProblemAbstract<T>[] allProblems =
+ (TestFieldProblemAbstract<T>[]) Array.newInstance(TestFieldProblemAbstract.class, 6);
+ allProblems[0] = new TestFieldProblem1<T>(field);
+ allProblems[1] = new TestFieldProblem2<T>(field);
+ allProblems[2] = new TestFieldProblem3<T>(field);
+ allProblems[3] = new TestFieldProblem4<T>(field);
+ allProblems[4] = new TestFieldProblem5<T>(field);
+ allProblems[5] = new TestFieldProblem6<T>(field);
+ for (TestFieldProblemAbstract<T> pb : allProblems) {
+
+ T previousValueError = null;
+ T previousTimeError = null;
+ for (int i = 4; i < 8; ++i) {
+
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(FastMath.pow(2.0, -i));
+
+ FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(field, step);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ FieldEventHandler<T>[] functions = pb.getEventsHandlers();
+ for (int l = 0; l < functions.length; ++l) {
+ integ.addEventHandler(functions[l],
+ Double.POSITIVE_INFINITY, 1.0e-6 * step.getReal(), 1000);
+ }
+ FieldODEStateAndDerivative<T> stop = integ.integrate(new FieldExpandableODE<T>(pb),
+ pb.getInitialState(),
+ pb.getFinalTime());
+ if (functions.length == 0) {
+ Assert.assertEquals(pb.getFinalTime().getReal(), stop.getTime().getReal(), 1.0e-10);
+ }
+
+ T valueError = handler.getMaximalValueError();
+ if (i > 4) {
+ Assert.assertTrue(valueError.subtract(previousValueError.abs()).getReal() < 0);
+ }
+ previousValueError = valueError;
+
+ T timeError = handler.getMaximalTimeError();
+ if (i > 4) {
+ Assert.assertTrue(timeError.subtract(previousTimeError.abs()).getReal() <= 0);
+ }
+ previousTimeError = timeError;
+
+ }
+
+ }
+
+ }
+
+ @Test
+ public void testSmallStep()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ doTestSmallStep(Decimal64Field.getInstance());
+ }
+
+ private <T extends RealFieldElement<T>> void doTestSmallStep(Field<T> field)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001);
+
+ FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(field, step);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertTrue(handler.getLastError().getReal() < 2.0e-4);
+ Assert.assertTrue(handler.getMaximalValueError().getReal() < 1.0e-3);
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
+ Assert.assertEquals("Euler", integ.getName());
+
+ }
+
+ @Test
+ public void testBigStep()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ doTestBigStep(Decimal64Field.getInstance());
+ }
+
+ private <T extends RealFieldElement<T>> void doTestBigStep(Field<T> field)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ TestFieldProblem1<T> pb = new TestFieldProblem1<T>(field);
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.2);
+
+ FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(field, step);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<T>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertTrue(handler.getLastError().getReal() > 0.01);
+ Assert.assertTrue(handler.getMaximalValueError().getReal() > 0.2);
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
+
+ }
+
+ @Test
+ public void testBackward()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ doTestBackward(Decimal64Field.getInstance());
+ }
+
+ private <T extends RealFieldElement<T>> void doTestBackward(Field<T> field)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ TestFieldProblem5<T> pb = new TestFieldProblem5<T>(field);
+ T step = pb.getFinalTime().subtract(pb.getInitialState().getTime()).multiply(0.001).abs();
+
+ FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(field, step);
+ TestFieldProblemHandler<T> handler = new TestFieldProblemHandler<T>(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(new FieldExpandableODE<>(pb), pb.getInitialState(), pb.getFinalTime());
+
+ Assert.assertTrue(handler.getLastError().getReal() < 0.45);
+ Assert.assertTrue(handler.getMaximalValueError().getReal() < 0.45);
+ Assert.assertEquals(0, handler.getMaximalTimeError().getReal(), 1.0e-12);
+ Assert.assertEquals("Euler", integ.getName());
+ }
+
+ @Test
+ public void testStepSize()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ doTestStepSize(Decimal64Field.getInstance());
+ }
+
+ private <T extends RealFieldElement<T>> void doTestStepSize(Field<T> field)
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ final T step = field.getZero().add(1.23456);
+ FieldFirstOrderIntegrator<T> integ = new EulerFieldIntegrator<T>(field, step);
+ integ.addStepHandler(new FieldStepHandler<T>() {
+ public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) {
+ if (! isLast) {
+ Assert.assertEquals(step.getReal(),
+ interpolator.getCurrentState().getTime().subtract(interpolator.getPreviousState().getTime()).getReal(),
+ 1.0e-12);
+ }
+ }
+ public void init(FieldODEStateAndDerivative<T> s0, T t) {
+ }
+ });
+ integ.integrate(new FieldExpandableODE<T>(new FieldFirstOrderDifferentialEquations<T>() {
+ public void init(T t0, T[] y0, T t) {
+ }
+ public T[] computeDerivatives(T t, T[] y) {
+ T[] dot = MathArrays.buildArray(t.getField(), 1);
+ dot[0] = t.getField().getOne();
+ return dot;
+ }
+ public int getDimension() {
+ return 1;
+ }
+ }), new FieldODEState<T>(field.getZero(), MathArrays.buildArray(field, 1)), field.getZero().add(5.0));
+ }
+
+}
[2/4] [math] One step towards immutable step interpolators.
Posted by lu...@apache.org.
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java
index 8a6b38e..8990dc5 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/GillFieldStepInterpolator.java
@@ -22,7 +22,6 @@ 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.FastMath;
-import org.apache.commons.math3.util.MathArrays;
/**
* This class implements a step interpolator for the Gill fourth
@@ -68,17 +67,13 @@ class GillFieldStepInterpolator<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
*/
GillFieldStepInterpolator(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.
@@ -98,6 +93,7 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
+ @SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
@@ -111,48 +107,30 @@ class GillFieldStepInterpolator<T extends RealFieldElement<T>>
final T coeffDot2 = cDot23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeffDot3 = cDot23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeffDot4 = theta.multiply(twoTheta.subtract(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)) {
- final T s = theta.multiply(h).divide(6.0);
- final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
- final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(6)).add(6));
- final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2);
- final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2);
- 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 yDot2 = yDotK[1][i];
- final T yDot3 = yDotK[2][i];
- final T yDot4 = yDotK[3][i];
- interpolatedState[i] = previousState[i].
- add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2)).
- add(coeff3.multiply(yDot3)).add(coeff4.multiply(yDot4));
- interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
- add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
- }
+ final T[] interpolatedState;
+ final T[] interpolatedDerivatives;
+
+ if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
+ final T s = theta.multiply(h).divide(6.0);
+ final T c23 = s.multiply(theta.multiply(6).subtract(fourTheta2));
+ final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(6)).add(6));
+ final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2);
+ final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2);
+ final T coeff4 = s.multiply(fourTheta2.subtract(theta.multiply(3)));
+ interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4);
} else {
- final T s = oneMinusThetaH.divide(6.0);
- final T c23 = s .multiply(twoTheta.add(2).subtract(fourTheta2));
+ final T s = oneMinusThetaH.divide(-6.0);
+ final T c23 = s.multiply(twoTheta.add(2).subtract(fourTheta2));
final T coeff1 = s.multiply(fourTheta2.subtract(theta.multiply(5)).add(1));
final T coeff2 = c23.multiply(ONE_MINUS_INV_SQRT_2);
final T coeff3 = c23.multiply(ONE_PLUS_INV_SQRT_2);
final T coeff4 = s.multiply(fourTheta2.add(theta).add(1));
- for (int i = 0; i < interpolatedState.length; ++i) {
- final T yDot1 = yDotK[0][i];
- final T yDot2 = yDotK[1][i];
- final T yDot3 = yDotK[2][i];
- final T yDot4 = yDotK[3][i];
- interpolatedState[i] = currentState[i].
- subtract(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2)).
- subtract(coeff3.multiply(yDot3)).subtract(coeff4.multiply(yDot4));
- interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
- add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
- }
+ interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, 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/HighamHall54FieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java
index 092393a..8ee90b7 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldIntegrator.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;
@@ -64,7 +63,7 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double scalAbsoluteTolerance,
final double scalRelativeTolerance) {
- super(field, METHOD_NAME, false,
+ super(field, METHOD_NAME, -1,
minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
e = MathArrays.buildArray(field, 7);
e[0] = fraction(-1, 20);
@@ -92,7 +91,7 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>>
final double minStep, final double maxStep,
final double[] vecAbsoluteTolerance,
final double[] vecRelativeTolerance) {
- super(field, METHOD_NAME, false,
+ super(field, METHOD_NAME, -1,
minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
e = MathArrays.buildArray(field, 7);
e[0] = fraction(-1, 20);
@@ -165,10 +164,8 @@ public class HighamHall54FieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected HighamHall54FieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new HighamHall54FieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new HighamHall54FieldStepInterpolator<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/HighamHall54FieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldStepInterpolator.java
index 3ee52a9..3f8499c 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/HighamHall54FieldStepInterpolator.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
@@ -38,17 +37,13 @@ class HighamHall54FieldStepInterpolator<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
*/
HighamHall54FieldStepInterpolator(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.
@@ -68,72 +63,44 @@ class HighamHall54FieldStepInterpolator<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 bDot0 = theta.multiply(theta.multiply(theta.multiply( -10.0 ).add( 16.0 )).add(-15.0 / 2.0)).add(1);
+ final T bDot1 = time.getField().getZero();
final T bDot2 = theta.multiply(theta.multiply(theta.multiply( 135.0 / 2.0).add(-729.0 / 8.0)).add(459.0 / 16.0));
final T bDot3 = theta.multiply(theta.multiply(theta.multiply(-120.0 ).add( 152.0 )).add(-44.0 ));
final T bDot4 = theta.multiply(theta.multiply(theta.multiply( 125.0 / 2.0).add(-625.0 / 8.0)).add(375.0 / 16.0));
final T bDot5 = theta.multiply( 5.0 / 8.0).multiply(theta.multiply(2).subtract(1));
- 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)) {
+ if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T hTheta = h.multiply(theta);
final T b0 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply( -5.0 / 2.0).add( 16.0 / 3.0)).add(-15.0 / 4.0)).add(1));
+ final T b1 = time.getField().getZero();
final T b2 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(135.0 / 8.0).add(-243.0 / 8.0)).add(459.0 / 32.0)));
final T b3 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(-30.0 ).add( 152.0 / 3.0)).add(-22.0 )));
final T b4 = hTheta.multiply(theta.multiply(theta.multiply(theta.multiply(125.0 / 8.0).add(-625.0 / 24.0)).add(375.0 / 32.0)));
final T b5 = hTheta.multiply(theta.multiply(theta.multiply( 5.0 / 12.0)).add( -5.0 / 16.0));
- for (int i = 0; i < interpolatedState.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];
- interpolatedState[i] = previousState[i].
- add(b0.multiply(yDot0)).
- add(b2.multiply(yDot2)).
- add(b3.multiply(yDot3)).
- add(b4.multiply(yDot4)).
- add(b5.multiply(yDot5));
- interpolatedDerivatives[i] = bDot0.multiply(yDot0).
- add(bDot2.multiply(yDot2)).
- add(bDot3.multiply(yDot3)).
- add(bDot4.multiply(yDot4)).
- add(bDot5.multiply(yDot5));
- }
+ interpolatedState = previousStateLinearCombination(b0, b1, b2, b3, b4, b5);
+ interpolatedDerivatives = derivativeLinearCombination(bDot0, bDot1, bDot2, bDot3, bDot4, bDot5);
} else {
final T theta2 = theta.multiply(theta);
final T b0 = h.multiply( theta.multiply(theta.multiply(theta.multiply(theta.multiply(-5.0 / 2.0).add( 16.0 / 3.0)).add( -15.0 / 4.0)).add( 1.0 )).add( -1.0 / 12.0));
+ final T b1 = time.getField().getZero();
final T b2 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( 135.0 / 8.0 ).add(-243.0 / 8.0)).add(459.0 / 32.0)).add( -27.0 / 32.0));
final T b3 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( -30.0 ).add( 152.0 / 3.0)).add(-22.0 )).add( 4.0 / 3.0));
final T b4 = h.multiply(theta2.multiply(theta.multiply(theta.multiply( 125.0 / 8.0 ).add(-625.0 / 24.0)).add(375.0 / 32.0)).add(-125.0 / 96.0));
final T b5 = h.multiply(theta2.multiply(theta.multiply( 5.0 / 12.0 ).add(-5.0 / 16.0)).add( -5.0 / 48.0));
- for (int i = 0; i < interpolatedState.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];
- interpolatedState[i] = currentState[i].
- add(b0.multiply(yDot0)).
- add(b2.multiply(yDot2)).
- add(b3.multiply(yDot3)).
- add(b4.multiply(yDot4)).
- add(b5.multiply(yDot5));
- interpolatedDerivatives[i] = bDot0.multiply(yDot0).
- add(bDot2.multiply(yDot2)).
- add(bDot3.multiply(yDot3)).
- add(bDot4.multiply(yDot4)).
- add(bDot5.multiply(yDot5));
- }
+ interpolatedState = currentStateLinearCombination(b0, b1, b2, b3, b4, b5);
+ interpolatedDerivatives = derivativeLinearCombination(bDot0, bDot1, bDot2, bDot3, bDot4, bDot5);
}
- 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/LutherFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java
index 363e042..c36e5dd 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldIntegrator.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;
@@ -138,10 +137,8 @@ public class LutherFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected LutherFieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new LutherFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new LutherFieldStepInterpolator<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/LutherFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java
index 8a0d407..6f643d5 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/LutherFieldStepInterpolator.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
@@ -83,17 +82,13 @@ class LutherFieldStepInterpolator<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
*/
LutherFieldStepInterpolator(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);
final T q = rkIntegrator.getField().getOne().multiply(21).sqrt();
c5a = q.multiply( -49).add( -49);
c5b = q.multiply( 287).add( 392);
@@ -143,6 +138,7 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
+ @SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
@@ -192,86 +188,42 @@ class LutherFieldStepInterpolator<T extends RealFieldElement<T>>
// At the end, we get the b_i as polynomials in theta.
final T coeffDot1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 ).add( -47 )).add( 36 )).add( -54 / 5.0)).add(1);
- // not really needed as it is zero: final T coeffDot2 = theta.getField().getZero();
+ final T coeffDot2 = theta.getField().getZero();
final T coeffDot3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 ).add(-608 / 3.0)).add( 320 / 3.0 )).add(-208 / 15.0));
final T coeffDot4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -567 / 5.0).add( 972 / 5.0)).add( -486 / 5.0 )).add( 324 / 25.0));
final T coeffDot5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(5)).add(c5b.divide(15))).add(c5c.divide(30))).add(c5d.divide(150)));
final T coeffDot6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c6a.divide(5)).add(c6b.divide(15))).add(c6c.divide(30))).add(c6d.divide(150)));
final T coeffDot7 = theta.multiply(theta.multiply(theta.multiply( 3 )).add( -3 )).add( 3 / 5.0);
- 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)) {
-
- final T coeff1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1);
- // not really needed as it is zero: final T coeff2 = theta.getField().getZero();
- final T coeff3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0));
- final T coeff4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0));
- final T coeff5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300)));
- final T coeff6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300)));
- final T coeff7 = theta.multiply(theta.multiply(theta.multiply( 3 / 4.0)).add( -1 )).add( 3 / 10.0);
- for (int i = 0; i < interpolatedState.length; ++i) {
- final T yDot1 = yDotK[0][i];
- // not really needed as associated coefficients are zero: final T yDot2 = yDotK[1][i];
- final T yDot3 = yDotK[2][i];
- final T yDot4 = yDotK[3][i];
- final T yDot5 = yDotK[4][i];
- final T yDot6 = yDotK[5][i];
- final T yDot7 = yDotK[6][i];
- interpolatedState[i] = previousState[i].
- add(theta.multiply(h).
- multiply( coeff1.multiply(yDot1).
- // not really needed as it is zero: add(coeff2.multiply(yDot2)).
- add(coeff3.multiply(yDot3)).
- add(coeff4.multiply(yDot4)).
- add(coeff5.multiply(yDot5)).
- add(coeff6.multiply(yDot6)).
- add(coeff7.multiply(yDot7))));
- interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).
- // not really needed as it is zero: add(coeffDot2.multiply(yDot2)).
- add(coeffDot3.multiply(yDot3)).
- add(coeffDot4.multiply(yDot4)).
- add(coeffDot5.multiply(yDot5)).
- add(coeffDot6.multiply(yDot6)).
- add(coeffDot7.multiply(yDot7));
- }
+ final T[] interpolatedState;
+ final T[] interpolatedDerivatives;
+
+ if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
+
+ final T s = theta.multiply(theta.multiply(h));
+ final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 21 / 5.0).add( -47 / 4.0)).add( 12 )).add( -27 / 5.0)).add(1);
+ final T coeff2 = s.getField().getZero();
+ final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 112 / 5.0).add(-152 / 3.0)).add( 320 / 9.0 )).add(-104 / 15.0));
+ final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply(-567 / 25.0).add( 243 / 5.0)).add( -162 / 5.0 )).add( 162 / 25.0));
+ final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c5b.divide(60))).add(c5c.divide(90))).add(c5d.divide(300)));
+ final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(c5a.divide(25)).add(c6b.divide(60))).add(c6c.divide(90))).add(c6d.divide(300)));
+ final T coeff7 = s.multiply(theta.multiply(theta.multiply( 3 / 4.0)).add( -1 )).add( 3 / 10.0);
+ interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
} else {
- final T coeff1 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add( -1 / 20.0);
- // not really needed as it is zero: final T coeff2 = theta.getField().getZero();
- final T coeff3 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0);
- final T coeff4 = theta.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0)));
- final T coeff5 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0);
- final T coeff6 = theta.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0);
- final T coeff7 = theta.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0);
- for (int i = 0; i < interpolatedState.length; ++i) {
- final T yDot1 = yDotK[0][i];
- // not really needed as associated coefficients are zero: final T yDot2 = yDotK[1][i];
- final T yDot3 = yDotK[2][i];
- final T yDot4 = yDotK[3][i];
- final T yDot5 = yDotK[4][i];
- final T yDot6 = yDotK[5][i];
- final T yDot7 = yDotK[6][i];
- interpolatedState[i] = currentState[i].
- add(oneMinusThetaH.
- multiply( coeff1.multiply(yDot1).
- // not really needed as it is zero: add(coeff2.multiply(yDot2)).
- add(coeff3.multiply(yDot3)).
- add(coeff4.multiply(yDot4)).
- add(coeff5.multiply(yDot5)).
- add(coeff6.multiply(yDot6)).
- add(coeff7.multiply(yDot7))));
- interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).
- // not really needed as it is zero: add(coeffDot2.multiply(yDot2)).
- add(coeffDot3.multiply(yDot3)).
- add(coeffDot4.multiply(yDot4)).
- add(coeffDot5.multiply(yDot5)).
- add(coeffDot6.multiply(yDot6)).
- add(coeffDot7.multiply(yDot7));
- }
+ final T s = oneMinusThetaH.multiply(theta);
+ final T coeff1 = s.multiply(theta.multiply(theta.multiply(theta.multiply( -21 / 5.0).add( 151 / 20.0)).add( -89 / 20.0)).add( 19 / 20.0)).add( -1 / 20.0);
+ final T coeff2 = s.getField().getZero();
+ final T coeff3 = s.multiply(theta.multiply(theta.multiply(theta.multiply(-112 / 5.0).add( 424 / 15.0)).add( -328 / 45.0)).add( -16 / 45.0)).add(-16 / 45.0);
+ final T coeff4 = s.multiply(theta.multiply(theta.multiply(theta.multiply( 567 / 25.0).add( -648 / 25.0)).add( 162 / 25.0)));
+ final T coeff5 = s.multiply(theta.multiply(theta.multiply(theta.multiply(d5a.divide(25)).add(d5b.divide(300))).add(d5c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0);
+ final T coeff6 = s.multiply(theta.multiply(theta.multiply(theta.multiply(d6a.divide(25)).add(d6b.divide(300))).add(d6c.divide(900))).add( -49 / 180.0)).add(-49 / 180.0);
+ final T coeff7 = s.multiply(theta.multiply(theta.multiply( -3 / 4.0 ).add( 1 / 4.0)).add( -1 / 20.0)).add( -1 / 20.0);
+ interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, coeff7);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4, coeffDot5, coeffDot6, coeffDot7);
}
- 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/MidpointFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldIntegrator.java
index f74fdba..3cc4c71 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldIntegrator.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 MidpointFieldIntegrator<T extends RealFieldElement<T>> extends Rung
/** {@inheritDoc} */
@Override
protected MidpointFieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new MidpointFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new MidpointFieldStepInterpolator<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/MidpointFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldStepInterpolator.java
index 6569c4c..943533a 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/MidpointFieldStepInterpolator.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 second order
@@ -54,17 +53,13 @@ class MidpointFieldStepInterpolator<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
*/
MidpointFieldStepInterpolator(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.
@@ -84,37 +79,30 @@ class MidpointFieldStepInterpolator<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 coeffDot2 = theta.multiply(2);
- final T coeffDot1 = time.getField().getOne().subtract(coeffDot2);
- final T[] interpolatedState = MathArrays.buildArray(theta.getField(), previousState.length);
- final T[] interpolatedDerivatives = MathArrays.buildArray(theta.getField(), previousState.length);
+ final T coeffDot2 = theta.multiply(2);
+ final T coeffDot1 = time.getField().getOne().subtract(coeffDot2);
+ final T[] interpolatedState;
+ final T[] interpolatedDerivatives;
- if ((previousState != null) && (theta.getReal() <= 0.5)) {
+ if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T coeff1 = theta.multiply(oneMinusThetaH);
final T coeff2 = theta.multiply(theta).multiply(h);
- for (int i = 0; i < previousState.length; ++i) {
- final T yDot1 = yDotK[0][i];
- final T yDot2 = yDotK[1][i];
- interpolatedState[i] = previousState[i].add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2));
- interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2));
- }
+ interpolatedState = previousStateLinearCombination(coeff1, coeff2);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2);
} else {
final T coeff1 = oneMinusThetaH.multiply(theta);
- final T coeff2 = oneMinusThetaH.multiply(theta.add(1));
- for (int i = 0; i < previousState.length; ++i) {
- final T yDot1 = yDotK[0][i];
- final T yDot2 = yDotK[1][i];
- interpolatedState[i] = currentState[i].add(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2));
- interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2));
- }
+ final T coeff2 = oneMinusThetaH.multiply(theta.add(1)).negate();
+ interpolatedState = currentStateLinearCombination(coeff1, coeff2);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2);
}
- 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/RungeKuttaFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java
index 61750f6..a511ac8 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldIntegrator.java
@@ -112,18 +112,11 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
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);
/** {@inheritDoc} */
@@ -147,7 +140,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
// 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
@@ -172,6 +165,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
interpolator.shift();
// first stage
+ y = equations.getMapper().mapState(stepStart);
yDotK[0] = equations.getMapper().mapDerivative(stepStart);
// next stages
@@ -202,6 +196,7 @@ public abstract class RungeKuttaFieldIntegrator<T extends RealFieldElement<T>>
final FieldODEStateAndDerivative<T> stateTmp = new FieldODEStateAndDerivative<T>(stepEnd, yTmp, yDotTmp);
// discrete events handling
+ 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/RungeKuttaFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java
index fcc0396..9a86d82 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/RungeKuttaFieldStepInterpolator.java
@@ -36,31 +36,23 @@ import org.apache.commons.math3.util.MathArrays;
abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
extends AbstractFieldStepInterpolator<T> {
- /** Previous state. */
- protected T[] previousState;
-
- /** Slopes at the intermediate points */
- protected T[][] yDotK;
-
/** Reference to the integrator. */
protected AbstractFieldIntegrator<T> integrator;
+ /** Slopes at the intermediate points. */
+ private T[][] yDotK;
+
/** 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
*/
protected RungeKuttaFieldStepInterpolator(final AbstractFieldIntegrator<T> rkIntegrator,
- final T[] y, final T[][] yDotArray, final boolean forward,
+ final boolean forward,
final FieldEquationsMapper<T> mapper) {
- super(y, forward, mapper);
- this.previousState = null;
- this.yDotK = yDotArray;
- this.integrator = rkIntegrator;
+ super(forward, mapper);
+ this.yDotK = null;
+ this.integrator = rkIntegrator;
}
/** Copy constructor.
@@ -84,18 +76,14 @@ abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
super(interpolator);
- if (interpolator.currentState != null) {
-
- previousState = interpolator.previousState.clone();
-
+ if (yDotK != null) {
yDotK = MathArrays.buildArray(interpolator.integrator.getField(),
- interpolator.yDotK.length, interpolator.yDotK[0].length);
- for (int k = 0; k < interpolator.yDotK.length; ++k) {
- System.arraycopy(interpolator.yDotK[k], 0, yDotK[k], 0, interpolator.yDotK[k].length);
+ interpolator.yDotK.length, -1);
+ for (int k = 0; k < yDotK.length; ++k) {
+ yDotK[k] = interpolator.yDotK[k].clone();
}
} else {
- previousState = null;
yDotK = null;
}
@@ -105,11 +93,56 @@ abstract class RungeKuttaFieldStepInterpolator<T extends RealFieldElement<T>>
}
- /** {@inheritDoc} */
- @Override
- public void shift() {
- previousState = currentState.clone();
- super.shift();
+ /** Store the slopes at the intermediate points.
+ * @param slopes slopes at the intermediate points
+ */
+ void setSlopes(final T[][] slopes) {
+ this.yDotK = slopes.clone();
+ }
+
+ /** Compute a state by linear combination added to previous state.
+ * @param coefficients coefficients to apply to the method staged derivatives
+ * @return combined state
+ */
+ @SafeVarargs
+ protected final T[] previousStateLinearCombination(final T ... coefficients) {
+ return combine(getPreviousState().getState(),
+ coefficients);
+ }
+
+ /** Compute a state by linear combination added to current state.
+ * @param coefficients coefficients to apply to the method staged derivatives
+ * @return combined state
+ */
+ @SuppressWarnings("unchecked")
+ protected T[] currentStateLinearCombination(final T ... coefficients) {
+ return combine(getCurrentState().getState(),
+ coefficients);
+ }
+
+ /** Compute a state derivative by linear combination.
+ * @param coefficients coefficients to apply to the method staged derivatives
+ * @return combined state
+ */
+ @SuppressWarnings("unchecked")
+ protected T[] derivativeLinearCombination(final T ... coefficients) {
+ return combine(MathArrays.buildArray(integrator.getField(), yDotK[0].length),
+ coefficients);
+ }
+
+ /** Linearly combine arrays.
+ * @param a array to add to
+ * @param coefficients coefficients to apply to the method staged derivatives
+ * @return a itself, as a conveniency for fluent API
+ */
+ @SuppressWarnings("unchecked")
+ private T[] combine(final T[] a, final T ... coefficients) {
+ for (int i = 0; i < a.length; ++i) {
+ for (int k = 0; k < coefficients.length; ++k) {
+ a[i] = a[i].add(coefficients[k].multiply(yDotK[k][i]));
+ }
+ }
+ return a;
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.java
index 6ba4b72..e0d5044 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldIntegrator.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;
@@ -100,10 +99,8 @@ public class ThreeEighthesFieldIntegrator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
@Override
protected ThreeEighthesFieldStepInterpolator<T>
- createInterpolator(final AbstractFieldIntegrator<T> rkIntegrator, final T[] y,
- final T[][] yDotArray, final boolean forward,
- final FieldEquationsMapper<T> mapper) {
- return new ThreeEighthesFieldStepInterpolator<T>(rkIntegrator, y, yDotArray, forward, mapper);
+ createInterpolator(final boolean forward, final FieldEquationsMapper<T> mapper) {
+ return new ThreeEighthesFieldStepInterpolator<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/ThreeEighthesFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java
index 3de88f6..39cdc1c 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesFieldStepInterpolator.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 3/8 fourth
@@ -64,17 +63,13 @@ class ThreeEighthesFieldStepInterpolator<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
*/
ThreeEighthesFieldStepInterpolator(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.
@@ -94,6 +89,7 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>>
/** {@inheritDoc} */
+ @SuppressWarnings("unchecked")
@Override
protected FieldODEStateAndDerivative<T> computeInterpolatedStateAndDerivatives(final FieldEquationsMapper<T> mapper,
final T time, final T theta,
@@ -103,51 +99,31 @@ class ThreeEighthesFieldStepInterpolator<T extends RealFieldElement<T>>
final T coeffDot1 = coeffDot3.multiply(theta.multiply(4).subtract(5)).add(1);
final T coeffDot2 = coeffDot3.multiply(theta.multiply(-6).add(5));
final T coeffDot4 = coeffDot3.multiply(theta.multiply(2).subtract(1));
- 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)) {
+ if (getGlobalPreviousState() != null && theta.getReal() <= 0.5) {
final T s = theta.multiply(h).divide(8);
final T fourTheta2 = theta.multiply(theta).multiply(4);
final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(15)).add(8));
final T coeff2 = s.multiply(theta.multiply(5).subtract(fourTheta2)).multiply(3);
final T coeff3 = s.multiply(theta).multiply(3);
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 yDot2 = yDotK[1][i];
- final T yDot3 = yDotK[2][i];
- final T yDot4 = yDotK[3][i];
- interpolatedState[i] = previousState[i].
- add(coeff1.multiply(yDot1)).add(coeff2.multiply(yDot2)).
- add(coeff3.multiply(yDot3)).add(coeff4.multiply(yDot4));
- interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
- add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
-
- }
+ interpolatedState = previousStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, coeffDot4);
} else {
- final T s = oneMinusThetaH.divide(8);
+ final T s = oneMinusThetaH.divide(-8);
final T fourTheta2 = theta.multiply(theta).multiply(4);
final T thetaPlus1 = theta.add(1);
final T coeff1 = s.multiply(fourTheta2.multiply(2).subtract(theta.multiply(7)).add(1));
final T coeff2 = s.multiply(thetaPlus1.subtract(fourTheta2)).multiply(3);
final T coeff3 = s.multiply(thetaPlus1).multiply(3);
final T coeff4 = s.multiply(thetaPlus1.add(fourTheta2));
- for (int i = 0; i < interpolatedState.length; ++i) {
- final T yDot1 = yDotK[0][i];
- final T yDot2 = yDotK[1][i];
- final T yDot3 = yDotK[2][i];
- final T yDot4 = yDotK[3][i];
- interpolatedState[i] = currentState[i].
- subtract(coeff1.multiply(yDot1)).subtract(coeff2.multiply(yDot2)).
- subtract(coeff3.multiply(yDot3)).subtract(coeff4.multiply(yDot4));
- interpolatedDerivatives[i] = coeffDot1.multiply(yDot1).add(coeffDot2.multiply(yDot2)).
- add(coeffDot3.multiply(yDot3)).add(coeffDot4.multiply(yDot4));
-
- }
+ interpolatedState = currentStateLinearCombination(coeff1, coeff2, coeff3, coeff4);
+ interpolatedDerivatives = derivativeLinearCombination(coeffDot1, coeffDot2, coeffDot3, 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/sampling/AbstractFieldStepInterpolator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java
index b450a77..e9d6f29 100644
--- a/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java
+++ b/src/main/java/org/apache/commons/math3/ode/sampling/AbstractFieldStepInterpolator.java
@@ -43,9 +43,6 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T
/** Current time step. */
protected T h;
- /** Current state. */
- protected T[] currentState;
-
/** Global previous state. */
private FieldODEStateAndDerivative<T> globalPreviousState;
@@ -68,18 +65,16 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T
private FieldEquationsMapper<T> mapper;
/** Simple constructor.
- * @param y reference to the integrator array holding the state at the end of the step
* @param isForward integration direction indicator
* @param equationsMapper mapper for ODE equations primary and secondary components
*/
- protected AbstractFieldStepInterpolator(final T[] y, final boolean isForward,
+ protected AbstractFieldStepInterpolator(final boolean isForward,
final FieldEquationsMapper<T> equationsMapper) {
globalPreviousState = null;
globalCurrentState = null;
softPreviousState = null;
softCurrentState = null;
h = null;
- currentState = y;
finalized = false;
this.forward = isForward;
this.mapper = equationsMapper;
@@ -109,17 +104,9 @@ public abstract class AbstractFieldStepInterpolator<T extends RealFieldElement<T
softPreviousState = interpolator.softPreviousState;
softCurrentState = interpolator.softCurrentState;
h = interpolator.h;
-
- if (interpolator.currentState == null) {
- currentState = null;
- mapper = null;
- } else {
- currentState = interpolator.currentState.clone();
- }
-
- finalized = interpolator.finalized;
- forward = interpolator.forward;
- mapper = interpolator.mapper;
+ finalized = interpolator.finalized;
+ forward = interpolator.forward;
+ mapper = interpolator.mapper;
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java
index f79fad0..278276b 100644
--- a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java
+++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem1.java
@@ -66,9 +66,10 @@ public class TestFieldProblem1<T extends RealFieldElement<T>>
@Override
public T[] computeTheoreticalState(T t) {
- final T[] y0 = getInitialState();
+ final FieldODEState<T> s0 = getInitialState();
+ final T[] y0 = s0.getState();
final T[] y = MathArrays.buildArray(getField(), getDimension());
- T c = getInitialTime().subtract(t).exp();
+ T c = s0.getTime().subtract(t).exp();
for (int i = 0; i < getDimension(); ++i) {
y[i] = c.multiply(y0[i]);
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java
index f84781a..80f0ea2 100644
--- a/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java
+++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblem5.java
@@ -35,7 +35,7 @@ public class TestFieldProblem5<T extends RealFieldElement<T>>
*/
public TestFieldProblem5(Field<T> field) {
super(field);
- setFinalConditions(getInitialTime().multiply(2).subtract(getFinalTime()));
+ setFinalConditions(getInitialState().getTime().multiply(2).subtract(getFinalTime()));
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java
index 09a10d8..02ff65f 100644
--- a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java
+++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemAbstract.java
@@ -109,20 +109,12 @@ public abstract class TestFieldProblemAbstract<T extends RealFieldElement<T>>
return n;
}
- /**
- * Get the initial time.
- * @return initial time
- */
- public T getInitialTime() {
- return t0;
- }
-
- /**
- * Get the initial state vector.
- * @return initial state vector
+ /**
+ * Get the initial state.
+ * @return initial state
*/
- public T[] getInitialState() {
- return y0;
+ public FieldODEState<T> getInitialState() {
+ return new FieldODEState<T>(t0, y0);
}
/**
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java
index d4e20e6..101f4fa 100644
--- a/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java
+++ b/src/test/java/org/apache/commons/math3/ode/TestFieldProblemHandler.java
@@ -74,7 +74,7 @@ public class TestFieldProblemHandler<T extends RealFieldElement<T>>
public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast) throws MaxCountExceededException {
T start = integrator.getCurrentStepStart().getTime();
- if (start.subtract(problem.getInitialTime()).divide(integrator.getCurrentSignedStepsize()).abs().getReal() > 0.001) {
+ if (start.subtract(problem.getInitialState().getTime()).divide(integrator.getCurrentSignedStepsize()).abs().getReal() > 0.001) {
// multistep integrators do not handle the first steps themselves
// so we have to make sure the integrator we look at has really started its work
if (expectedStepStart != null) {
http://git-wip-us.apache.org/repos/asf/commons-math/blob/07e2b944/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java
index 54f0330..9ec73d2 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerIntegratorTest.java
@@ -58,8 +58,8 @@ public class EulerIntegratorTest {
MaxCountExceededException, NoBracketingException {
for (TestProblemAbstract pb : new TestProblemAbstract[] {
- new TestProblem1(), new TestProblem2(), new TestProblem3(),
- new TestProblem4(), new TestProblem5(), new TestProblem6()
+ //new TestProblem1(), new TestProblem2(), new TestProblem3(),
+ new TestProblem4()//, new TestProblem5(), new TestProblem6()
}) {
double previousValueError = Double.NaN;
@@ -94,6 +94,7 @@ public class EulerIntegratorTest {
}
previousTimeError = timeError;
+ System.exit(1);
}
}