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 <l...@apache.org> Authored: Thu Dec 24 22:49:59 2015 +0100 Committer: Luc Maisonobe <l...@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); } }
