Author: luc
Date: Sun May 31 22:02:30 2009
New Revision: 780513
URL: http://svn.apache.org/viewvc?rev=780513&view=rev
Log:
added a step interpolator for multistep methods in Nordsieck form
Added:
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
- copied, changed from r777113,
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
(with props)
Copied:
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
(from r777113,
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java)
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java?p2=commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java&p1=commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java&r1=777113&r2=780513&rev=780513&view=diff
==============================================================================
---
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java
(original)
+++
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
Sun May 31 22:02:30 2009
@@ -15,59 +15,44 @@
* limitations under the License.
*/
-package org.apache.commons.math.ode.nonstiff;
+package org.apache.commons.math.ode.sampling;
import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
+import java.util.Arrays;
-import org.apache.commons.math.fraction.Fraction;
+import org.apache.commons.math.MathRuntimeException;
+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.ode.DerivativeException;
-import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
-import org.apache.commons.math.ode.sampling.MultistepStepInterpolator;
-import org.apache.commons.math.ode.sampling.StepInterpolator;
+import org.apache.commons.math.ode.nonstiff.AdamsIntegrator;
/**
- * This class implements an interpolator for Adams-Bashforth multiple steps.
+ * This class implements an interpolator for integrators using Nordsieck
representation.
*
- * <p>This interpolator computes dense output inside the last few
- * steps computed. The interpolation equation is consistent with the
- * integration scheme, it is based on a kind of <em>rollback</em> of the
- * integration from step end to interpolation date:
- * <pre>
- * y(t<sub>n</sub> + theta h) = y (t<sub>n</sub> + h) -
∫<sub>t<sub>n</sub> + theta h</sub><sup>t<sub>n</sub> + h</sup>p(t)dt
- * </pre>
- * where theta belongs to [0 ; 1] and p(t) is the interpolation polynomial
based on
- * the derivatives at previous steps f<sub>n-k+1</sub>, f<sub>n-k+2</sub> ...
- * f<sub>n</sub> and f<sub>n</sub>.</p>
+ * <p>This interpolator computes dense output around the current point.
+ * The interpolation equation is based on Taylor series formulas.
*
- * @see AdamsBashforthIntegrator
+ * @see AdamsIntegrator
* @version $Revision$ $Date$
* @since 2.0
*/
-class AdamsBashforthStepInterpolator extends MultistepStepInterpolator {
+public class NordsieckStepInterpolator extends AbstractStepInterpolator {
/** Serializable version identifier */
private static final long serialVersionUID = -7179861704951334960L;
- /** Neville's interpolation array. */
- private double[] neville;
+ /** Step size used in the first scaled derivative and Nordsieck vector. */
+ private double scalingH;
- /** Integration rollback array. */
- private double[] rollback;
+ /** First scaled derivative. */
+ private double[] scaled;
- /** γ array. */
- private double[] gamma;
-
- /** Backward differences array. */
- private int[][] bdArray;
-
- /** Original non-truncated step end time. */
- private double nonTruncatedEnd;
-
- /** Original non-truncated step size. */
- private double nonTruncatedH;
+ /** Nordsieck vector. */
+ private RealMatrix nordsieck;
/** Simple constructor.
* This constructor builds an instance that is not usable yet, the
@@ -76,7 +61,7 @@
* constructor is used only in order to delay the initialization in
* some cases.
*/
- public AdamsBashforthStepInterpolator() {
+ public NordsieckStepInterpolator() {
}
/** Copy constructor.
@@ -84,214 +69,200 @@
* copy: its arrays are separated from the original arrays of the
* instance
*/
- public AdamsBashforthStepInterpolator(final AdamsBashforthStepInterpolator
interpolator) {
+ public NordsieckStepInterpolator(final NordsieckStepInterpolator
interpolator) {
super(interpolator);
- nonTruncatedEnd = interpolator.nonTruncatedEnd;
- nonTruncatedH = interpolator.nonTruncatedH;
+ scalingH = interpolator.scalingH;
+ if (interpolator.scaled != null) {
+ scaled = interpolator.scaled.clone();
+ }
+ if (interpolator.nordsieck != null) {
+ nordsieck = interpolator.nordsieck.copy();
+ }
}
/** {...@inheritdoc} */
@Override
protected StepInterpolator doCopy() {
- return new AdamsBashforthStepInterpolator(this);
- }
-
- /** {...@inheritdoc} */
- @Override
- protected void initializeCoefficients() {
-
- neville = new double[previousF.length];
- rollback = new double[previousF.length];
-
- bdArray =
AdamsBashforthIntegrator.computeBackwardDifferencesArray(previousF.length);
-
- Fraction[] fGamma =
AdamsBashforthIntegrator.computeGammaArray(previousF.length);
- gamma = new double[fGamma.length];
- for (int i = 0; i < fGamma.length; ++i) {
- gamma[i] = fGamma[i].doubleValue();
- }
-
+ return new NordsieckStepInterpolator(this);
}
- /** {...@inheritdoc} */
+ /** 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
+ * the end of the step
+ * @param forward integration direction indicator
+ */
@Override
- public void storeTime(final double t) {
- nonTruncatedEnd = t;
- nonTruncatedH = nonTruncatedEnd - previousTime;
- super.storeTime(t);
+ public void reinitialize(final double[] y, final boolean forward) {
+ super.reinitialize(y, forward);
}
- /** Truncate a step.
- * <p>Truncating a step is necessary when an event is triggered
- * before the nominal end of the step.</p>
- * @param truncatedEndTime end time of truncated step
+ /** 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 scalingH step size used in the scaled and nordsieck arrays
+ * @param scaled reference to the integrator array holding the first
+ * scaled derivative
+ * @param nordsieck reference to the integrator matrix holding the
+ * nordsieck vector
*/
- void truncateStep(final double truncatedEndTime) {
- currentTime = truncatedEndTime;
- h = currentTime - previousTime;
+ public void reinitialize(final double scalingH, final double[] scaled,
+ final RealMatrix nordsieck) {
+ this.scalingH = scalingH;
+ this.scaled = scaled;
+ this.nordsieck = nordsieck;
}
- /** {...@inheritdoc} */
+ /** Store the current step time.
+ * @param t current time
+ */
@Override
- public void setInterpolatedTime(final double time)
- throws DerivativeException {
- interpolatedTime = time;
- final double oneMinusThetaH = nonTruncatedEnd - interpolatedTime;
- final double theta = (nonTruncatedH == 0) ?
- 0 : (nonTruncatedH - oneMinusThetaH) /
nonTruncatedH;
- computeInterpolatedState(theta, oneMinusThetaH);
+ public void storeTime(final double t) {
+ currentTime = t;
+ h = currentTime - previousTime;
+ interpolatedTime = t;
+ computeInterpolatedState(1.0, 0.0);
}
/** {...@inheritdoc} */
@Override
protected void computeInterpolatedState(final double theta, final double
oneMinusThetaH) {
- interpolateDerivatives();
- interpolateState(theta);
+ final double x = theta * h;
+ nordsieck.walkInOptimizedOrder(new StateEstimator(x, x / scalingH));
}
- /** Interpolate the derivatives.
- * <p>The Adams method is based on a polynomial interpolation of the
- * derivatives based on the preceding steps. So the interpolation of
- * the derivatives here is strictly equivalent: it is a simple polynomial
- * interpolation.</p>
- */
- private void interpolateDerivatives() {
-
- for (int i = 0; i < interpolatedDerivatives.length; ++i) {
+ /** State estimator. */
+ private class StateEstimator implements RealMatrixPreservingVisitor {
- // initialize the Neville's interpolation algorithm
- for (int k = 0; k < previousF.length; ++k) {
- neville[k] = previousF[k][i];
+ /** Scaling factor for derivative. */
+ private final double scale;
+
+ /** First order power. */
+ private final double lowPower;
+
+ /** High order powers. */
+ private final double[] highPowers;
+
+ /** Simple constructor.
+ * @param scale scaling factor for derivative
+ * @param theta normalized interpolation abscissa within the step
+ */
+ public StateEstimator(final double scale, final double theta) {
+ this.scale = scale;
+ lowPower = theta;
+ highPowers = new double[nordsieck.getRowDimension()];
+ double thetaN = theta;
+ for (int i = 0; i < highPowers.length; ++i) {
+ thetaN *= theta;
+ highPowers[i] = thetaN;
}
+ }
- // combine the contributions of each points
- for (int l = 1; l < neville.length; ++l) {
- for (int m = neville.length - 1; m >= l; --m) {
- final double xm = previousT[m];
- final double xmMl = previousT[m - l];
- neville[m] = ((interpolatedTime - xm) * neville[m-1] +
- (xmMl - interpolatedTime) * neville[m]) /
(xmMl - xm);
- }
+ /** {...@inheritdoc} */
+ public void start(int rows, int columns,
+ int startRow, int endRow, int startColumn, int
endColumn) {
+ Arrays.fill(interpolatedState, 0.0);
+ Arrays.fill(interpolatedDerivatives, 0.0);
+ }
+
+ /** {...@inheritdoc} */
+ public void visit(int row, int column, double value) {
+ final double d = value * highPowers[row];
+ interpolatedState[column] += d;
+ interpolatedDerivatives[column] += (row + 2) * d;
+ }
+
+ /** {...@inheritdoc} */
+ public double end() {
+ for (int j = 0; j < currentState.length; ++j) {
+ interpolatedState[j] += currentState[j] + scaled[j] * lowPower;
+ interpolatedDerivatives[j] =
+ (interpolatedDerivatives[j] + scaled[j] * lowPower) /
scale;
}
-
- // the interpolation polynomial value is in the array last element
- interpolatedDerivatives[i] = neville[neville.length - 1];
-
+ return 0;
}
}
- /** Interpolate the state.
- * <p>The Adams method is based on a polynomial interpolation of the
- * derivatives based on the preceding steps. The polynomial model is
- * integrated analytically throughout the last step. Using the notations
- * found in the second edition of the first volume (Nonstiff Problems)
- * of the reference book by Hairer, Norsett and Wanner: <i>Solving Ordinary
- * Differential Equations</i> (Springer-Verlag, ISBN 3-540-56670-8), this
- * process leads to the following expression:</p>
- * <pre>
- * y<sub>n+1</sub> = y<sub>n</sub> +
- * h × ∑<sub>j=0</sub><sup>j=k-1</sup>
γ<sub>j</sub>∇<sup>j</sup>f<sub>n</sub>
- * </pre>
- * <p>In the previous expression, the γ<sub>j</sub> terms are the
- * ones that result from the analytical integration, and can be computed
form
- * the binomial coefficients C<sub>j</sub><sup>-s</sup>:</p>
- * <p>
- * γ<sub>j</sub> =
(-1)<sup>j</sup>∫<sub>0</sub><sup>1</sup>C<sub>j</sub><sup>-s</sup>ds
- * </p>
- * <p>In order to interpolate the state in a manner that is consistent
with the
- * integration scheme, we simply subtract from the current state (at the
end of the step)
- * the integral computed from interpolation time to step end time.</p>
- * <p>
- * η<sub>j</sub>(θ)=
- *
(-1)<sup>j</sup>∫<sub>θ</sub><sup>1</sup>C<sub>j</sub><sup>-s</sup>ds
- * </p>
- * The method described in the Hairer, Norsett and Wanner book to compute
γ<sub>j</sub>
- * is easily extended to compute γ<sub>j</sub>(θ)=
- *
(-1)<sup>j</sup>∫<sub>0</sub><sup>θ</sup>C<sub>j</sub><sup>-s</sup>ds.
From this,
- * we can compute η<sub>j</sub>(θ) =
γ<sub>j</sub>-γ<sub>j</sub>(θ).
- * The first few values are:</p>
- * <table>
- *
<tr><td>j</td><td>γ<sub>j</sub></td><td>γ<sub>j</sub>(θ)</td><td>η<sub>j</sub>(θ)</td></tr>
- * <tr><td>0</td><td>1</td><td></td>θ<td>1-θ</td></tr>
- *
<tr><td>1</td><td>1/2</td><td></td>θ<sup>2</sup>/2<td>(1-θ<sup>2</sup>)/2</td></tr>
- *
<tr><td>2</td><td>5/12</td><td></td>(3θ<sup>2</sup>+2θ<sup>3</sup>)/12<td>(5-3θ<sup>2</sup>-2θ<sup>3</sup>)/12</td></tr>
- * </table>
- * <p>
- * The η<sub>j</sub>(θ) functions appear to be polynomial ones.
As expected,
- * we see that η<sub>j</sub>(1)= 0. The recurrence relation derived for
- * γ<sub>j</sub>(θ) is:
- * </p>
- * <p>
- * &sum<sub>j=0</sub><sup>j=m</sup>γ<sub>j</sub>(θ)/(m+1-j) =
- * 1/(m+1)! ∏<sub>k=0</sub><sup>k=m</sup>(θ+k)
- * </p>
- * @param theta location of the interpolation point within the last step
- */
- private void interpolateState(final double theta) {
+ /** {...@inheritdoc} */
+ @Override
+ public void writeExternal(final ObjectOutput out)
+ throws IOException {
- // compute the integrals to remove from the final state
- computeRollback(previousT.length - 1, theta);
+ // save the state of the base class
+ writeBaseExternal(out);
- // remove these integrals from the final state
- for (int j = 0; j < interpolatedState.length; ++j) {
- double sum = 0;
- for (int l = 0; l < previousT.length; ++l) {
- sum += rollback[l] * previousF[l][j];
+ // save the local attributes
+ final int n = (currentState == null) ? -1 : currentState.length;
+ if (scaled == null) {
+ out.writeBoolean(false);
+ } else {
+ out.writeBoolean(true);
+ for (int j = 0; j < n; ++j) {
+ out.writeDouble(scaled[j]);
+ }
+ }
+
+ if (nordsieck == null) {
+ out.writeBoolean(false);
+ } else {
+ out.writeBoolean(true);
+ final int rows = nordsieck.getRowDimension();
+ out.writeInt(rows);
+ for (int i = 0; i < rows; ++i) {
+ for (int j = 0; j < n; ++j) {
+ out.writeDouble(nordsieck.getEntry(i, j));
+ }
}
- interpolatedState[j] = currentState[j] - h * sum;
}
}
- /** Compute the rollback coefficients.
- * @param order order of the integration method
- * @param theta current value for theta
- */
- private void computeRollback(final int order, final double theta) {
+ /** {...@inheritdoc} */
+ @Override
+ public void readExternal(final ObjectInput in)
+ throws IOException {
+
+ // read the base class
+ final double t = readBaseExternal(in);
- // compute the gamma(theta) values from the recurrence relation
- double product = theta;
- rollback[0] = theta;
- for (int i = 1; i < order; ++i) {
- product *= (i + theta) / (i + 1);
- double g = product;
- for (int j = 1; j <= i; ++j) {
- g -= rollback[i - j] / (j + 1);
+ // read the local attributes
+ final int n = (currentState == null) ? -1 : currentState.length;
+ final boolean hasScaled = in.readBoolean();
+ if (hasScaled) {
+ scaled = new double[n];
+ for (int j = 0; j < n; ++j) {
+ scaled[j] = in.readDouble();
}
- rollback[i] = g;
+ } else {
+ scaled = null;
}
- // subtract it from gamma to get eta(theta)
- for (int i = 0; i < order; ++i) {
- rollback[i] -= gamma[i];
+ final boolean hasNordsieck = in.readBoolean();
+ if (hasNordsieck) {
+ final int rows = in.readInt();
+ final double[][] nData = new double[rows][n];
+ for (int i = 0; i < rows; ++i) {
+ final double[] nI = nData[i];
+ for (int j = 0; j < n; ++j) {
+ nI[j] = in.readDouble();
+ }
+ }
+ nordsieck = new RealMatrixImpl(nData, false);
+ } else {
+ nordsieck = null;
}
- // combine the eta integrals with the backward differences array
- // to get the rollback coefficients
- for (int i = 0; i < order; ++i) {
- double f = 0;
- for (int j = i; j <= order; ++j) {
- f -= rollback[j] * bdArray[j][i];
+ try {
+ if (hasScaled && hasNordsieck) {
+ // we can now set the interpolated time and state
+ setInterpolatedTime(t);
}
- rollback[i] = f;
+ } catch (DerivativeException e) {
+ throw MathRuntimeException.createIOException(e);
}
}
- /** {...@inheritdoc} */
- @Override
- public void writeExternal(final ObjectOutput out)
- throws IOException {
- super.writeExternal(out);
- out.writeDouble(nonTruncatedEnd);
- }
-
- /** {...@inheritdoc} */
- @Override
- public void readExternal(final ObjectInput in)
- throws IOException {
- nonTruncatedEnd = in.readDouble();
- }
-
}
Added:
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java?rev=780513&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
(added)
+++
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
Sun May 31 22:02:30 2009
@@ -0,0 +1,94 @@
+/*
+ * 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.sampling;
+
+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 java.util.Random;
+
+import org.apache.commons.math.ode.ContinuousOutputModel;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.IntegratorException;
+import org.apache.commons.math.ode.nonstiff.AdamsIntegrator;
+import org.apache.commons.math.ode.nonstiff.TestProblem1;
+import org.apache.commons.math.ode.nonstiff.TestProblem3;
+import org.junit.Test;
+
+public class NordsieckStepInterpolatorTest {
+
+ @Test
+ public void derivativesConsistency()
+ throws DerivativeException, IntegratorException {
+ TestProblem3 pb = new TestProblem3();
+ double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+ AdamsIntegrator integ = new AdamsIntegrator(5, false, step);
+ StepInterpolatorTestUtils.checkDerivativesConsistency(integ, pb,
1.0e-10);
+ }
+
+ @Test
+ public void serialization()
+ throws DerivativeException, IntegratorException,
+ IOException, ClassNotFoundException {
+
+ TestProblem1 pb = new TestProblem1();
+ double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.001;
+ AdamsIntegrator integ = new AdamsIntegrator(5, false, step);
+ integ.addStepHandler(new ContinuousOutputModel());
+ integ.integrate(pb,
+ pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
+
+ ByteArrayOutputStream bos = new ByteArrayOutputStream();
+ ObjectOutputStream oos = new ObjectOutputStream(bos);
+ for (StepHandler handler : integ.getStepHandlers()) {
+ oos.writeObject(handler);
+ }
+
+ assertTrue(bos.size () > 148000);
+ assertTrue(bos.size () < 149000);
+
+ ByteArrayInputStream bis = new
ByteArrayInputStream(bos.toByteArray());
+ ObjectInputStream ois = new ObjectInputStream(bis);
+ ContinuousOutputModel cm = (ContinuousOutputModel) ois.readObject();
+
+ Random random = new Random(347588535632l);
+ double maxError = 0.0;
+ for (int i = 0; i < 1000; ++i) {
+ double r = random.nextDouble();
+ double time = r * pb.getInitialTime() + (1.0 - r) *
pb.getFinalTime();
+ cm.setInterpolatedTime(time);
+ double[] interpolatedY = cm.getInterpolatedState ();
+ double[] theoreticalY = pb.computeTheoreticalState(time);
+ double dx = interpolatedY[0] - theoreticalY[0];
+ double dy = interpolatedY[1] - theoreticalY[1];
+ double error = dx * dx + dy * dy;
+ if (error > maxError) {
+ maxError = error;
+ }
+ }
+
+ assertTrue(maxError < 1.0e-6);
+
+ }
+
+}
Propchange:
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
