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/01 00:07:27 UTC
svn commit: r780514 - in /commons/proper/math/trunk/src:
java/org/apache/commons/math/ode/nonstiff/
test/org/apache/commons/math/ode/nonstiff/
Author: luc
Date: Sun May 31 22:07:26 2009
New Revision: 780514
URL: http://svn.apache.org/viewvc?rev=780514&view=rev
Log:
completely redesigned Adams-Bashforth and Adams-Moulton integrators
they are now provided by a single class and use the Nordsieck form
with higher derivatives at current step instead of classical form
with only first derivatives but at several steps.
The implementation is simpler for both the integrators and the step
interpolator and it allows a future enhancement with adaptive stepsize.
Added:
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java
- copied, changed from r780292, commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsIntegratorTest.java
- copied, changed from r780292, commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
Removed:
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/AdamsBashforthStepInterpolator.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/nonstiff/AdamsMoultonStepInterpolator.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/nonstiff/StepInterpolatorAbstractTest.java
Copied: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java (from r780292, 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/AdamsIntegrator.java?p2=commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java&p1=commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java&r1=780292&r2=780514&rev=780514&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/AdamsIntegrator.java Sun May 31 22:07:26 2009
@@ -17,17 +17,36 @@
package org.apache.commons.math.ode.nonstiff;
-import org.apache.commons.math.fraction.Fraction;
+import java.io.IOException;
+import java.io.ObjectInputStream;
+import java.io.ObjectOutputStream;
+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.FieldMatrix;
+import org.apache.commons.math.linear.FieldMatrixImpl;
+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.RealMatrixImpl;
+import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
+import org.apache.commons.math.linear.decomposition.FieldLUDecompositionImpl;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.MultistepIntegrator;
import org.apache.commons.math.ode.events.CombinedEventsManager;
+import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
import org.apache.commons.math.ode.sampling.StepHandler;
/**
- * This class implements explicit Adams-Bashforth integrators for Ordinary
+ * This class implements explicit Adams-Bashforth and Adams-Moulton integrators for Ordinary
* Differential Equations.
*
* <p>Adams-Bashforth (in fact due to Adams alone) methods are explicit
@@ -37,58 +56,192 @@
* steps one wants to use for computing the next value, different formulas
* are available:</p>
* <ul>
- * <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + h f<sub>n</sub></li>
- * <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + h (3f<sub>n</sub>-f<sub>n-1</sub>)/2</li>
- * <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + h (23f<sub>n</sub>-16f<sub>n-1</sub>+5f<sub>n-2</sub>)/12</li>
- * <li>k = 4: y<sub>n+1</sub> = y<sub>n</sub> + h (55f<sub>n</sub>-59f<sub>n-1</sub>+37f<sub>n-2</sub>-9f<sub>n-3)/24</sub></li>
+ * <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + h y'<sub>n</sub></li>
+ * <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + h (3y'<sub>n</sub>-y'<sub>n-1</sub>)/2</li>
+ * <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + h (23y'<sub>n</sub>-16y'<sub>n-1</sub>+5y'<sub>n-2</sub>)/12</li>
+ * <li>k = 4: y<sub>n+1</sub> = y<sub>n</sub> + h (55y'<sub>n</sub>-59y'<sub>n-1</sub>+37y'<sub>n-2</sub>-9y'<sub>n-3)/24</sub></li>
* <li>...</li>
* </ul>
*
- * <p>A k-steps Adams-Bashforth method is of order k. There is no limit to the
- * value of k.</p>
+ * <p>A k-steps Adams-Bashforth method is of order k. There is no theoretical limit to the
+ * value of k, but due to an implementation limitation k must be greater than 1.</p>
*
- * <p>These methods are explicit: f<sub>n+1</sub> is not used to compute
- * y<sub>n+1</sub>. More accurate implicit Adams methods exist: the
- * Adams-Moulton methods (which are also due to Adams alone). They are
- * provided by the {@link AdamsMoultonIntegrator AdamsMoultonIntegrator} class.</p>
+ * <p>Adams-Moulton (also due to Adams alone) methods are implicit
+ * multistep ODE solvers witch fixed stepsize. The value of state vector
+ * at step n+1 is a simple combination of the value at step n and of the
+ * derivatives at steps n+1, n, n-1 ... Depending on the number k of previous
+ * steps one wants to use for computing the next value, different formulas
+ * are available:</p>
+ * <ul>
+ * <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + h y'<sub>n+1</sub></li>
+ * <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + h (y'<sub>n+1</sub>+y'<sub>n</sub>)/2</li>
+ * <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + h (5y'<sub>n+1</sub>+8y'<sub>n</sub>-y'<sub>n-1</sub>)/12</li>
+ * <li>k = 4: y<sub>n+1</sub> = y<sub>n</sub> + h (9y'<sub>n+1</sub>+19y'<sub>n</sub>-5y'<sub>n-1</sub>+y'<sub>n-2)/24</sub></li>
+ * <li>...</li>
+ * </ul>
+ *
+ * <p>A k-steps Adams-Moulton method is of order k+1. There is no theoretical limit to the
+ * value of k, but due to an implementation limitation k must be greater than 1.</p>
+ *
+ * <h3>Implementation details</h3>
+ *
+ * <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>The definitions above use the classical representation with several previous first
+ * derivatives. Lets define
+ * <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). With these definitions,
+ * Adams-Bashforth methods can be written:
+ * <ul>
+ * <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n)</li>
+ * <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + 3/2 s<sub>1</sub>(n) + [ -1/2 ] q<sub>n</sub></li>
+ * <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + 23/12 s<sub>1</sub>(n) + [ -16/12 5/12 ] q<sub>n</sub></li>
+ * <li>k = 4: y<sub>n+1</sub> = y<sub>n</sub> + 55/24 s<sub>1</sub>(n) + [ -59/24 37/24 -9/24 ] q<sub>n</sub></li>
+ * <li>...</li>
+ * </ul>and Adams-Moulton methods can be written:
+ * <ul>
+ * <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + s<sub>1</sub>(n+1)</li>
+ * <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + 1/2 s<sub>1</sub>(n+1) + [ 1/2 ] q<sub>n+1</sub></li>
+ * <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + 5/12 s<sub>1</sub>(n+1) + [ 8/12 -1/12 ] q<sub>n+1</sub></li>
+ * <li>k = 4: y<sub>n+1</sub> = y<sub>n</sub> + 9/24 s<sub>1</sub>(n+1) + [ 19/24 -5/24 1/24 ] q<sub>n+1</sub></li>
+ * <li>...</li>
+ * </ul></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 (q<sub>n</sub> for Adams-Bashforth or q<sub>n+1</sub> for Adams-Moulton)
+ * and Nordsieck vector
+ * <pre>
+ * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<supT</sup>
+ * </pre>
+ * (here again we omit the k index in the notation for clarity). 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>This class implements the Adams-Bashforth and Adams-Moulton method using the Nordsieck vector
+ * (i.e. y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) rather than the classical representation.
+ * Using the Nordsieck vector has several advantages:
+ * <ul>
+ * <li>it leverages Adams-Bashforth and Adams-Moulton methods as in this representation
+ * they share most of their coefficients and most of their implementation,</li>
+ * <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>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>
+ * If the method is an Adams-Moulton method, the following additional correction is performed:
+ * <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) - S<sub>1</sub>(n+1)) P<sup>-1</sup> u</li>
+ * </ul>
+ * where the upper case Y<sub>n+1</sub>, S<sub>n+1</sub> and R<sub>n+1</sub> represent the
+ * corrected states whereas the lower case y<sub>n+1</sub>, s<sub>n+1</sub> and r<sub>n+1</sub>
+ * represent the predicted states.</p>
+ *
+ * <p>The P<sup>-1</sup>u vector and the P<sup>-1</sup> A P matrix do not depend on the state,
+ * they are precomputed once for all.</p>
*
- * @see AdamsMoultonIntegrator
* @version $Revision$ $Date$
* @since 2.0
*/
-public class AdamsBashforthIntegrator extends MultistepIntegrator {
+public class AdamsIntegrator extends MultistepIntegrator {
/** Serializable version identifier. */
- private static final long serialVersionUID = 1676381657635800870L;
+ private static final long serialVersionUID = -5893911062100008922L;
- /** Integrator method name. */
- private static final String METHOD_NAME = "Adams-Bashforth";
+ /** Cache for already computed coefficients. */
+ private static final Map<Integer, CachedCoefficients> cache =
+ new HashMap<Integer, CachedCoefficients>();
- /** Coefficients for the current method. */
- private final double[] coeffs;
+ /** No correction integrator method name. */
+ private static final String NO_CORRECTION_METHOD_NAME = "Adams-Bashforth";
+
+ /** Correction integrator method name. */
+ private static final String CORRECTION_METHOD_NAME = "Adams-Moulton";
+
+ /** Correction indicator (to choose between Adams-Bashforth and Adams-Moulton). */
+ private final boolean withCorrection;
+
+ /** Coefficients of the method. */
+ private final transient CachedCoefficients coefficients;
/** Integration step. */
private final double step;
/**
- * Build an Adams-Bashforth integrator with the given order and step size.
- * @param order order of the method (must be strictly positive)
+ * Build an Adams-Bashforth or 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 withCorrection if true apply Adams-Moulton correction at end of
+ * step, otherwise use only Adams-Bashforth prediction
* @param step integration step size
+ * @exception IllegalArgumentException if order is 1 or less
*/
- public AdamsBashforthIntegrator(final int order, final double step) {
-
- super(METHOD_NAME, order, new AdamsBashforthStepInterpolator());
+ public AdamsIntegrator(final int order, final boolean withCorrection,
+ final double step)
+ throws IllegalArgumentException {
+
+ super(withCorrection ? CORRECTION_METHOD_NAME : NO_CORRECTION_METHOD_NAME,
+ order, new NordsieckStepInterpolator());
+ if (order <= 1) {
+ throw MathRuntimeException.createIllegalArgumentException(
+ "{0} is supported only for orders 2 or more",
+ getName());
+ }
+ this.withCorrection = withCorrection;
- // compute the integration coefficients
- int[][] bdArray = computeBackwardDifferencesArray(order);
- Fraction[] gamma = computeGammaArray(order);
- coeffs = new double[order];
- for (int i = 0; i < order; ++i) {
- Fraction f = Fraction.ZERO;
- for (int j = i; j < order; ++j) {
- f = f.add(gamma[j].multiply(new Fraction(bdArray[j][i], 1)));
+ // 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);
}
- coeffs[i] = f.doubleValue();
+ coefficients = coeff;
}
this.step = Math.abs(step);
@@ -96,23 +249,25 @@
}
/** {@inheritDoc} */
- public double integrate(FirstOrderDifferentialEquations equations,
- double t0, double[] y0, double t, double[] y)
+ public double integrate(final FirstOrderDifferentialEquations equations,
+ final double t0, final double[] y0,
+ final double t, final double[] y)
throws DerivativeException, IntegratorException {
+ final int n = y0.length;
sanityChecks(equations, t0, y0, t, y);
final boolean forward = (t > t0);
// initialize working arrays
if (y != y0) {
- System.arraycopy(y0, 0, y, 0, y0.length);
+ System.arraycopy(y0, 0, y, 0, n);
}
final double[] yTmp = new double[y0.length];
// set up an interpolator sharing the integrator arrays
- final AdamsBashforthStepInterpolator interpolator =
- (AdamsBashforthStepInterpolator) prototype.copy();
- interpolator.reinitialize(yTmp, previousT, previousF, forward);
+ final NordsieckStepInterpolator interpolator =
+ (NordsieckStepInterpolator) prototype.copy();
+ interpolator.reinitialize(yTmp, forward);
// set up integration control objects
stepStart = t0;
@@ -129,6 +284,12 @@
return stopTime;
}
stepStart = previousT[0];
+ System.arraycopy(y, 0, yTmp, 0, n);
+
+ // convert to Nordsieck representation
+ double[] scaled = convertToNordsieckLow();
+ RealMatrix nordsieck = convertToNordsieckHigh(scaled);
+ interpolator.reinitialize(stepSize, scaled, nordsieck);
interpolator.storeTime(stepStart);
boolean lastStep = false;
@@ -137,29 +298,39 @@
// shift all data
interpolator.shift();
- // estimate the state at the end of the step
- for (int j = 0; j < y0.length; ++j) {
- double sum = 0;
- for (int l = 0; l < coeffs.length; ++l) {
- sum += coeffs[l] * previousF[l][j];
+ if (withCorrection) {
+
+ // evaluate derivative at predicted state
+ final double stepEnd = stepStart + stepSize;
+ final double[] f0 = previousF[0];
+ previousT[0] = stepEnd;
+ equations.computeDerivatives(stepEnd, interpolator.getInterpolatedState(), f0);
+
+ // update Nordsieck vector
+ 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];
}
- yTmp[j] = y[j] + stepSize * sum;
+ nordsieck.walkInOptimizedOrder(new NordsieckUpdater(scaled, end, coefficients.c1));
+ scaled = end;
+
+ // update interpolator
+ nordsieck.walkInOptimizedOrder(new Corrector(y, scaled, yTmp));
+ interpolator.reinitialize(stepSize, scaled, nordsieck);
+
}
// discrete events handling
interpolator.storeTime(stepStart + stepSize);
- final boolean truncated;
if (manager.evaluateStep(interpolator)) {
- truncated = true;
- interpolator.truncateStep(manager.getEventTime());
- } else {
- truncated = false;
+ stepSize = manager.getEventTime() - stepStart;
}
// the step has been accepted (may have been truncated)
- final double nextStep = interpolator.getCurrentTime();
+ final double nextStep = stepStart + stepSize;
interpolator.setInterpolatedTime(nextStep);
- System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
+ System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, n);
manager.stepAccepted(nextStep, y);
lastStep = manager.stop();
@@ -183,25 +354,31 @@
}
stepStart = previousT[0];
- } else {
+ // convert to Nordsieck representation
+ scaled = convertToNordsieckLow();
+ nordsieck = convertToNordsieckHigh(scaled);
- if (truncated) {
- // the step has been truncated, we need to adjust the previous steps
- for (int i = 1; i < previousF.length; ++i) {
- previousT[i] = stepStart - i * stepSize;
- interpolator.setInterpolatedTime(previousT[i]);
- System.arraycopy(interpolator.getInterpolatedState(), 0,
- previousF[i], 0, y0.length);
- }
- } else {
- rotatePreviousSteps();
- }
+ } else {
// evaluate differential equations for next step
+ final double[] f0 = previousF[0];
previousT[0] = stepStart;
- equations.computeDerivatives(stepStart, y, previousF[0]);
+ equations.computeDerivatives(stepStart, y, f0);
+ if (!withCorrection) {
+ 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;
}
+
+ System.arraycopy(y, 0, yTmp, 0, n);
+ interpolator.reinitialize(stepSize, scaled, nordsieck);
+
}
}
@@ -213,74 +390,253 @@
}
- /** Get the coefficients of the method.
- * <p>The coefficients are the c<sub>i</sub> terms in the following formula:</p>
- * <pre>
- * y<sub>n+1</sub> = y<sub>n</sub> + h × ∑<sub>i=0</sub><sup>i=k-1</sup> c<sub>i</sub>f<sub>n-i</sub></li>
- * </pre>
- * @return a copy of the coefficients of the method
+ /** 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
*/
- public double[] getCoeffs() {
- return coeffs.clone();
+ 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 RealMatrixImpl(multistep, false));
+
}
- /** Compute the backward differences coefficients array.
- * <p>This is quite similar to the Pascal triangle containing the
- * binomial coefficients, except for an additional (-1)<sup>i</sup> sign.
- * We use a straightforward approach here, since we don't expect this to
- * be run too many times with too high k. It is based on the recurrence
- * relations:</p>
+ /** Corrector for current state in Adams-Moulton method.
+ * <p>
+ * This visitor implements the Taylor series formula:
* <pre>
- * ∇<sup>0</sup> f<sub>n</sub> = f<sub>n</sub>
- * ∇<sup>i+1</sup> f<sub>n</sub> = ∇<sup>i</sup>f<sub>n</sub> - ∇<sup>i</sup>f<sub>n-1</sub>
+ * 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>
* </pre>
- * @param order order of the integration method
- * @return the coefficients array for backward differences
+ * </p>
*/
- static int[][] computeBackwardDifferencesArray(final int order) {
+ private static class Corrector implements RealMatrixPreservingVisitor {
+
+ /** Previous state. */
+ private final double[] previous;
+
+ /** Current scaled first derivative. */
+ private final double[] scaled;
- // create the array
- int[][] bdArray = new int[order][];
+ /** Placeholder where to put the recomputed current state. */
+ private final double[] corrected;
+
+ /** Simple constructor.
+ * @param previous previous state
+ * @param scaled current scaled first derivative
+ * @param corrected placeholder where to put the corrected current state
+ */
+ public Corrector(final double[] previous, final double[] scaled, final double[] corrected) {
+ this.previous = previous;
+ this.scaled = scaled;
+ this.corrected = corrected;
+ }
- // recurrence initialization
- bdArray[0] = new int[] { 1 };
+ /** {@inheritDoc} */
+ public void start(int rows, int columns,
+ int startRow, int endRow, int startColumn, int endColumn) {
+ Arrays.fill(corrected, 0.0);
+ }
- // fill up array using recurrence relation
- for (int i = 1; i < order; ++i) {
- bdArray[i] = new int[i + 1];
- bdArray[i][0] = 1;
- for (int j = 0; j < i - 1; ++j) {
- bdArray[i][j + 1] = bdArray[i - 1][j + 1] - bdArray[i - 1][j];
+ /** {@inheritDoc} */
+ public void visit(int row, int column, double value)
+ throws MatrixVisitorException {
+ if ((row & 0x1) == 0) {
+ corrected[column] -= value;
+ } else {
+ corrected[column] += value;
}
- bdArray[i][i] = -bdArray[i - 1][i - 1];
}
- return bdArray;
+ /** {@inheritDoc} */
+ public double end() {
+ for (int i = 0; i < corrected.length; ++i) {
+ corrected[i] += previous[i] + scaled[i];
+ }
+ 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]);
+ }
}
- /** Compute the gamma coefficients.
- * @param order order of the integration method
- * @return gamma coefficients array
- */
- static Fraction[] computeGammaArray(final int order) {
+ /** Cache for already computed coefficients. */
+ private static class CachedCoefficients {
- // create the array
- Fraction[] gammaArray = new Fraction[order];
+ /** Transformer between multistep and Nordsieck representations. */
+ private final RealMatrix msToN;
- // recurrence initialization
- gammaArray[0] = Fraction.ONE;
+ /** Update coefficients of the higher order derivatives wrt y'', y''' ... */
+ private final RealMatrix msUpdate;
- // fill up array using recurrence relation
- for (int i = 1; i < order; ++i) {
- Fraction gamma = Fraction.ONE;
- for (int j = 1; j <= i; ++j) {
- gamma = gamma.subtract(gammaArray[i - j].multiply(new Fraction(1, j + 1)));
+ /** 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];
}
- gammaArray[i] = gamma;
+ shiftedP[0] = new BigFraction[order - 1];
+ Arrays.fill(shiftedP[0], BigFraction.ZERO);
+ FieldMatrix<BigFraction> bigMSupdate =
+ bigMStoN.multiply(new FieldMatrixImpl<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 FieldMatrixImpl<BigFraction>(pData, false);
+
}
- return gammaArray;
+ }
+
+ /** 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<AdamsIntegrator> cl = AdamsIntegrator.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) {
+ throw new IOException(nsfe);
+ } catch (IllegalAccessException iae) {
+ throw new IOException(iae);
+ }
}
Copied: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsIntegratorTest.java (from r780292, 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/AdamsIntegratorTest.java?p2=commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsIntegratorTest.java&p1=commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java&r1=780292&r2=780514&rev=780514&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/AdamsIntegratorTest.java Sun May 31 22:07:26 2009
@@ -17,208 +17,298 @@
package org.apache.commons.math.ode.nonstiff;
-import junit.framework.Test;
-import junit.framework.TestCase;
-import junit.framework.TestSuite;
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertTrue;
+
+import java.io.ByteArrayInputStream;
+import java.io.ByteArrayOutputStream;
+import java.io.IOException;
+import java.io.ObjectInputStream;
+import java.io.ObjectOutputStream;
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.MultistepIntegrator;
import org.apache.commons.math.ode.events.EventHandler;
+import org.junit.Test;
-public class AdamsBashforthIntegratorTest
- extends TestCase {
+public class AdamsIntegratorTest {
- public AdamsBashforthIntegratorTest(String name) {
- super(name);
- }
-
- public void testCoefficients() {
-
- double[] coeffs1 = new AdamsBashforthIntegrator(1, 0.01).getCoeffs();
- assertEquals(1, coeffs1.length);
- assertEquals(1.0, coeffs1[0], 1.0e-16);
-
- double[] coeffs2 = new AdamsBashforthIntegrator(2, 0.01).getCoeffs();
- assertEquals(2, coeffs2.length);
- assertEquals( 3.0 / 2.0, coeffs2[0], 1.0e-16);
- assertEquals(-1.0 / 2.0, coeffs2[1], 1.0e-16);
-
- double[] coeffs3 = new AdamsBashforthIntegrator(3, 0.01).getCoeffs();
- assertEquals(3, coeffs3.length);
- assertEquals( 23.0 / 12.0, coeffs3[0], 1.0e-16);
- assertEquals(-16.0 / 12.0, coeffs3[1], 1.0e-16);
- assertEquals( 5.0 / 12.0, coeffs3[2], 1.0e-16);
-
- double[] coeffs4 = new AdamsBashforthIntegrator(4, 0.01).getCoeffs();
- assertEquals(4, coeffs4.length);
- assertEquals( 55.0 / 24.0, coeffs4[0], 1.0e-16);
- assertEquals(-59.0 / 24.0, coeffs4[1], 1.0e-16);
- assertEquals( 37.0 / 24.0, coeffs4[2], 1.0e-16);
- assertEquals( -9.0 / 24.0, coeffs4[3], 1.0e-16);
-
- double[] coeffs5 = new AdamsBashforthIntegrator(5, 0.01).getCoeffs();
- assertEquals(5, coeffs5.length);
- assertEquals( 1901.0 / 720.0, coeffs5[0], 1.0e-16);
- assertEquals(-2774.0 / 720.0, coeffs5[1], 1.0e-16);
- assertEquals( 2616.0 / 720.0, coeffs5[2], 1.0e-16);
- assertEquals(-1274.0 / 720.0, coeffs5[3], 1.0e-16);
- assertEquals( 251.0 / 720.0, coeffs5[4], 1.0e-16);
-
- double[] coeffs6 = new AdamsBashforthIntegrator(6, 0.01).getCoeffs();
- assertEquals(6, coeffs6.length);
- assertEquals( 4277.0 / 1440.0, coeffs6[0], 1.0e-16);
- assertEquals(-7923.0 / 1440.0, coeffs6[1], 1.0e-16);
- assertEquals( 9982.0 / 1440.0, coeffs6[2], 1.0e-16);
- assertEquals(-7298.0 / 1440.0, coeffs6[3], 1.0e-16);
- assertEquals( 2877.0 / 1440.0, coeffs6[4], 1.0e-16);
- assertEquals( -475.0 / 1440.0, coeffs6[5], 1.0e-16);
-
- double[] coeffs7 = new AdamsBashforthIntegrator(7, 0.01).getCoeffs();
- assertEquals(7, coeffs7.length);
- assertEquals( 198721.0 / 60480.0, coeffs7[0], 1.0e-16);
- assertEquals(-447288.0 / 60480.0, coeffs7[1], 1.0e-16);
- assertEquals( 705549.0 / 60480.0, coeffs7[2], 1.0e-16);
- assertEquals(-688256.0 / 60480.0, coeffs7[3], 1.0e-16);
- assertEquals( 407139.0 / 60480.0, coeffs7[4], 1.0e-16);
- assertEquals(-134472.0 / 60480.0, coeffs7[5], 1.0e-16);
- assertEquals( 19087.0 / 60480.0, coeffs7[6], 1.0e-16);
-
- double[] coeffs8 = new AdamsBashforthIntegrator(8, 0.01).getCoeffs();
- assertEquals(8, coeffs8.length);
- assertEquals( 434241.0 / 120960.0, coeffs8[0], 1.0e-16);
- assertEquals(-1152169.0 / 120960.0, coeffs8[1], 1.0e-16);
- assertEquals( 2183877.0 / 120960.0, coeffs8[2], 1.0e-16);
- assertEquals(-2664477.0 / 120960.0, coeffs8[3], 1.0e-16);
- assertEquals( 2102243.0 / 120960.0, coeffs8[4], 1.0e-16);
- assertEquals(-1041723.0 / 120960.0, coeffs8[5], 1.0e-16);
- assertEquals( 295767.0 / 120960.0, coeffs8[6], 1.0e-16);
- assertEquals( -36799.0 / 120960.0, coeffs8[7], 1.0e-16);
-
- double[] coeffs9 = new AdamsBashforthIntegrator(9, 0.01).getCoeffs();
- assertEquals(9, coeffs9.length);
- assertEquals( 14097247.0 / 3628800.0, coeffs9[0], 1.0e-16);
- assertEquals( -43125206.0 / 3628800.0, coeffs9[1], 1.0e-16);
- assertEquals( 95476786.0 / 3628800.0, coeffs9[2], 1.0e-16);
- assertEquals(-139855262.0 / 3628800.0, coeffs9[3], 1.0e-16);
- assertEquals( 137968480.0 / 3628800.0, coeffs9[4], 1.0e-16);
- assertEquals( -91172642.0 / 3628800.0, coeffs9[5], 1.0e-16);
- assertEquals( 38833486.0 / 3628800.0, coeffs9[6], 1.0e-16);
- assertEquals( -9664106.0 / 3628800.0, coeffs9[7], 1.0e-16);
- assertEquals( 1070017.0 / 3628800.0, coeffs9[8], 1.0e-16);
-
- }
-
- public void testDimensionCheck() {
- try {
- TestProblem1 pb = new TestProblem1();
- new AdamsBashforthIntegrator(3, 0.01).integrate(pb,
+ @Test(expected=IntegratorException.class)
+ public void dimensionCheckBashforth() throws DerivativeException, IntegratorException {
+ TestProblem1 pb = new TestProblem1();
+ new AdamsIntegrator(3, false, 0.01).integrate(pb,
0.0, new double[pb.getDimension()+10],
1.0, new double[pb.getDimension()+10]);
- fail("an exception should have been thrown");
- } catch(DerivativeException de) {
- fail("wrong exception caught");
- } catch(IntegratorException ie) {
}
- }
-
- public void testDecreasingSteps()
- throws DerivativeException, IntegratorException {
- TestProblemAbstract[] problems = TestProblemFactory.getProblems();
- for (int k = 0; k < problems.length; ++k) {
+ @Test
+ public void decreasingStepsBashforth() 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 AdamsIntegrator(5, false, 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;
- 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);
+ @Test
+ public void smallStepBashforth() throws DerivativeException, IntegratorException {
+
+ TestProblem1 pb = new TestProblem1();
+ double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+
+ FirstOrderIntegrator integ = new AdamsIntegrator(3, false, 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-6 * 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);
- }
+ integ.integrate(pb,
+ pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
+
+ assertTrue(handler.getLastError() < 2.0e-9);
+ assertTrue(handler.getMaximalValueError() < 9.0e-9);
+ assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
+ assertEquals("Adams-Bashforth", integ.getName());
+
+ }
+
+ @Test
+ public void bigStepBashforth() throws DerivativeException, IntegratorException {
+
+ TestProblem1 pb = new TestProblem1();
+ double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
+
+ FirstOrderIntegrator integ = new AdamsIntegrator(3, false, 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-14);
+
+ }
+
+ @Test
+ public void backwardBashforth() throws DerivativeException, IntegratorException {
+
+ TestProblem5 pb = new TestProblem5();
+ double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+
+ FirstOrderIntegrator integ = new AdamsIntegrator(5, false, step);
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
- double error = handler.getMaximalValueError();
- if (i > 6) {
- assertTrue(error < Math.abs(previousError));
+ assertTrue(handler.getLastError() < 8.0e-11);
+ assertTrue(handler.getMaximalValueError() < 8.0e-11);
+ assertEquals(0, handler.getMaximalTimeError(), 1.0e-15);
+ assertEquals("Adams-Bashforth", integ.getName());
+ }
+
+ @Test
+ public void polynomialBashforth() throws DerivativeException, IntegratorException {
+ TestProblem6 pb = new TestProblem6();
+ double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.02;
+
+ for (int order = 2; order < 9; ++order) {
+ MultistepIntegrator integ = new AdamsIntegrator(order, false, step);
+ integ.setStarterIntegrator(new DormandPrince853Integrator(1.0e-3 * step, 1.0e3 * step,
+ 1.0e-5, 1.0e-5));
+ 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);
+ } else {
+ assertTrue(handler.getMaximalValueError() < 7.0e-12);
+ }
}
- previousError = error;
- }
+ }
+
+ @Test
+ public void serializationBashforth()
+ throws IntegratorException, DerivativeException,
+ IOException, ClassNotFoundException {
+
+ TestProblem6 pb = new TestProblem6();
+ double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.01;
+
+ ByteArrayOutputStream bos = new ByteArrayOutputStream();
+ ObjectOutputStream oos = new ObjectOutputStream(bos);
+ oos.writeObject(new AdamsIntegrator(8, false, step));
+ assertTrue(bos.size() > 3000);
+ assertTrue(bos.size() < 3100);
+
+ ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
+ ObjectInputStream ois = new ObjectInputStream(bis);
+ FirstOrderIntegrator integ = (AdamsIntegrator) ois.readObject();
+ assertEquals("Adams-Bashforth", integ.getName());
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
+ assertTrue(handler.getMaximalValueError() < 7.0e-13);
+
+ }
+ @Test(expected=IntegratorException.class)
+ public void dimensionCheckMoulton()
+ throws DerivativeException, IntegratorException {
+ TestProblem1 pb = new TestProblem1();
+ new AdamsIntegrator(3, true, 0.01).integrate(pb,
+ 0.0, new double[pb.getDimension()+10],
+ 1.0, new double[pb.getDimension()+10]);
}
- }
+ @Test
+ public void decreasingStepsMoulton()
+ 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 AdamsIntegrator(5, true, 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;
+
+ }
+
+ }
+
+ }
- public void testSmallStep()
- throws DerivativeException, IntegratorException {
+ @Test
+ public void smallStepMoulton()
+ throws DerivativeException, IntegratorException {
- TestProblem1 pb = new TestProblem1();
- double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+ TestProblem1 pb = new TestProblem1();
+ double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
- 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()]);
+ FirstOrderIntegrator integ = new AdamsIntegrator(3, true, 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() < 1.0e-14);
+ assertTrue(handler.getMaximalValueError() < 2.0e-17);
+ assertEquals(0, handler.getMaximalTimeError(), 1.0e-15);
+ assertEquals("Adams-Moulton", integ.getName());
- assertTrue(handler.getLastError() < 2.0e-9);
- assertTrue(handler.getMaximalValueError() < 3.0e-8);
- assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
- assertEquals("Adams-Bashforth", integ.getName());
+ }
- }
+ @Test
+ public void bigStepMoulton()
+ throws DerivativeException, IntegratorException {
- public void testBigStep()
- throws DerivativeException, IntegratorException {
+ TestProblem1 pb = new TestProblem1();
+ double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
- TestProblem1 pb = new TestProblem1();
- double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
+ FirstOrderIntegrator integ = new AdamsIntegrator(3, true, step);
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(pb,
+ pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
- 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.getMaximalValueError() > 6.0e-6);
- assertTrue(handler.getLastError() > 0.05);
- assertTrue(handler.getMaximalValueError() > 0.1);
- assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
+ }
- }
+ @Test
+ public void backwardMoulton()
+ throws DerivativeException, IntegratorException {
- public void testBackward()
- throws DerivativeException, IntegratorException {
+ TestProblem5 pb = new TestProblem5();
+ double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001;
- TestProblem5 pb = new TestProblem5();
- double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+ FirstOrderIntegrator integ = new AdamsIntegrator(5, true, step);
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
- FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, 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() < 1.0e-15);
+ assertTrue(handler.getMaximalValueError() < 3.0e-16);
+ assertEquals(0, handler.getMaximalTimeError(), 1.0e-15);
+ assertEquals("Adams-Moulton", integ.getName());
+ }
- assertTrue(handler.getLastError() < 8.0e-11);
- assertTrue(handler.getMaximalValueError() < 8.0e-11);
- assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
- assertEquals("Adams-Bashforth", integ.getName());
- }
+ @Test
+ public void polynomialMoulton()
+ throws DerivativeException, IntegratorException {
+ TestProblem6 pb = new TestProblem6();
+ double step = Math.abs(pb.getFinalTime() - pb.getInitialTime()) * 0.02;
+
+ for (int order = 2; order < 9; ++order) {
+ MultistepIntegrator integ = new AdamsIntegrator(order, true, step);
+ integ.setStarterIntegrator(new DormandPrince853Integrator(1.0e-3 * step, 1.0e3 * step,
+ 1.0e-5, 1.0e-5));
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
+ assertTrue(handler.getMaximalValueError() < 2.0e-13);
+ }
- public static Test suite() {
- return new TestSuite(AdamsBashforthIntegratorTest.class);
- }
+ }
}