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 2011/04/29 19:15:59 UTC
svn commit: r1097891 - in /commons/proper/math/trunk/src:
main/java/org/apache/commons/math/exception/util/
main/java/org/apache/commons/math/ode/
main/java/org/apache/commons/math/ode/nonstiff/
main/resources/META-INF/localization/ site/xdoc/ test/jav...
Author: luc
Date: Fri Apr 29 17:15:59 2011
New Revision: 1097891
URL: http://svn.apache.org/viewvc?rev=1097891&view=rev
Log:
Fixed initialization of multistep ODE integrators. Relying on the interpolation model of the starter integrator inside only one step was wrong. The model may have a too low order to compute high degrees derivatives in the Nordsieck vector. Now we use several steps and use only grid points instead of interpolated points.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.java
commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java Fri Apr 29 17:15:59 2011
@@ -128,7 +128,7 @@ public enum LocalizedFormats implements
DIMENSION("dimension ({0})"), /* keep */
INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE("sample contains {0} observed points, at least {1} are required"),
INSUFFICIENT_ROWS_AND_COLUMNS("insufficient data: only {0} rows and {1} columns."),
- INTEGRATION_METHOD_NEEDS_AT_LEAST_ONE_PREVIOUS_POINT("{0} method needs at least one previous point"),
+ INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS("{0} method needs at least two previous points"),
INTERNAL_ERROR("internal error, please fill a bug report at {0}"),
INVALID_BRACKETING_PARAMETERS("invalid bracketing parameters: lower bound={0}, initial={1}, upper bound={2}"),
INVALID_INTERVAL_INITIAL_VALUE_PARAMETERS("invalid interval, initial value parameters: lower={0}, initial={1}, upper={2}"),
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MultistepIntegrator.java Fri Apr 29 17:15:59 2011
@@ -21,7 +21,6 @@ import org.apache.commons.math.MathRunti
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
-import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math.ode.sampling.StepHandler;
@@ -37,7 +36,7 @@ import org.apache.commons.math.util.Fast
* 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
+ * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
* </pre></p>
* <p>Rather than storing several previous steps separately, this implementation uses
* the Nordsieck vector with higher degrees scaled derivatives all taken at the same
@@ -48,7 +47,7 @@ import org.apache.commons.math.util.Fast
* (we omit the k index in the notation for clarity)</p>
* <p>
* Multistep integrators with Nordsieck representation are highly sensitive to
- * large step changes because when the step is multiplied by a factor a, the
+ * large step changes because when the step is multiplied by factor a, the
* k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
* and the last components are the least accurate ones. The default max growth
* factor is therefore set to a quite low value: 2<sup>1/order</sup>.
@@ -65,7 +64,7 @@ public abstract class MultistepIntegrato
protected double[] scaled;
/** Nordsieck matrix of the higher scaled derivatives.
- * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y(k))</p>
+ * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
*/
protected Array2DRowRealMatrix nordsieck;
@@ -114,9 +113,9 @@ public abstract class MultistepIntegrato
super(name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);
- if (nSteps <= 0) {
+ if (nSteps <= 1) {
throw MathRuntimeException.createIllegalArgumentException(
- LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_ONE_PREVIOUS_POINT,
+ LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
name);
}
@@ -219,17 +218,14 @@ public abstract class MultistepIntegrato
starter.clearStepHandlers();
// set up one specific step handler to extract initial Nordsieck vector
- starter.addStepHandler(new NordsieckInitializer(y0.length));
+ starter.addStepHandler(new NordsieckInitializer(nSteps, y0.length));
// start integration, expecting a InitializationCompletedMarkerException
try {
starter.integrate(new CountingDifferentialEquations(y0.length),
t0, y0, t, new double[y0.length]);
- } catch (MathUserException mue) {
- if (!(mue instanceof InitializationCompletedMarkerException)) {
- // this is not the expected nominal interruption of the start integrator
- throw mue;
- }
+ } catch (InitializationCompletedMarkerException icme) {
+ // this is the expected nominal interruption of the start integrator
}
// remove the specific step handler
@@ -238,13 +234,16 @@ public abstract class MultistepIntegrato
}
/** 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 scaled derivatives at step start
- */
- protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
- final double[][] multistep);
+ * @param h step size to use for scaling
+ * @param t first steps times
+ * @param y first steps states
+ * @param yDot first steps derivatives
+ * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
+ * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
+ */
+ protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
+ final double[][] y,
+ final double[][] yDot);
/** Get the minimal reduction factor for stepsize control.
* @return minimal reduction factor
@@ -299,57 +298,88 @@ public abstract class MultistepIntegrato
/** Transformer used to convert the first step to Nordsieck representation. */
public static interface NordsieckTransformer {
/** 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
+ * @param h step size to use for scaling
+ * @param t first steps times
+ * @param y first steps states
+ * @param yDot first steps derivatives
+ * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
+ * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
*/
- RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);
+ Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
+ final double[][] y,
+ final double[][] yDot);
}
/** Specialized step handler storing the first step. */
private class NordsieckInitializer implements StepHandler {
- /** Problem dimension. */
- private final int n;
+ /** Steps counter. */
+ int count;
+
+ /** First steps times. */
+ final double[] t;
+
+ /** First steps states. */
+ final double[][] y;
+
+ /** First steps derivatives. */
+ final double[][] yDot;
/** Simple constructor.
+ * @param nSteps number of steps of the multistep method (excluding the one being computed)
* @param n problem dimension
*/
- public NordsieckInitializer(final int n) {
- this.n = n;
+ public NordsieckInitializer(final int nSteps, final int n) {
+ this.count = 0;
+ this.t = new double[nSteps];
+ this.y = new double[nSteps][n];
+ this.yDot = new double[nSteps][n];
}
/** {@inheritDoc} */
- public void handleStep(StepInterpolator interpolator, boolean isLast)
- throws MathUserException {
+ public void handleStep(StepInterpolator interpolator, boolean isLast) {
final double prev = interpolator.getPreviousTime();
final double curr = interpolator.getCurrentTime();
- stepStart = prev;
- stepSize = (curr - prev) / (nSteps + 1);
- // compute the first scaled derivative
- interpolator.setInterpolatedTime(prev);
- scaled = interpolator.getInterpolatedDerivatives().clone();
- for (int j = 0; j < n; ++j) {
- scaled[j] *= stepSize;
+ if (count == 0) {
+ // first step, we need to store also the beginning of the step
+ interpolator.setInterpolatedTime(prev);
+ t[0] = prev;
+ System.arraycopy(interpolator.getInterpolatedState(), 0,
+ y[0], 0, y[0].length);
+ System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
+ yDot[0], 0, yDot[0].length);
}
- // compute the high order scaled derivatives
- final double[][] multistep = new double[nSteps][];
- for (int i = 1; i <= nSteps; ++i) {
- interpolator.setInterpolatedTime(prev + stepSize * i);
- final double[] msI = interpolator.getInterpolatedDerivatives().clone();
- for (int j = 0; j < n; ++j) {
- msI[j] *= stepSize;
+ // store the end of the step
+ ++count;
+ interpolator.setInterpolatedTime(curr);
+ t[count] = curr;
+ System.arraycopy(interpolator.getInterpolatedState(), 0,
+ y[count], 0, y[count].length);
+ System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
+ yDot[count], 0, yDot[count].length);
+
+ if (count == t.length - 1) {
+
+ // this was the last step we needed, we can compute the derivatives
+ stepStart = t[0];
+ stepSize = (t[t.length - 1] - t[0]) / (t.length - 1);
+
+ // first scaled derivative
+ scaled = yDot[0].clone();
+ for (int j = 0; j < scaled.length; ++j) {
+ scaled[j] *= stepSize;
}
- multistep[i - 1] = msI;
- }
- nordsieck = initializeHighOrderDerivatives(scaled, multistep);
- // stop the integrator after the first step has been handled
- throw new InitializationCompletedMarkerException();
+ // higher order derivatives
+ nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
+
+ // stop the integrator now that all needed steps have been handled
+ throw new InitializationCompletedMarkerException();
+
+ }
}
@@ -367,10 +397,10 @@ public abstract class MultistepIntegrato
/** Marker exception used ONLY to stop the starter integrator after first step. */
private static class InitializationCompletedMarkerException
- extends MathUserException {
+ extends MathRuntimeException {
/** Serializable version identifier. */
- private static final long serialVersionUID = -4105805787353488365L;
+ private static final long serialVersionUID = -1914085471038046418L;
/** Simple constructor. */
public InitializationCompletedMarkerException() {
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Fri Apr 29 17:15:59 2011
@@ -56,7 +56,7 @@ import org.apache.commons.math.util.Fast
* 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
+ * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
* </pre></p>
*
* <p>The definitions above use the classical representation with several previous first
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java Fri Apr 29 17:15:59 2011
@@ -35,7 +35,7 @@ public abstract class AdamsIntegrator ex
private final AdamsNordsieckTransformer transformer;
/**
- * Build an Adams integrator with the given order and step control prameters.
+ * Build an Adams integrator with the given order and step control parameters.
* @param name name of the method
* @param nSteps number of steps of the method excluding the one being computed
* @param order order of the method
@@ -93,9 +93,10 @@ public abstract class AdamsIntegrator ex
/** {@inheritDoc} */
@Override
- protected Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
- final double[][] multistep) {
- return transformer.initializeHighOrderDerivatives(first, multistep);
+ protected Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
+ final double[][] y,
+ final double[][] yDot) {
+ return transformer.initializeHighOrderDerivatives(h, t, y, yDot);
}
/** Update the high order scaled derivatives for Adams integrators (phase 1).
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Fri Apr 29 17:15:59 2011
@@ -62,7 +62,7 @@ import org.apache.commons.math.util.Fast
* 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
+ * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
* </pre></p>
*
* <p>The definitions above use the classical representation with several previous first
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.java?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsNordsieckTransformer.java Fri Apr 29 17:15:59 2011
@@ -24,14 +24,15 @@ 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.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.QRDecomposition;
+import org.apache.commons.math.linear.QRDecompositionImpl;
/** Transformer to Nordsieck vectors for Adams integrators.
- * <p>This class i used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
+ * <p>This class is used by {@link AdamsBashforthIntegrator Adams-Bashforth} and
* {@link AdamsMoultonIntegrator Adams-Moulton} integrators to convert between
* classical representation with several previous first derivatives and Nordsieck
* representation with higher order scaled derivatives.</p>
@@ -42,7 +43,7 @@ import org.apache.commons.math.linear.Ma
* 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
+ * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
* </pre></p>
*
* <p>With the previous definition, the classical representation of multistep methods
@@ -136,9 +137,6 @@ public class AdamsNordsieckTransformer {
private static final Map<Integer, AdamsNordsieckTransformer> CACHE =
new HashMap<Integer, AdamsNordsieckTransformer>();
- /** Initialization matrix for the higher order derivatives wrt y'', y''' ... */
- private final Array2DRowRealMatrix initialization;
-
/** Update matrix for the higher order derivatives h<sup>2</sup>/2y'', h<sup>3</sup>/6 y''' ... */
private final Array2DRowRealMatrix update;
@@ -173,19 +171,7 @@ public class AdamsNordsieckTransformer {
FieldMatrix<BigFraction> bigMSupdate =
pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false));
- // initialization coefficients, computed from a R 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> bigRInverse =
- new FieldLUDecompositionImpl<BigFraction>(bigP).getSolver().getInverse();
-
// convert coefficients to double
- initialization = MatrixUtils.bigFractionMatrixToRealMatrix(bigRInverse);
update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate);
c1 = new double[nSteps];
for (int i = 0; i < nSteps; ++i) {
@@ -251,21 +237,53 @@ public class AdamsNordsieckTransformer {
}
- /** 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 Array2DRowRealMatrix 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];
+ /** {@inheritDoc} */
+ public Array2DRowRealMatrix initializeHighOrderDerivatives(final double h, final double[] t,
+ final double[][] y,
+ final double[][] yDot) {
+
+ // using Taylor series with di = ti - t0, we get:
+ // y(ti) - y(t0) - di y'(t0) = di^2 / h^2 s2 + ... + di^k / h^k sk + O(h^(k+1))
+ // y'(ti) - y'(t0) = 2 di / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^k)
+ // we write these relations for i = 1 to i= n-1 as a set of 2(n-1) linear
+ // equations depending on the Nordsieck vector [s2 ... sk]
+ final double[][] a = new double[2 * (y.length - 1)][c1.length];
+ final double[][] b = new double[2 * (y.length - 1)][y[0].length];
+ final double[] y0 = y[0];
+ final double[] yDot0 = yDot[0];
+ for (int i = 1; i < y.length; ++i) {
+
+ final double di = t[i] - t[0];
+ final double ratio = di / h;
+ double dikM1Ohk = 1 / h;
+
+ // linear coefficients of equations
+ // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0)
+ final double[] aI = a[2 * i - 2];
+ final double[] aDotI = a[2 * i - 1];
+ for (int j = 0; j < aI.length; ++j) {
+ dikM1Ohk *= ratio;
+ aI[j] = di * dikM1Ohk;
+ aDotI[j] = (j + 2) * dikM1Ohk;
}
+
+ // expected value of the previous equations
+ final double[] yI = y[i];
+ final double[] yDotI = yDot[i];
+ final double[] bI = b[2 * i - 2];
+ final double[] bDotI = b[2 * i - 1];
+ for (int j = 0; j < yI.length; ++j) {
+ bI[j] = yI[j] - y0[j] - di * yDot0[j];
+ bDotI[j] = yDotI[j] - yDot0[j];
+ }
+
}
- return initialization.multiply(new Array2DRowRealMatrix(multistep, false));
+
+ // solve the rectangular system in the least square sense
+ // to get the best estimate of the Nordsieck vector [s2 ... sk]
+ QRDecomposition decomposition = new QRDecompositionImpl(new Array2DRowRealMatrix(a, false));
+ return new Array2DRowRealMatrix(decomposition.getSolver().solve(b), false);
+
}
/** Update the high order scaled derivatives for Adams integrators (phase 1).
Modified: commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties (original)
+++ commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties Fri Apr 29 17:15:59 2011
@@ -100,7 +100,7 @@ INSUFFICIENT_DIMENSION = dimension {0} i
DIMENSION = dimension ({0})
INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE = l''\u00e9chantillon ne contient que {0} points alors qu''au moins {1} sont n\u00e9cessaires
INSUFFICIENT_ROWS_AND_COLUMNS = donn\u00e9es insuffisantes : seulement {0} lignes et {1} colonnes.
-INTEGRATION_METHOD_NEEDS_AT_LEAST_ONE_PREVIOUS_POINT = la m\u00e9thode {0} n\u00e9cessite au moins un point pr\u00e9c\u00e9dent
+INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS = la m\u00e9thode {0} n\u00e9cessite au moins deux points pr\u00e9c\u00e9dents
INTERNAL_ERROR = erreur interne, veuillez signaler l''erreur \u00e0 {0}
INVALID_BRACKETING_PARAMETERS = param\u00e8tres d''encadrement invalides : borne inf\u00e9rieure = {0}, valeur initiale = {1}, borne sup\u00e9rieure = {2}
INVALID_INTERVAL_INITIAL_VALUE_PARAMETERS = param\u00e8tres de l''intervalle initial invalides : borne inf = {0}, valeur initiale = {1}, borne sup = {2}
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Apr 29 17:15:59 2011
@@ -52,6 +52,12 @@ The <action> type attribute can be add,u
If the output is not quite correct, check for invisible trailing spaces!
-->
<release version="3.0" date="TBD" description="TBD">
+ <action dev="luc" type="fix" >
+ Fixed initialization of multistep ODE integrators. Relying on the interpolation model
+ of the starter integrator inside only one step was wrong. The model may have a too
+ low order to compute high degrees derivatives in the Nordsieck vector. Now we use several
+ steps and use only grid points instead of interpolated points.
+ </action>
<action dev="luc" type="add" issue="MATH-564">
Added solve methods using double[][] to linear algebra decomposition solvers.
</action>
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java Fri Apr 29 17:15:59 2011
@@ -82,11 +82,11 @@ public class AdamsBashforthIntegratorTes
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- // the 31 and 36 factors are only valid for this test
+ // the 50 and 300 factors are only valid for this test
// and has been obtained from trial and error
// there is no general relation between local and global errors
- Assert.assertTrue(handler.getMaximalValueError() > (31.0 * scalAbsoluteTolerance));
- Assert.assertTrue(handler.getMaximalValueError() < (36.0 * scalAbsoluteTolerance));
+ Assert.assertTrue(handler.getMaximalValueError() > (50.0 * scalAbsoluteTolerance));
+ Assert.assertTrue(handler.getMaximalValueError() < (300.0 * scalAbsoluteTolerance));
Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
int calls = pb.getCalls();
@@ -126,8 +126,8 @@ public class AdamsBashforthIntegratorTes
integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- Assert.assertTrue(handler.getLastError() < 1.0e-8);
- Assert.assertTrue(handler.getMaximalValueError() < 1.0e-8);
+ Assert.assertTrue(handler.getLastError() < 1.5e-8);
+ Assert.assertTrue(handler.getMaximalValueError() < 1.5e-8);
Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
Assert.assertEquals("Adams-Bashforth", integ.getName());
}
@@ -137,17 +137,17 @@ public class AdamsBashforthIntegratorTes
TestProblem6 pb = new TestProblem6();
double range = FastMath.abs(pb.getFinalTime() - pb.getInitialTime());
- for (int nSteps = 1; nSteps < 8; ++nSteps) {
+ for (int nSteps = 2; nSteps < 8; ++nSteps) {
AdamsBashforthIntegrator integ =
- new AdamsBashforthIntegrator(nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-10, 1.0e-10);
+ new AdamsBashforthIntegrator(nSteps, 1.0e-6 * range, 0.1 * range, 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 (nSteps < 4) {
- Assert.assertTrue(integ.getEvaluations() > 150);
+ Assert.assertTrue(handler.getMaximalValueError() > 1.0e-03);
} else {
- Assert.assertTrue(integ.getEvaluations() < 70);
+ Assert.assertTrue(handler.getMaximalValueError() < 4.0e-12);
}
}
Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java?rev=1097891&r1=1097890&r2=1097891&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java Fri Apr 29 17:15:59 2011
@@ -82,11 +82,11 @@ public class AdamsMoultonIntegratorTest
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- // the 0.15 and 3.0 factors are only valid for this test
+ // the 0.5 and 11.0 factors are only valid for this test
// and has been obtained from trial and error
// there is no general relation between local and global errors
- Assert.assertTrue(handler.getMaximalValueError() > (0.15 * scalAbsoluteTolerance));
- Assert.assertTrue(handler.getMaximalValueError() < (3.0 * scalAbsoluteTolerance));
+ Assert.assertTrue(handler.getMaximalValueError() > ( 0.5 * scalAbsoluteTolerance));
+ Assert.assertTrue(handler.getMaximalValueError() < (11.0 * scalAbsoluteTolerance));
Assert.assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
int calls = pb.getCalls();
@@ -137,17 +137,17 @@ public class AdamsMoultonIntegratorTest
TestProblem6 pb = new TestProblem6();
double range = FastMath.abs(pb.getFinalTime() - pb.getInitialTime());
- for (int nSteps = 1; nSteps < 7; ++nSteps) {
+ for (int nSteps = 2; nSteps < 8; ++nSteps) {
AdamsMoultonIntegrator integ =
- new AdamsMoultonIntegrator(nSteps, 1.0e-6 * range, 0.1 * range, 1.0e-9, 1.0e-9);
+ new AdamsMoultonIntegrator(nSteps, 1.0e-6 * range, 0.1 * range, 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 (nSteps < 4) {
- Assert.assertTrue(integ.getEvaluations() > 140);
+ Assert.assertTrue(handler.getMaximalValueError() > 7.0e-04);
} else {
- Assert.assertTrue(integ.getEvaluations() < 90);
+ Assert.assertTrue(handler.getMaximalValueError() < 3.0e-13);
}
}