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/24 23:06:09 UTC
[math] Fixed stability issues with Adams integrators.
Repository: commons-math
Updated Branches:
refs/heads/MATH_3_X bf317177d -> 6a01c7b2d
Fixed stability issues with Adams integrators.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/6a01c7b2
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/6a01c7b2
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/6a01c7b2
Branch: refs/heads/MATH_3_X
Commit: 6a01c7b2df5b8e61d491b7ff2ce8b0bdb9409cb2
Parents: bf31717
Author: Luc Maisonobe <lu...@apache.org>
Authored: Thu Dec 24 22:49:59 2015 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Thu Dec 24 23:05:21 2015 +0100
----------------------------------------------------------------------
src/changes/changes.xml | 6 +
.../commons/math3/ode/MultistepIntegrator.java | 30 ++-
.../ode/nonstiff/AdamsBashforthIntegrator.java | 123 +++++++----
.../ode/nonstiff/AdamsMoultonIntegrator.java | 14 +-
.../ode/nonstiff/AdamsNordsieckTransformer.java | 100 +++++----
.../commons/math3/ode/TestProblemHandler.java | 40 ++--
.../nonstiff/AdamsBashforthIntegratorTest.java | 203 +++++++++++++++++--
.../nonstiff/AdamsMoultonIntegratorTest.java | 191 ++++++++++++++++-
.../nonstiff/AdamsNordsieckTransformerTest.java | 120 +++++++++++
...ClassicalRungeKuttaStepInterpolatorTest.java | 2 +-
.../DormandPrince54StepInterpolatorTest.java | 2 +-
.../DormandPrince853StepInterpolatorTest.java | 2 +-
.../ode/nonstiff/EulerStepInterpolatorTest.java | 2 +-
.../ode/nonstiff/GillStepInterpolatorTest.java | 2 +-
.../GraggBulirschStoerStepInterpolatorTest.java | 2 +-
.../HighamHall54StepInterpolatorTest.java | 2 +-
.../nonstiff/LutherStepInterpolatorTest.java | 2 +-
.../nonstiff/MidpointStepInterpolatorTest.java | 2 +-
.../ThreeEighthesStepInterpolatorTest.java | 2 +-
.../sampling/NordsieckStepInterpolatorTest.java | 8 +-
.../ode/sampling/StepInterpolatorTestUtils.java | 8 +-
21 files changed, 708 insertions(+), 155 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index e5f9246..fcf3bad 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -51,6 +51,12 @@ If the output is not quite correct, check for invisible trailing spaces!
</properties>
<body>
<release version="3.6" date="XXXX-XX-XX" description="">
+ <action dev="luc" type="fix">
+ Fixed stability issues with Adams-Bashforth and Adams-Moulton ODE integrators.
+ The integrators did not estimate the local error properly and were sometimes
+ stuck attempting to reduce indefinitely the step size as they thought the
+ error was too high.
+ </action>
<action dev="luc" type="add" issue="MATH-1288">
Added a field-based version of Ordinary Differential Equations framework.
This allows integrating ode that refer to RealField elements instead of
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java b/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java
index 17baed7..f82ccbb 100644
--- a/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java
@@ -223,7 +223,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
starter.clearStepHandlers();
// set up one specific step handler to extract initial Nordsieck vector
- starter.addStepHandler(new NordsieckInitializer(nSteps, y0.length));
+ starter.addStepHandler(new NordsieckInitializer((nSteps + 3) / 2, y0.length));
// start integration, expecting a InitializationCompletedMarkerException
try {
@@ -313,6 +313,13 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
this.safety = safety;
}
+ /** Get the number of steps of the multistep method (excluding the one being computed).
+ * @return number of steps of the multistep method (excluding the one being computed)
+ */
+ public int getNSteps() {
+ return nSteps;
+ }
+
/** Compute step grow/shrink factor according to normalized error.
* @param error normalized error of the current step
* @return grow/shrink factor for next step
@@ -321,7 +328,10 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
return FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
}
- /** Transformer used to convert the first step to Nordsieck representation. */
+ /** Transformer used to convert the first step to Nordsieck representation.
+ * @deprecated as of 3.6 this unused interface is deprecated
+ */
+ @Deprecated
public interface NordsieckTransformer {
/** Initialize the high order scaled derivatives at step start.
* @param h step size to use for scaling
@@ -352,14 +362,14 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
private final double[][] yDot;
/** Simple constructor.
- * @param nSteps number of steps of the multistep method (excluding the one being computed)
+ * @param nbStartPoints number of start points (including the initial point)
* @param n problem dimension
*/
- NordsieckInitializer(final int nSteps, final int n) {
+ NordsieckInitializer(final int nbStartPoints, final int n) {
this.count = 0;
- this.t = new double[nSteps];
- this.y = new double[nSteps][n];
- this.yDot = new double[nSteps][n];
+ this.t = new double[nbStartPoints];
+ this.y = new double[nbStartPoints][n];
+ this.yDot = new double[nbStartPoints][n];
}
/** {@inheritDoc} */
@@ -370,7 +380,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
final double curr = interpolator.getCurrentTime();
if (count == 0) {
- // first step, we need to store also the beginning of the step
+ // first step, we need to store also the point at the beginning of the step
interpolator.setInterpolatedTime(prev);
t[0] = prev;
final ExpandableStatefulODE expandable = getExpandable();
@@ -385,7 +395,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
}
}
- // store the end of the step
+ // store the point at the end of the step
++count;
interpolator.setInterpolatedTime(curr);
t[count] = curr;
@@ -403,7 +413,7 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
if (count == t.length - 1) {
- // this was the last step we needed, we can compute the derivatives
+ // this was the last point we needed, we can compute the derivatives
stepStart = t[0];
stepSize = (t[t.length - 1] - t[0]) / (t.length - 1);
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java
index fd80d85..af4a435 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java
@@ -22,6 +22,7 @@ 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.linear.Array2DRowRealMatrix;
+import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.ode.EquationsMapper;
import org.apache.commons.math3.ode.ExpandableStatefulODE;
import org.apache.commons.math3.ode.sampling.NordsieckStepInterpolator;
@@ -90,7 +91,7 @@ import org.apache.commons.math3.util.FastMath;
* computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
* for degree k polynomials.
* <pre>
- * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + ∑<sub>j</sub> j (-i)<sup>j-1</sup> s<sub>j</sub>(n)
+ * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + ∑<sub>j>0</sub> (j+1) (-i)<sup>j</sup> s<sub>j+1</sub>(n)
* </pre>
* The previous formula can be used with several values for i to compute the transform between
* classical representation and Nordsieck vector. The transform between r<sub>n</sub>
@@ -99,7 +100,8 @@ import org.apache.commons.math3.util.FastMath;
* q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
* </pre>
* where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built
- * with the j (-i)<sup>j-1</sup> terms:
+ * with the (j+1) (-i)<sup>j</sup> terms with i being the row number starting from 1 and j being
+ * the column number starting from 1:
* <pre>
* [ -2 3 -4 5 ... ]
* [ -4 12 -32 80 ... ]
@@ -188,6 +190,48 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
vecAbsoluteTolerance, vecRelativeTolerance);
}
+ /** Estimate error.
+ * <p>
+ * Error is estimated by interpolating back to previous state using
+ * the state Taylor expansion and comparing to real previous state.
+ * </p>
+ * @param previousState state vector at step start
+ * @param predictedState predicted state vector at step end
+ * @param predictedScaled predicted value of the scaled derivatives at step end
+ * @param predictedNordsieck predicted value of the Nordsieck vector at step end
+ * @return estimated normalized local discretization error
+ */
+ private double errorEstimation(final double[] previousState,
+ final double[] predictedState,
+ final double[] predictedScaled,
+ final RealMatrix predictedNordsieck) {
+
+ double error = 0;
+ for (int i = 0; i < mainSetDimension; ++i) {
+ final double yScale = FastMath.abs(predictedState[i]);
+ final double tol = (vecAbsoluteTolerance == null) ?
+ (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
+ (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
+
+ // apply Taylor formula from high order to low order,
+ // for the sake of numerical accuracy
+ double variation = 0;
+ int sign = predictedNordsieck.getRowDimension() % 2 == 0 ? -1 : 1;
+ for (int k = predictedNordsieck.getRowDimension() - 1; k >= 0; --k) {
+ variation += sign * predictedNordsieck.getEntry(k, i);
+ sign = -sign;
+ }
+ variation -= predictedScaled[i];
+
+ final double ratio = (predictedState[i] - previousState[i] + variation) / tol;
+ error += ratio * ratio;
+
+ }
+
+ return FastMath.sqrt(error / mainSetDimension);
+
+ }
+
/** {@inheritDoc} */
@Override
public void integrate(final ExpandableStatefulODE equations, final double t)
@@ -199,8 +243,7 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
final boolean forward = t > equations.getTime();
// initialize working arrays
- final double[] y0 = equations.getCompleteState();
- final double[] y = y0.clone();
+ final double[] y = equations.getCompleteState();
final double[] yDot = new double[y.length];
// set up an interpolator sharing the integrator arrays
@@ -209,13 +252,12 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
equations.getPrimaryMapper(), equations.getSecondaryMappers());
// set up integration control objects
- initIntegration(equations.getTime(), y0, t);
+ initIntegration(equations.getTime(), y, t);
// compute the initial Nordsieck vector using the configured starter integrator
start(equations.getTime(), y, t);
interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
interpolator.storeTime(stepStart);
- final int lastRow = nordsieck.getRowDimension() - 1;
// reuse the step that was chosen by the starter integrator
double hNew = stepSize;
@@ -225,62 +267,57 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
isLastStep = false;
do {
+ interpolator.shift();
+ final double[] predictedY = new double[y.length];
+ final double[] predictedScaled = new double[y.length];
+ Array2DRowRealMatrix predictedNordsieck = null;
double error = 10;
while (error >= 1.0) {
- stepSize = hNew;
-
- // evaluate error using the last term of the Taylor expansion
- error = 0;
- for (int i = 0; i < mainSetDimension; ++i) {
- final double yScale = FastMath.abs(y[i]);
- final double tol = (vecAbsoluteTolerance == null) ?
- (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
- (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
- final double ratio = nordsieck.getEntry(lastRow, i) / tol;
- error += ratio * ratio;
+ // predict a first estimate of the state at step end
+ final double stepEnd = stepStart + hNew;
+ interpolator.storeTime(stepEnd);
+ final ExpandableStatefulODE expandable = getExpandable();
+ final EquationsMapper primary = expandable.getPrimaryMapper();
+ primary.insertEquationData(interpolator.getInterpolatedState(), predictedY);
+ int index = 0;
+ for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
+ secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), predictedY);
+ ++index;
}
- error = FastMath.sqrt(error / mainSetDimension);
+
+ // evaluate the derivative
+ computeDerivatives(stepEnd, predictedY, yDot);
+
+ // predict Nordsieck vector at step end
+ for (int j = 0; j < predictedScaled.length; ++j) {
+ predictedScaled[j] = hNew * yDot[j];
+ }
+ predictedNordsieck = updateHighOrderDerivativesPhase1(nordsieck);
+ updateHighOrderDerivativesPhase2(scaled, predictedScaled, predictedNordsieck);
+
+ // evaluate error
+ error = errorEstimation(y, predictedY, predictedScaled, predictedNordsieck);
if (error >= 1.0) {
// reject the step and attempt to reduce error by stepsize control
final double factor = computeStepGrowShrinkFactor(error);
- hNew = filterStep(stepSize * factor, forward, false);
+ hNew = filterStep(hNew * factor, forward, false);
interpolator.rescale(hNew);
}
}
- // predict a first estimate of the state at step end
+ stepSize = hNew;
final double stepEnd = stepStart + stepSize;
- interpolator.shift();
- interpolator.setInterpolatedTime(stepEnd);
- final ExpandableStatefulODE expandable = getExpandable();
- final EquationsMapper primary = expandable.getPrimaryMapper();
- primary.insertEquationData(interpolator.getInterpolatedState(), y);
- int index = 0;
- for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
- secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y);
- ++index;
- }
-
- // evaluate the derivative
- computeDerivatives(stepEnd, y, yDot);
-
- // update Nordsieck vector
- final double[] predictedScaled = new double[y0.length];
- for (int j = 0; j < y0.length; ++j) {
- predictedScaled[j] = stepSize * yDot[j];
- }
- final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
- updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
- interpolator.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
+ interpolator.reinitialize(stepEnd, stepSize, predictedScaled, predictedNordsieck);
// discrete events handling
interpolator.storeTime(stepEnd);
+ System.arraycopy(predictedY, 0, y, 0, y.length);
stepStart = acceptStep(interpolator, y, yDot, t);
scaled = predictedScaled;
- nordsieck = nordsieckTmp;
+ nordsieck = predictedNordsieck;
interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
if (!isLastStep) {
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java
index c12c19b..69d2469 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java
@@ -96,7 +96,7 @@ import org.apache.commons.math3.util.FastMath;
* computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
* for degree k polynomials.
* <pre>
- * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + ∑<sub>j</sub> j (-i)<sup>j-1</sup> s<sub>j</sub>(n)
+ * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + ∑<sub>j>0</sub> (j+1) (-i)<sup>j</sup> s<sub>j+1</sub>(n)
* </pre>
* The previous formula can be used with several values for i to compute the transform between
* classical representation and Nordsieck vector. The transform between r<sub>n</sub>
@@ -105,7 +105,8 @@ import org.apache.commons.math3.util.FastMath;
* q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
* </pre>
* where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built
- * with the j (-i)<sup>j-1</sup> terms:
+ * with the (j+1) (-i)<sup>j</sup> terms with i being the row number starting from 1 and j being
+ * the column number starting from 1:
* <pre>
* [ -2 3 -4 5 ... ]
* [ -4 12 -32 80 ... ]
@@ -204,7 +205,6 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
vecAbsoluteTolerance, vecRelativeTolerance);
}
-
/** {@inheritDoc} */
@Override
public void integrate(final ExpandableStatefulODE equations,final double t)
@@ -405,10 +405,10 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
after[i] += previous[i] + scaled[i];
if (i < mainSetDimension) {
final double yScale = FastMath.max(FastMath.abs(previous[i]), FastMath.abs(after[i]));
- final double tol = (vecAbsoluteTolerance == null) ?
- (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
- (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
- final double ratio = (after[i] - before[i]) / tol;
+ final double tol = (vecAbsoluteTolerance == null) ?
+ (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
+ (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
+ final double ratio = (after[i] - before[i]) / tol; // (corrected-predicted)/tol
error += ratio * ratio;
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformer.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformer.java b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformer.java
index 68f66c6..1c54b17 100644
--- a/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformer.java
+++ b/src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformer.java
@@ -68,7 +68,7 @@ import org.apache.commons.math3.linear.RealMatrix;
* computed from s<sub>1</sub>(n), s<sub>2</sub>(n) ... s<sub>k</sub>(n), the formula being exact
* for degree k polynomials.
* <pre>
- * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + ∑<sub>j>1</sub> j (-i)<sup>j-1</sup> s<sub>j</sub>(n)
+ * s<sub>1</sub>(n-i) = s<sub>1</sub>(n) + ∑<sub>j>0</sub> (j+1) (-i)<sup>j</sup> s<sub>j+1</sub>(n)
* </pre>
* The previous formula can be used with several values for i to compute the transform between
* classical representation and Nordsieck vector at step end. The transform between r<sub>n</sub>
@@ -77,7 +77,8 @@ import org.apache.commons.math3.linear.RealMatrix;
* q<sub>n</sub> = s<sub>1</sub>(n) u + P r<sub>n</sub>
* </pre>
* where u is the [ 1 1 ... 1 ]<sup>T</sup> vector and P is the (k-1)×(k-1) matrix built
- * with the j (-i)<sup>j-1</sup> terms:
+ * with the (j+1) (-i)<sup>j</sup> terms with i being the row number starting from 1 and j being
+ * the column number starting from 1:
* <pre>
* [ -2 3 -4 5 ... ]
* [ -4 12 -32 80 ... ]
@@ -137,27 +138,28 @@ public class AdamsNordsieckTransformer {
private static final Map<Integer, AdamsNordsieckTransformer> CACHE =
new HashMap<Integer, AdamsNordsieckTransformer>();
- /** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */
+ /** Update matrix for the higher order derivatives h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ... */
private final Array2DRowRealMatrix update;
/** Update coefficients of the higher order derivatives wrt y'. */
private final double[] c1;
/** Simple constructor.
- * @param nSteps number of steps of the multistep method
+ * @param n number of steps of the multistep method
* (excluding the one being computed)
*/
- private AdamsNordsieckTransformer(final int nSteps) {
+ private AdamsNordsieckTransformer(final int n) {
+
+ final int rows = n - 1;
// compute exact coefficients
- FieldMatrix<BigFraction> bigP = buildP(nSteps);
+ FieldMatrix<BigFraction> bigP = buildP(rows);
FieldDecompositionSolver<BigFraction> pSolver =
new FieldLUDecomposition<BigFraction>(bigP).getSolver();
- BigFraction[] u = new BigFraction[nSteps];
+ BigFraction[] u = new BigFraction[rows];
Arrays.fill(u, BigFraction.ONE);
- BigFraction[] bigC1 = pSolver
- .solve(new ArrayFieldVector<BigFraction>(u, false)).toArray();
+ BigFraction[] bigC1 = pSolver.solve(new ArrayFieldVector<BigFraction>(u, false)).toArray();
// update coefficients are computed by combining transform from
// Nordsieck to multistep, then shifting rows to represent step advance
@@ -167,15 +169,15 @@ public class AdamsNordsieckTransformer {
// shift rows
shiftedP[i] = shiftedP[i - 1];
}
- shiftedP[0] = new BigFraction[nSteps];
+ shiftedP[0] = new BigFraction[rows];
Arrays.fill(shiftedP[0], BigFraction.ZERO);
FieldMatrix<BigFraction> bigMSupdate =
pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
// convert coefficients to double
update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
- c1 = new double[nSteps];
- for (int i = 0; i < nSteps; ++i) {
+ c1 = new double[rows];
+ for (int i = 0; i < rows; ++i) {
c1[i] = bigC1[i].doubleValue();
}
@@ -201,13 +203,17 @@ public class AdamsNordsieckTransformer {
* (excluding the one being computed).
* @return number of steps of the method
* (excluding the one being computed)
+ * @deprecated as of 3.6, this method is not used anymore
*/
+ @Deprecated
public int getNSteps() {
return c1.length;
}
/** Build the P matrix.
- * <p>The P matrix general terms are shifted j (-i)<sup>j-1</sup> terms:
+ * <p>The P matrix general terms are shifted (j+1) (-i)<sup>j</sup> terms
+ * with i being the row number starting from 1 and j being the column
+ * number starting from 1:
* <pre>
* [ -2 3 -4 5 ... ]
* [ -4 12 -32 80 ... ]
@@ -215,21 +221,20 @@ public class AdamsNordsieckTransformer {
* [ -8 48 -256 1280 ... ]
* [ ... ]
* </pre></p>
- * @param nSteps number of steps of the multistep method
- * (excluding the one being computed)
+ * @param rows number of rows of the matrix
* @return P matrix
*/
- private FieldMatrix<BigFraction> buildP(final int nSteps) {
+ private FieldMatrix<BigFraction> buildP(final int rows) {
- final BigFraction[][] pData = new BigFraction[nSteps][nSteps];
+ final BigFraction[][] pData = new BigFraction[rows][rows];
- for (int i = 0; i < pData.length; ++i) {
+ for (int i = 1; i <= pData.length; ++i) {
// build the P matrix elements from Taylor series formulas
- final BigFraction[] pI = pData[i];
- final int factor = -(i + 1);
+ final BigFraction[] pI = pData[i - 1];
+ final int factor = -i;
int aj = factor;
- for (int j = 0; j < pI.length; ++j) {
- pI[j] = new BigFraction(aj * (j + 2));
+ for (int j = 1; j <= pI.length; ++j) {
+ pI[j - 1] = new BigFraction(aj * (j + 1));
aj *= factor;
}
}
@@ -243,20 +248,25 @@ public class AdamsNordsieckTransformer {
* @param t first steps times
* @param y first steps states
* @param yDot first steps derivatives
- * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
+ * @return Nordieck vector at start of first step (h<sup>2</sup>/2 y''<sub>n</sub>,
* h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
*/
+
public Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
final double[][] y,
final double[][] yDot) {
// using Taylor series with di = ti - t0, we get:
- // y(ti) - y(t0) - di y'(t0) = di^2 / h^2 s2 + ... + di^k / h^k sk + O(h^(k+1))
- // y'(ti) - y'(t0) = 2 di / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^k)
- // we write these relations for i = 1 to i= n-1 as a set of 2(n-1) linear
- // equations depending on the Nordsieck vector [s2 ... sk]
- final double[][] a = new double[2 * (y.length - 1)][c1.length];
- final double[][] b = new double[2 * (y.length - 1)][y[0].length];
+ // y(ti) - y(t0) - di y'(t0) = di^2 / h^2 s2 + ... + di^k / h^k sk + O(h^k)
+ // y'(ti) - y'(t0) = 2 di / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^(k-1))
+ // we write these relations for i = 1 to i= 1+n/2 as a set of n + 2 linear
+ // equations depending on the Nordsieck vector [s2 ... sk rk], so s2 to sk correspond
+ // to the appropriately truncated Taylor expansion, and rk is the Taylor remainder.
+ // The goal is to have s2 to sk as accurate as possible considering the fact the sum is
+ // truncated and we don't want the error terms to be included in s2 ... sk, so we need
+ // to solve also for the remainder
+ final double[][] a = new double[c1.length + 1][c1.length + 1];
+ final double[][] b = new double[c1.length + 1][y[0].length];
final double[] y0 = y[0];
final double[] yDot0 = yDot[0];
for (int i = 1; i < y.length; ++i) {
@@ -268,31 +278,43 @@ public class AdamsNordsieckTransformer {
// linear coefficients of equations
// y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0)
final double[] aI = a[2 * i - 2];
- final double[] aDotI = a[2 * i - 1];
+ final double[] aDotI = (2 * i - 1) < a.length ? a[2 * i - 1] : null;
for (int j = 0; j < aI.length; ++j) {
dikM1Ohk *= ratio;
aI[j] = di * dikM1Ohk;
- aDotI[j] = (j + 2) * dikM1Ohk;
+ if (aDotI != null) {
+ aDotI[j] = (j + 2) * dikM1Ohk;
+ }
}
// expected value of the previous equations
final double[] yI = y[i];
final double[] yDotI = yDot[i];
final double[] bI = b[2 * i - 2];
- final double[] bDotI = b[2 * i - 1];
+ final double[] bDotI = (2 * i - 1) < b.length ? b[2 * i - 1] : null;
for (int j = 0; j < yI.length; ++j) {
bI[j] = yI[j] - y0[j] - di * yDot0[j];
- bDotI[j] = yDotI[j] - yDot0[j];
+ if (bDotI != null) {
+ bDotI[j] = yDotI[j] - yDot0[j];
+ }
}
}
- // solve the rectangular system in the least square sense
- // to get the best estimate of the Nordsieck vector [s2 ... sk]
- QRDecomposition decomposition;
- decomposition = new QRDecomposition(new Array2DRowRealMatrix(a, false));
- RealMatrix x = decomposition.getSolver().solve(new Array2DRowRealMatrix(b, false));
- return new Array2DRowRealMatrix(x.getData(), false);
+ // solve the linear system to get the best estimate of the Nordsieck vector [s2 ... sk],
+ // with the additional terms s(k+1) and c grabbing the parts after the truncated Taylor expansion
+ final QRDecomposition decomposition = new QRDecomposition(new Array2DRowRealMatrix(a, false));
+ final RealMatrix x = decomposition.getSolver().solve(new Array2DRowRealMatrix(b, false));
+
+ // extract just the Nordsieck vector [s2 ... sk]
+ final Array2DRowRealMatrix truncatedX = new Array2DRowRealMatrix(x.getRowDimension() - 1, x.getColumnDimension());
+ for (int i = 0; i < truncatedX.getRowDimension(); ++i) {
+ for (int j = 0; j < truncatedX.getColumnDimension(); ++j) {
+ truncatedX.setEntry(i, j, x.getEntry(i, j));
+ }
+ }
+ return truncatedX;
+
}
/** Update the high order scaled derivatives for Adams integrators (phase 1).
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/TestProblemHandler.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/TestProblemHandler.java b/src/test/java/org/apache/commons/math3/ode/TestProblemHandler.java
index 40ff404..05d83b6 100644
--- a/src/test/java/org/apache/commons/math3/ode/TestProblemHandler.java
+++ b/src/test/java/org/apache/commons/math3/ode/TestProblemHandler.java
@@ -87,35 +87,36 @@ public class TestProblemHandler
expectedStepStart = start + integrator.getCurrentSignedStepsize();
}
+
double pT = interpolator.getPreviousTime();
double cT = interpolator.getCurrentTime();
double[] errorScale = problem.getErrorScale();
// store the error at the last step
if (isLast) {
- double[] interpolatedY = interpolator.getInterpolatedState();
- double[] theoreticalY = problem.computeTheoreticalState(cT);
- for (int i = 0; i < interpolatedY.length; ++i) {
- double error = FastMath.abs(interpolatedY[i] - theoreticalY[i]);
- lastError = FastMath.max(error, lastError);
- }
- lastTime = cT;
+ double[] interpolatedY = interpolator.getInterpolatedState();
+ double[] theoreticalY = problem.computeTheoreticalState(cT);
+ for (int i = 0; i < interpolatedY.length; ++i) {
+ double error = FastMath.abs(interpolatedY[i] - theoreticalY[i]);
+ lastError = FastMath.max(error, lastError);
+ }
+ lastTime = cT;
}
-
// walk through the step
for (int k = 0; k <= 20; ++k) {
- double time = pT + (k * (cT - pT)) / 20;
- interpolator.setInterpolatedTime(time);
- double[] interpolatedY = interpolator.getInterpolatedState();
- double[] theoreticalY = problem.computeTheoreticalState(interpolator.getInterpolatedTime());
+ double time = pT + (k * (cT - pT)) / 20;
+ interpolator.setInterpolatedTime(time);
+ double[] interpolatedY = interpolator.getInterpolatedState();
+ double[] theoreticalY = problem.computeTheoreticalState(interpolator.getInterpolatedTime());
- // update the errors
- for (int i = 0; i < interpolatedY.length; ++i) {
- double error = errorScale[i] * FastMath.abs(interpolatedY[i] - theoreticalY[i]);
- maxValueError = FastMath.max(error, maxValueError);
- }
+ // update the errors
+ for (int i = 0; i < interpolatedY.length; ++i) {
+ double error = errorScale[i] * FastMath.abs(interpolatedY[i] - theoreticalY[i]);
+ maxValueError = FastMath.max(error, maxValueError);
+ }
}
+
}
/**
@@ -134,6 +135,11 @@ public class TestProblemHandler
return maxTimeError;
}
+
+ public int getCalls() {
+ return problem.getCalls();
+ }
+
/**
* Get the error at the end of the integration.
* @return error at the end of the integration
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegratorTest.java
index f5a96a6..ba7409a 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegratorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegratorTest.java
@@ -18,15 +18,29 @@
package org.apache.commons.math3.ode.nonstiff;
+import java.io.File;
+import java.io.IOException;
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+import java.io.PrintStream;
+import java.util.Locale;
+
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
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.AbstractIntegrator;
+import org.apache.commons.math3.ode.ExpandableStatefulODE;
import org.apache.commons.math3.ode.FirstOrderIntegrator;
import org.apache.commons.math3.ode.TestProblem1;
import org.apache.commons.math3.ode.TestProblem5;
import org.apache.commons.math3.ode.TestProblem6;
+import org.apache.commons.math3.ode.TestProblemAbstract;
import org.apache.commons.math3.ode.TestProblemHandler;
+import org.apache.commons.math3.ode.sampling.StepHandler;
+import org.apache.commons.math3.ode.sampling.StepInterpolator;
+import org.apache.commons.math3.util.CombinatoricsUtils;
import org.apache.commons.math3.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
@@ -64,8 +78,7 @@ public class AdamsBashforthIntegratorTest {
}
@Test
- public void testIncreasingTolerance() throws DimensionMismatchException, NumberIsTooSmallException, MaxCountExceededException, NoBracketingException
- {
+ public void testIncreasingTolerance() throws DimensionMismatchException, NumberIsTooSmallException, MaxCountExceededException, NoBracketingException {
int previousCalls = Integer.MAX_VALUE;
for (int i = -12; i < -5; ++i) {
@@ -84,12 +97,11 @@ public class AdamsBashforthIntegratorTest {
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- // the 50 and 300 factors are only valid for this test
+ // the 8 and 122 factors are only valid for this test
// and has been obtained from trial and error
- // there is no general relation between local and global errors
- Assert.assertTrue(handler.getMaximalValueError() > (50.0 * scalAbsoluteTolerance));
- Assert.assertTrue(handler.getMaximalValueError() < (300.0 * scalAbsoluteTolerance));
- Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
+ // there are no general relationship between local and global errors
+ Assert.assertTrue(handler.getMaximalValueError() > ( 8 * scalAbsoluteTolerance));
+ Assert.assertTrue(handler.getMaximalValueError() < (122 * scalAbsoluteTolerance));
int calls = pb.getCalls();
Assert.assertEquals(integ.getEvaluations(), calls);
@@ -122,14 +134,15 @@ public class AdamsBashforthIntegratorTest {
TestProblem5 pb = new TestProblem5();
double range = FastMath.abs(pb.getFinalTime() - pb.getInitialTime());
- FirstOrderIntegrator integ = new AdamsBashforthIntegrator(4, 0, range, 1.0e-12, 1.0e-12);
+ AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(4, 0, range, 1.0e-12, 1.0e-12);
+ integ.setStarterIntegrator(new PerfectStarter(pb, (integ.getNSteps() + 5) / 2));
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- Assert.assertTrue(handler.getLastError() < 1.5e-8);
- Assert.assertTrue(handler.getMaximalValueError() < 1.5e-8);
+ Assert.assertEquals(0.0, handler.getLastError(), 4.3e-8);
+ Assert.assertEquals(0.0, handler.getMaximalValueError(), 4.3e-8);
Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
Assert.assertEquals("Adams-Bashforth", integ.getName());
}
@@ -141,16 +154,178 @@ public class AdamsBashforthIntegratorTest {
for (int nSteps = 2; nSteps < 8; ++nSteps) {
AdamsBashforthIntegrator integ =
- new AdamsBashforthIntegrator(nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-5, 1.0e-5);
+ new AdamsBashforthIntegrator(nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-4, 1.0e-4);
+ integ.setStarterIntegrator(new PerfectStarter(pb, nSteps));
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- if (nSteps < 4) {
- Assert.assertTrue(handler.getMaximalValueError() > 1.0e-03);
+ if (nSteps < 5) {
+ Assert.assertTrue(handler.getMaximalValueError() > 0.005);
} else {
- Assert.assertTrue(handler.getMaximalValueError() < 4.0e-12);
+ Assert.assertTrue(handler.getMaximalValueError() < 5.0e-10);
+ }
+ }
+
+ }
+
+ private static class PerfectStarter extends AbstractIntegrator {
+
+ private final PerfectInterpolator interpolator;
+ private final int nbSteps;
+
+ public PerfectStarter(final TestProblemAbstract problem, final int nbSteps) {
+ this.interpolator = new PerfectInterpolator(problem);
+ this.nbSteps = nbSteps;
+ }
+
+ public void integrate(ExpandableStatefulODE equations, double t) {
+ double tStart = equations.getTime() + 0.01 * (t - equations.getTime());
+ for (int i = 0; i < nbSteps; ++i) {
+ double tK = ((nbSteps - 1 - (i + 1)) * equations.getTime() + (i + 1) * tStart) / (nbSteps - 1);
+ interpolator.setPreviousTime(interpolator.getCurrentTime());
+ interpolator.setCurrentTime(tK);
+ interpolator.setInterpolatedTime(tK);
+ for (StepHandler handler : getStepHandlers()) {
+ handler.handleStep(interpolator, i == nbSteps - 1);
+ }
+ }
+ }
+
+ }
+
+ private static class PerfectInterpolator implements StepInterpolator {
+ private final TestProblemAbstract problem;
+ private double previousTime;
+ private double currentTime;
+ private double interpolatedTime;
+
+ public PerfectInterpolator(final TestProblemAbstract problem) {
+ this.problem = problem;
+ this.previousTime = problem.getInitialTime();
+ this.currentTime = problem.getInitialTime();
+ this.interpolatedTime = problem.getInitialTime();
+ }
+
+ public void readExternal(ObjectInput arg0) {
+ }
+
+ public void writeExternal(ObjectOutput arg0) {
+ }
+
+ public double getPreviousTime() {
+ return previousTime;
+ }
+
+ public void setPreviousTime(double time) {
+ previousTime = time;
+ }
+
+ public double getCurrentTime() {
+ return currentTime;
+ }
+
+ public void setCurrentTime(double time) {
+ currentTime = time;
+ }
+
+ public double getInterpolatedTime() {
+ return interpolatedTime;
+ }
+
+ public void setInterpolatedTime(double time) {
+ interpolatedTime = time;
+ }
+
+ public double[] getInterpolatedState() {
+ return problem.computeTheoreticalState(interpolatedTime);
+ }
+
+ public double[] getInterpolatedDerivatives() {
+ double[] y = problem.computeTheoreticalState(interpolatedTime);
+ double[] yDot = new double[y.length];
+ problem.computeDerivatives(interpolatedTime, y, yDot);
+ return yDot;
+ }
+
+ public double[] getInterpolatedSecondaryState(int index) {
+ return null;
+ }
+
+ public double[] getInterpolatedSecondaryDerivatives(int index) {
+ return null;
+ }
+
+ public boolean isForward() {
+ return problem.getFinalTime() > problem.getInitialTime();
+ }
+
+ public StepInterpolator copy() {
+ return this;
+ }
+
+ }
+
+ @Test
+ public void testTmp() throws IOException, DimensionMismatchException, NumberIsTooSmallException, MaxCountExceededException, NoBracketingException {
+ final PrintStream out = new PrintStream(new File(new File(System.getProperty("user.home")), "x.dat"));
+ final int n = 4;
+ final AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(n, 1.0e-12, 1.0e3, 1.0e-10, 1.0e-12);
+ final SinT4 sinT4 = new SinT4();
+ integ.setStarterIntegrator(new PerfectStarter(sinT4, n));
+ final StepHandler handler = new StepHandler() {
+ double previousError = 0;
+ public void init(double t0, double[] y0, double t) {
}
+ public void handleStep(StepInterpolator interpolator, boolean isLast) {
+ final double t = interpolator.getCurrentTime();
+ final double h = t - interpolator.getPreviousTime();
+ final double y = interpolator.getInterpolatedState()[0];
+ final double error = sinT4.computeTheoreticalState(t)[0] - y;
+ out.format(Locale.US, "%10.8f %16.9e %10.8f %16.9e",
+ t, h, y, error - previousError);
+ for (int i = 1; i < n + 1; ++i) {
+ out.format(Locale.US, " %16.9e", sinT4.scaledDerivative(t, h, i));
+ }
+ out.format(Locale.US, "%n");
+ previousError = error;
+ }
+ };
+ integ.addStepHandler(handler);
+ final ExpandableStatefulODE expandableODE = new ExpandableStatefulODE(sinT4);
+ final double t0 = 0.75;
+ final double h0 = 1.0 / 128;
+ expandableODE.setTime(t0);
+ expandableODE.setPrimaryState(sinT4.computeTheoreticalState(t0));
+ integ.setInitialStepSize(h0);
+ integ.integrate(expandableODE, 0.9);
+ out.close();
+ }
+
+ private static class SinT4 extends TestProblemAbstract {
+
+ @Override
+ public int getDimension() {
+ return 1;
+ }
+
+ @Override
+ public void doComputeDerivatives(double t, double[] y, double[] yDot) {
+ // compute the derivatives
+ double t3 = FastMath.pow(t, 3);
+ yDot[0] = 4 * t3 * FastMath.cos(t3 * t);
+ }
+
+ public double[] computeTheoreticalState(final double t) {
+ return new double[] { FastMath.sin(FastMath.pow(t, 4)) };
+ }
+
+ public DerivativeStructure valueDS(final double t, final int n) {
+ return new DerivativeStructure(1, n, 0, t).pow(4).sin();
+ }
+
+ public double scaledDerivative(final double t, final double h, final int n) {
+ return FastMath.pow(h, n) * valueDS(t, n).getPartialDerivative(n) / CombinatoricsUtils.factorial(n);
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegratorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegratorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegratorTest.java
index 712eba4..ebab362 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegratorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegratorTest.java
@@ -18,15 +18,29 @@
package org.apache.commons.math3.ode.nonstiff;
+import java.io.File;
+import java.io.IOException;
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+import java.io.PrintStream;
+import java.util.Locale;
+
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
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.AbstractIntegrator;
+import org.apache.commons.math3.ode.ExpandableStatefulODE;
import org.apache.commons.math3.ode.FirstOrderIntegrator;
import org.apache.commons.math3.ode.TestProblem1;
import org.apache.commons.math3.ode.TestProblem5;
import org.apache.commons.math3.ode.TestProblem6;
+import org.apache.commons.math3.ode.TestProblemAbstract;
import org.apache.commons.math3.ode.TestProblemHandler;
+import org.apache.commons.math3.ode.sampling.StepHandler;
+import org.apache.commons.math3.ode.sampling.StepInterpolator;
+import org.apache.commons.math3.util.CombinatoricsUtils;
import org.apache.commons.math3.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
@@ -89,11 +103,11 @@ public class AdamsMoultonIntegratorTest {
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- // the 0.5 and 11.0 factors are only valid for this test
+ // the 0.45 and 8.69 factors are only valid for this test
// and has been obtained from trial and error
// there is no general relation between local and global errors
- Assert.assertTrue(handler.getMaximalValueError() > ( 0.5 * scalAbsoluteTolerance));
- Assert.assertTrue(handler.getMaximalValueError() < (11.0 * scalAbsoluteTolerance));
+ Assert.assertTrue(handler.getMaximalValueError() > (0.45 * scalAbsoluteTolerance));
+ Assert.assertTrue(handler.getMaximalValueError() < (8.69 * scalAbsoluteTolerance));
Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
int calls = pb.getCalls();
@@ -137,8 +151,8 @@ public class AdamsMoultonIntegratorTest {
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- Assert.assertTrue(handler.getLastError() < 1.0e-9);
- Assert.assertTrue(handler.getMaximalValueError() < 1.0e-9);
+ Assert.assertTrue(handler.getLastError() < 3.0e-9);
+ Assert.assertTrue(handler.getMaximalValueError() < 3.0e-9);
Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
Assert.assertEquals("Adams-Moulton", integ.getName());
}
@@ -153,15 +167,176 @@ public class AdamsMoultonIntegratorTest {
for (int nSteps = 2; nSteps < 8; ++nSteps) {
AdamsMoultonIntegrator integ =
new AdamsMoultonIntegrator(nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-5, 1.0e-5);
+ integ.setStarterIntegrator(new PerfectStarter(pb, nSteps));
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- if (nSteps < 4) {
- Assert.assertTrue(handler.getMaximalValueError() > 7.0e-04);
+ if (nSteps < 5) {
+ Assert.assertTrue(handler.getMaximalValueError() > 2.2e-05);
} else {
- Assert.assertTrue(handler.getMaximalValueError() < 3.0e-13);
+ Assert.assertTrue(handler.getMaximalValueError() < 1.1e-11);
+ }
+ }
+
+ }
+
+ private static class PerfectStarter extends AbstractIntegrator {
+
+ private final PerfectInterpolator interpolator;
+ private final int nbSteps;
+
+ public PerfectStarter(final TestProblemAbstract problem, final int nbSteps) {
+ this.interpolator = new PerfectInterpolator(problem);
+ this.nbSteps = nbSteps;
+ }
+
+ public void integrate(ExpandableStatefulODE equations, double t) {
+ double tStart = equations.getTime() + 0.01 * (t - equations.getTime());
+ for (int i = 0; i < nbSteps; ++i) {
+ double tK = ((nbSteps - 1 - (i + 1)) * equations.getTime() + (i + 1) * tStart) / (nbSteps - 1);
+ interpolator.setPreviousTime(interpolator.getCurrentTime());
+ interpolator.setCurrentTime(tK);
+ interpolator.setInterpolatedTime(tK);
+ for (StepHandler handler : getStepHandlers()) {
+ handler.handleStep(interpolator, i == nbSteps - 1);
+ }
+ }
+ }
+
+ }
+
+ private static class PerfectInterpolator implements StepInterpolator {
+ private final TestProblemAbstract problem;
+ private double previousTime;
+ private double currentTime;
+ private double interpolatedTime;
+
+ public PerfectInterpolator(final TestProblemAbstract problem) {
+ this.problem = problem;
+ this.previousTime = problem.getInitialTime();
+ this.currentTime = problem.getInitialTime();
+ this.interpolatedTime = problem.getInitialTime();
+ }
+
+ public void readExternal(ObjectInput arg0) {
+ }
+
+ public void writeExternal(ObjectOutput arg0) {
+ }
+
+ public double getPreviousTime() {
+ return previousTime;
+ }
+
+ public void setPreviousTime(double time) {
+ previousTime = time;
+ }
+
+ public double getCurrentTime() {
+ return currentTime;
+ }
+
+ public void setCurrentTime(double time) {
+ currentTime = time;
+ }
+
+ public double getInterpolatedTime() {
+ return interpolatedTime;
+ }
+
+ public void setInterpolatedTime(double time) {
+ interpolatedTime = time;
+ }
+
+ public double[] getInterpolatedState() {
+ return problem.computeTheoreticalState(interpolatedTime);
+ }
+
+ public double[] getInterpolatedDerivatives() {
+ double[] y = problem.computeTheoreticalState(interpolatedTime);
+ double[] yDot = new double[y.length];
+ problem.computeDerivatives(interpolatedTime, y, yDot);
+ return yDot;
+ }
+
+ public double[] getInterpolatedSecondaryState(int index) {
+ return null;
+ }
+
+ public double[] getInterpolatedSecondaryDerivatives(int index) {
+ return null;
+ }
+
+ public boolean isForward() {
+ return problem.getFinalTime() > problem.getInitialTime();
+ }
+
+ public StepInterpolator copy() {
+ return this;
+ }
+
+ }
+
+ @Test
+ public void testTmp() throws IOException, DimensionMismatchException, NumberIsTooSmallException, MaxCountExceededException, NoBracketingException {
+ final PrintStream out = new PrintStream(new File(new File(System.getProperty("user.home")), "x.dat"));
+ final int n = 4;
+ final AdamsMoultonIntegrator integ = new AdamsMoultonIntegrator(n, 1.0e-12, 1.0e3, 1.0e-10, 1.0e-12);
+ final SinT4 sinT4 = new SinT4();
+ final StepHandler handler = new StepHandler() {
+ double previousError = 0;
+ public void init(double t0, double[] y0, double t) {
}
+ public void handleStep(StepInterpolator interpolator, boolean isLast) {
+ final double t = interpolator.getCurrentTime();
+ final double h = t - interpolator.getPreviousTime();
+ final double y = interpolator.getInterpolatedState()[0];
+ final double error = y - sinT4.computeTheoreticalState(t)[0];
+ out.format(Locale.US, "%10.8f %16.9e %10.8f %16.9e",
+ t, h, y, error - previousError);
+ for (int i = 1; i < n + 2; ++i) {
+ out.format(Locale.US, " %16.9e", sinT4.scaledDerivative(t, h, i));
+ }
+ out.format(Locale.US, "%n");
+ previousError = error;
+ }
+ };
+ integ.addStepHandler(handler);
+ final ExpandableStatefulODE expandableODE = new ExpandableStatefulODE(sinT4);
+ final double t0 = 0.75;
+ final double h0 = 1.0 / 128;
+ expandableODE.setTime(t0);
+ expandableODE.setPrimaryState(sinT4.computeTheoreticalState(t0));
+ integ.setInitialStepSize(h0);
+ integ.integrate(expandableODE, 0.9);
+ out.close();
+ }
+
+ private static class SinT4 extends TestProblemAbstract {
+
+ @Override
+ public int getDimension() {
+ return 1;
+ }
+
+ @Override
+ public void doComputeDerivatives(double t, double[] y, double[] yDot) {
+ // compute the derivatives
+ double t3 = FastMath.pow(t, 3);
+ yDot[0] = 4 * t3 * FastMath.cos(t3 * t);
+ }
+
+ public double[] computeTheoreticalState(final double t) {
+ return new double[] { FastMath.sin(FastMath.pow(t, 4)) };
+ }
+
+ public DerivativeStructure valueDS(final double t, final int n) {
+ return new DerivativeStructure(1, n, 0, t).pow(4).sin();
+ }
+
+ public double scaledDerivative(final double t, final double h, final int n) {
+ return FastMath.pow(h, n) * valueDS(t, n).getPartialDerivative(n) / CombinatoricsUtils.factorial(n);
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformerTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformerTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformerTest.java
new file mode 100644
index 0000000..42e8363
--- /dev/null
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/AdamsNordsieckTransformerTest.java
@@ -0,0 +1,120 @@
+/*
+ * 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 org.apache.commons.math3.analysis.UnivariateFunction;
+import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
+import org.apache.commons.math3.linear.Array2DRowRealMatrix;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class AdamsNordsieckTransformerTest {
+
+ @Test
+ public void testPolynomialExtraDerivative() {
+ checkNordsieckStart(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }),
+ 5, 0.0, 0.125, 3.2e-16);
+ }
+
+ @Test
+ public void testPolynomialRegular() {
+ checkNordsieckStart(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }),
+ 4, 0.0, 0.125, 3.1e-16);
+ }
+
+ @Test
+ public void testPolynomialMissingLastDerivative() {
+ // this test intentionally uses not enough start points,
+ // the Nordsieck vector is therefore not expected to match the exact scaled derivatives
+ checkNordsieckStart(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }),
+ 3, 0.0, 0.125, 1.6e-4);
+ }
+
+ @Test
+ public void testTransformExact() {
+ // a 5 steps transformer handles a degree 5 polynomial exactly
+ // the Nordsieck vector holds the full information about the function
+ // transforming the vector from t0 to t0+h or recomputing it from scratch
+ // at t0+h yields the same result
+ checkTransform(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }), 5, 2.567e-15);
+ }
+
+ @Test
+ public void testTransformInexact() {
+ // a 4 steps transformer cannot handle a degree 5 polynomial exactly
+ // the Nordsieck vector lacks some high degree information about the function
+ // transforming the vector from t0 to t0+h or recomputing it from scratch
+ // at t0+h yields different results
+ checkTransform(new PolynomialFunction(new double[] { 6, 5, 4, 3, 2, 1 }), 4, 5.658e-4);
+ }
+
+ private void checkNordsieckStart(final PolynomialFunction polynomial, final int nbSteps, final double t0,
+ final double h, final double epsilon) {
+
+ final AdamsNordsieckTransformer transformer = AdamsNordsieckTransformer.getInstance(nbSteps);
+ PolynomialFunction derivative = polynomial.polynomialDerivative();
+ final Array2DRowRealMatrix nordsieck = start(transformer, nbSteps, t0, h, polynomial, derivative);
+
+ Assert.assertEquals(nbSteps - 1, nordsieck.getRowDimension());
+ double coeff = h;
+ for (int i = 0; i < nordsieck.getRowDimension(); ++i) {
+ coeff *= h / (i + 2);
+ derivative = derivative.polynomialDerivative();
+ Assert.assertEquals(derivative.value(t0) * coeff, nordsieck.getEntry(i, 0), epsilon);
+ }
+
+ }
+
+ private void checkTransform(final PolynomialFunction polynomial, final int nbSteps, final double expectedError) {
+
+ final AdamsNordsieckTransformer transformer = AdamsNordsieckTransformer.getInstance(nbSteps);
+ final PolynomialFunction derivative = polynomial.polynomialDerivative();
+
+ final double t0 = 0.0;
+ final double h = 0.125;
+ final Array2DRowRealMatrix n0 = start(transformer, nbSteps, t0, h, polynomial, derivative);
+ final Array2DRowRealMatrix n1 = transformer.updateHighOrderDerivativesPhase1(n0);
+ transformer.updateHighOrderDerivativesPhase2(new double[] { h * derivative.value(t0) },
+ new double[] { h * derivative.value(t0 + h) },
+ n1);
+ final Array2DRowRealMatrix n2 = start(transformer, nbSteps, t0 + h, h, polynomial, derivative);
+
+ Assert.assertEquals(expectedError, n2.subtract(n1).getNorm(), expectedError * 0.001);
+
+ }
+
+ private Array2DRowRealMatrix start(final AdamsNordsieckTransformer transformer, final int nbSteps,
+ final double t0, final double h,
+ final UnivariateFunction f0, final UnivariateFunction f1) {
+
+ final int nbStartPoints = (nbSteps + 3) / 2;
+ final double[] t = new double[nbStartPoints];
+ final double[][] y = new double[nbStartPoints][1];
+ final double[][] yDot = new double[nbStartPoints][1];
+ for (int i = 0; i < nbStartPoints; ++i) {
+ t[i] = t0 + i * h;
+ y[i][0] = f0.value(t[i]);
+ yDot[i][0] = f1.value(t[i]);
+ }
+
+ return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
+
+ }
+
+}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java
index 6f618cd..94aab6b 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java
@@ -45,7 +45,7 @@ public class ClassicalRungeKuttaStepInterpolatorTest {
TestProblem3 pb = new TestProblem3();
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
ClassicalRungeKuttaIntegrator integ = new ClassicalRungeKuttaIntegrator(step);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 6.6e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolatorTest.java
index 2b3b8e7..4f41853 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince54StepInterpolatorTest.java
@@ -52,7 +52,7 @@ public class DormandPrince54StepInterpolatorTest {
DormandPrince54Integrator integ = new DormandPrince54Integrator(minStep, maxStep,
scalAbsoluteTolerance,
scalRelativeTolerance);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 3.6e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853StepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853StepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853StepInterpolatorTest.java
index e6b8fd5..499272a 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853StepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/DormandPrince853StepInterpolatorTest.java
@@ -52,7 +52,7 @@ public class DormandPrince853StepInterpolatorTest {
DormandPrince853Integrator integ = new DormandPrince853Integrator(minStep, maxStep,
scalAbsoluteTolerance,
scalRelativeTolerance);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 1.8e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerStepInterpolatorTest.java
index de99547..87a5466 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/EulerStepInterpolatorTest.java
@@ -136,7 +136,7 @@ public class EulerStepInterpolatorTest {
TestProblem3 pb = new TestProblem3();
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
EulerIntegrator integ = new EulerIntegrator(step);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 5.1e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/GillStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/GillStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/GillStepInterpolatorTest.java
index d441ac0..8d3ba7e 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/GillStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/GillStepInterpolatorTest.java
@@ -45,7 +45,7 @@ public class GillStepInterpolatorTest {
TestProblem3 pb = new TestProblem3();
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
GillIntegrator integ = new GillIntegrator(step);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 6.6e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/GraggBulirschStoerStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/GraggBulirschStoerStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/GraggBulirschStoerStepInterpolatorTest.java
index e71e5c0..6c163c5 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/GraggBulirschStoerStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/GraggBulirschStoerStepInterpolatorTest.java
@@ -53,7 +53,7 @@ public class GraggBulirschStoerStepInterpolatorTest {
GraggBulirschStoerIntegrator integ =
new GraggBulirschStoerIntegrator(minStep, maxStep,
absTolerance, relTolerance);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-8);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 5.9e-10);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/HighamHall54StepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/HighamHall54StepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/HighamHall54StepInterpolatorTest.java
index 7962e94..1242e68 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/HighamHall54StepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/HighamHall54StepInterpolatorTest.java
@@ -52,7 +52,7 @@ public class HighamHall54StepInterpolatorTest {
HighamHall54Integrator integ = new HighamHall54Integrator(minStep, maxStep,
scalAbsoluteTolerance,
scalRelativeTolerance);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.1e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 4.8e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
index d29fe20..de137b9 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/LutherStepInterpolatorTest.java
@@ -45,7 +45,7 @@ public class LutherStepInterpolatorTest {
TestProblem3 pb = new TestProblem3();
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
LutherIntegrator integ = new LutherIntegrator(step);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 6.5e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/MidpointStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/MidpointStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/MidpointStepInterpolatorTest.java
index 0f32602..f2d4460 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/MidpointStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/MidpointStepInterpolatorTest.java
@@ -46,7 +46,7 @@ public class MidpointStepInterpolatorTest {
TestProblem3 pb = new TestProblem3();
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
MidpointIntegrator integ = new MidpointIntegrator(step);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 6.6e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java
index 965bf6b..b47cf37 100644
--- a/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java
@@ -45,7 +45,7 @@ public class ThreeEighthesStepInterpolatorTest {
TestProblem3 pb = new TestProblem3();
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
ThreeEighthesIntegrator integ = new ThreeEighthesIntegrator(step);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 1.0e-10);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.01, 6.6e-12);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/sampling/NordsieckStepInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/sampling/NordsieckStepInterpolatorTest.java b/src/test/java/org/apache/commons/math3/ode/sampling/NordsieckStepInterpolatorTest.java
index f11df71..0f00fc6 100644
--- a/src/test/java/org/apache/commons/math3/ode/sampling/NordsieckStepInterpolatorTest.java
+++ b/src/test/java/org/apache/commons/math3/ode/sampling/NordsieckStepInterpolatorTest.java
@@ -44,7 +44,7 @@ public class NordsieckStepInterpolatorTest {
MaxCountExceededException, NoBracketingException {
TestProblem3 pb = new TestProblem3();
AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(4, 0.0, 1.0, 1.0e-10, 1.0e-10);
- StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 5e-9);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb, 0.05, 2.8e-9);
}
@Test
@@ -66,8 +66,8 @@ public class NordsieckStepInterpolatorTest {
oos.writeObject(handler);
}
- Assert.assertTrue(bos.size () > 25500);
- Assert.assertTrue(bos.size () < 26500);
+ Assert.assertTrue(bos.size() > 47000);
+ Assert.assertTrue(bos.size() < 48000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);
@@ -79,7 +79,7 @@ public class NordsieckStepInterpolatorTest {
double r = random.nextDouble();
double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
cm.setInterpolatedTime(time);
- double[] interpolatedY = cm.getInterpolatedState ();
+ double[] interpolatedY = cm.getInterpolatedState();
double[] theoreticalY = pb.computeTheoreticalState(time);
double dx = interpolatedY[0] - theoreticalY[0];
double dy = interpolatedY[1] - theoreticalY[1];
http://git-wip-us.apache.org/repos/asf/commons-math/blob/6a01c7b2/src/test/java/org/apache/commons/math3/ode/sampling/StepInterpolatorTestUtils.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/ode/sampling/StepInterpolatorTestUtils.java b/src/test/java/org/apache/commons/math3/ode/sampling/StepInterpolatorTestUtils.java
index 9450df2..184f3e0 100644
--- a/src/test/java/org/apache/commons/math3/ode/sampling/StepInterpolatorTestUtils.java
+++ b/src/test/java/org/apache/commons/math3/ode/sampling/StepInterpolatorTestUtils.java
@@ -35,6 +35,7 @@ public class StepInterpolatorTestUtils {
public static void checkDerivativesConsistency(final FirstOrderIntegrator integrator,
final TestProblemAbstract problem,
+ final double finiteDifferencesRatio,
final double threshold)
throws DimensionMismatchException, NumberIsTooSmallException,
MaxCountExceededException, NoBracketingException {
@@ -43,8 +44,9 @@ public class StepInterpolatorTestUtils {
public void handleStep(StepInterpolator interpolator, boolean isLast)
throws MaxCountExceededException {
- final double h = 0.001 * (interpolator.getCurrentTime() - interpolator.getPreviousTime());
- final double t = interpolator.getCurrentTime() - 300 * h;
+ final double dt = interpolator.getCurrentTime() - interpolator.getPreviousTime();
+ final double h = finiteDifferencesRatio * dt;
+ final double t = interpolator.getCurrentTime() - 0.3 * dt;
if (FastMath.abs(h) < 10 * FastMath.ulp(t)) {
return;
@@ -75,7 +77,7 @@ public class StepInterpolatorTestUtils {
32 * (yP3h[i] - yM3h[i]) +
-168 * (yP2h[i] - yM2h[i]) +
672 * (yP1h[i] - yM1h[i])) / (840 * h);
- Assert.assertEquals(approYDot, yDot[i], threshold);
+ Assert.assertEquals("" + (approYDot - yDot[i]), approYDot, yDot[i], threshold);
}
}