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 2008/07/14 23:13:23 UTC
svn commit: r676737 [1/2] - in /commons/proper/math/branches/MATH_2_0/src:
java/org/apache/commons/math/ode/nonstiff/ site/xdoc/
test/org/apache/commons/math/ode/nonstiff/
Author: luc
Date: Mon Jul 14 14:13:23 2008
New Revision: 676737
URL: http://svn.apache.org/viewvc?rev=676737&view=rev
Log:
New ODE integrators have been added: the explicit Adams-Bashforth and implicit
Adams-Moulton multistep methods. These methods support customizable starter
integrators and support discrete events even during the start phase.
All these methods provide the same rich features has the existing ones:
continuous output, step handlers, discrete events, G-stop ...
Added:
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (with props)
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java (with props)
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (with props)
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java (with props)
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java (with props)
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java (with props)
commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java (with props)
commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java (with props)
Modified:
commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml
Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,285 @@
+/*
+ * 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.nonstiff;
+
+import org.apache.commons.math.fraction.Fraction;
+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.events.CombinedEventsManager;
+import org.apache.commons.math.ode.sampling.StepHandler;
+
+
+/**
+ * This class implements explicit Adams-Bashforth integrators for Ordinary
+ * Differential Equations.
+ *
+ * <p>Adams-Bashforth (in fact due to Adams alone) methods are explicit
+ * 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, n-1, n-2 ... 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 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>...</li>
+ * </ul>
+ *
+ * <p>A k-steps Adams-Bashforth method is of order k. There is no limit to the
+ * value of k.</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>
+ *
+ * @see AdamsMoultonIntegrator
+ * @see BDFIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class AdamsBashforthIntegrator extends MultistepIntegrator {
+
+ /** Serializable version identifier. */
+ private static final long serialVersionUID = 1676381657635800870L;
+
+ /** Integrator method name. */
+ private static final String METHOD_NAME = "Adams-Bashforth";
+
+ /** Coefficients for the current method. */
+ private final double[] coeffs;
+
+ /** 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)
+ * @param step integration step size
+ */
+ public AdamsBashforthIntegrator(final int order, final double step) {
+
+ super(METHOD_NAME, order, new AdamsBashforthStepInterpolator());
+
+ // 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)));
+ }
+ coeffs[i] = f.doubleValue();
+ }
+
+ this.step = step;
+
+ }
+
+ /** {@inheritDoc} */
+ public double integrate(FirstOrderDifferentialEquations equations,
+ double t0, double[] y0, double t, double[] y)
+ throws DerivativeException, IntegratorException {
+
+ 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);
+ }
+ 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);
+
+ // set up integration control objects
+ stepStart = t0;
+ stepSize = step;
+ for (StepHandler handler : stepHandlers) {
+ handler.reset();
+ }
+ CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
+
+ // compute the first few steps using the configured starter integrator
+ double stopTime =
+ start(previousF.length, stepSize, manager, equations, stepStart, y);
+ if (Double.isNaN(previousT[0])) {
+ return stopTime;
+ }
+ stepStart = previousT[0];
+ interpolator.storeTime(stepStart);
+
+ boolean lastStep = false;
+ while (!lastStep) {
+
+ // 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];
+ }
+ yTmp[j] = y[j] + stepSize * sum;
+ }
+
+ // discrete events handling
+ interpolator.storeTime(stepStart + stepSize);
+ final boolean truncated;
+ if (manager.evaluateStep(interpolator)) {
+ truncated = true;
+ interpolator.truncateStep(manager.getEventTime());
+ } else {
+ truncated = false;
+ }
+
+ // the step has been accepted (may have been truncated)
+ final double nextStep = interpolator.getCurrentTime();
+ interpolator.setInterpolatedTime(nextStep);
+ System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
+ manager.stepAccepted(nextStep, y);
+ lastStep = manager.stop();
+
+ // provide the step data to the step handler
+ for (StepHandler handler : stepHandlers) {
+ handler.handleStep(interpolator, lastStep);
+ }
+ stepStart = nextStep;
+
+ if (!lastStep) {
+ // prepare next step
+
+ if (manager.reset(stepStart, y)) {
+
+ // some events handler has triggered changes that
+ // invalidate the derivatives, we need to restart from scratch
+ stopTime =
+ start(previousF.length, stepSize, manager, equations, stepStart, y);
+ if (Double.isNaN(previousT[0])) {
+ return stopTime;
+ }
+ stepStart = previousT[0];
+
+ } else {
+
+ 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();
+ }
+
+ // evaluate differential equations for next step
+ previousT[0] = stepStart;
+ equations.computeDerivatives(stepStart, y, previousF[0]);
+
+ }
+ }
+
+ }
+
+ stopTime = stepStart;
+ stepStart = Double.NaN;
+ stepSize = Double.NaN;
+ return stopTime;
+
+ }
+
+ /** 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
+ */
+ public double[] getCoeffs() {
+ return coeffs.clone();
+ }
+
+ /** Compute the backward differences coefficients array.
+ * <p>This is quite similar to the Pascal triangle, except for a
+ * (-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>
+ * <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>
+ * </pre>
+ * @param order order of the integration method
+ */
+ static int[][] computeBackwardDifferencesArray(final int order) {
+
+ // create the array
+ int[][] bdArray = new int[order][];
+
+ // recurrence initialization
+ bdArray[0] = new int[] { 1 };
+
+ // 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];
+ }
+ bdArray[i][i] = -bdArray[i - 1][i - 1];
+ }
+
+ return bdArray;
+
+ }
+
+ /** Compute the gamma coefficients.
+ * @param order order of the integration method
+ * @return gamma coefficients array
+ */
+ static Fraction[] computeGammaArray(final int order) {
+
+ // create the array
+ Fraction[] gammaArray = new Fraction[order];
+
+ // recurrence initialization
+ gammaArray[0] = Fraction.ONE;
+
+ // 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)));
+ }
+ gammaArray[i] = gamma;
+ }
+
+ return gammaArray;
+
+ }
+
+}
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,288 @@
+/*
+ * 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.nonstiff;
+
+import java.io.IOException;
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+
+import org.apache.commons.math.fraction.Fraction;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
+import org.apache.commons.math.ode.sampling.StepInterpolator;
+
+/**
+ * This class implements an interpolator for Adams-Bashforth multiple steps.
+ *
+ * <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>
+ *
+ * @see AdamsBashforthIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+
+class AdamsBashforthStepInterpolator extends MultistepStepInterpolator {
+
+ /** Serializable version identifier */
+ private static final long serialVersionUID = -7179861704951334960L;
+
+ /** Neville's interpolation array. */
+ private double[] neville;
+
+ /** Integration rollback array. */
+ private double[] rollback;
+
+ /** γ 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;
+
+ /** Simple constructor.
+ * This constructor builds an instance that is not usable yet, the
+ * {@link AbstractStepInterpolator#reinitialize} method should be called
+ * before using the instance in order to initialize the internal arrays. This
+ * constructor is used only in order to delay the initialization in
+ * some cases.
+ */
+ public AdamsBashforthStepInterpolator() {
+ }
+
+ /** Copy constructor.
+ * @param interpolator interpolator to copy from. The copy is a deep
+ * copy: its arrays are separated from the original arrays of the
+ * instance
+ */
+ public AdamsBashforthStepInterpolator(final AdamsBashforthStepInterpolator interpolator) {
+ super(interpolator);
+ nonTruncatedEnd = interpolator.nonTruncatedEnd;
+ nonTruncatedH = interpolator.nonTruncatedH;
+ }
+
+ /** {@inheritDoc} */
+ protected StepInterpolator doCopy() {
+ return new AdamsBashforthStepInterpolator(this);
+ }
+
+ /** {@inheritDoc} */
+ 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();
+ }
+
+ }
+
+ /** {@inheritDoc} */
+ public void storeTime(final double t) {
+ nonTruncatedEnd = t;
+ nonTruncatedH = nonTruncatedEnd - previousTime;
+ super.storeTime(t);
+ }
+
+ /** Truncate a step.
+ * <p>Truncating a step is necessary when an event is triggered
+ * before the nominal end of the step.</p>
+ */
+ void truncateStep(final double truncatedEndTime) {
+ currentTime = truncatedEndTime;
+ h = currentTime - previousTime;
+ }
+
+ /** {@inheritDoc} */
+ 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);
+ }
+
+ /** {@inheritDoc} */
+ protected void computeInterpolatedState(final double theta, final double oneMinusThetaH) {
+ interpolateDerivatives();
+ interpolateState(theta);
+ }
+
+ /** 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) {
+
+ // initialize the Neville's interpolation algorithm
+ for (int k = 0; k < previousF.length; ++k) {
+ neville[k] = previousF[k][i];
+ }
+
+ // 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);
+ }
+ }
+
+ // the interpolation polynomial value is in the array last element
+ interpolatedDerivatives[i] = neville[neville.length - 1];
+
+ }
+
+ }
+
+ /** 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) {
+
+ // compute the integrals to remove from the final state
+ computeRollback(previousT.length - 1, theta);
+
+ // 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];
+ }
+ 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) {
+
+ // 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);
+ }
+ rollback[i] = g;
+ }
+
+ // subtract it from gamma to get eta(theta)
+ for (int i = 0; i < order; ++i) {
+ rollback[i] -= gamma[i];
+ }
+
+ // 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];
+ }
+ rollback[i] = f;
+ }
+
+ }
+
+ /** {@inheritDoc} */
+ public void writeExternal(final ObjectOutput out)
+ throws IOException {
+ super.writeExternal(out);
+ out.writeDouble(nonTruncatedEnd);
+ }
+
+ /** {@inheritDoc} */
+ public void readExternal(final ObjectInput in)
+ throws IOException {
+ nonTruncatedEnd = in.readDouble();
+ }
+
+}
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,299 @@
+/*
+ * 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.nonstiff;
+
+import org.apache.commons.math.fraction.Fraction;
+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.events.CombinedEventsManager;
+import org.apache.commons.math.ode.sampling.StepHandler;
+
+
+/**
+ * This class implements implicit Adams-Moulton integrators for Ordinary
+ * Differential Equations.
+ *
+ * <p>Adams-Moulton (in fact 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 = 0: y<sub>n+1</sub> = y<sub>n</sub> + h f<sub>n+1</sub></li>
+ * <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + h (f<sub>n+1</sub>+f<sub>n</sub>)/2</li>
+ * <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + h (5f<sub>n+1</sub>+8f<sub>n</sub>-f<sub>n-1</sub>)/12</li>
+ * <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + h (9f<sub>n+1</sub>+19f<sub>n</sub>-5f<sub>n-1</sub>+f<sub>n-2)/24</sub></li>
+ * <li>...</li>
+ * </ul>
+ *
+ * <p>The coefficients are computed (and cached) as needed, so their are no
+ * theoretical limitations to the number of steps</p>
+ *
+ * <p>A k-steps Adams-Moulton method is of order k+1. There is no limit to the
+ * value of k.</p>
+ *
+ * <p>These methods are implicit: f<sub>n+1</sub> is used to compute
+ * y<sub>n+1</sub>. Simpler explicit Adams methods exist: the
+ * Adams-Bashforth methods (which are also due to Adams alone). They are
+ * provided by the {@link AdamsBashforthIntegrator AdamsBashforthIntegrator} class.</p>
+ *
+ * @see AdamsBashforthIntegrator
+ * @see BDFIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public class AdamsMoultonIntegrator extends MultistepIntegrator {
+
+ /** Serializable version identifier. */
+ private static final long serialVersionUID = 4990335331377040417L;
+
+ /** Integrator method name. */
+ private static final String METHOD_NAME = "Adams-Moulton";
+
+ /** Coefficients for the predictor phase of the method. */
+ private final double[] predictorCoeffs;
+
+ /** Coefficients for the corrector phase of the method. */
+ private final double[] correctorCoeffs;
+
+ /** Integration step. */
+ private final double step;
+
+ /**
+ * Build an Adams-Moulton integrator with the given order and step size.
+ * @param order order of the method (must be strictly positive)
+ * @param step integration step size
+ */
+ public AdamsMoultonIntegrator(final int order, final double step) {
+
+ super(METHOD_NAME, order + 1, new AdamsMoultonStepInterpolator());
+
+ // compute the integration coefficients
+ int[][] bdArray = AdamsBashforthIntegrator.computeBackwardDifferencesArray(order + 1);
+
+ Fraction[] gamma = AdamsBashforthIntegrator.computeGammaArray(order);
+ predictorCoeffs = new double[order];
+ for (int i = 0; i < order; ++i) {
+ Fraction fPredictor = Fraction.ZERO;
+ for (int j = i; j < order; ++j) {
+ Fraction f = new Fraction(bdArray[j][i], 1);
+ fPredictor = fPredictor.add(gamma[j].multiply(f));
+ }
+ predictorCoeffs[i] = fPredictor.doubleValue();
+ }
+
+ Fraction[] gammaStar = computeGammaStarArray(order);
+ correctorCoeffs = new double[order + 1];
+ for (int i = 0; i <= order; ++i) {
+ Fraction fCorrector = Fraction.ZERO;
+ for (int j = i; j <= order; ++j) {
+ Fraction f = new Fraction(bdArray[j][i], 1);
+ fCorrector = fCorrector.add(gammaStar[j].multiply(f));
+ }
+ correctorCoeffs[i] = fCorrector.doubleValue();
+ }
+
+ this.step = step;
+
+ }
+
+ /** {@inheritDoc} */
+ public double integrate(FirstOrderDifferentialEquations equations,
+ double t0, double[] y0, double t, double[] y)
+ throws DerivativeException, IntegratorException {
+
+ 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);
+ }
+ final double[] yTmp = new double[y0.length];
+
+ // set up an interpolator sharing the integrator arrays
+ final AdamsMoultonStepInterpolator interpolator =
+ (AdamsMoultonStepInterpolator) prototype.copy();
+ interpolator.reinitialize(yTmp, previousT, previousF, forward);
+
+ // set up integration control objects
+ stepStart = t0;
+ stepSize = step;
+ for (StepHandler handler : stepHandlers) {
+ handler.reset();
+ }
+ CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
+
+ // compute the first few steps using the configured starter integrator
+ double stopTime =
+ start(previousF.length - 1, stepSize, manager, equations, stepStart, y);
+ if (Double.isNaN(previousT[0])) {
+ return stopTime;
+ }
+ stepStart = previousT[0];
+ rotatePreviousSteps();
+ previousF[0] = new double[y0.length];
+ interpolator.storeTime(stepStart);
+
+ boolean lastStep = false;
+ while (!lastStep) {
+
+ // shift all data
+ interpolator.shift();
+
+ // predict state at end of step
+ for (int j = 0; j < y0.length; ++j) {
+ double sum = 0;
+ for (int l = 0; l < predictorCoeffs.length; ++l) {
+ sum += predictorCoeffs[l] * previousF[l+1][j];
+ }
+ yTmp[j] = y[j] + stepSize * sum;
+ }
+
+ // evaluate the derivatives
+ final double stepEnd = stepStart + stepSize;
+ equations.computeDerivatives(stepEnd, yTmp, previousF[0]);
+
+ // apply corrector
+ for (int j = 0; j < y0.length; ++j) {
+ double sum = 0;
+ for (int l = 0; l < correctorCoeffs.length; ++l) {
+ sum += correctorCoeffs[l] * previousF[l][j];
+ }
+ yTmp[j] = y[j] + stepSize * sum;
+ }
+
+ // discrete events handling
+ interpolator.storeTime(stepEnd);
+ final boolean truncated;
+ if (manager.evaluateStep(interpolator)) {
+ truncated = true;
+ interpolator.truncateStep(manager.getEventTime());
+ } else {
+ truncated = false;
+ }
+
+ // the step has been accepted (may have been truncated)
+ final double nextStep = interpolator.getCurrentTime();
+ interpolator.setInterpolatedTime(nextStep);
+ System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
+ manager.stepAccepted(nextStep, y);
+ lastStep = manager.stop();
+
+ // provide the step data to the step handler
+ for (StepHandler handler : stepHandlers) {
+ handler.handleStep(interpolator, lastStep);
+ }
+ stepStart = nextStep;
+
+ if (!lastStep) {
+ // prepare next step
+
+ if (manager.reset(stepStart, y)) {
+
+ // some events handler has triggered changes that
+ // invalidate the derivatives, we need to restart from scratch
+ stopTime =
+ start(previousF.length - 1, stepSize, manager, equations, stepStart, y);
+ if (Double.isNaN(previousT[0])) {
+ return stopTime;
+ }
+ stepStart = previousT[0];
+ rotatePreviousSteps();
+ previousF[0] = new double[y0.length];
+
+ } else {
+
+ 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();
+ }
+
+ // evaluate differential equations for next step
+ previousT[0] = stepStart;
+ equations.computeDerivatives(stepStart, y, previousF[0]);
+
+ }
+ }
+
+ }
+
+ stopTime = stepStart;
+ stepStart = Double.NaN;
+ stepSize = Double.NaN;
+ return stopTime;
+
+ }
+
+ /** Get the coefficients of the prdictor phase 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
+ */
+ public double[] getPredictorCoeffs() {
+ return predictorCoeffs.clone();
+ }
+
+ /** Get the coefficients of the corrector phase 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</sup> c<sub>i</sub>f<sub>n-i</sub></li>
+ * </pre>
+ * @return a copy of the coefficients of the method
+ */
+ public double[] getCorrectorCoeffs() {
+ return correctorCoeffs.clone();
+ }
+
+ /** Compute the gamma star coefficients.
+ * @param order order of the integration method
+ * @return gamma star coefficients array
+ */
+ static Fraction[] computeGammaStarArray(final int order) {
+
+ // create the array
+ Fraction[] gammaStarArray = new Fraction[order + 1];
+
+ // recurrence initialization
+ gammaStarArray[0] = Fraction.ONE;
+
+ // fill up array using recurrence relation
+ for (int i = 1; i <= order; ++i) {
+ Fraction gammaStar = Fraction.ZERO;
+ for (int j = 1; j <= i; ++j) {
+ gammaStar = gammaStar.subtract(gammaStarArray[i - j].multiply(new Fraction(1, j + 1)));
+ }
+ gammaStarArray[i] = gammaStar;
+ }
+
+ return gammaStarArray;
+
+ }
+
+}
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,289 @@
+/*
+ * 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.nonstiff;
+
+import java.io.IOException;
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+
+import org.apache.commons.math.fraction.Fraction;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
+import org.apache.commons.math.ode.sampling.StepInterpolator;
+
+/**
+ * This class implements an interpolator for Adams-Moulton multiple steps.
+ *
+ * <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>, f<sub>n</sub> and f<sub>n+1</sub>.</p>
+ *
+ * @see AdamsMoultonIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+
+class AdamsMoultonStepInterpolator extends MultistepStepInterpolator {
+
+ /** Serializable version identifier */
+ private static final long serialVersionUID = 735568489801241899L;
+
+ /** Neville's interpolation array. */
+ private double[] neville;
+
+ /** Integration rollback array. */
+ private double[] rollback;
+
+ /** γ star array. */
+ private double[] gammaStar;
+
+ /** Backward differences array. */
+ private int[][] bdArray;
+
+ /** Original non-truncated step end time. */
+ private double nonTruncatedEnd;
+
+ /** Original non-truncated step size. */
+ private double nonTruncatedH;
+
+ /** Simple constructor.
+ * This constructor builds an instance that is not usable yet, the
+ * {@link AbstractStepInterpolator#reinitialize} method should be called
+ * before using the instance in order to initialize the internal arrays. This
+ * constructor is used only in order to delay the initialization in
+ * some cases.
+ */
+ public AdamsMoultonStepInterpolator() {
+ }
+
+ /** Copy constructor.
+ * @param interpolator interpolator to copy from. The copy is a deep
+ * copy: its arrays are separated from the original arrays of the
+ * instance
+ */
+ public AdamsMoultonStepInterpolator(final AdamsMoultonStepInterpolator interpolator) {
+ super(interpolator);
+ nonTruncatedEnd = interpolator.nonTruncatedEnd;
+ nonTruncatedH = interpolator.nonTruncatedH;
+ }
+
+ /** {@inheritDoc} */
+ protected StepInterpolator doCopy() {
+ return new AdamsMoultonStepInterpolator(this);
+ }
+
+ /** {@inheritDoc} */
+ protected void initializeCoefficients() {
+
+ neville = new double[previousF.length];
+ rollback = new double[previousF.length];
+
+ bdArray = AdamsBashforthIntegrator.computeBackwardDifferencesArray(previousF.length);
+
+ Fraction[] fGammaStar = AdamsMoultonIntegrator.computeGammaStarArray(previousF.length);
+ gammaStar = new double[fGammaStar.length];
+ for (int i = 0; i < fGammaStar.length; ++i) {
+ gammaStar[i] = fGammaStar[i].doubleValue();
+ }
+
+ }
+
+ /** {@inheritDoc} */
+ public void storeTime(final double t) {
+ nonTruncatedEnd = t;
+ nonTruncatedH = nonTruncatedEnd - previousTime;
+ super.storeTime(t);
+ }
+
+ /** Truncate a step.
+ * <p>Truncating a step is necessary when an event is triggered
+ * before the nominal end of the step.</p>
+ */
+ void truncateStep(final double truncatedEndTime) {
+ currentTime = truncatedEndTime;
+ h = currentTime - previousTime;
+ }
+
+ /** {@inheritDoc} */
+ 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);
+ }
+
+ /** {@inheritDoc} */
+ protected void computeInterpolatedState(final double theta, final double oneMinusThetaH) {
+ interpolateDerivatives();
+ interpolateState(theta);
+ }
+
+ /** 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) {
+
+ // initialize the Neville's interpolation algorithm
+ for (int k = 0; k < previousF.length; ++k) {
+ neville[k] = previousF[k][i];
+ }
+
+ // 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);
+ }
+ }
+
+ // the interpolation polynomial value is in the array last element
+ interpolatedDerivatives[i] = neville[neville.length - 1];
+
+ }
+
+ }
+
+ /** 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</sup> γ<sub>j</sub><sup>*</sup>∇<sup>j</sup>f<sub>n+1</sub>
+ * </pre>
+ * <p>In the previous expression, the γ<sub>j</sub><sup>*</sup> 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><sup>*</sup> = (-1)<sup>j</sup>∫<sub>0</sub><sup>1</sup>C<sub>j</sub><sup>1-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><sup>*</sup>(θ)=
+ * (-1)<sup>j</sup>∫<sub>θ</sub><sup>1</sup>C<sub>j</sub><sup>1-s</sup>ds
+ * </p>
+ * The method described in the Hairer, Norsett and Wanner book to compute γ<sub>j</sub><sup>*</sup>
+ * is easily extended to compute γ<sub>j</sub><sup>*</sup>(θ)=
+ * (-1)<sup>j</sup>∫<sub>0</sub><sup>θ</sup>C<sub>j</sub><sup>1-s</sup>ds. From this,
+ * we can compute η<sub>j</sub><sup>*</sup>(θ) =
+ * γ<sub>j</sub><sup>*</sup>-γ<sub>j</sub><sup>*</sup>(θ).
+ * The first few values are:</p>
+ * <table>
+ * <tr><td>j</td><td>γ<sub>j</sub><sup>*</sup></td><td>γ<sub>j</sub><sup>*</sup>(θ)</td><td>η<sub>j</sub><sup>*</sup>(θ)</td></tr>
+ * <tr><td>0</td><td>1</td><td>θ</td><td>1-θ</td></tr>
+ * <tr><td>1</td><td>-1/2</td><td>(θ<sup>2</sup>-2θ)/2</td><td>(-1+2θ-θ<sup>2</sup>)/2</td></tr>
+ * <tr><td>2</td><td>-1/12</td><td>(2θ<sup>3</sup>-3θ<sup>2</sup>)/12</td><td>(-1+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><sup>*</sup>(θ)/(m+1-j) =
+ * 1/(m+1)! ∏<sub>k=0</sub><sup>k=m</sup>(θ+k-1)
+ * </p>
+ * @param theta location of the interpolation point within the last step
+ */
+ private void interpolateState(final double theta) {
+
+ // compute the integrals to remove from the final state
+ computeRollback(previousT.length - 1, theta);
+
+ // 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];
+ }
+ 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) {
+
+ // compute the gamma star(theta) values from the recurrence relation
+ double product = theta - 1;
+ rollback[0] = theta;
+ for (int i = 1; i <= order; ++i) {
+ product *= (i - 1 + theta) / (i + 1);
+ double gStar = product;
+ for (int j = 1; j <= i; ++j) {
+ gStar -= rollback[i - j] / (j + 1);
+ }
+ rollback[i] = gStar;
+ }
+
+ // subtract it from gamma star to get eta star(theta)
+ for (int i = 0; i <= order; ++i) {
+ rollback[i] -= gammaStar[i];
+ }
+
+ // combine the eta star 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];
+ }
+ rollback[i] = f;
+ }
+
+ }
+
+ /** {@inheritDoc} */
+ public void writeExternal(final ObjectOutput out)
+ throws IOException {
+ super.writeExternal(out);
+ out.writeDouble(nonTruncatedEnd);
+ }
+
+ /** {@inheritDoc} */
+ public void readExternal(final ObjectInput in)
+ throws IOException {
+ nonTruncatedEnd = in.readDouble();
+ }
+
+}
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,325 @@
+/*
+ * 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.nonstiff;
+
+import java.util.Arrays;
+
+import org.apache.commons.math.ode.AbstractIntegrator;
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.FirstOrderIntegrator;
+import org.apache.commons.math.ode.IntegratorException;
+import org.apache.commons.math.ode.ODEIntegrator;
+import org.apache.commons.math.ode.events.CombinedEventsManager;
+import org.apache.commons.math.ode.events.EventException;
+import org.apache.commons.math.ode.events.EventHandler;
+import org.apache.commons.math.ode.events.EventState;
+import org.apache.commons.math.ode.sampling.FixedStepHandler;
+import org.apache.commons.math.ode.sampling.StepHandler;
+import org.apache.commons.math.ode.sampling.StepInterpolator;
+import org.apache.commons.math.ode.sampling.StepNormalizer;
+
+/**
+ * This class is the base class for multistep integrators for Ordinary
+ * Differential Equations.
+ *
+ * @see AdamsBashforthIntegrator
+ * @see AdamsMoultonIntegrator
+ * @see BDFIntegrator
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+public abstract class MultistepIntegrator extends AbstractIntegrator {
+
+ /** Starter integrator. */
+ private FirstOrderIntegrator starter;
+
+ /** Previous steps times. */
+ protected double[] previousT;
+
+ /** Previous steps derivatives. */
+ protected double[][] previousF;
+
+ /** Time of last detected reset. */
+ private double resetTime;
+
+ /** Prototype of the step interpolator. */
+ protected MultistepStepInterpolator prototype;
+
+ /**
+ * Build a multistep integrator with the given number of steps.
+ * <p>The default starter integrator is set to the {@link
+ * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
+ * some defaults settings.</p>
+ * @param name name of the method
+ * @param k number of steps of the multistep method
+ * (including the one being computed)
+ * @param prototype prototype of the step interpolator to use
+ */
+ protected MultistepIntegrator(final String name, final int k,
+ final MultistepStepInterpolator prototype) {
+ super(name);
+ starter = new DormandPrince853Integrator(1.0e-6, 1.0e6, 1.0e-5, 1.0e-6);
+ previousT = new double[k];
+ previousF = new double[k][];
+ this.prototype = prototype;
+ }
+
+ /**
+ * Get the starter integrator.
+ * @return starter integrator
+ */
+ public ODEIntegrator getStarterIntegrator() {
+ return starter;
+ }
+
+ /**
+ * Set the starter integrator.
+ * <p>The various step and event handlers for this starter integrator
+ * will be managed automatically by the multi-step integrator. Any
+ * user configuration for these elements will be cleared before use.</p>
+ * @param starter starter integrator
+ */
+ public void setStarterIntegrator(FirstOrderIntegrator starter) {
+ this.starter = starter;
+ }
+
+ /** Start the integration.
+ * <p>This method computes the first few steps of the multistep method,
+ * using the underlying starter integrator, ensuring the returned steps
+ * all belong to the same smooth range.</p>
+ * <p>In order to ensure smoothness, the start phase is automatically
+ * restarted when a state or derivative reset is triggered by the
+ * registered events handlers before this start phase is completed. As
+ * an example, consider integrating a differential equation from t=0
+ * to t=100 with a 4 steps method and step size equal to 0.2. If an event
+ * resets the state at t=0.5, the start phase will not end at t=0.7 with
+ * steps at [0.0, 0.2, 0.4, 0.6] but instead will end at t=1.1 with steps
+ * at [0.5, 0.7, 0.9, 1.1].</p>
+ * <p>A side effect of the need for smoothness is that an ODE triggering
+ * short period regular resets will remain in the start phase throughout
+ * the integration range if the step size or the number of steps to store
+ * are too large.</p>
+ * <p>If the start phase ends prematurely (because of some triggered event
+ * for example), then the time of latest previous steps will be set to
+ * <code>Double.NaN</code>.</p>
+ * @param n number of steps to store
+ * @param h signed step size to use for the first steps
+ * @param manager discrete events manager to use
+ * @param equations differential equations to integrate
+ * @param t0 initial time
+ * @param y state vector: contains the initial value of the state vector at t0,
+ * will be used to put the state vector at each successful step and hence
+ * contains the final value at the end of the start phase
+ * @return time of the end of the start phase
+ * @throws IntegratorException if the integrator cannot perform integration
+ * @throws DerivativeException this exception is propagated to the caller if
+ * the underlying user function triggers one
+ */
+ protected double start(final int n, final double h,
+ final CombinedEventsManager manager,
+ final FirstOrderDifferentialEquations equations,
+ final double t0, final double[] y)
+ throws DerivativeException, IntegratorException {
+
+ // clear the first steps
+ Arrays.fill(previousT, Double.NaN);
+ Arrays.fill(previousF, null);
+
+ // configure the event handlers
+ starter.clearEventHandlers();
+ for (EventState state : manager.getEventsStates()) {
+ starter.addEventHandler(new ResetCheckingWrapper(state.getEventHandler()),
+ state.getMaxCheckInterval(),
+ state.getConvergence(), state.getMaxIterationCount());
+ }
+
+ // configure the step handlers
+ starter.clearStepHandlers();
+ for (final StepHandler handler : stepHandlers) {
+ // add the user defined step handlers, filtering out the isLast indicator
+ starter.addStepHandler(new FilteringWrapper(handler));
+ }
+
+ // add one specific step handler to store the first steps
+ final StoringStepHandler store = new StoringStepHandler(n);
+ starter.addStepHandler(new StepNormalizer(h, store));
+
+ // integrate over the first few steps, ensuring no intermediate reset occurs
+ double t = t0;
+ double stopTime = Double.NaN;
+ do {
+ resetTime = Double.NaN;
+ store.restart();
+ // we overshoot by 1/10000 step the end to make sure we get don't miss the last point
+ stopTime = starter.integrate(equations, t, y, t + (n - 0.9999) * h, y);
+ if (!Double.isNaN(resetTime)) {
+ // there was an intermediate reset, we restart
+ t = resetTime;
+ }
+ } while (!Double.isNaN(resetTime));
+
+ // clear configuration
+ starter.clearEventHandlers();
+ starter.clearStepHandlers();
+
+ if (store.getFinalState() != null) {
+ System.arraycopy(store.getFinalState(), 0, y, 0, y.length);
+ }
+ return stopTime;
+
+ }
+
+ /** Rotate the previous steps arrays.
+ */
+ protected void rotatePreviousSteps() {
+ final double[] rolled = previousF[previousT.length - 1];
+ for (int k = previousF.length - 1; k > 0; --k) {
+ previousT[k] = previousT[k - 1];
+ previousF[k] = previousF[k - 1];
+ }
+ previousF[0] = rolled;
+ }
+
+ /** Event handler wrapper to check if state or derivatives have been reset. */
+ private class ResetCheckingWrapper implements EventHandler {
+
+ /** Serializable version identifier. */
+ private static final long serialVersionUID = 4922660285376467937L;
+
+ /** Wrapped event handler. */
+ private final EventHandler handler;
+
+ /** Build a new instance.
+ * @param handler event handler to wrap
+ */
+ public ResetCheckingWrapper(final EventHandler handler) {
+ this.handler = handler;
+ }
+
+ /** {@inheritDoc} */
+ public int eventOccurred(double t, double[] y) throws EventException {
+ final int action = handler.eventOccurred(t, y);
+ if ((action == RESET_DERIVATIVES) || (action == RESET_STATE)) {
+ // a singularity has been encountered
+ // we need to restart the start phase
+ resetTime = t;
+ return STOP;
+ }
+ return action;
+ }
+
+ /** {@inheritDoc} */
+ public double g(double t, double[] y) throws EventException {
+ return handler.g(t, y);
+ }
+
+ /** {@inheritDoc} */
+ public void resetState(double t, double[] y) throws EventException {
+ handler.resetState(t, y);
+ }
+
+ }
+
+ /** Step handler wrapper filtering out the isLast indicator. */
+ private class FilteringWrapper implements StepHandler {
+
+ /** Serializable version identifier. */
+ private static final long serialVersionUID = 4607975253344802232L;
+
+ /** Wrapped step handler. */
+ private final StepHandler handler;
+
+ /** Build a new instance.
+ * @param handler step handler to wrap
+ */
+ public FilteringWrapper(final StepHandler handler) {
+ this.handler = handler;
+ }
+
+ /** {@inheritDoc} */
+ public void handleStep(StepInterpolator interpolator, boolean isLast)
+ throws DerivativeException {
+ // we force the isLast indicator to false EXCEPT if some event handler triggered a stop
+ handler.handleStep(interpolator, eventsHandlersManager.stop());
+ }
+
+ /** {@inheritDoc} */
+ public boolean requiresDenseOutput() {
+ return handler.requiresDenseOutput();
+ }
+
+ /** {@inheritDoc} */
+ public void reset() {
+ handler.reset();
+ }
+
+ }
+
+ /** Specialized step handler storing the first few steps. */
+ private class StoringStepHandler implements FixedStepHandler {
+
+ /** Serializable version identifier. */
+ private static final long serialVersionUID = 4592974435520688797L;
+
+ /** Number of steps to store. */
+ private final int n;
+
+ /** Counter for already stored steps. */
+ private int count;
+
+ /** Final state. */
+ private double[] finalState;
+
+ /** Build a new instance.
+ * @param number of steps to store
+ */
+ public StoringStepHandler(final int n) {
+ this.n = n;
+ restart();
+ }
+
+ /** Restart storage.
+ */
+ public void restart() {
+ count = 0;
+ finalState = null;
+ }
+
+ /** Get the final state.
+ * @return final state
+ */
+ public double[] getFinalState() {
+ return finalState;
+ }
+
+ /** {@inheritDoc} */
+ public void handleStep(final double t, final double[] y, final double[] yDot,
+ final boolean isLast) {
+ if (count++ < n) {
+ previousT[n - count] = t;
+ previousF[n - count] = yDot.clone();
+ if (count == n) {
+ finalState = y.clone();
+ }
+ }
+ }
+
+ }
+
+}
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Added: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,165 @@
+/*
+ * 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.nonstiff;
+
+import java.io.IOException;
+import java.io.ObjectInput;
+import java.io.ObjectOutput;
+
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
+
+/** This class represents an interpolator over the last step during an
+ * ODE integration for multistep integrators.
+ *
+ * @see MultistepIntegrator
+ *
+ * @version $Revision$ $Date$
+ * @since 2.0
+ */
+
+abstract class MultistepStepInterpolator
+ extends AbstractStepInterpolator {
+
+ /** Previous steps times. */
+ protected double[] previousT;
+
+ /** Previous steps derivatives. */
+ protected double[][] previousF;
+
+ /** Simple constructor.
+ * This constructor builds an instance that is not usable yet, the
+ * {@link #reinitialize} method should be called before using the
+ * instance in order to initialize the internal arrays. This
+ * constructor is used only in order to delay the initialization in
+ * some cases. The {@link MultistepIntegrator} classe uses the
+ * prototyping design pattern to create the step interpolators by
+ * cloning an uninitialized model and latter initializing the copy.
+ */
+ protected MultistepStepInterpolator() {
+ previousT = null;
+ previousF = null;
+ }
+
+ /** Copy constructor.
+
+ * <p>The copied interpolator should have been finalized before the
+ * copy, otherwise the copy will not be able to perform correctly any
+ * interpolation and will throw a {@link NullPointerException}
+ * later. Since we don't want this constructor to throw the
+ * exceptions finalization may involve and since we don't want this
+ * method to modify the state of the copied interpolator,
+ * finalization is <strong>not</strong> done automatically, it
+ * remains under user control.</p>
+
+ * <p>The copy is a deep copy: its arrays are separated from the
+ * original arrays of the instance.</p>
+
+ * @param interpolator interpolator to copy from.
+
+ */
+ public MultistepStepInterpolator(final MultistepStepInterpolator interpolator) {
+
+ super(interpolator);
+
+ if (interpolator.currentState != null) {
+ previousT = interpolator.previousT.clone();
+ previousF = new double[interpolator.previousF.length][];
+ for (int k = 0; k < interpolator.previousF.length; ++k) {
+ previousF[k] = interpolator.previousF[k].clone();
+ }
+ initializeCoefficients();
+ } else {
+ previousT = null;
+ previousF = null;
+ }
+
+ }
+
+ /** Reinitialize the instance
+ * @param y reference to the integrator array holding the state at
+ * the end of the step
+ * @param previousT reference to the integrator array holding the times
+ * of the previous steps
+ * @param previousF reference to the integrator array holding the
+ * previous slopes
+ * @param forward integration direction indicator
+ */
+ public void reinitialize(final double[] y,
+ final double[] previousT, final double[][] previousF,
+ final boolean forward) {
+ reinitialize(y, forward);
+ this.previousT = previousT;
+ this.previousF = previousF;
+ initializeCoefficients();
+ }
+
+ /** Initialize the coefficients arrays.
+ */
+ protected abstract void initializeCoefficients();
+
+ /** {@inheritDoc} */
+ public void writeExternal(final ObjectOutput out)
+ throws IOException {
+
+ // save the state of the base class
+ writeBaseExternal(out);
+
+ // save the local attributes
+ out.writeInt(previousT.length);
+ for (int k = 0; k < previousF.length; ++k) {
+ out.writeDouble(previousT[k]);
+ for (int i = 0; i < currentState.length; ++i) {
+ out.writeDouble(previousF[k][i]);
+ }
+ }
+
+ }
+
+ /** {@inheritDoc} */
+ public void readExternal(final ObjectInput in)
+ throws IOException {
+
+ // read the base class
+ final double t = readBaseExternal(in);
+
+ // read the local attributes
+ final int kMax = in.readInt();
+ previousT = new double[kMax];
+ previousF = new double[kMax][];
+ for (int k = 0; k < kMax; ++k) {
+ previousT[k] = in.readDouble();
+ previousF[k] = new double[currentState.length];
+ for (int i = 0; i < currentState.length; ++i) {
+ previousF[k][i] = in.readDouble();
+ }
+ }
+
+ // initialize the coefficients
+ initializeCoefficients();
+
+ try {
+ // we can now set the interpolated time and state
+ setInterpolatedTime(t);
+ } catch (DerivativeException e) {
+ throw new IOException(e.getMessage());
+ }
+
+ }
+
+}
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision
Modified: commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml?rev=676737&r1=676736&r2=676737&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/branches/MATH_2_0/src/site/xdoc/changes.xml Mon Jul 14 14:13:23 2008
@@ -39,6 +39,13 @@
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
+ <action dev="luc" type="add">
+ New ODE integrators have been added: the explicit Adams-Bashforth and implicit
+ Adams-Moulton multistep methods. These methods support customizable starter
+ integrators and support discrete events even during the start phase.
+ All these methods provide the same rich features has the existing ones:
+ continuous output, step handlers, discrete events, G-stop ...
+ </action>
<action dev="luc" type="fix" issue="MATH-214" >
Replaced size adjustment of all steps of fixed steps Runge-Kutta integrators by
a truncation of the last step only.
Added: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java?rev=676737&view=auto
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java (added)
+++ commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java Mon Jul 14 14:13:23 2008
@@ -0,0 +1,206 @@
+/*
+ * 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.nonstiff;
+
+import junit.framework.Test;
+import junit.framework.TestCase;
+import junit.framework.TestSuite;
+
+import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.FirstOrderIntegrator;
+import org.apache.commons.math.ode.IntegratorException;
+import org.apache.commons.math.ode.events.EventHandler;
+
+public class AdamsBashforthIntegratorTest
+ extends TestCase {
+
+ 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,
+ 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) {
+
+ double previousError = Double.NaN;
+ for (int i = 6; i < 10; ++i) {
+
+ TestProblemAbstract pb = (TestProblemAbstract) problems[k].clone();
+ double step = (pb.getFinalTime() - pb.getInitialTime()) * Math.pow(2.0, -i);
+
+ FirstOrderIntegrator integ = new AdamsBashforthIntegrator(5, step);
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ EventHandler[] functions = pb.getEventsHandlers();
+ for (int l = 0; l < functions.length; ++l) {
+ integ.addEventHandler(functions[l],
+ Double.POSITIVE_INFINITY, 1.0e-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);
+ }
+
+ double error = handler.getMaximalValueError();
+ if (i > 6) {
+ assertTrue(error < Math.abs(previousError));
+ }
+ previousError = error;
+
+ }
+
+ }
+
+ }
+
+ public void testSmallStep()
+ throws DerivativeException, IntegratorException {
+
+ 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()]);
+
+ assertTrue(handler.getLastError() < 2.0e-9);
+ assertTrue(handler.getMaximalValueError() < 3.0e-8);
+ assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
+ assertEquals("Adams-Bashforth", integ.getName());
+
+ }
+
+ public void testBigStep()
+ throws DerivativeException, IntegratorException {
+
+ TestProblem1 pb = new TestProblem1();
+ double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
+
+ FirstOrderIntegrator integ = new AdamsBashforthIntegrator(3, step);
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.integrate(pb,
+ pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
+
+ assertTrue(handler.getLastError() > 0.05);
+ assertTrue(handler.getMaximalValueError() > 0.1);
+ assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
+
+ }
+
+ public static Test suite() {
+ return new TestSuite(AdamsBashforthIntegratorTest.class);
+ }
+
+}
Propchange: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange: commons/proper/math/branches/MATH_2_0/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
------------------------------------------------------------------------------
svn:keywords = Author Date Id Revision