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:
+ *A k-steps Adams-Bashforth method is of order k. There is no limit to the + * value of k.
+ * + *These methods are explicit: fn+1 is not used to compute + * yn+1. 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.
+ * + * @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. + *The coefficients are the ci terms in the following formula:
+ *+ * yn+1 = yn + h × ∑i=0i=k-1 cifn-i + *+ * @return a copy of the coefficients of the method + */ + public double[] getCoeffs() { + return coeffs.clone(); + } + + /** Compute the backward differences coefficients array. + *
This is quite similar to the Pascal triangle, except for a + * (-1)i 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:
+ *+ * ∇0 fn = fn + * ∇i+1 fn = ∇ifn - ∇ifn-1 + *+ * @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; + + } + +} diff --git a/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java b/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java new file mode 100644 index 000000000..0d46419dc --- /dev/null +++ b/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java @@ -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. + * + *
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 rollback of the + * integration from step end to interpolation date: + *
+ * y(tn + theta h) = y (tn + h) - ∫tn + theta htn + hp(t)dt + *+ * where theta belongs to [0 ; 1] and p(t) is the interpolation polynomial based on + * the derivatives at previous steps fn-k+1, fn-k+2 ... + * fn and fn. + * + * @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. + *
Truncating a step is necessary when an event is triggered + * before the nominal end of the step.
+ */ + 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. + *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.
+ */ + 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. + *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: Solving Ordinary + * Differential Equations (Springer-Verlag, ISBN 3-540-56670-8), this + * process leads to the following expression:
+ *+ * yn+1 = yn + + * h × ∑j=0j=k-1 γj∇jfn + *+ *
In the previous expression, the γj terms are the + * ones that result from the analytical integration, and can be computed form + * the binomial coefficients Cj-s:
+ *+ * γj = (-1)j∫01Cj-sds + *
+ *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.
+ *+ * ηj(θ)= + * (-1)j∫θ1Cj-sds + *
+ * The method described in the Hairer, Norsett and Wanner book to compute γj + * is easily extended to compute γj(θ)= + * (-1)j∫0θCj-sds. From this, + * we can compute ηj(θ) = γj-γj(θ). + * The first few values are: + *| j | γj | γj(θ) | ηj(θ) |
| 0 | 1 | θ | 1-θ |
| 1 | 1/2 | θ2/2 | (1-θ2)/2 |
| 2 | 5/12 | (3θ2+2θ3)/12 | (5-3θ2-2θ3)/12 |
+ * The ηj(θ) functions appear to be polynomial ones. As expected, + * we see that ηj(1)= 0. The recurrence relation derived for + * γj(θ) is: + *
+ *+ * &sumj=0j=mγj(θ)/(m+1-j) = + * 1/(m+1)! ∏k=0k=m(θ+k) + *
+ * @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(); + } + +} diff --git a/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java b/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java new file mode 100644 index 000000000..1e0351385 --- /dev/null +++ b/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java @@ -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. + * + *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:
+ *The coefficients are computed (and cached) as needed, so their are no + * theoretical limitations to the number of steps
+ * + *A k-steps Adams-Moulton method is of order k+1. There is no limit to the + * value of k.
+ * + *These methods are implicit: fn+1 is used to compute + * yn+1. 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.
+ * + * @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. + *The coefficients are the ci terms in the following formula:
+ *+ * yn+1 = yn + h × ∑i=0i=k-1 cifn-i + *+ * @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. + *
The coefficients are the ci terms in the following formula:
+ *+ * yn+1 = yn + h × ∑i=0i=k cifn-i + *+ * @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; + + } + +} diff --git a/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java b/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java new file mode 100644 index 000000000..35ed7d98a --- /dev/null +++ b/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java @@ -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. + * + *
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 rollback of the + * integration from step end to interpolation date: + *
+ * y(tn + theta h) = y (tn + h) - ∫tn + theta htn + hp(t)dt + *+ * where theta belongs to [0 ; 1] and p(t) is the interpolation polynomial based on + * the derivatives at previous steps fn-k+1, fn-k+2 ... + * fn, fn and fn+1. + * + * @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. + *
Truncating a step is necessary when an event is triggered + * before the nominal end of the step.
+ */ + 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. + *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.
+ */ + 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. + *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: Solving Ordinary + * Differential Equations (Springer-Verlag, ISBN 3-540-56670-8), this + * process leads to the following expression:
+ *+ * yn+1 = yn + + * h × ∑j=0j=k γj*∇jfn+1 + *+ *
In the previous expression, the γj* terms are the + * ones that result from the analytical integration, and can be computed form + * the binomial coefficients Cj-s:
+ *+ * γj* = (-1)j∫01Cj1-sds + *
+ *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.
+ *+ * ηj*(θ)= + * (-1)j∫θ1Cj1-sds + *
+ * The method described in the Hairer, Norsett and Wanner book to compute γj* + * is easily extended to compute γj*(θ)= + * (-1)j∫0θCj1-sds. From this, + * we can compute ηj*(θ) = + * γj*-γj*(θ). + * The first few values are: + *| j | γj* | γj*(θ) | ηj*(θ) |
| 0 | 1 | θ | 1-θ |
| 1 | -1/2 | (θ2-2θ)/2 | (-1+2θ-θ2)/2 |
| 2 | -1/12 | (2θ3-3θ2)/12 | (-1+3θ2-2θ3)/12 |
+ * The ηj(θ) functions appear to be polynomial ones. As expected, + * we see that ηj(1)= 0. The recurrence relation derived for + * γj(θ) is: + *
+ *+ * &sumj=0j=mγj*(θ)/(m+1-j) = + * 1/(m+1)! ∏k=0k=m(θ+k-1) + *
+ * @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(); + } + +} diff --git a/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java b/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java new file mode 100644 index 000000000..4991e8a9d --- /dev/null +++ b/src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java @@ -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. + *The default starter integrator is set to the {@link + * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with + * some defaults settings.
+ * @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. + *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.
+ * @param starter starter integrator + */ + public void setStarterIntegrator(FirstOrderIntegrator starter) { + this.starter = starter; + } + + /** Start the integration. + *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.
+ *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].
+ *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.
+ *If the start phase ends prematurely (because of some triggered event
+ * for example), then the time of latest previous steps will be set to
+ * Double.NaN.
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 not done automatically, it + * remains under user control.
+ + *The copy is a deep copy: its arrays are separated from the + * original arrays of the instance.
+ + * @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()); + } + + } + +} diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index d6a5c251c..c6e85fc67 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -39,6 +39,13 @@ The