You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2009/06/20 20:15:54 UTC
svn commit: r786877 [1/2] - in /commons/proper/math/trunk/src:
java/org/apache/commons/math/ode/ java/org/apache/commons/math/ode/nonstiff/
java/org/apache/commons/math/ode/sampling/
test/org/apache/commons/math/ode/nonstiff/ test/org/apache/commons/ma...
Author: luc
Date: Sat Jun 20 18:15:54 2009
New Revision: 786877
URL: http://svn.apache.org/viewvc?rev=786877&view=rev
Log:
changed the Adams-Bashforth and Adams-Moulton multistep integrators to adaptive stepsize.
this was made possible thanks to the Nordsieck representation of integration state
Added:
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java (with props)
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
Added: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java?rev=786877&view=auto
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java (added)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java Sat Jun 20 18:15:54 2009
@@ -0,0 +1,326 @@
+/*
+ * 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.math.ode;
+
+import java.util.Arrays;
+import java.util.HashMap;
+import java.util.Map;
+
+import org.apache.commons.math.fraction.BigFraction;
+import org.apache.commons.math.linear.Array2DRowFieldMatrix;
+import org.apache.commons.math.linear.Array2DRowRealMatrix;
+import org.apache.commons.math.linear.DefaultFieldMatrixChangingVisitor;
+import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
+import org.apache.commons.math.linear.FieldDecompositionSolver;
+import org.apache.commons.math.linear.FieldLUDecompositionImpl;
+import org.apache.commons.math.linear.FieldMatrix;
+import org.apache.commons.math.linear.MatrixUtils;
+import org.apache.commons.math.linear.MatrixVisitorException;
+import org.apache.commons.math.linear.RealMatrix;
+import org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator;
+import org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator;
+
+/** Transformer for Nordsieck vectors.
+ * <p>This class i used by {@link MultistepIntegrator multistep integrators}
+ * to convert between classical representation with several previous first
+ * derivatives and Nordsieck representation with higher order scaled derivatives.</p>
+ *
+ * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
+ * <pre>
+ * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
+ * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
+ * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
+ * ...
+ * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup> derivative
+ * </pre></p>
+ *
+ * <p>With the previous definition, the classical representation of multistep methods
+ * uses first derivatives only, i.e. it handles y<sub>n</sub>, s<sub>1</sub>(n) and
+ * q<sub>n</sub> where q<sub>n</sub> is defined as:
+ * <pre>
+ * q<sub>n</sub> = [ s<sub>1</sub>(n-1) s<sub>1</sub>(n-2) ... s<sub>1</sub>(n-(k-1)) ]<sup>T</sup>
+ * </pre>
+ * (we omit the k index in the notation for clarity).</p>
+ *
+ * <p>Another possible representation uses the Nordsieck vector with
+ * higher degrees scaled derivatives all taken at the same step, i.e it handles y<sub>n</sub>,
+ * s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
+ * <pre>
+ * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
+ * </pre>
+ * (here again we omit the k index in the notation for clarity)
+ * </p>
+ *
+ * <p>Taylor series formulas show that for any index offset i, s<sub>1</sub>(n-i) can be
+ * 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)
+ * </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>
+ * and q<sub>n</sub> resulting from the Taylor series formulas above is:
+ * <pre>
+ * 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:
+ * <pre>
+ * [ -2 3 -4 5 ... ]
+ * [ -4 12 -32 80 ... ]
+ * P = [ -6 27 -108 405 ... ]
+ * [ -8 48 -256 1280 ... ]
+ * [ ... ]
+ * </pre></p>
+ *
+ * <p>Changing -i into +i in the formula above can be used to compute a similar transform between
+ * classical representation and Nordsieck vector at step start. The resulting Q matrix is simply
+ * the absolute value of matrix P.</p>
+ *
+ * <p>Using the Nordsieck vector has several advantages:
+ * <ul>
+ * <li>it greatly simplifies step interpolation as the interpolator mainly applies
+ * Taylor series formulas,</li>
+ * <li>it simplifies step changes that occur when discrete events that truncate
+ * the step are triggered,</li>
+ * <li>it allows to extend the methods in order to support adaptive stepsize (not implemented yet).</li>
+ * </ul></p>
+ *
+ * <p>For {@link AdamsBashforthIntegrator Adams-Bashforth} method, the Nordsieck vector
+ * at step n+1 is computed from the Nordsieck vector at step n as follows:
+ * <ul>
+ * <li>y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
+ * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
+ * <li>r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub></li>
+ * </ul>
+ * where A is a rows shifting matrix (the lower left part is an identity matrix):
+ * <pre>
+ * [ 0 0 ... 0 0 | 0 ]
+ * [ ---------------+---]
+ * [ 1 0 ... 0 0 | 0 ]
+ * A = [ 0 1 ... 0 0 | 0 ]
+ * [ ... | 0 ]
+ * [ 0 0 ... 1 0 | 0 ]
+ * [ 0 0 ... 0 1 | 0 ]
+ * </pre></p>
+ *
+ * <p>For {@link AdamsMoultonIntegrator Adams-Moulton} method, the predicted Nordsieck vector
+ * at step n+1 is computed from the Nordsieck vector at step n as follows:
+ * <ul>
+ * <li>Y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n) + u<sup>T</sup> r<sub>n</sub></li>
+ * <li>S<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, Y<sub>n+1</sub>)</li>
+ * <li>R<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub></li>
+ * </ul>
+ * where A is a rows shifting matrix (the lower left part is an identity matrix):
+ * <pre>
+ * [ 0 0 ... 0 0 | 0 ]
+ * [ ---------------+---]
+ * [ 1 0 ... 0 0 | 0 ]
+ * A = [ 0 1 ... 0 0 | 0 ]
+ * [ ... | 0 ]
+ * [ 0 0 ... 1 0 | 0 ]
+ * [ 0 0 ... 0 1 | 0 ]
+ * </pre>
+ * From this predicted vector, the corrected vector is computed as follows:
+ * <ul>
+ * <li>y<sub>n+1</sub> = y<sub>n</sub> + S<sub>1</sub>(n+1) + [ -1 +1 -1 +1 ... ±1 ] r<sub>n+1</sub></li>
+ * <li>s<sub>1</sub>(n+1) = h f(t<sub>n+1</sub>, y<sub>n+1</sub>)</li>
+ * <li>r<sub>n+1</sub> = R<sub>n+1</sub> + (s<sub>1</sub>(n+1) - S<sub>1</sub>(n+1)) P<sup>-1</sup> u</li>
+ * </ul>
+ * where the upper case Y<sub>n+1</sub>, S<sub>1</sub>(n+1) and R<sub>n+1</sub> represent the
+ * predicted states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
+ * represent the corrected states.</p>
+ *
+ * <p>We observe that both methods use similar update formulas. In both cases a P<sup>-1</sup>u
+ * vector and a P<sup>-1</sup> A P matrix are used that do not depend on the state,
+ * they only depend on k. This class handles these transformations.</p>
+ *
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class NordsieckTransformer {
+
+ /** Cache for already computed coefficients. */
+ private static final Map<Integer, NordsieckTransformer> cache =
+ new HashMap<Integer, NordsieckTransformer>();
+
+ /** Initialization matrix for the higher order derivatives wrt y'', y''' ... */
+ private final RealMatrix initialization;
+
+ /** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */
+ private final RealMatrix 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
+ * (including the one being computed)
+ */
+ private NordsieckTransformer(final int nSteps) {
+
+ // compute exact coefficients
+ FieldMatrix<BigFraction> bigP = buildP(nSteps);
+ FieldDecompositionSolver<BigFraction> pSolver =
+ new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver();
+ BigFraction[] u = new BigFraction[nSteps - 1];
+ Arrays.fill(u, BigFraction.ONE);
+ BigFraction[] bigC1 = pSolver.solve(u);
+
+ // update coefficients are computed by combining transform from
+ // Nordsieck to multistep, then shifting rows to represent step advance
+ // then applying inverse transform
+ BigFraction[][] shiftedP = bigP.getData();
+ for (int i = shiftedP.length - 1; i > 0; --i) {
+ // shift rows
+ shiftedP[i] = shiftedP[i - 1];
+ }
+ shiftedP[0] = new BigFraction[nSteps - 1];
+ Arrays.fill(shiftedP[0], BigFraction.ZERO);
+ FieldMatrix<BigFraction> bigMSupdate =
+ pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
+
+ // initialization coefficients, computed from a Q matrix = abs(P)
+ bigP.walkInOptimizedOrder(new DefaultFieldMatrixChangingVisitor<BigFraction>(BigFraction.ZERO) {
+ /** {@inheritDoc} */
+ @Override
+ public BigFraction visit(int row, int column, BigFraction value) {
+ return ((column & 0x1) == 0x1) ? value : value.negate();
+ }
+ });
+ FieldMatrix<BigFraction> bigQInverse =
+ new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse();
+
+ // convert coefficients to double
+ initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigQInverse);
+ update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
+ c1 = new double[nSteps - 1];
+ for (int i = 0; i < nSteps - 1; ++i) {
+ c1[i] = bigC1[i].doubleValue();
+ }
+
+ }
+
+ /** Get the Nordsieck transformer for a given number of steps.
+ * @param nSteps number of steps of the multistep method
+ * (including the one being computed)
+ * @return Nordsieck transformer for the specified number of steps
+ */
+ public static NordsieckTransformer getInstance(final int nSteps) {
+ synchronized(cache) {
+ NordsieckTransformer t = cache.get(nSteps);
+ if (t == null) {
+ t = new NordsieckTransformer(nSteps);
+ cache.put(nSteps, t);
+ }
+ return t;
+ }
+ }
+
+ /** Build the P matrix transforming multistep to Nordsieck.
+ * <p>
+ * Multistep representation uses y(k), s<sub>1</sub>(k), s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)).
+ * Nordsieck representation uses y(k), s<sub>1</sub>(k), s<sub>2</sub>(k) ... s<sub>n</sub>(k).
+ * The two representations share their two first components y(k) and
+ * s<sub>1</sub>(k). The P matrix is used to transform the remaining ones:
+ * <pre>
+ * [ s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)]<sup>T</sup> = s<sub>1</sub>(k) [1 ... 1]<sup>T</sup> + P [s<sub>2</sub>(k) ... s<sub>n</sub>(k)]<sup>T</sup>
+ * </pre>
+ * </p>
+ * @param nSteps number of steps of the multistep method
+ * (including the one being computed)
+ * @return P matrix
+ */
+ private FieldMatrix<BigFraction> buildP(final int nSteps) {
+
+ final BigFraction[][] pData = new BigFraction[nSteps - 1][nSteps - 1];
+
+ for (int i = 0; i < pData.length; ++i) {
+ // build the P matrix elements from Taylor series formulas
+ final BigFraction[] pI = pData[i];
+ final int factor = -(i + 1);
+ int aj = factor;
+ for (int j = 0; j < pI.length; ++j) {
+ pI[j] = new BigFraction(aj * (j + 2));
+ aj *= factor;
+ }
+ }
+
+ return new Array2DRowFieldMatrix<BigFraction>(pData, false);
+
+ }
+
+ /** Initialize the high order scaled derivatives at step start.
+ * @param first first scaled derivative at step start
+ * @param multistep scaled derivatives after step start (hy'1, ..., hy'k-1)
+ * will be modified
+ * @return high order derivatives at step start
+ */
+ public RealMatrix initializeHighOrderDerivatives(final double[] first,
+ final double[][] multistep) {
+ for (int i = 0; i < multistep.length; ++i) {
+ final double[] msI = multistep[i];
+ for (int j = 0; j < first.length; ++j) {
+ msI[j] -= first[j];
+ }
+ }
+ return initialization.multiply(new Array2DRowRealMatrix(multistep, false));
+ }
+
+ /** Update the high order scaled derivatives (phase 1).
+ * <p>The complete update of high order derivatives has a form similar to:
+ * <pre>
+ * r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub>
+ * </pre>
+ * this method computes the P<sup>-1</sup> A P r<sub>n</sub> part.</p>
+ * @param highOrder high order scaled derivatives
+ * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
+ * @return updated high order derivatives
+ * @see #updateHighOrderDerivativesPhase2(double[], double[], RealMatrix)
+ */
+ public RealMatrix updateHighOrderDerivativesPhase1(final RealMatrix highOrder) {
+ return update.multiply(highOrder);
+ }
+
+ /** Update the high order scaled derivatives (phase 2).
+ * <p>The complete update of high order derivatives has a form similar to:
+ * <pre>
+ * r<sub>n+1</sub> = (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u + P<sup>-1</sup> A P r<sub>n</sub>
+ * </pre>
+ * this method computes the (s<sub>1</sub>(n) - s<sub>1</sub>(n+1)) P<sup>-1</sup> u part.</p>
+ * <p>Phase 1 of the update must already have been performed.</p>
+ * @param start first order scaled derivatives at step start
+ * @param end first order scaled derivatives at step end
+ * @param highOrder high order scaled derivatives, will be modified
+ * (h<sup>2</sup>/2 y'', ... h<sup>k</sup>/k! y(k))
+ * @see #updateHighOrderDerivativesPhase1(RealMatrix)
+ */
+ public void updateHighOrderDerivativesPhase2(final double[] start,
+ final double[] end,
+ final RealMatrix highOrder) {
+ highOrder.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
+ /** {@inheritDoc} */
+ @Override
+ public double visit(int row, int column, double value)
+ throws MatrixVisitorException {
+ return value + c1[row] * (start[column] - end[column]);
+ }
+ });
+ }
+
+}
Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/NordsieckTransformer.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=786877&r1=786876&r2=786877&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Sat Jun 20 18:15:54 2009
@@ -17,25 +17,9 @@
package org.apache.commons.math.ode.nonstiff;
-import java.io.IOException;
-import java.io.ObjectInputStream;
-import java.io.ObjectOutputStream;
import java.io.Serializable;
-import java.lang.reflect.Field;
-import java.util.Arrays;
-import java.util.HashMap;
-import java.util.Map;
-
-import org.apache.commons.math.MathRuntimeException;
-import org.apache.commons.math.fraction.BigFraction;
-import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
-import org.apache.commons.math.linear.FieldLUDecompositionImpl;
-import org.apache.commons.math.linear.FieldMatrix;
-import org.apache.commons.math.linear.Array2DRowFieldMatrix;
-import org.apache.commons.math.linear.MatrixUtils;
-import org.apache.commons.math.linear.MatrixVisitorException;
+
import org.apache.commons.math.linear.RealMatrix;
-import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
@@ -161,48 +145,50 @@
/** Serializable version identifier. */
private static final long serialVersionUID = 67792782787082199L;
- /** Cache for already computed coefficients. */
- private static final Map<Integer, CachedCoefficients> cache =
- new HashMap<Integer, CachedCoefficients>();
-
- /** Coefficients of the method. */
- private final transient CachedCoefficients coefficients;
-
- /** Integration step. */
- private final double step;
-
/**
* Build an Adams-Bashforth with the given order and step size.
* @param order order of the method (must be greater than 1: due to
* an implementation limitation the order 1 method is not supported)
- * @param step integration step size
+ * @param minStep minimal step (must be positive even for backward
+ * integration), the last step can be smaller than this
+ * @param maxStep maximal step (must be positive even for backward
+ * integration)
+ * @param scalAbsoluteTolerance allowed absolute error
+ * @param scalRelativeTolerance allowed relative error
* @exception IllegalArgumentException if order is 1 or less
*/
- public AdamsBashforthIntegrator(final int order, final double step)
+ public AdamsBashforthIntegrator(final int order,
+ final double minStep, final double maxStep,
+ final double scalAbsoluteTolerance,
+ final double scalRelativeTolerance)
throws IllegalArgumentException {
+ super("Adams-Bashforth", order, order, minStep, maxStep,
+ scalAbsoluteTolerance, scalRelativeTolerance);
+ }
- super("Adams-Bashforth", order);
- if (order <= 1) {
- throw MathRuntimeException.createIllegalArgumentException(
- "{0} is supported only for orders 2 or more",
- getName());
- }
-
- // cache the coefficients for each order, to avoid recomputing them
- synchronized(cache) {
- CachedCoefficients coeff = cache.get(order);
- if (coeff == null) {
- coeff = new CachedCoefficients(order);
- cache.put(order, coeff);
- }
- coefficients = coeff;
- }
-
- this.step = Math.abs(step);
-
+ /**
+ * Build an Adams-Bashforth with the given order and step size.
+ * @param order order of the method (must be greater than 1: due to
+ * an implementation limitation the order 1 method is not supported)
+ * @param minStep minimal step (must be positive even for backward
+ * integration), the last step can be smaller than this
+ * @param maxStep maximal step (must be positive even for backward
+ * integration)
+ * @param vecAbsoluteTolerance allowed absolute error
+ * @param vecRelativeTolerance allowed relative error
+ * @exception IllegalArgumentException if order is 1 or less
+ */
+ public AdamsBashforthIntegrator(final int order,
+ final double minStep, final double maxStep,
+ final double[] vecAbsoluteTolerance,
+ final double[] vecRelativeTolerance)
+ throws IllegalArgumentException {
+ super("Adams-Bashforth", order, order, minStep, maxStep,
+ vecAbsoluteTolerance, vecRelativeTolerance);
}
/** {@inheritDoc} */
+ @Override
public double integrate(final FirstOrderDifferentialEquations equations,
final double t0, final double[] y0,
final double t, final double[] y)
@@ -218,66 +204,115 @@
if (y != y0) {
System.arraycopy(y0, 0, y, 0, n);
}
+ final double[] yDot = new double[n];
+ final double[] yTmp = new double[y0.length];
// set up an interpolator sharing the integrator arrays
final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
interpolator.reinitialize(y, forward);
+ final NordsieckStepInterpolator interpolatorTmp = new NordsieckStepInterpolator();
+ interpolatorTmp.reinitialize(yTmp, forward);
// set up integration control objects
- stepStart = t0;
- stepSize = forward ? step : -step;
for (StepHandler handler : stepHandlers) {
handler.reset();
}
CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
- // compute the first few steps using the configured starter integrator
- double stopTime = start(previousF.length, stepSize, manager, stepStart, y);
- if (Double.isNaN(previousT[0])) {
- return stopTime;
- }
- stepStart = previousT[0];
-
- // convert to Nordsieck representation
- double[] scaled = convertToNordsieckLow();
- RealMatrix nordsieck = convertToNordsieckHigh(scaled);
+ // compute the initial Nordsieck vector using the configured starter integrator
+ start(t0, 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;
+ interpolator.rescale(hNew);
+
boolean lastStep = false;
while (!lastStep) {
// shift all data
interpolator.shift();
- // discrete events handling
- interpolator.storeTime(stepStart + stepSize);
- if (manager.evaluateStep(interpolator)) {
- stepSize = manager.getEventTime() - stepStart;
+ double error = 0;
+ for (boolean loop = true; loop;) {
+
+ stepSize = hNew;
+
+ // evaluate error using the last term of the Taylor expansion
+ error = 0;
+ for (int i = 0; i < y0.length; ++i) {
+ final double yScale = Math.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;
+ }
+ error = Math.sqrt(error / y0.length);
+
+ if (error <= 1.0) {
+
+ // predict a first estimate of the state at step end
+ final double stepEnd = stepStart + stepSize;
+ interpolator.setInterpolatedTime(stepEnd);
+ System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length);
+
+ // evaluate the derivative
+ computeDerivatives(stepEnd, yTmp, 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 RealMatrix nordsieckTmp =
+ transformer.updateHighOrderDerivativesPhase1(nordsieck);
+ transformer.updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
+
+ // discrete events handling
+ interpolatorTmp.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
+ interpolatorTmp.storeTime(stepStart);
+ interpolatorTmp.shift();
+ interpolatorTmp.storeTime(stepEnd);
+ if (manager.evaluateStep(interpolatorTmp)) {
+ final double dt = manager.getEventTime() - stepStart;
+ if (Math.abs(dt) <= Math.ulp(stepStart)) {
+ // rejecting the step would lead to a too small next step, we accept it
+ loop = false;
+ } else {
+ // reject the step to match exactly the next switch time
+ hNew = dt;
+ interpolator.rescale(hNew);
+ }
+ } else {
+ // accept the step
+ scaled = predictedScaled;
+ nordsieck = nordsieckTmp;
+ interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
+ loop = false;
+ }
+
+ } else {
+ // reject the step and attempt to reduce error by stepsize control
+ final double factor = computeStepGrowShrinkFactor(error);
+ hNew = filterStep(stepSize * factor, forward, false);
+ interpolator.rescale(hNew);
+ }
+
}
// the step has been accepted (may have been truncated)
final double nextStep = stepStart + stepSize;
+ System.arraycopy(yTmp, 0, y, 0, n);
interpolator.storeTime(nextStep);
- System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, n);
manager.stepAccepted(nextStep, y);
lastStep = manager.stop();
- // update the Nordsieck vector
- final double[] f0 = previousF[0];
- previousT[0] = nextStep;
- computeDerivatives(nextStep, y, f0);
- nordsieck = coefficients.msUpdate.multiply(nordsieck);
- final double[] end = new double[y0.length];
- for (int j = 0; j < y0.length; ++j) {
- end[j] = stepSize * f0[j];
- }
- nordsieck.walkInOptimizedOrder(new NordsieckUpdater(scaled, end, coefficients.c1));
- scaled = end;
- interpolator.reinitialize(nextStep, stepSize, scaled, nordsieck);
-
// provide the step data to the step handler
for (StepHandler handler : stepHandlers) {
+ interpolator.setInterpolatedTime(nextStep);
handler.handleStep(interpolator, lastStep);
}
stepStart = nextStep;
@@ -286,225 +321,33 @@
// some events handler has triggered changes that
// invalidate the derivatives, we need to restart from scratch
- stopTime = start(previousF.length, stepSize, manager, stepStart, y);
- if (Double.isNaN(previousT[0])) {
- return stopTime;
- }
- stepStart = previousT[0];
-
- // convert to Nordsieck representation
- scaled = convertToNordsieckLow();
- nordsieck = convertToNordsieckHigh(scaled);
+ start(stepStart, y, t);
interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
}
+ if (! lastStep) {
+ // in some rare cases we may get here with stepSize = 0, for example
+ // when an event occurs at integration start, reducing the first step
+ // to zero; we have to reset the step to some safe non zero value
+ stepSize = filterStep(stepSize, forward, true);
+
+ // stepsize control for next step
+ final double factor = computeStepGrowShrinkFactor(error);
+ final double scaledH = stepSize * factor;
+ final double nextT = stepStart + scaledH;
+ final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
+ hNew = filterStep(scaledH, forward, nextIsLast);
+ interpolator.rescale(hNew);
+ }
+
}
- stopTime = stepStart;
+ final double stopTime = stepStart;
stepStart = Double.NaN;
stepSize = Double.NaN;
return stopTime;
}
- /** Convert the multistep representation after a restart to Nordsieck representation.
- * @return first scaled derivative
- */
- private double[] convertToNordsieckLow() {
-
- final double[] f0 = previousF[0];
- final double[] scaled = new double[f0.length];
- for (int j = 0; j < f0.length; ++j) {
- scaled[j] = stepSize * f0[j];
- }
- return scaled;
-
- }
-
- /** Convert the multistep representation after a restart to Nordsieck representation.
- * @param scaled first scaled derivative
- * @return Nordsieck matrix of the higher scaled derivatives
- */
- private RealMatrix convertToNordsieckHigh(final double[] scaled) {
-
- final double[] f0 = previousF[0];
- final double[][] multistep = new double[coefficients.msToN.getColumnDimension()][f0.length];
- for (int i = 0; i < multistep.length; ++i) {
- final double[] msI = multistep[i];
- final double[] fI = previousF[i + 1];
- for (int j = 0; j < f0.length; ++j) {
- msI[j] = stepSize * fI[j] - scaled[j];
- }
- }
-
- return coefficients.msToN.multiply(new Array2DRowRealMatrix(multistep, false));
-
- }
-
- /** Updater for Nordsieck vector. */
- private static class NordsieckUpdater extends DefaultRealMatrixChangingVisitor {
-
- /** Scaled first derivative at step start. */
- private final double[] start;
-
- /** Scaled first derivative at step end. */
- private final double[] end;
-
- /** Update coefficients. */
- private final double[] c1;
-
- /** Simple constructor.
- * @param start scaled first derivative at step start
- * @param end scaled first derivative at step end
- * @param c1 update coefficients
- */
- public NordsieckUpdater(final double[] start, final double[] end,
- final double[] c1) {
- this.start = start;
- this.end = end;
- this.c1 = c1;
- }
-
- /** {@inheritDoc} */
- @Override
- public double visit(int row, int column, double value)
- throws MatrixVisitorException {
- return value + c1[row] * (start[column] - end[column]);
- }
-
- }
-
- /** Cache for already computed coefficients. */
- private static class CachedCoefficients {
-
- /** Transformer between multistep and Nordsieck representations. */
- private final RealMatrix msToN;
-
- /** Update coefficients of the higher order derivatives wrt y'', y''' ... */
- private final RealMatrix msUpdate;
-
- /** Update coefficients of the higher order derivatives wrt y'. */
- private final double[] c1;
-
- /** Simple constructor.
- * @param order order of the method (must be greater than 1: due to
- * an implementation limitation the order 1 method is not supported)
- */
- public CachedCoefficients(int order) {
-
- // compute exact coefficients
- FieldMatrix<BigFraction> bigNtoMS = buildP(order);
- FieldMatrix<BigFraction> bigMStoN =
- new FieldLUDecompositionImpl<BigFraction>(bigNtoMS).getSolver().getInverse();
- BigFraction[] u = new BigFraction[order - 1];
- Arrays.fill(u, BigFraction.ONE);
- BigFraction[] bigC1 = bigMStoN.operate(u);
-
- // update coefficients are computed by combining transform from
- // Nordsieck to multistep, then shifting rows to represent step advance
- // then applying inverse transform
- BigFraction[][] shiftedP = bigNtoMS.getData();
- for (int i = shiftedP.length - 1; i > 0; --i) {
- // shift rows
- shiftedP[i] = shiftedP[i - 1];
- }
- shiftedP[0] = new BigFraction[order - 1];
- Arrays.fill(shiftedP[0], BigFraction.ZERO);
- FieldMatrix<BigFraction> bigMSupdate =
- bigMStoN.multiply(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
-
- // convert coefficients to double
- msToN = MatrixUtils.bigFractionMatrixToRealMatrix(bigMStoN);
- msUpdate = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
- c1 = new double[order - 1];
- for (int i = 0; i < order - 1; ++i) {
- c1[i] = bigC1[i].doubleValue();
- }
-
- }
-
- /** Build the P matrix transforming multistep to Nordsieck.
- * <p>
- * <p>
- * Multistep representation uses y(k), s<sub>1</sub>(k), s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)).
- * Nordsieck representation uses y(k), s<sub>1</sub>(k), s<sub>2</sub>(k) ... s<sub>n</sub>(k).
- * The two representations share their two first components y(k) and
- * s<sub>1</sub>(k). The P matrix is used to transform the remaining ones:
- * <pre>
- * [ s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)]<sup>T</sup> = s<sub>1</sub>(k) [1 ... 1]<sup>T</sup> + P [s<sub>2</sub>(k) ... s<sub>n</sub>(k)]<sup>T</sup>
- * </pre>
- * </p>
- * @param order order of the method (must be strictly positive)
- * @return P matrix
- */
- private static FieldMatrix<BigFraction> buildP(final int order) {
-
- final BigFraction[][] pData = new BigFraction[order - 1][order - 1];
-
- for (int i = 0; i < pData.length; ++i) {
- // build the P matrix elements from Taylor series formulas
- final BigFraction[] pI = pData[i];
- final int factor = -(i + 1);
- int aj = factor;
- for (int j = 0; j < pI.length; ++j) {
- pI[j] = new BigFraction(aj * (j + 2));
- aj *= factor;
- }
- }
-
- return new Array2DRowFieldMatrix<BigFraction>(pData, false);
-
- }
-
- }
-
- /** Serialize the instance.
- * @param oos stream where object should be written
- * @throws IOException if object cannot be written to stream
- */
- private void writeObject(ObjectOutputStream oos)
- throws IOException {
- oos.defaultWriteObject();
- oos.writeInt(coefficients.msToN.getRowDimension() + 1);
- }
-
- /** Deserialize the instance.
- * @param ois stream from which the object should be read
- * @throws ClassNotFoundException if a class in the stream cannot be found
- * @throws IOException if object cannot be read from the stream
- */
- private void readObject(ObjectInputStream ois)
- throws ClassNotFoundException, IOException {
- try {
-
- ois.defaultReadObject();
- final int order = ois.readInt();
-
- final Class<AdamsBashforthIntegrator> cl = AdamsBashforthIntegrator.class;
- final Field f = cl.getDeclaredField("coefficients");
- f.setAccessible(true);
-
- // cache the coefficients for each order, to avoid recomputing them
- synchronized(cache) {
- CachedCoefficients coeff = cache.get(order);
- if (coeff == null) {
- coeff = new CachedCoefficients(order);
- cache.put(order, coeff);
- }
- f.set(this, coeff);
- }
-
- } catch (NoSuchFieldException nsfe) {
- IOException ioe = new IOException();
- ioe.initCause(nsfe);
- throw ioe;
- } catch (IllegalAccessException iae) {
- IOException ioe = new IOException();
- ioe.initCause(iae);
- throw ioe;
- }
-
- }
-
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=786877&r1=786876&r2=786877&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Sat Jun 20 18:15:54 2009
@@ -17,25 +17,11 @@
package org.apache.commons.math.ode.nonstiff;
-import java.io.IOException;
-import java.io.ObjectInputStream;
-import java.io.ObjectOutputStream;
import java.io.Serializable;
-import java.lang.reflect.Field;
import java.util.Arrays;
-import java.util.HashMap;
-import java.util.Map;
-import org.apache.commons.math.MathRuntimeException;
-import org.apache.commons.math.fraction.BigFraction;
-import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
-import org.apache.commons.math.linear.FieldLUDecompositionImpl;
-import org.apache.commons.math.linear.FieldMatrix;
-import org.apache.commons.math.linear.Array2DRowFieldMatrix;
-import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.MatrixVisitorException;
import org.apache.commons.math.linear.RealMatrix;
-import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
@@ -177,49 +163,51 @@
*/
private static final long serialVersionUID = 3624292432281962886L;
- /** Cache for already computed coefficients. */
- private static final Map<Integer, CachedCoefficients> cache =
- new HashMap<Integer, CachedCoefficients>();
-
- /** Coefficients of the method. */
- private final transient CachedCoefficients coefficients;
-
- /** Integration step. */
- private final double step;
-
/**
* Build an Adams-Moulton integrator with the given order and step size.
* @param order order of the method (must be greater than 1: due to
* an implementation limitation the order 1 method is not supported)
- * @param step integration step size
+ * @param minStep minimal step (must be positive even for backward
+ * integration), the last step can be smaller than this
+ * @param maxStep maximal step (must be positive even for backward
+ * integration)
+ * @param scalAbsoluteTolerance allowed absolute error
+ * @param scalRelativeTolerance allowed relative error
* @exception IllegalArgumentException if order is 1 or less
*/
- public AdamsMoultonIntegrator(final int order, final double step)
+ public AdamsMoultonIntegrator(final int order,
+ final double minStep, final double maxStep,
+ final double scalAbsoluteTolerance,
+ final double scalRelativeTolerance)
throws IllegalArgumentException {
+ super("Adams-Moulton", order, order, minStep, maxStep,
+ scalAbsoluteTolerance, scalRelativeTolerance);
+ }
- super("Adams-Moulton", order);
- if (order <= 1) {
- throw MathRuntimeException.createIllegalArgumentException(
- "{0} is supported only for orders 2 or more",
- getName());
- }
-
- // cache the coefficients for each order, to avoid recomputing them
- synchronized(cache) {
- CachedCoefficients coeff = cache.get(order);
- if (coeff == null) {
- coeff = new CachedCoefficients(order);
- cache.put(order, coeff);
- }
- coefficients = coeff;
- }
-
- this.step = Math.abs(step);
-
+ /**
+ * Build an Adams-Moulton integrator with the given order and step size.
+ * @param order order of the method (must be greater than 1: due to
+ * an implementation limitation the order 1 method is not supported)
+ * @param minStep minimal step (must be positive even for backward
+ * integration), the last step can be smaller than this
+ * @param maxStep maximal step (must be positive even for backward
+ * integration)
+ * @param vecAbsoluteTolerance allowed absolute error
+ * @param vecRelativeTolerance allowed relative error
+ * @exception IllegalArgumentException if order is 1 or less
+ */
+ public AdamsMoultonIntegrator(final int order,
+ final double minStep, final double maxStep,
+ final double[] vecAbsoluteTolerance,
+ final double[] vecRelativeTolerance)
+ throws IllegalArgumentException {
+ super("Adams-Moulton", order, order, minStep, maxStep,
+ vecAbsoluteTolerance, vecRelativeTolerance);
}
/** {@inheritDoc} */
+ @Override
public double integrate(final FirstOrderDifferentialEquations equations,
final double t0, final double[] y0,
final double t, final double[] y)
@@ -235,6 +223,7 @@
if (y != y0) {
System.arraycopy(y0, 0, y, 0, n);
}
+ final double[] yDot = new double[y0.length];
final double[] yTmp = new double[y0.length];
// set up two interpolators sharing the integrator arrays
@@ -244,34 +233,27 @@
interpolatorTmp.reinitialize(yTmp, forward);
// set up integration control objects
- stepStart = t0;
- stepSize = forward ? step : -step;
for (StepHandler handler : stepHandlers) {
handler.reset();
}
CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
- // compute the first few steps using the configured starter integrator
- double stopTime = start(previousF.length, stepSize, manager, stepStart, y);
- if (Double.isNaN(previousT[0])) {
- return stopTime;
- }
- stepStart = previousT[0];
- double hNew = 0;
- // convert to Nordsieck representation
- double[] scaled = convertToNordsieckLow();
- RealMatrix nordsieck = convertToNordsieckHigh(scaled);
+ // compute the initial Nordsieck vector using the configured starter integrator
+ start(t0, y, t);
interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
interpolator.storeTime(stepStart);
+ double hNew = stepSize;
+ interpolator.rescale(hNew);
+
boolean lastStep = false;
while (!lastStep) {
// shift all data
interpolator.shift();
- hNew = forward ? step : -step;
+ double error = 0;
for (boolean loop = true; loop;) {
stepSize = hNew;
@@ -282,51 +264,60 @@
System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length);
// evaluate a first estimate of the derivative (first E in the PECE sequence)
- final double[] f0 = previousF[0];
- previousT[0] = stepEnd;
- computeDerivatives(stepEnd, yTmp, f0);
+ computeDerivatives(stepEnd, yTmp, yDot);
// update Nordsieck vector
- final RealMatrix nordsieckTmp = coefficients.msUpdate.multiply(nordsieck);
final double[] predictedScaled = new double[y0.length];
for (int j = 0; j < y0.length; ++j) {
- predictedScaled[j] = stepSize * f0[j];
+ predictedScaled[j] = stepSize * yDot[j];
}
- nordsieckTmp.walkInOptimizedOrder(new NordsieckUpdater(scaled, predictedScaled, coefficients.c1));
+ final RealMatrix nordsieckTmp =
+ transformer.updateHighOrderDerivativesPhase1(nordsieck);
+ transformer.updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
// apply correction (C in the PECE sequence)
- nordsieckTmp.walkInOptimizedOrder(new Corrector(y, predictedScaled, yTmp));
+ error = nordsieckTmp.walkInOptimizedOrder(new Corrector(y, predictedScaled, yTmp));
- // evaluate a final estimate of the derivative (second E in the PECE sequence)
- computeDerivatives(stepEnd, yTmp, f0);
+ if (error <= 1.0) {
- // update Nordsieck vector
- final double[] correctedScaled = new double[y0.length];
- for (int j = 0; j < y0.length; ++j) {
- correctedScaled[j] = stepSize * f0[j];
- }
- nordsieckTmp.walkInOptimizedOrder(new NordsieckUpdater(predictedScaled, correctedScaled, coefficients.c1));
+ // evaluate a final estimate of the derivative (second E in the PECE sequence)
+ computeDerivatives(stepEnd, yTmp, yDot);
- // discrete events handling
- interpolatorTmp.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
- interpolatorTmp.storeTime(stepStart);
- interpolatorTmp.shift();
- interpolatorTmp.storeTime(stepEnd);
- if (manager.evaluateStep(interpolatorTmp)) {
- final double dt = manager.getEventTime() - stepStart;
- if (Math.abs(dt) <= Math.ulp(stepStart)) {
- // rejecting the step would lead to a too small next step, we accept it
- loop = false;
+ // update Nordsieck vector
+ final double[] correctedScaled = new double[y0.length];
+ for (int j = 0; j < y0.length; ++j) {
+ correctedScaled[j] = stepSize * yDot[j];
+ }
+ transformer.updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
+
+ // discrete events handling
+ interpolatorTmp.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
+ interpolatorTmp.storeTime(stepStart);
+ interpolatorTmp.shift();
+ interpolatorTmp.storeTime(stepEnd);
+ if (manager.evaluateStep(interpolatorTmp)) {
+ final double dt = manager.getEventTime() - stepStart;
+ if (Math.abs(dt) <= Math.ulp(stepStart)) {
+ // rejecting the step would lead to a too small next step, we accept it
+ loop = false;
+ } else {
+ // reject the step to match exactly the next switch time
+ hNew = dt;
+ interpolator.rescale(hNew);
+ }
} else {
- // reject the step to match exactly the next switch time
- hNew = dt;
+ // accept the step
+ scaled = correctedScaled;
+ nordsieck = nordsieckTmp;
+ interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
+ loop = false;
}
+
} else {
- // accept the step
- scaled = correctedScaled;
- nordsieck = nordsieckTmp;
- interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
- loop = false;
+ // reject the step and attempt to reduce error by stepsize control
+ final double factor = computeStepGrowShrinkFactor(error);
+ hNew = filterStep(stepSize * factor, forward, false);
+ interpolator.rescale(hNew);
}
}
@@ -339,8 +330,8 @@
lastStep = manager.stop();
// provide the step data to the step handler
- interpolator.setInterpolatedTime(nextStep);
for (StepHandler handler : stepHandlers) {
+ interpolator.setInterpolatedTime(nextStep);
handler.handleStep(interpolator, lastStep);
}
stepStart = nextStep;
@@ -349,63 +340,35 @@
// some events handler has triggered changes that
// invalidate the derivatives, we need to restart from scratch
- stopTime = start(previousF.length, stepSize, manager, stepStart, y);
- if (Double.isNaN(previousT[0])) {
- return stopTime;
- }
- stepStart = previousT[0];
-
- // convert to Nordsieck representation
- scaled = convertToNordsieckLow();
- nordsieck = convertToNordsieckHigh(scaled);
+ start(stepStart, y, t);
interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
- interpolator.storeTime(stepStart);
}
+ if (! lastStep) {
+ // in some rare cases we may get here with stepSize = 0, for example
+ // when an event occurs at integration start, reducing the first step
+ // to zero; we have to reset the step to some safe non zero value
+ stepSize = filterStep(stepSize, forward, true);
+
+ // stepsize control for next step
+ final double factor = computeStepGrowShrinkFactor(error);
+ final double scaledH = stepSize * factor;
+ final double nextT = stepStart + scaledH;
+ final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
+ hNew = filterStep(scaledH, forward, nextIsLast);
+ interpolator.rescale(hNew);
+ }
+
}
- stopTime = stepStart;
+ final double stopTime = stepStart;
stepStart = Double.NaN;
stepSize = Double.NaN;
return stopTime;
}
- /** Convert the multistep representation after a restart to Nordsieck representation.
- * @return first scaled derivative
- */
- private double[] convertToNordsieckLow() {
-
- final double[] f0 = previousF[0];
- final double[] scaled = new double[f0.length];
- for (int j = 0; j < f0.length; ++j) {
- scaled[j] = stepSize * f0[j];
- }
- return scaled;
-
- }
-
- /** Convert the multistep representation after a restart to Nordsieck representation.
- * @param scaled first scaled derivative
- * @return Nordsieck matrix of the higher scaled derivatives
- */
- private RealMatrix convertToNordsieckHigh(final double[] scaled) {
-
- final double[] f0 = previousF[0];
- final double[][] multistep = new double[coefficients.msToN.getColumnDimension()][f0.length];
- for (int i = 0; i < multistep.length; ++i) {
- final double[] msI = multistep[i];
- final double[] fI = previousF[i + 1];
- for (int j = 0; j < f0.length; ++j) {
- msI[j] = stepSize * fI[j] - scaled[j];
- }
- }
-
- return coefficients.msToN.multiply(new Array2DRowRealMatrix(multistep, false));
-
- }
-
/** Corrector for current state in Adams-Moulton method.
* <p>
* This visitor implements the Taylor series formula:
@@ -414,7 +377,7 @@
* </pre>
* </p>
*/
- private static class Corrector implements RealMatrixPreservingVisitor {
+ private class Corrector implements RealMatrixPreservingVisitor {
/** Previous state. */
private final double[] previous;
@@ -422,214 +385,65 @@
/** Current scaled first derivative. */
private final double[] scaled;
- /** Placeholder where to put the recomputed current state. */
- private final double[] corrected;
+ /** Current state before correction. */
+ private final double[] before;
+
+ /** Current state after correction. */
+ private final double[] after;
/** Simple constructor.
* @param previous previous state
* @param scaled current scaled first derivative
- * @param corrected placeholder where to put the corrected current state
+ * @param state state to correct (will be overwritten after visit)
*/
- public Corrector(final double[] previous, final double[] scaled, final double[] corrected) {
+ public Corrector(final double[] previous, final double[] scaled, final double[] state) {
this.previous = previous;
- this.scaled = scaled;
- this.corrected = corrected;
+ this.scaled = scaled;
+ this.after = state;
+ this.before = state.clone();
}
/** {@inheritDoc} */
public void start(int rows, int columns,
int startRow, int endRow, int startColumn, int endColumn) {
- Arrays.fill(corrected, 0.0);
+ Arrays.fill(after, 0.0);
}
/** {@inheritDoc} */
public void visit(int row, int column, double value)
throws MatrixVisitorException {
if ((row & 0x1) == 0) {
- corrected[column] -= value;
+ after[column] -= value;
} else {
- corrected[column] += value;
- }
- }
-
- /** {@inheritDoc} */
- public double end() {
- for (int i = 0; i < corrected.length; ++i) {
- corrected[i] += previous[i] + scaled[i];
+ after[column] += value;
}
- return 0;
- }
- }
-
- /** Updater for Nordsieck vector. */
- private static class NordsieckUpdater extends DefaultRealMatrixChangingVisitor {
-
- /** Scaled first derivative at step start. */
- private final double[] start;
-
- /** Scaled first derivative at step end. */
- private final double[] end;
-
- /** Update coefficients. */
- private final double[] c1;
-
- /** Simple constructor.
- * @param start scaled first derivative at step start
- * @param end scaled first derivative at step end
- * @param c1 update coefficients
- */
- public NordsieckUpdater(final double[] start, final double[] end,
- final double[] c1) {
- this.start = start;
- this.end = end;
- this.c1 = c1;
- }
-
- /** {@inheritDoc} */
- @Override
- public double visit(int row, int column, double value)
- throws MatrixVisitorException {
- return value + c1[row] * (start[column] - end[column]);
}
- }
-
- /** Cache for already computed coefficients.
- * @param <impements>*/
- private static class CachedCoefficients implements Serializable {
-
/**
- * Serialization UID
+ * End visiting te Nordsieck vector.
+ * <p>The correction is used to control stepsize. So its amplitude is
+ * considered to be an error, which must be normalized according to
+ * error control settings. If the normalized value is greater than 1,
+ * the correction was too large and the step must be rejected.</p>
+ * @return the normalized correction, if greater than 1, the step
+ * must be rejected
*/
- private static final long serialVersionUID = -8464316300182136812L;
-
- /** Transformer between multistep and Nordsieck representations. */
- private final RealMatrix msToN;
-
- /** Update coefficients of the higher order derivatives wrt y'', y''' ... */
- private final RealMatrix msUpdate;
-
- /** Update coefficients of the higher order derivatives wrt y'. */
- private final double[] c1;
-
- /** Simple constructor.
- * @param order order of the method (must be greater than 1: due to
- * an implementation limitation the order 1 method is not supported)
- */
- public CachedCoefficients(int order) {
-
- // compute exact coefficients
- FieldMatrix<BigFraction> bigNtoMS = buildP(order);
- FieldMatrix<BigFraction> bigMStoN =
- new FieldLUDecompositionImpl<BigFraction>(bigNtoMS).getSolver().getInverse();
- BigFraction[] u = new BigFraction[order - 1];
- Arrays.fill(u, BigFraction.ONE);
- BigFraction[] bigC1 = bigMStoN.operate(u);
-
- // update coefficients are computed by combining transform from
- // Nordsieck to multistep, then shifting rows to represent step advance
- // then applying inverse transform
- BigFraction[][] shiftedP = bigNtoMS.getData();
- for (int i = shiftedP.length - 1; i > 0; --i) {
- // shift rows
- shiftedP[i] = shiftedP[i - 1];
- }
- shiftedP[0] = new BigFraction[order - 1];
- Arrays.fill(shiftedP[0], BigFraction.ZERO);
- FieldMatrix<BigFraction> bigMSupdate =
- bigMStoN.multiply(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
-
- // convert coefficients to double
- msToN = MatrixUtils.bigFractionMatrixToRealMatrix(bigMStoN);
- msUpdate = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
- c1 = new double[order - 1];
- for (int i = 0; i < order - 1; ++i) {
- c1[i] = bigC1[i].doubleValue();
- }
-
- }
-
- /** Build the P matrix transforming multistep to Nordsieck.
- * <p>
- * <p>
- * Multistep representation uses y(k), s<sub>1</sub>(k), s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)).
- * Nordsieck representation uses y(k), s<sub>1</sub>(k), s<sub>2</sub>(k) ... s<sub>n</sub>(k).
- * The two representations share their two first components y(k) and
- * s<sub>1</sub>(k). The P matrix is used to transform the remaining ones:
- * <pre>
- * [ s<sub>1</sub>(k-1) ... s<sub>1</sub>(k-(n-1)]<sup>T</sup> = s<sub>1</sub>(k) [1 ... 1]<sup>T</sup> + P [s<sub>2</sub>(k) ... s<sub>n</sub>(k)]<sup>T</sup>
- * </pre>
- * </p>
- * @param order order of the method (must be strictly positive)
- * @return P matrix
- */
- private static FieldMatrix<BigFraction> buildP(final int order) {
-
- final BigFraction[][] pData = new BigFraction[order - 1][order - 1];
+ public double end() {
- for (int i = 0; i < pData.length; ++i) {
- // build the P matrix elements from Taylor series formulas
- final BigFraction[] pI = pData[i];
- final int factor = -(i + 1);
- int aj = factor;
- for (int j = 0; j < pI.length; ++j) {
- pI[j] = new BigFraction(aj * (j + 2));
- aj *= factor;
- }
+ double error = 0;
+ for (int i = 0; i < after.length; ++i) {
+ after[i] += previous[i] + scaled[i];
+ final double yScale = Math.max(Math.abs(previous[i]), Math.abs(after[i]));
+ final double tol = (vecAbsoluteTolerance == null) ?
+ (scalAbsoluteTolerance + scalRelativeTolerance * yScale) :
+ (vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * yScale);
+ final double ratio = (after[i] - before[i]) / tol;
+ error += ratio * ratio;
}
- return new Array2DRowFieldMatrix<BigFraction>(pData, false);
+ return Math.sqrt(error / after.length);
}
-
- }
-
- /** Serialize the instance.
- * @param oos stream where object should be written
- * @throws IOException if object cannot be written to stream
- */
- private void writeObject(ObjectOutputStream oos)
- throws IOException {
- oos.defaultWriteObject();
- oos.writeInt(coefficients.msToN.getRowDimension() + 1);
- }
-
- /** Deserialize the instance.
- * @param ois stream from which the object should be read
- * @throws ClassNotFoundException if a class in the stream cannot be found
- * @throws IOException if object cannot be read from the stream
- */
- private void readObject(ObjectInputStream ois)
- throws ClassNotFoundException, IOException {
- try {
-
- ois.defaultReadObject();
- final int order = ois.readInt();
-
- final Class<AdamsMoultonIntegrator> cl = AdamsMoultonIntegrator.class;
- final Field f = cl.getDeclaredField("coefficients");
- f.setAccessible(true);
-
- // cache the coefficients for each order, to avoid recomputing them
- synchronized(cache) {
- CachedCoefficients coeff = cache.get(order);
- if (coeff == null) {
- coeff = new CachedCoefficients(order);
- cache.put(order, coeff);
- }
- f.set(this, coeff);
- }
-
- } catch (NoSuchFieldException nsfe) {
- IOException ioe = new IOException();
- ioe.initCause(nsfe);
- throw ioe;
- } catch (IllegalAccessException iae) {
- IOException ioe = new IOException();
- ioe.initCause(iae);
- throw ioe;
- }
-
}
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java?rev=786877&r1=786876&r2=786877&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java Sat Jun 20 18:15:54 2009
@@ -22,8 +22,9 @@
import java.io.ObjectOutput;
import java.util.Arrays;
-import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
+import org.apache.commons.math.linear.DefaultRealMatrixChangingVisitor;
+import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
/**
@@ -93,7 +94,7 @@
return new NordsieckStepInterpolator(this);
}
- /** Reinitialize the instance
+ /** Reinitialize the instance.
* <p>Beware that all arrays <em>must</em> be references to integrator
* arrays, in order to ensure proper update without copy.</p>
* @param y reference to the integrator array holding the state at
@@ -105,7 +106,7 @@
super.reinitialize(y, forward);
}
- /** Reinitialize the instance
+ /** Reinitialize the instance.
* <p>Beware that all arrays <em>must</em> be references to integrator
* arrays, in order to ensure proper update without copy.</p>
* @param referenceTime time at which all arrays are defined
@@ -127,6 +128,20 @@
}
+ /** Rescale the instance.
+ * <p>Since the scaled and Nordiseck arrays are shared with the caller,
+ * this method has the side effect of rescaling this arrays in the caller too.</p>
+ * @param scalingH new step size to use in the scaled and nordsieck arrays
+ */
+ public void rescale(final double scalingH) {
+ final double ratio = scalingH / this.scalingH;
+ for (int i = 0; i < scaled.length; ++i) {
+ scaled[i] *= ratio;
+ }
+ nordsieck.walkInOptimizedOrder(new Rescaler(ratio));
+ this.scalingH = scalingH;
+ }
+
/** {@inheritDoc} */
@Override
protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
@@ -151,7 +166,7 @@
* @param theta normalized interpolation abscissa within the step
*/
public StateEstimator(final double scale, final double theta) {
- this.scale = scale;
+ this.scale = scale;
lowPower = theta;
highPowers = new double[nordsieck.getRowDimension()];
double thetaN = theta;
@@ -187,6 +202,32 @@
}
+ /** Visitor rescaling the Nordsieck vector. */
+ private class Rescaler extends DefaultRealMatrixChangingVisitor {
+
+ /** Powers of the rescaling ratio. */
+ private final double[] powers;
+
+ /** Simple constructor.
+ * @param ratio rescaling ratio
+ */
+ public Rescaler(final double ratio) {
+ powers = new double[nordsieck.getRowDimension()];
+ double f = ratio;
+ for (int i = 0; i < powers.length; ++i) {
+ f *= ratio;
+ powers[i] = f;
+ }
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double visit(final int row, final int column, final double value) {
+ return value * powers[row];
+ }
+
+ }
+
/** {@inheritDoc} */
@Override
public void writeExternal(final ObjectOutput out)
Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java?rev=786877&r1=786876&r2=786877&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java Sat Jun 20 18:15:54 2009
@@ -29,7 +29,6 @@
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.EventHandler;
import org.junit.Test;
public class AdamsBashforthIntegratorTest {
@@ -37,70 +36,67 @@
@Test(expected=IntegratorException.class)
public void dimensionCheck() throws DerivativeException, IntegratorException {
TestProblem1 pb = new TestProblem1();
- new AdamsBashforthIntegrator(3, 0.01).integrate(pb,
- 0.0, new double[pb.getDimension()+10],
- 1.0, new double[pb.getDimension()+10]);
+ FirstOrderIntegrator integ =
+ new AdamsBashforthIntegrator(3, 0.0, 1.0, 1.0e-10, 1.0e-10);
+ integ.integrate(pb,
+ 0.0, new double[pb.getDimension()+10],
+ 1.0, new double[pb.getDimension()+10]);
}
- @Test
- public void decreasingSteps() throws DerivativeException, IntegratorException {
-
- TestProblemAbstract[] problems = TestProblemFactory.getProblems();
- for (int k = 0; k < problems.length; ++k) {
-
- double previousError = Double.NaN;
- for (int i = 6; i < 10; ++i) {
-
- TestProblemAbstract pb = (TestProblemAbstract) problems[k].clone();
- double step = (pb.getFinalTime() - pb.getInitialTime()) * Math.pow(2.0, -i);
-
- FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, step);
- TestProblemHandler handler = new TestProblemHandler(pb, integ);
- integ.addStepHandler(handler);
- EventHandler[] functions = pb.getEventsHandlers();
- for (int l = 0; l < functions.length; ++l) {
- integ.addEventHandler(functions[l],
- Double.POSITIVE_INFINITY, 1.0e-3 * step, 1000);
- }
- double stopTime = integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
- pb.getFinalTime(), new double[pb.getDimension()]);
- if (functions.length == 0) {
- assertEquals(pb.getFinalTime(), stopTime, 1.0e-10);
- }
-
- double error = handler.getMaximalValueError();
- if ((i > 6) && !(pb instanceof TestProblem4) && !(pb instanceof TestProblem6)) {
- assertTrue(error <= Math.abs(1.05 * previousError));
- }
- previousError = error;
+ @Test(expected=IntegratorException.class)
+ public void testMinStep() throws DerivativeException, IntegratorException {
- }
-
- }
+ TestProblem1 pb = new TestProblem1();
+ double minStep = 0.1 * (pb.getFinalTime() - pb.getInitialTime());
+ double maxStep = pb.getFinalTime() - pb.getInitialTime();
+ double[] vecAbsoluteTolerance = { 1.0e-15, 1.0e-16 };
+ double[] vecRelativeTolerance = { 1.0e-15, 1.0e-16 };
+
+ FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, minStep, maxStep,
+ vecAbsoluteTolerance,
+ vecRelativeTolerance);
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(pb,
+ pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
}
@Test
- public void smallStep() throws DerivativeException, IntegratorException {
-
- TestProblem1 pb = new TestProblem1();
- double range = pb.getFinalTime() - pb.getInitialTime();
- double step = range * 0.001;
+ public void testIncreasingTolerance()
+ throws DerivativeException, IntegratorException {
+
+ int previousCalls = Integer.MAX_VALUE;
+ for (int i = -12; i < -2; ++i) {
+ TestProblem1 pb = new TestProblem1();
+ double minStep = 0;
+ double maxStep = pb.getFinalTime() - pb.getInitialTime();
+ double scalAbsoluteTolerance = Math.pow(10.0, i);
+ double scalRelativeTolerance = 0.01 * scalAbsoluteTolerance;
+
+ FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, minStep, maxStep,
+ scalAbsoluteTolerance,
+ scalRelativeTolerance);
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(pb,
+ pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
- AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(3, step);
- integ.setStarterIntegrator(new DormandPrince853Integrator(0, range, 1.0e-12, 1.0e-12));
- TestProblemHandler handler = new TestProblemHandler(pb, integ);
- integ.addStepHandler(handler);
- integ.integrate(pb,
- pb.getInitialTime(), pb.getInitialState(),
- pb.getFinalTime(), new double[pb.getDimension()]);
+ // the 28 and 42 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
+ assertTrue(handler.getMaximalValueError() > (28.0 * scalAbsoluteTolerance));
+ assertTrue(handler.getMaximalValueError() < (42.0 * scalAbsoluteTolerance));
+ assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
+
+ int calls = pb.getCalls();
+ assertEquals(integ.getEvaluations(), calls);
+ assertTrue(calls <= previousCalls);
+ previousCalls = calls;
- assertTrue(handler.getLastError() < 2.0e-9);
- assertTrue(handler.getMaximalValueError() < 9.0e-9);
- assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
- assertEquals("Adams-Bashforth", integ.getName());
- assertTrue(integ.getEvaluations() > 1000);
- assertEquals(Integer.MAX_VALUE, integ.getMaxEvaluations());
+ }
}
@@ -109,13 +105,11 @@
TestProblem1 pb = new TestProblem1();
double range = pb.getFinalTime() - pb.getInitialTime();
- double step = range * 0.001;
- AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(3, step);
- integ.setStarterIntegrator(new DormandPrince853Integrator(0, range, 1.0e-12, 1.0e-12));
+ AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(3, 0, range, 1.0e-12, 1.0e-12);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
- integ.setMaxEvaluations(1000);
+ integ.setMaxEvaluations(650);
integ.integrate(pb,
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
@@ -123,38 +117,19 @@
}
@Test
- public void bigStep() throws DerivativeException, IntegratorException {
-
- TestProblem1 pb = new TestProblem1();
- double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
-
- FirstOrderIntegrator integ = new AdamsBashforthIntegrator(3, step);
- TestProblemHandler handler = new TestProblemHandler(pb, integ);
- integ.addStepHandler(handler);
- integ.integrate(pb,
- pb.getInitialTime(), pb.getInitialState(),
- pb.getFinalTime(), new double[pb.getDimension()]);
-
- assertTrue(handler.getLastError() > 0.06);
- assertTrue(handler.getMaximalValueError() > 0.06);
- assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
-
- }
-
- @Test
public void backward() throws DerivativeException, IntegratorException {
TestProblem5 pb = new TestProblem5();
- double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+ double range = Math.abs(pb.getFinalTime() - pb.getInitialTime());
- FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, step);
+ FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, 0, range, 1.0e-12, 1.0e-12);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- assertTrue(handler.getLastError() < 8.0e-11);
- assertTrue(handler.getMaximalValueError() < 8.0e-11);
+ assertTrue(handler.getLastError() < 1.0e-8);
+ assertTrue(handler.getMaximalValueError() < 1.0e-8);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
assertEquals("Adams-Bashforth", integ.getName());
}
@@ -162,20 +137,19 @@
@Test
public void polynomial() throws DerivativeException, IntegratorException {
TestProblem6 pb = new TestProblem6();
- double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.02;
+ double range = Math.abs(pb.getFinalTime() - pb.getInitialTime());
for (int order = 2; order < 9; ++order) {
- AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(order, step);
- integ.setStarterIntegrator(new DormandPrince853Integrator(1.0e-3 * step, 1.0e3 * step,
- 1.0e-5, 1.0e-5));
+ AdamsBashforthIntegrator integ =
+ new AdamsBashforthIntegrator(order, 1.0e-6 * range, 0.1 * range, 1.0e-10, 1.0e-10);
TestProblemHandler handler = new TestProblemHandler(pb, integ);
integ.addStepHandler(handler);
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
if (order < 5) {
- assertTrue(handler.getMaximalValueError() > 1.0e-5);
+ assertTrue(integ.getEvaluations() > 160);
} else {
- assertTrue(handler.getMaximalValueError() < 7.0e-12);
+ assertTrue(integ.getEvaluations() < 70);
}
}
@@ -187,12 +161,14 @@
IOException, ClassNotFoundException {
TestProblem6 pb = new TestProblem6();
- double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.01;
+ double range = Math.abs(pb.getFinalTime() - pb.getInitialTime());
ByteArrayOutputStream bos = new ByteArrayOutputStream();
ObjectOutputStream oos = new ObjectOutputStream(bos);
- oos.writeObject(new AdamsBashforthIntegrator(8, step));
-
+ oos.writeObject(new AdamsBashforthIntegrator(5, 0, range, 1.0e-12, 1.0e-12));
+ assertTrue(bos.size() > 2800);
+ assertTrue(bos.size() < 3000);
+
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);
FirstOrderIntegrator integ = (AdamsBashforthIntegrator) ois.readObject();
@@ -201,7 +177,7 @@
integ.addStepHandler(handler);
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- assertTrue(handler.getMaximalValueError() < 7.0e-13);
+ assertTrue(handler.getMaximalValueError() < 2.0e-11);
}