Return-Path: Delivered-To: apmail-commons-commits-archive@locus.apache.org Received: (qmail 67986 invoked from network); 14 Jul 2008 21:13:49 -0000 Received: from hermes.apache.org (HELO mail.apache.org) (140.211.11.2) by minotaur.apache.org with SMTP; 14 Jul 2008 21:13:49 -0000 Received: (qmail 35269 invoked by uid 500); 14 Jul 2008 21:13:47 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 35207 invoked by uid 500); 14 Jul 2008 21:13:46 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 35191 invoked by uid 99); 14 Jul 2008 21:13:46 -0000 Received: from athena.apache.org (HELO athena.apache.org) (140.211.11.136) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 14 Jul 2008 14:13:46 -0700 X-ASF-Spam-Status: No, hits=-2000.0 required=10.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 14 Jul 2008 21:13:00 +0000 Received: by eris.apache.org (Postfix, from userid 65534) id 1A58723889F1; Mon, 14 Jul 2008 14:13:24 -0700 (PDT) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: 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/ Date: Mon, 14 Jul 2008 21:13:23 -0000 To: commits@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.8 Message-Id: <20080714211324.1A58723889F1@eris.apache.org> X-Virus-Checked: Checked by ClamAV on apache.org 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. + * + *

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:

+ *
    + *
  • k = 1: yn+1 = yn + h fn
    • + *
  • k = 2: yn+1 = yn + h (3fn-fn-1)/2
    • + *
  • k = 3: yn+1 = yn + h (23fn-16fn-1+5fn-2)/12
    • + *
  • k = 4: yn+1 = yn + h (55fn-59fn-1+37fn-2-9fn-3)/24
    • + *
  • ...
    • + *
+ * + *

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; + + } + +} 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. + * + *

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:

+ * + * + * θ + * θ2/2 + * (3θ2+2θ3)/12 + *
jγjγj(θ)ηj(θ)
011-θ
11/2(1-θ2)/2
25/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(); + } + +} 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. + * + *

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:

+ *
    + *
  • k = 0: yn+1 = yn + h fn+1
    • + *
  • k = 1: yn+1 = yn + h (fn+1+fn)/2
    • + *
  • k = 2: yn+1 = yn + h (5fn+1+8fn-fn-1)/12
    • + *
  • k = 3: yn+1 = yn + h (9fn+1+19fn-5fn-1+fn-2)/24
    • + *
  • ...
    • + *
+ * + *

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; + + } + +} 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. + * + *

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*(θ)
01θ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(); + } + +} 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. + *

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.

+ * @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. + + *

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()); + } + + } + +} 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 @@ + + 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 ... + 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