Fixed stability issues with Adams integrators.
This commit is contained in:
Luc Maisonobe
parent
9a87c766bb
commit
bf803b119b
|
|
@ -54,6 +54,12 @@ If the output is not quite correct, check for invisible trailing spaces!
|
|||
</release>
|
||||
|
||||
<release version="4.0" 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="erans" type="fix" issue="MATH-1301">
|
||||
"JDKRandomGenerator": Method "nextInt(int)" now throws a "NotStrictlyPositiveException".
|
||||
</action>
|
||||
|
|
|
|||
|
|
@ -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 {
|
||||
|
|
@ -315,6 +315,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
|
||||
|
|
@ -323,7 +330,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
|
||||
|
|
@ -354,14 +364,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} */
|
||||
|
|
@ -373,7 +383,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();
|
||||
|
|
@ -388,7 +398,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;
|
||||
|
|
@ -406,7 +416,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);
|
||||
|
||||
|
|
|
|||
|
|
@ -22,6 +22,7 @@ import org.apache.commons.math4.exception.MaxCountExceededException;
|
|||
import org.apache.commons.math4.exception.NoBracketingException;
|
||||
import org.apache.commons.math4.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math4.linear.Array2DRowRealMatrix;
|
||||
import org.apache.commons.math4.linear.RealMatrix;
|
||||
import org.apache.commons.math4.ode.EquationsMapper;
|
||||
import org.apache.commons.math4.ode.ExpandableStatefulODE;
|
||||
import org.apache.commons.math4.ode.sampling.NordsieckStepInterpolator;
|
||||
|
|
@ -90,7 +91,7 @@ import org.apache.commons.math4.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.math4.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) {
|
||||
|
|
|
|||
|
|
@ -96,7 +96,7 @@ import org.apache.commons.math4.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.math4.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)
|
||||
|
|
@ -407,10 +407,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;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -68,7 +68,7 @@ import org.apache.commons.math4.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.math4.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).
|
||||
|
|
|
|||
|
|
@ -88,35 +88,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);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -135,6 +136,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
|
||||
|
|
|
|||
|
|
@ -18,16 +18,23 @@
|
|||
package org.apache.commons.math4.ode.nonstiff;
|
||||
|
||||
|
||||
import java.io.ObjectInput;
|
||||
import java.io.ObjectOutput;
|
||||
|
||||
import org.apache.commons.math4.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math4.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math4.exception.NoBracketingException;
|
||||
import org.apache.commons.math4.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math4.ode.AbstractIntegrator;
|
||||
import org.apache.commons.math4.ode.ExpandableStatefulODE;
|
||||
import org.apache.commons.math4.ode.FirstOrderIntegrator;
|
||||
import org.apache.commons.math4.ode.TestProblem1;
|
||||
import org.apache.commons.math4.ode.TestProblem5;
|
||||
import org.apache.commons.math4.ode.TestProblem6;
|
||||
import org.apache.commons.math4.ode.TestProblemAbstract;
|
||||
import org.apache.commons.math4.ode.TestProblemHandler;
|
||||
import org.apache.commons.math4.ode.nonstiff.AdamsBashforthIntegrator;
|
||||
import org.apache.commons.math4.ode.sampling.StepHandler;
|
||||
import org.apache.commons.math4.ode.sampling.StepInterpolator;
|
||||
import org.apache.commons.math4.util.FastMath;
|
||||
import org.junit.Assert;
|
||||
import org.junit.Test;
|
||||
|
|
@ -65,8 +72,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) {
|
||||
|
|
@ -85,12 +91,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);
|
||||
|
|
@ -123,14 +128,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());
|
||||
}
|
||||
|
|
@ -142,18 +148,116 @@ 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;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -18,16 +18,24 @@
|
|||
package org.apache.commons.math4.ode.nonstiff;
|
||||
|
||||
|
||||
import java.io.ObjectInput;
|
||||
import java.io.ObjectOutput;
|
||||
|
||||
import org.apache.commons.math4.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math4.exception.MaxCountExceededException;
|
||||
import org.apache.commons.math4.exception.NoBracketingException;
|
||||
import org.apache.commons.math4.exception.NumberIsTooSmallException;
|
||||
import org.apache.commons.math4.ode.AbstractIntegrator;
|
||||
import org.apache.commons.math4.ode.ExpandableStatefulODE;
|
||||
import org.apache.commons.math4.ode.FirstOrderIntegrator;
|
||||
import org.apache.commons.math4.ode.TestProblem1;
|
||||
import org.apache.commons.math4.ode.TestProblem5;
|
||||
import org.apache.commons.math4.ode.TestProblem6;
|
||||
import org.apache.commons.math4.ode.TestProblemAbstract;
|
||||
import org.apache.commons.math4.ode.TestProblemHandler;
|
||||
import org.apache.commons.math4.ode.nonstiff.AdamsMoultonIntegrator;
|
||||
import org.apache.commons.math4.ode.sampling.StepHandler;
|
||||
import org.apache.commons.math4.ode.sampling.StepInterpolator;
|
||||
import org.apache.commons.math4.util.FastMath;
|
||||
import org.junit.Assert;
|
||||
import org.junit.Test;
|
||||
|
|
@ -90,11 +98,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();
|
||||
|
|
@ -138,8 +146,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());
|
||||
}
|
||||
|
|
@ -154,17 +162,115 @@ 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;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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.math4.ode.nonstiff;
|
||||
|
||||
|
||||
import org.apache.commons.math4.analysis.UnivariateFunction;
|
||||
import org.apache.commons.math4.analysis.polynomials.PolynomialFunction;
|
||||
import org.apache.commons.math4.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);
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -46,7 +46,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
|
||||
|
|
|
|||
|
|
@ -53,7 +53,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
|
||||
|
|
|
|||
|
|
@ -53,7 +53,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
|
||||
|
|
|
|||
|
|
@ -140,7 +140,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
|
||||
|
|
|
|||
|
|
@ -46,7 +46,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
|
||||
|
|
|
|||
|
|
@ -54,7 +54,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
|
||||
|
|
|
|||
|
|
@ -53,7 +53,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
|
||||
|
|
|
|||
|
|
@ -46,7 +46,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
|
||||
|
|
|
|||
|
|
@ -47,7 +47,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
|
||||
|
|
|
|||
|
|
@ -46,7 +46,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
|
||||
|
|
|
|||
|
|
@ -45,7 +45,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
|
||||
|
|
@ -67,8 +67,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);
|
||||
|
|
@ -80,7 +80,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];
|
||||
|
|
|
|||
|
|
@ -32,6 +32,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 {
|
||||
|
|
@ -40,8 +41,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;
|
||||
|
|
@ -72,7 +74,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);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue