Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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 ...


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_0@676737 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-07-14 21:13:23 +00:00
parent 4fb3ee4e70
commit 8365d33be2
9 changed files with 2080 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

285
src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Normal file
View File

@ -0,0 +1,285 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.fraction.Fraction;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.sampling.StepHandler;
/**
* This class implements explicit Adams-Bashforth integrators for Ordinary
* Differential Equations.
*
* <p>Adams-Bashforth (in fact due to Adams alone) methods are explicit
* multistep ODE solvers witch fixed stepsize. The value of state vector
* at step n+1 is a simple combination of the value at step n and of the
* derivatives at steps n, n-1, n-2 ... Depending on the number k of previous
* steps one wants to use for computing the next value, different formulas
* are available:</p>
* <ul>
* <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + h f<sub>n</sub></li>
* <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + h (3f<sub>n</sub>-f<sub>n-1</sub>)/2</li>
* <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + h (23f<sub>n</sub>-16f<sub>n-1</sub>+5f<sub>n-2</sub>)/12</li>
* <li>k = 4: y<sub>n+1</sub> = y<sub>n</sub> + h (55f<sub>n</sub>-59f<sub>n-1</sub>+37f<sub>n-2</sub>-9f<sub>n-3)/24</sub></li>
* <li>...</li>
* </ul>
*
* <p>A k-steps Adams-Bashforth method is of order k. There is no limit to the
* value of k.</p>
*
* <p>These methods are explicit: f<sub>n+1</sub> is not used to compute
* y<sub>n+1</sub>. More accurate implicit Adams methods exist: the
* Adams-Moulton methods (which are also due to Adams alone). They are
* provided by the {@link AdamsMoultonIntegrator AdamsMoultonIntegrator} class.</p>
*
* @see AdamsMoultonIntegrator
* @see BDFIntegrator
* @version $Revision$ $Date$
* @since 2.0
*/
public class AdamsBashforthIntegrator extends MultistepIntegrator {
/** Serializable version identifier. */
private static final long serialVersionUID = 1676381657635800870L;
/** Integrator method name. */
private static final String METHOD_NAME = "Adams-Bashforth";
/** Coefficients for the current method. */
private final double[] coeffs;
/** Integration step. */
private final double step;
/**
* Build an Adams-Bashforth integrator with the given order and step size.
* @param order order of the method (must be strictly positive)
* @param step integration step size
*/
public AdamsBashforthIntegrator(final int order, final double step) {
super(METHOD_NAME, order, new AdamsBashforthStepInterpolator());
// compute the integration coefficients
int[][] bdArray = computeBackwardDifferencesArray(order);
Fraction[] gamma = computeGammaArray(order);
coeffs = new double[order];
for (int i = 0; i < order; ++i) {
Fraction f = Fraction.ZERO;
for (int j = i; j < order; ++j) {
f = f.add(gamma[j].multiply(new Fraction(bdArray[j][i], 1)));
}
coeffs[i] = f.doubleValue();
}
this.step = step;
}
/** {@inheritDoc} */
public double integrate(FirstOrderDifferentialEquations equations,
double t0, double[] y0, double t, double[] y)
throws DerivativeException, IntegratorException {
sanityChecks(equations, t0, y0, t, y);
final boolean forward = (t > t0);
// initialize working arrays
if (y != y0) {
System.arraycopy(y0, 0, y, 0, y0.length);
}
final double[] yTmp = new double[y0.length];
// set up an interpolator sharing the integrator arrays
final AdamsBashforthStepInterpolator interpolator =
(AdamsBashforthStepInterpolator) prototype.copy();
interpolator.reinitialize(yTmp, previousT, previousF, forward);
// set up integration control objects
stepStart = t0;
stepSize = step;
for (StepHandler handler : stepHandlers) {
handler.reset();
}
CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
// compute the first few steps using the configured starter integrator
double stopTime =
start(previousF.length, stepSize, manager, equations, stepStart, y);
if (Double.isNaN(previousT[0])) {
return stopTime;
}
stepStart = previousT[0];
interpolator.storeTime(stepStart);
boolean lastStep = false;
while (!lastStep) {
// shift all data
interpolator.shift();
// estimate the state at the end of the step
for (int j = 0; j < y0.length; ++j) {
double sum = 0;
for (int l = 0; l < coeffs.length; ++l) {
sum += coeffs[l] * previousF[l][j];
}
yTmp[j] = y[j] + stepSize * sum;
}
// discrete events handling
interpolator.storeTime(stepStart + stepSize);
final boolean truncated;
if (manager.evaluateStep(interpolator)) {
truncated = true;
interpolator.truncateStep(manager.getEventTime());
} else {
truncated = false;
}
// the step has been accepted (may have been truncated)
final double nextStep = interpolator.getCurrentTime();
interpolator.setInterpolatedTime(nextStep);
System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
manager.stepAccepted(nextStep, y);
lastStep = manager.stop();
// provide the step data to the step handler
for (StepHandler handler : stepHandlers) {
handler.handleStep(interpolator, lastStep);
}
stepStart = nextStep;
if (!lastStep) {
// prepare next step
if (manager.reset(stepStart, y)) {
// some events handler has triggered changes that
// invalidate the derivatives, we need to restart from scratch
stopTime =
start(previousF.length, stepSize, manager, equations, stepStart, y);
if (Double.isNaN(previousT[0])) {
return stopTime;
}
stepStart = previousT[0];
} else {
if (truncated) {
// the step has been truncated, we need to adjust the previous steps
for (int i = 1; i < previousF.length; ++i) {
previousT[i] = stepStart - i * stepSize;
interpolator.setInterpolatedTime(previousT[i]);
System.arraycopy(interpolator.getInterpolatedState(), 0,
previousF[i], 0, y0.length);
}
} else {
rotatePreviousSteps();
}
// evaluate differential equations for next step
previousT[0] = stepStart;
equations.computeDerivatives(stepStart, y, previousF[0]);
}
}
}
stopTime = stepStart;
stepStart = Double.NaN;
stepSize = Double.NaN;
return stopTime;
}
/** Get the coefficients of the method.
* <p>The coefficients are the c<sub>i</sub> terms in the following formula:</p>
* <pre>
* y<sub>n+1</sub> = y<sub>n</sub> + h &times; &sum;<sub>i=0</sub><sup>i=k-1</sup> c<sub>i</sub>f<sub>n-i</sub></li>
* </pre>
* @return a copy of the coefficients of the method
*/
public double[] getCoeffs() {
return coeffs.clone();
}
/** Compute the backward differences coefficients array.
* <p>This is quite similar to the Pascal triangle, except for a
* (-1)<sup>i</sup> sign. We use a straightforward approach here,
* since we don't expect this to be run too many times with too
* high k. It is based on the recurrence relations:</p>
* <pre>
* &nabla;<sup>0</sup> f<sub>n</sub> = f<sub>n</sub>
* &nabla;<sup>i+1</sup> f<sub>n</sub> = &nabla;<sup>i</sup>f<sub>n</sub> - &nabla;<sup>i</sup>f<sub>n-1</sub>
* </pre>
* @param order order of the integration method
*/
static int[][] computeBackwardDifferencesArray(final int order) {
// create the array
int[][] bdArray = new int[order][];
// recurrence initialization
bdArray[0] = new int[] { 1 };
// fill up array using recurrence relation
for (int i = 1; i < order; ++i) {
bdArray[i] = new int[i + 1];
bdArray[i][0] = 1;
for (int j = 0; j < i - 1; ++j) {
bdArray[i][j + 1] = bdArray[i - 1][j + 1] - bdArray[i - 1][j];
}
bdArray[i][i] = -bdArray[i - 1][i - 1];
}
return bdArray;
}
/** Compute the gamma coefficients.
* @param order order of the integration method
* @return gamma coefficients array
*/
static Fraction[] computeGammaArray(final int order) {
// create the array
Fraction[] gammaArray = new Fraction[order];
// recurrence initialization
gammaArray[0] = Fraction.ONE;
// fill up array using recurrence relation
for (int i = 1; i < order; ++i) {
Fraction gamma = Fraction.ONE;
for (int j = 1; j <= i; ++j) {
gamma = gamma.subtract(gammaArray[i - j].multiply(new Fraction(1, j + 1)));
}
gammaArray[i] = gamma;
}
return gammaArray;
}
}

288
src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthStepInterpolator.java Normal file
View File

@ -0,0 +1,288 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.ode.nonstiff;
import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math.fraction.Fraction;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.StepInterpolator;
/**
* This class implements an interpolator for Adams-Bashforth multiple steps.
*
* <p>This interpolator computes dense output inside the last few
* steps computed. The interpolation equation is consistent with the
* integration scheme, it is based on a kind of <em>rollback</em> of the
* integration from step end to interpolation date:
* <pre>
* y(t<sub>n</sub> + theta h) = y (t<sub>n</sub> + h) - &int;<sub>t<sub>n</sub> + theta h</sub><sup>t<sub>n</sub> + h</sup>p(t)dt
* </pre>
* where theta belongs to [0 ; 1] and p(t) is the interpolation polynomial based on
* the derivatives at previous steps f<sub>n-k+1</sub>, f<sub>n-k+2</sub> ...
* f<sub>n</sub> and f<sub>n</sub>.</p>
*
* @see AdamsBashforthIntegrator
* @version $Revision$ $Date$
* @since 2.0
*/
class AdamsBashforthStepInterpolator extends MultistepStepInterpolator {
/** Serializable version identifier */
private static final long serialVersionUID = -7179861704951334960L;
/** Neville's interpolation array. */
private double[] neville;
/** Integration rollback array. */
private double[] rollback;
/** &gamma; array. */
private double[] gamma;
/** Backward differences array. */
private int[][] bdArray;
/** Original non-truncated step end time. */
private double nonTruncatedEnd;
/** Original non-truncated step size. */
private double nonTruncatedH;
/** Simple constructor.
* This constructor builds an instance that is not usable yet, the
* {@link AbstractStepInterpolator#reinitialize} method should be called
* before using the instance in order to initialize the internal arrays. This
* constructor is used only in order to delay the initialization in
* some cases.
*/
public AdamsBashforthStepInterpolator() {
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
public AdamsBashforthStepInterpolator(final AdamsBashforthStepInterpolator interpolator) {
super(interpolator);
nonTruncatedEnd = interpolator.nonTruncatedEnd;
nonTruncatedH = interpolator.nonTruncatedH;
}
/** {@inheritDoc} */
protected StepInterpolator doCopy() {
return new AdamsBashforthStepInterpolator(this);
}
/** {@inheritDoc} */
protected void initializeCoefficients() {
neville = new double[previousF.length];
rollback = new double[previousF.length];
bdArray = AdamsBashforthIntegrator.computeBackwardDifferencesArray(previousF.length);
Fraction[] fGamma = AdamsBashforthIntegrator.computeGammaArray(previousF.length);
gamma = new double[fGamma.length];
for (int i = 0; i < fGamma.length; ++i) {
gamma[i] = fGamma[i].doubleValue();
}
}
/** {@inheritDoc} */
public void storeTime(final double t) {
nonTruncatedEnd = t;
nonTruncatedH = nonTruncatedEnd - previousTime;
super.storeTime(t);
}
/** Truncate a step.
* <p>Truncating a step is necessary when an event is triggered
* before the nominal end of the step.</p>
*/
void truncateStep(final double truncatedEndTime) {
currentTime = truncatedEndTime;
h = currentTime - previousTime;
}
/** {@inheritDoc} */
public void setInterpolatedTime(final double time)
throws DerivativeException {
interpolatedTime = time;
final double oneMinusThetaH = nonTruncatedEnd - interpolatedTime;
final double theta = (nonTruncatedH == 0) ?
0 : (nonTruncatedH - oneMinusThetaH) / nonTruncatedH;
computeInterpolatedState(theta, oneMinusThetaH);
}
/** {@inheritDoc} */
protected void computeInterpolatedState(final double theta, final double oneMinusThetaH) {
interpolateDerivatives();
interpolateState(theta);
}
/** Interpolate the derivatives.
* <p>The Adams method is based on a polynomial interpolation of the
* derivatives based on the preceding steps. So the interpolation of
* the derivatives here is strictly equivalent: it is a simple polynomial
* interpolation.</p>
*/
private void interpolateDerivatives() {
for (int i = 0; i < interpolatedDerivatives.length; ++i) {
// initialize the Neville's interpolation algorithm
for (int k = 0; k < previousF.length; ++k) {
neville[k] = previousF[k][i];
}
// combine the contributions of each points
for (int l = 1; l < neville.length; ++l) {
for (int m = neville.length - 1; m >= l; --m) {
final double xm = previousT[m];
final double xmMl = previousT[m - l];
neville[m] = ((interpolatedTime - xm) * neville[m-1] +
(xmMl - interpolatedTime) * neville[m]) / (xmMl - xm);
}
}
// the interpolation polynomial value is in the array last element
interpolatedDerivatives[i] = neville[neville.length - 1];
}
}
/** Interpolate the state.
* <p>The Adams method is based on a polynomial interpolation of the
* derivatives based on the preceding steps. The polynomial model is
* integrated analytically throughout the last step. Using the notations
* found in the second edition of the first volume (Nonstiff Problems)
* of the reference book by Hairer, Norsett and Wanner: <i>Solving Ordinary
* Differential Equations</i> (Springer-Verlag, ISBN 3-540-56670-8), this
* process leads to the following expression:</p>
* <pre>
* y<sub>n+1</sub> = y<sub>n</sub> +
* h &times; &sum;<sub>j=0</sub><sup>j=k-1</sup> &gamma;<sub>j</sub>&nabla;<sup>j</sup>f<sub>n</sub>
* </pre>
* <p>In the previous expression, the &gamma;<sub>j</sub> terms are the
* ones that result from the analytical integration, and can be computed form
* the binomial coefficients C<sub>j</sub><sup>-s</sup>:</p>
* <p>
* &gamma;<sub>j</sub> = (-1)<sup>j</sup>&int;<sub>0</sub><sup>1</sup>C<sub>j</sub><sup>-s</sup>ds
* </p>
* <p>In order to interpolate the state in a manner that is consistent with the
* integration scheme, we simply subtract from the current state (at the end of the step)
* the integral computed from interpolation time to step end time.</p>
* <p>
* &eta;<sub>j</sub>(&theta;)=
* (-1)<sup>j</sup>&int;<sub>&theta;</sub><sup>1</sup>C<sub>j</sub><sup>-s</sup>ds
* </p>
* The method described in the Hairer, Norsett and Wanner book to compute &gamma;<sub>j</sub>
* is easily extended to compute &gamma;<sub>j</sub>(&theta;)=
* (-1)<sup>j</sup>&int;<sub>0</sub><sup>&theta;</sup>C<sub>j</sub><sup>-s</sup>ds. From this,
* we can compute &eta;<sub>j</sub>(&theta;) = &gamma;<sub>j</sub>-&gamma;<sub>j</sub>(&theta;).
* The first few values are:</p>
* <table>
* <tr><td>j</td><td>&gamma;<sub>j</sub></td><td>&gamma;<sub>j</sub>(&theta;)</td><td>&eta;<sub>j</sub>(&theta;)</td></tr>
* <tr><td>0</td><td>1</td><td></td>&theta;<td>1-&theta;</td></tr>
* <tr><td>1</td><td>1/2</td><td></td>&theta;<sup>2</sup>/2<td>(1-&theta;<sup>2</sup>)/2</td></tr>
* <tr><td>2</td><td>5/12</td><td></td>(3&theta;<sup>2</sup>+2&theta;<sup>3</sup>)/12<td>(5-3&theta;<sup>2</sup>-2&theta;<sup>3</sup>)/12</td></tr>
* </table>
* <p>
* The &eta;<sub>j</sub>(&theta;) functions appear to be polynomial ones. As expected,
* we see that &eta;<sub>j</sub>(1)= 0. The recurrence relation derived for
* &gamma;<sub>j</sub>(&theta;) is:
* </p>
* <p>
* &sum<sub>j=0</sub><sup>j=m</sup>&gamma;<sub>j</sub>(&theta;)/(m+1-j) =
* 1/(m+1)! &prod;<sub>k=0</sub><sup>k=m</sup>(&theta;+k)
* </p>
* @param theta location of the interpolation point within the last step
*/
private void interpolateState(final double theta) {
// compute the integrals to remove from the final state
computeRollback(previousT.length - 1, theta);
// remove these integrals from the final state
for (int j = 0; j < interpolatedState.length; ++j) {
double sum = 0;
for (int l = 0; l < previousT.length; ++l) {
sum += rollback[l] * previousF[l][j];
}
interpolatedState[j] = currentState[j] - h * sum;
}
}
/** Compute the rollback coefficients.
* @param order order of the integration method
* @param theta current value for theta
*/
private void computeRollback(final int order, final double theta) {
// compute the gamma(theta) values from the recurrence relation
double product = theta;
rollback[0] = theta;
for (int i = 1; i < order; ++i) {
product *= (i + theta) / (i + 1);
double g = product;
for (int j = 1; j <= i; ++j) {
g -= rollback[i - j] / (j + 1);
}
rollback[i] = g;
}
// subtract it from gamma to get eta(theta)
for (int i = 0; i < order; ++i) {
rollback[i] -= gamma[i];
}
// combine the eta integrals with the backward differences array
// to get the rollback coefficients
for (int i = 0; i < order; ++i) {
double f = 0;
for (int j = i; j <= order; ++j) {
f -= rollback[j] * bdArray[j][i];
}
rollback[i] = f;
}
}
/** {@inheritDoc} */
public void writeExternal(final ObjectOutput out)
throws IOException {
super.writeExternal(out);
out.writeDouble(nonTruncatedEnd);
}
/** {@inheritDoc} */
public void readExternal(final ObjectInput in)
throws IOException {
nonTruncatedEnd = in.readDouble();
}
}

299
src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Normal file
View File

@ -0,0 +1,299 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.fraction.Fraction;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.sampling.StepHandler;
/**
* This class implements implicit Adams-Moulton integrators for Ordinary
* Differential Equations.
*
* <p>Adams-Moulton (in fact due to Adams alone) methods are implicit
* multistep ODE solvers witch fixed stepsize. The value of state vector
* at step n+1 is a simple combination of the value at step n and of the
* derivatives at steps n+1, n, n-1 ... Depending on the number k of previous
* steps one wants to use for computing the next value, different formulas
* are available:</p>
* <ul>
* <li>k = 0: y<sub>n+1</sub> = y<sub>n</sub> + h f<sub>n+1</sub></li>
* <li>k = 1: y<sub>n+1</sub> = y<sub>n</sub> + h (f<sub>n+1</sub>+f<sub>n</sub>)/2</li>
* <li>k = 2: y<sub>n+1</sub> = y<sub>n</sub> + h (5f<sub>n+1</sub>+8f<sub>n</sub>-f<sub>n-1</sub>)/12</li>
* <li>k = 3: y<sub>n+1</sub> = y<sub>n</sub> + h (9f<sub>n+1</sub>+19f<sub>n</sub>-5f<sub>n-1</sub>+f<sub>n-2)/24</sub></li>
* <li>...</li>
* </ul>
*
* <p>The coefficients are computed (and cached) as needed, so their are no
* theoretical limitations to the number of steps</p>
*
* <p>A k-steps Adams-Moulton method is of order k+1. There is no limit to the
* value of k.</p>
*
* <p>These methods are implicit: f<sub>n+1</sub> is used to compute
* y<sub>n+1</sub>. Simpler explicit Adams methods exist: the
* Adams-Bashforth methods (which are also due to Adams alone). They are
* provided by the {@link AdamsBashforthIntegrator AdamsBashforthIntegrator} class.</p>
*
* @see AdamsBashforthIntegrator
* @see BDFIntegrator
* @version $Revision$ $Date$
* @since 2.0
*/
public class AdamsMoultonIntegrator extends MultistepIntegrator {
/** Serializable version identifier. */
private static final long serialVersionUID = 4990335331377040417L;
/** Integrator method name. */
private static final String METHOD_NAME = "Adams-Moulton";
/** Coefficients for the predictor phase of the method. */
private final double[] predictorCoeffs;
/** Coefficients for the corrector phase of the method. */
private final double[] correctorCoeffs;
/** Integration step. */
private final double step;
/**
* Build an Adams-Moulton integrator with the given order and step size.
* @param order order of the method (must be strictly positive)
* @param step integration step size
*/
public AdamsMoultonIntegrator(final int order, final double step) {
super(METHOD_NAME, order + 1, new AdamsMoultonStepInterpolator());
// compute the integration coefficients
int[][] bdArray = AdamsBashforthIntegrator.computeBackwardDifferencesArray(order + 1);
Fraction[] gamma = AdamsBashforthIntegrator.computeGammaArray(order);
predictorCoeffs = new double[order];
for (int i = 0; i < order; ++i) {
Fraction fPredictor = Fraction.ZERO;
for (int j = i; j < order; ++j) {
Fraction f = new Fraction(bdArray[j][i], 1);
fPredictor = fPredictor.add(gamma[j].multiply(f));
}
predictorCoeffs[i] = fPredictor.doubleValue();
}
Fraction[] gammaStar = computeGammaStarArray(order);
correctorCoeffs = new double[order + 1];
for (int i = 0; i <= order; ++i) {
Fraction fCorrector = Fraction.ZERO;
for (int j = i; j <= order; ++j) {
Fraction f = new Fraction(bdArray[j][i], 1);
fCorrector = fCorrector.add(gammaStar[j].multiply(f));
}
correctorCoeffs[i] = fCorrector.doubleValue();
}
this.step = step;
}
/** {@inheritDoc} */
public double integrate(FirstOrderDifferentialEquations equations,
double t0, double[] y0, double t, double[] y)
throws DerivativeException, IntegratorException {
sanityChecks(equations, t0, y0, t, y);
final boolean forward = (t > t0);
// initialize working arrays
if (y != y0) {
System.arraycopy(y0, 0, y, 0, y0.length);
}
final double[] yTmp = new double[y0.length];
// set up an interpolator sharing the integrator arrays
final AdamsMoultonStepInterpolator interpolator =
(AdamsMoultonStepInterpolator) prototype.copy();
interpolator.reinitialize(yTmp, previousT, previousF, forward);
// set up integration control objects
stepStart = t0;
stepSize = step;
for (StepHandler handler : stepHandlers) {
handler.reset();
}
CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
// compute the first few steps using the configured starter integrator
double stopTime =
start(previousF.length - 1, stepSize, manager, equations, stepStart, y);
if (Double.isNaN(previousT[0])) {
return stopTime;
}
stepStart = previousT[0];
rotatePreviousSteps();
previousF[0] = new double[y0.length];
interpolator.storeTime(stepStart);
boolean lastStep = false;
while (!lastStep) {
// shift all data
interpolator.shift();
// predict state at end of step
for (int j = 0; j < y0.length; ++j) {
double sum = 0;
for (int l = 0; l < predictorCoeffs.length; ++l) {
sum += predictorCoeffs[l] * previousF[l+1][j];
}
yTmp[j] = y[j] + stepSize * sum;
}
// evaluate the derivatives
final double stepEnd = stepStart + stepSize;
equations.computeDerivatives(stepEnd, yTmp, previousF[0]);
// apply corrector
for (int j = 0; j < y0.length; ++j) {
double sum = 0;
for (int l = 0; l < correctorCoeffs.length; ++l) {
sum += correctorCoeffs[l] * previousF[l][j];
}
yTmp[j] = y[j] + stepSize * sum;
}
// discrete events handling
interpolator.storeTime(stepEnd);
final boolean truncated;
if (manager.evaluateStep(interpolator)) {
truncated = true;
interpolator.truncateStep(manager.getEventTime());
} else {
truncated = false;
}
// the step has been accepted (may have been truncated)
final double nextStep = interpolator.getCurrentTime();
interpolator.setInterpolatedTime(nextStep);
System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
manager.stepAccepted(nextStep, y);
lastStep = manager.stop();
// provide the step data to the step handler
for (StepHandler handler : stepHandlers) {
handler.handleStep(interpolator, lastStep);
}
stepStart = nextStep;
if (!lastStep) {
// prepare next step
if (manager.reset(stepStart, y)) {
// some events handler has triggered changes that
// invalidate the derivatives, we need to restart from scratch
stopTime =
start(previousF.length - 1, stepSize, manager, equations, stepStart, y);
if (Double.isNaN(previousT[0])) {
return stopTime;
}
stepStart = previousT[0];
rotatePreviousSteps();
previousF[0] = new double[y0.length];
} else {
if (truncated) {
// the step has been truncated, we need to adjust the previous steps
for (int i = 1; i < previousF.length; ++i) {
previousT[i] = stepStart - i * stepSize;
interpolator.setInterpolatedTime(previousT[i]);
System.arraycopy(interpolator.getInterpolatedState(), 0,
previousF[i], 0, y0.length);
}
} else {
rotatePreviousSteps();
}
// evaluate differential equations for next step
previousT[0] = stepStart;
equations.computeDerivatives(stepStart, y, previousF[0]);
}
}
}
stopTime = stepStart;
stepStart = Double.NaN;
stepSize = Double.NaN;
return stopTime;
}
/** Get the coefficients of the prdictor phase of the method.
* <p>The coefficients are the c<sub>i</sub> terms in the following formula:</p>
* <pre>
* y<sub>n+1</sub> = y<sub>n</sub> + h &times; &sum;<sub>i=0</sub><sup>i=k-1</sup> c<sub>i</sub>f<sub>n-i</sub></li>
* </pre>
* @return a copy of the coefficients of the method
*/
public double[] getPredictorCoeffs() {
return predictorCoeffs.clone();
}
/** Get the coefficients of the corrector phase of the method.
* <p>The coefficients are the c<sub>i</sub> terms in the following formula:</p>
* <pre>
* y<sub>n+1</sub> = y<sub>n</sub> + h &times; &sum;<sub>i=0</sub><sup>i=k</sup> c<sub>i</sub>f<sub>n-i</sub></li>
* </pre>
* @return a copy of the coefficients of the method
*/
public double[] getCorrectorCoeffs() {
return correctorCoeffs.clone();
}
/** Compute the gamma star coefficients.
* @param order order of the integration method
* @return gamma star coefficients array
*/
static Fraction[] computeGammaStarArray(final int order) {
// create the array
Fraction[] gammaStarArray = new Fraction[order + 1];
// recurrence initialization
gammaStarArray[0] = Fraction.ONE;
// fill up array using recurrence relation
for (int i = 1; i <= order; ++i) {
Fraction gammaStar = Fraction.ZERO;
for (int j = 1; j <= i; ++j) {
gammaStar = gammaStar.subtract(gammaStarArray[i - j].multiply(new Fraction(1, j + 1)));
}
gammaStarArray[i] = gammaStar;
}
return gammaStarArray;
}
}

289
src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonStepInterpolator.java Normal file
View File

@ -0,0 +1,289 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.ode.nonstiff;
import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math.fraction.Fraction;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.StepInterpolator;
/**
* This class implements an interpolator for Adams-Moulton multiple steps.
*
* <p>This interpolator computes dense output inside the last few
* steps computed. The interpolation equation is consistent with the
* integration scheme, it is based on a kind of <em>rollback</em> of the
* integration from step end to interpolation date:
* <pre>
* y(t<sub>n</sub> + theta h) = y (t<sub>n</sub> + h) - &int;<sub>t<sub>n</sub> + theta h</sub><sup>t<sub>n</sub> + h</sup>p(t)dt
* </pre>
* where theta belongs to [0 ; 1] and p(t) is the interpolation polynomial based on
* the derivatives at previous steps f<sub>n-k+1</sub>, f<sub>n-k+2</sub> ...
* f<sub>n</sub>, f<sub>n</sub> and f<sub>n+1</sub>.</p>
*
* @see AdamsMoultonIntegrator
* @version $Revision$ $Date$
* @since 2.0
*/
class AdamsMoultonStepInterpolator extends MultistepStepInterpolator {
/** Serializable version identifier */
private static final long serialVersionUID = 735568489801241899L;
/** Neville's interpolation array. */
private double[] neville;
/** Integration rollback array. */
private double[] rollback;
/** &gamma; star array. */
private double[] gammaStar;
/** Backward differences array. */
private int[][] bdArray;
/** Original non-truncated step end time. */
private double nonTruncatedEnd;
/** Original non-truncated step size. */
private double nonTruncatedH;
/** Simple constructor.
* This constructor builds an instance that is not usable yet, the
* {@link AbstractStepInterpolator#reinitialize} method should be called
* before using the instance in order to initialize the internal arrays. This
* constructor is used only in order to delay the initialization in
* some cases.
*/
public AdamsMoultonStepInterpolator() {
}
/** Copy constructor.
* @param interpolator interpolator to copy from. The copy is a deep
* copy: its arrays are separated from the original arrays of the
* instance
*/
public AdamsMoultonStepInterpolator(final AdamsMoultonStepInterpolator interpolator) {
super(interpolator);
nonTruncatedEnd = interpolator.nonTruncatedEnd;
nonTruncatedH = interpolator.nonTruncatedH;
}
/** {@inheritDoc} */
protected StepInterpolator doCopy() {
return new AdamsMoultonStepInterpolator(this);
}
/** {@inheritDoc} */
protected void initializeCoefficients() {
neville = new double[previousF.length];
rollback = new double[previousF.length];
bdArray = AdamsBashforthIntegrator.computeBackwardDifferencesArray(previousF.length);
Fraction[] fGammaStar = AdamsMoultonIntegrator.computeGammaStarArray(previousF.length);
gammaStar = new double[fGammaStar.length];
for (int i = 0; i < fGammaStar.length; ++i) {
gammaStar[i] = fGammaStar[i].doubleValue();
}
}
/** {@inheritDoc} */
public void storeTime(final double t) {
nonTruncatedEnd = t;
nonTruncatedH = nonTruncatedEnd - previousTime;
super.storeTime(t);
}
/** Truncate a step.
* <p>Truncating a step is necessary when an event is triggered
* before the nominal end of the step.</p>
*/
void truncateStep(final double truncatedEndTime) {
currentTime = truncatedEndTime;
h = currentTime - previousTime;
}
/** {@inheritDoc} */
public void setInterpolatedTime(final double time)
throws DerivativeException {
interpolatedTime = time;
final double oneMinusThetaH = nonTruncatedEnd - interpolatedTime;
final double theta = (nonTruncatedH == 0) ?
0 : (nonTruncatedH - oneMinusThetaH) / nonTruncatedH;
computeInterpolatedState(theta, oneMinusThetaH);
}
/** {@inheritDoc} */
protected void computeInterpolatedState(final double theta, final double oneMinusThetaH) {
interpolateDerivatives();
interpolateState(theta);
}
/** Interpolate the derivatives.
* <p>The Adams method is based on a polynomial interpolation of the
* derivatives based on the preceding steps. So the interpolation of
* the derivatives here is strictly equivalent: it is a simple polynomial
* interpolation.</p>
*/
private void interpolateDerivatives() {
for (int i = 0; i < interpolatedDerivatives.length; ++i) {
// initialize the Neville's interpolation algorithm
for (int k = 0; k < previousF.length; ++k) {
neville[k] = previousF[k][i];
}
// combine the contributions of each points
for (int l = 1; l < neville.length; ++l) {
for (int m = neville.length - 1; m >= l; --m) {
final double xm = previousT[m];
final double xmMl = previousT[m - l];
neville[m] = ((interpolatedTime - xm) * neville[m-1] +
(xmMl - interpolatedTime) * neville[m]) / (xmMl - xm);
}
}
// the interpolation polynomial value is in the array last element
interpolatedDerivatives[i] = neville[neville.length - 1];
}
}
/** Interpolate the state.
* <p>The Adams method is based on a polynomial interpolation of the
* derivatives based on the preceding steps. The polynomial model is
* integrated analytically throughout the last step. Using the notations
* found in the second edition of the first volume (Nonstiff Problems)
* of the reference book by Hairer, Norsett and Wanner: <i>Solving Ordinary
* Differential Equations</i> (Springer-Verlag, ISBN 3-540-56670-8), this
* process leads to the following expression:</p>
* <pre>
* y<sub>n+1</sub> = y<sub>n</sub> +
* h &times; &sum;<sub>j=0</sub><sup>j=k</sup> &gamma;<sub>j</sub><sup>*</sup>&nabla;<sup>j</sup>f<sub>n+1</sub>
* </pre>
* <p>In the previous expression, the &gamma;<sub>j</sub><sup>*</sup> terms are the
* ones that result from the analytical integration, and can be computed form
* the binomial coefficients C<sub>j</sub><sup>-s</sup>:</p>
* <p>
* &gamma;<sub>j</sub><sup>*</sup> = (-1)<sup>j</sup>&int;<sub>0</sub><sup>1</sup>C<sub>j</sub><sup>1-s</sup>ds
* </p>
* <p>In order to interpolate the state in a manner that is consistent with the
* integration scheme, we simply subtract from the current state (at the end of the step)
* the integral computed from interpolation time to step end time.</p>
* <p>
* &eta;<sub>j</sub><sup>*</sup>(&theta;)=
* (-1)<sup>j</sup>&int;<sub>&theta;</sub><sup>1</sup>C<sub>j</sub><sup>1-s</sup>ds
* </p>
* The method described in the Hairer, Norsett and Wanner book to compute &gamma;<sub>j</sub><sup>*</sup>
* is easily extended to compute &gamma;<sub>j</sub><sup>*</sup>(&theta;)=
* (-1)<sup>j</sup>&int;<sub>0</sub><sup>&theta;</sup>C<sub>j</sub><sup>1-s</sup>ds. From this,
* we can compute &eta;<sub>j</sub><sup>*</sup>(&theta;) =
* &gamma;<sub>j</sub><sup>*</sup>-&gamma;<sub>j</sub><sup>*</sup>(&theta;).
* The first few values are:</p>
* <table>
* <tr><td>j</td><td>&gamma;<sub>j</sub><sup>*</sup></td><td>&gamma;<sub>j</sub><sup>*</sup>(&theta;)</td><td>&eta;<sub>j</sub><sup>*</sup>(&theta;)</td></tr>
* <tr><td>0</td><td>1</td><td>&theta;</td><td>1-&theta;</td></tr>
* <tr><td>1</td><td>-1/2</td><td>(&theta;<sup>2</sup>-2&theta;)/2</td><td>(-1+2&theta;-&theta;<sup>2</sup>)/2</td></tr>
* <tr><td>2</td><td>-1/12</td><td>(2&theta;<sup>3</sup>-3&theta;<sup>2</sup>)/12</td><td>(-1+3&theta;<sup>2</sup>-2&theta;<sup>3</sup>)/12</td></tr>
* </table>
* <p>
* The &eta;<sub>j</sub>(&theta;) functions appear to be polynomial ones. As expected,
* we see that &eta;<sub>j</sub>(1)= 0. The recurrence relation derived for
* &gamma;<sub>j</sub>(&theta;) is:
* </p>
* <p>
* &sum<sub>j=0</sub><sup>j=m</sup>&gamma;<sub>j</sub><sup>*</sup>(&theta;)/(m+1-j) =
* 1/(m+1)! &prod;<sub>k=0</sub><sup>k=m</sup>(&theta;+k-1)
* </p>
* @param theta location of the interpolation point within the last step
*/
private void interpolateState(final double theta) {
// compute the integrals to remove from the final state
computeRollback(previousT.length - 1, theta);
// remove these integrals from the final state
for (int j = 0; j < interpolatedState.length; ++j) {
double sum = 0;
for (int l = 0; l < previousT.length; ++l) {
sum += rollback[l] * previousF[l][j];
}
interpolatedState[j] = currentState[j] - h * sum;
}
}
/** Compute the rollback coefficients.
* @param order order of the integration method
* @param theta current value for theta
*/
private void computeRollback(final int order, final double theta) {
// compute the gamma star(theta) values from the recurrence relation
double product = theta - 1;
rollback[0] = theta;
for (int i = 1; i <= order; ++i) {
product *= (i - 1 + theta) / (i + 1);
double gStar = product;
for (int j = 1; j <= i; ++j) {
gStar -= rollback[i - j] / (j + 1);
}
rollback[i] = gStar;
}
// subtract it from gamma star to get eta star(theta)
for (int i = 0; i <= order; ++i) {
rollback[i] -= gammaStar[i];
}
// combine the eta star integrals with the backward differences array
// to get the rollback coefficients
for (int i = 0; i <= order; ++i) {
double f = 0;
for (int j = i; j <= order; ++j) {
f -= rollback[j] * bdArray[j][i];
}
rollback[i] = f;
}
}
/** {@inheritDoc} */
public void writeExternal(final ObjectOutput out)
throws IOException {
super.writeExternal(out);
out.writeDouble(nonTruncatedEnd);
}
/** {@inheritDoc} */
public void readExternal(final ObjectInput in)
throws IOException {
nonTruncatedEnd = in.readDouble();
}
}

325
src/java/org/apache/commons/math/ode/nonstiff/MultistepIntegrator.java Normal file
View File

@ -0,0 +1,325 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.ode.nonstiff;
import java.util.Arrays;
import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.ODEIntegrator;
import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.events.EventException;
import org.apache.commons.math.ode.events.EventHandler;
import org.apache.commons.math.ode.events.EventState;
import org.apache.commons.math.ode.sampling.FixedStepHandler;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.ode.sampling.StepInterpolator;
import org.apache.commons.math.ode.sampling.StepNormalizer;
/**
* This class is the base class for multistep integrators for Ordinary
* Differential Equations.
*
* @see AdamsBashforthIntegrator
* @see AdamsMoultonIntegrator
* @see BDFIntegrator
* @version $Revision$ $Date$
* @since 2.0
*/
public abstract class MultistepIntegrator extends AbstractIntegrator {
/** Starter integrator. */
private FirstOrderIntegrator starter;
/** Previous steps times. */
protected double[] previousT;
/** Previous steps derivatives. */
protected double[][] previousF;
/** Time of last detected reset. */
private double resetTime;
/** Prototype of the step interpolator. */
protected MultistepStepInterpolator prototype;
/**
* Build a multistep integrator with the given number of steps.
* <p>The default starter integrator is set to the {@link
* DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
* some defaults settings.</p>
* @param name name of the method
* @param k number of steps of the multistep method
* (including the one being computed)
* @param prototype prototype of the step interpolator to use
*/
protected MultistepIntegrator(final String name, final int k,
final MultistepStepInterpolator prototype) {
super(name);
starter = new DormandPrince853Integrator(1.0e-6, 1.0e6, 1.0e-5, 1.0e-6);
previousT = new double[k];
previousF = new double[k][];
this.prototype = prototype;
}
/**
* Get the starter integrator.
* @return starter integrator
*/
public ODEIntegrator getStarterIntegrator() {
return starter;
}
/**
* Set the starter integrator.
* <p>The various step and event handlers for this starter integrator
* will be managed automatically by the multi-step integrator. Any
* user configuration for these elements will be cleared before use.</p>
* @param starter starter integrator
*/
public void setStarterIntegrator(FirstOrderIntegrator starter) {
this.starter = starter;
}
/** Start the integration.
* <p>This method computes the first few steps of the multistep method,
* using the underlying starter integrator, ensuring the returned steps
* all belong to the same smooth range.</p>
* <p>In order to ensure smoothness, the start phase is automatically
* restarted when a state or derivative reset is triggered by the
* registered events handlers before this start phase is completed. As
* an example, consider integrating a differential equation from t=0
* to t=100 with a 4 steps method and step size equal to 0.2. If an event
* resets the state at t=0.5, the start phase will not end at t=0.7 with
* steps at [0.0, 0.2, 0.4, 0.6] but instead will end at t=1.1 with steps
* at [0.5, 0.7, 0.9, 1.1].</p>
* <p>A side effect of the need for smoothness is that an ODE triggering
* short period regular resets will remain in the start phase throughout
* the integration range if the step size or the number of steps to store
* are too large.</p>
* <p>If the start phase ends prematurely (because of some triggered event
* for example), then the time of latest previous steps will be set to
* <code>Double.NaN</code>.</p>
* @param n number of steps to store
* @param h signed step size to use for the first steps
* @param manager discrete events manager to use
* @param equations differential equations to integrate
* @param t0 initial time
* @param y state vector: contains the initial value of the state vector at t0,
* will be used to put the state vector at each successful step and hence
* contains the final value at the end of the start phase
* @return time of the end of the start phase
* @throws IntegratorException if the integrator cannot perform integration
* @throws DerivativeException this exception is propagated to the caller if
* the underlying user function triggers one
*/
protected double start(final int n, final double h,
final CombinedEventsManager manager,
final FirstOrderDifferentialEquations equations,
final double t0, final double[] y)
throws DerivativeException, IntegratorException {
// clear the first steps
Arrays.fill(previousT, Double.NaN);
Arrays.fill(previousF, null);
// configure the event handlers
starter.clearEventHandlers();
for (EventState state : manager.getEventsStates()) {
starter.addEventHandler(new ResetCheckingWrapper(state.getEventHandler()),
state.getMaxCheckInterval(),
state.getConvergence(), state.getMaxIterationCount());
}
// configure the step handlers
starter.clearStepHandlers();
for (final StepHandler handler : stepHandlers) {
// add the user defined step handlers, filtering out the isLast indicator
starter.addStepHandler(new FilteringWrapper(handler));
}
// add one specific step handler to store the first steps
final StoringStepHandler store = new StoringStepHandler(n);
starter.addStepHandler(new StepNormalizer(h, store));
// integrate over the first few steps, ensuring no intermediate reset occurs
double t = t0;
double stopTime = Double.NaN;
do {
resetTime = Double.NaN;
store.restart();
// we overshoot by 1/10000 step the end to make sure we get don't miss the last point
stopTime = starter.integrate(equations, t, y, t + (n - 0.9999) * h, y);
if (!Double.isNaN(resetTime)) {
// there was an intermediate reset, we restart
t = resetTime;
}
} while (!Double.isNaN(resetTime));
// clear configuration
starter.clearEventHandlers();
starter.clearStepHandlers();
if (store.getFinalState() != null) {
System.arraycopy(store.getFinalState(), 0, y, 0, y.length);
}
return stopTime;
}
/** Rotate the previous steps arrays.
*/
protected void rotatePreviousSteps() {
final double[] rolled = previousF[previousT.length - 1];
for (int k = previousF.length - 1; k > 0; --k) {
previousT[k] = previousT[k - 1];
previousF[k] = previousF[k - 1];
}
previousF[0] = rolled;
}
/** Event handler wrapper to check if state or derivatives have been reset. */
private class ResetCheckingWrapper implements EventHandler {
/** Serializable version identifier. */
private static final long serialVersionUID = 4922660285376467937L;
/** Wrapped event handler. */
private final EventHandler handler;
/** Build a new instance.
* @param handler event handler to wrap
*/
public ResetCheckingWrapper(final EventHandler handler) {
this.handler = handler;
}
/** {@inheritDoc} */
public int eventOccurred(double t, double[] y) throws EventException {
final int action = handler.eventOccurred(t, y);
if ((action == RESET_DERIVATIVES) || (action == RESET_STATE)) {
// a singularity has been encountered
// we need to restart the start phase
resetTime = t;
return STOP;
}
return action;
}
/** {@inheritDoc} */
public double g(double t, double[] y) throws EventException {
return handler.g(t, y);
}
/** {@inheritDoc} */
public void resetState(double t, double[] y) throws EventException {
handler.resetState(t, y);
}
}
/** Step handler wrapper filtering out the isLast indicator. */
private class FilteringWrapper implements StepHandler {
/** Serializable version identifier. */
private static final long serialVersionUID = 4607975253344802232L;
/** Wrapped step handler. */
private final StepHandler handler;
/** Build a new instance.
* @param handler step handler to wrap
*/
public FilteringWrapper(final StepHandler handler) {
this.handler = handler;
}
/** {@inheritDoc} */
public void handleStep(StepInterpolator interpolator, boolean isLast)
throws DerivativeException {
// we force the isLast indicator to false EXCEPT if some event handler triggered a stop
handler.handleStep(interpolator, eventsHandlersManager.stop());
}
/** {@inheritDoc} */
public boolean requiresDenseOutput() {
return handler.requiresDenseOutput();
}
/** {@inheritDoc} */
public void reset() {
handler.reset();
}
}
/** Specialized step handler storing the first few steps. */
private class StoringStepHandler implements FixedStepHandler {
/** Serializable version identifier. */
private static final long serialVersionUID = 4592974435520688797L;
/** Number of steps to store. */
private final int n;
/** Counter for already stored steps. */
private int count;
/** Final state. */
private double[] finalState;
/** Build a new instance.
* @param number of steps to store
*/
public StoringStepHandler(final int n) {
this.n = n;
restart();
}
/** Restart storage.
*/
public void restart() {
count = 0;
finalState = null;
}
/** Get the final state.
* @return final state
*/
public double[] getFinalState() {
return finalState;
}
/** {@inheritDoc} */
public void handleStep(final double t, final double[] y, final double[] yDot,
final boolean isLast) {
if (count++ < n) {
previousT[n - count] = t;
previousF[n - count] = yDot.clone();
if (count == n) {
finalState = y.clone();
}
}
}
}
}

165
src/java/org/apache/commons/math/ode/nonstiff/MultistepStepInterpolator.java Normal file
View File

@ -0,0 +1,165 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.ode.nonstiff;
import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
/** This class represents an interpolator over the last step during an
* ODE integration for multistep integrators.
*
* @see MultistepIntegrator
*
* @version $Revision$ $Date$
* @since 2.0
*/
abstract class MultistepStepInterpolator
extends AbstractStepInterpolator {
/** Previous steps times. */
protected double[] previousT;
/** Previous steps derivatives. */
protected double[][] previousF;
/** Simple constructor.
* This constructor builds an instance that is not usable yet, the
* {@link #reinitialize} method should be called before using the
* instance in order to initialize the internal arrays. This
* constructor is used only in order to delay the initialization in
* some cases. The {@link MultistepIntegrator} classe uses the
* prototyping design pattern to create the step interpolators by
* cloning an uninitialized model and latter initializing the copy.
*/
protected MultistepStepInterpolator() {
previousT = null;
previousF = null;
}
/** Copy constructor.
* <p>The copied interpolator should have been finalized before the
* copy, otherwise the copy will not be able to perform correctly any
* interpolation and will throw a {@link NullPointerException}
* later. Since we don't want this constructor to throw the
* exceptions finalization may involve and since we don't want this
* method to modify the state of the copied interpolator,
* finalization is <strong>not</strong> done automatically, it
* remains under user control.</p>
* <p>The copy is a deep copy: its arrays are separated from the
* original arrays of the instance.</p>
* @param interpolator interpolator to copy from.
*/
public MultistepStepInterpolator(final MultistepStepInterpolator interpolator) {
super(interpolator);
if (interpolator.currentState != null) {
previousT = interpolator.previousT.clone();
previousF = new double[interpolator.previousF.length][];
for (int k = 0; k < interpolator.previousF.length; ++k) {
previousF[k] = interpolator.previousF[k].clone();
}
initializeCoefficients();
} else {
previousT = null;
previousF = null;
}
}
/** Reinitialize the instance
* @param y reference to the integrator array holding the state at
* the end of the step
* @param previousT reference to the integrator array holding the times
* of the previous steps
* @param previousF reference to the integrator array holding the
* previous slopes
* @param forward integration direction indicator
*/
public void reinitialize(final double[] y,
final double[] previousT, final double[][] previousF,
final boolean forward) {
reinitialize(y, forward);
this.previousT = previousT;
this.previousF = previousF;
initializeCoefficients();
}
/** Initialize the coefficients arrays.
*/
protected abstract void initializeCoefficients();
/** {@inheritDoc} */
public void writeExternal(final ObjectOutput out)
throws IOException {
// save the state of the base class
writeBaseExternal(out);
// save the local attributes
out.writeInt(previousT.length);
for (int k = 0; k < previousF.length; ++k) {
out.writeDouble(previousT[k]);
for (int i = 0; i < currentState.length; ++i) {
out.writeDouble(previousF[k][i]);
}
}
}
/** {@inheritDoc} */
public void readExternal(final ObjectInput in)
throws IOException {
// read the base class
final double t = readBaseExternal(in);
// read the local attributes
final int kMax = in.readInt();
previousT = new double[kMax];
previousF = new double[kMax][];
for (int k = 0; k < kMax; ++k) {
previousT[k] = in.readDouble();
previousF[k] = new double[currentState.length];
for (int i = 0; i < currentState.length; ++i) {
previousF[k][i] = in.readDouble();
}
}
// initialize the coefficients
initializeCoefficients();
try {
// we can now set the interpolated time and state
setInterpolatedTime(t);
} catch (DerivativeException e) {
throw new IOException(e.getMessage());
}
}
}

7
src/site/xdoc/changes.xml
View File

@ -39,6 +39,13 @@ The <action> type attribute can be add,update,fix,remove.
</properties>
<body>
<release version="2.0" date="TBD" description="TBD">
<action dev="luc" type="add">
New ODE integrators have been added: the explicit Adams-Bashforth and implicit
Adams-Moulton multistep methods. These methods support customizable starter
integrators and support discrete events even during the start phase.
All these methods provide the same rich features has the existing ones:
continuous output, step handlers, discrete events, G-stop ...
</action>
<action dev="luc" type="fix" issue="MATH-214" >
Replaced size adjustment of all steps of fixed steps Runge-Kutta integrators by
a truncation of the last step only.

206
src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java Normal file
View File

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

216
src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java Normal file
View File

@ -0,0 +1,216 @@
/*
* 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.ode.DerivativeException;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.events.EventHandler;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
public class AdamsMoultonIntegratorTest
extends TestCase {
public AdamsMoultonIntegratorTest(String name) {
super(name);
}
public void testPredictorCoefficients() {
for (int order = 1; order < 10; ++order) {
double[] moulton = new AdamsMoultonIntegrator(order, 0.01).getPredictorCoeffs();
double[] bashforth = new AdamsBashforthIntegrator(order, 0.01).getCoeffs();
assertEquals(bashforth.length, moulton.length);
for (int i = 0; i < moulton.length; ++i) {
assertEquals(bashforth[i], moulton[i], 1.0e-16);
}
}
}
public void testCorrectorCoefficients() {
double[] coeffs1 = new AdamsMoultonIntegrator(1, 0.01).getCorrectorCoeffs();
assertEquals(2, coeffs1.length);
assertEquals(1.0 / 2.0, coeffs1[0], 1.0e-16);
assertEquals(1.0 / 2.0, coeffs1[1], 1.0e-16);
double[] coeffs2 = new AdamsMoultonIntegrator(2, 0.01).getCorrectorCoeffs();
assertEquals(3, coeffs2.length);
assertEquals( 5.0 / 12.0, coeffs2[0], 1.0e-16);
assertEquals( 8.0 / 12.0, coeffs2[1], 1.0e-16);
assertEquals(-1.0 / 12.0, coeffs2[2], 1.0e-16);
double[] coeffs3 = new AdamsMoultonIntegrator(3, 0.01).getCorrectorCoeffs();
assertEquals(4, coeffs3.length);
assertEquals( 9.0 / 24.0, coeffs3[0], 1.0e-16);
assertEquals(19.0 / 24.0, coeffs3[1], 1.0e-16);
assertEquals(-5.0 / 24.0, coeffs3[2], 1.0e-16);
assertEquals( 1.0 / 24.0, coeffs3[3], 1.0e-16);
double[] coeffs4 = new AdamsMoultonIntegrator(4, 0.01).getCorrectorCoeffs();
assertEquals(5, coeffs4.length);
assertEquals( 251.0 / 720.0, coeffs4[0], 1.0e-16);
assertEquals( 646.0 / 720.0, coeffs4[1], 1.0e-16);
assertEquals(-264.0 / 720.0, coeffs4[2], 1.0e-16);
assertEquals( 106.0 / 720.0, coeffs4[3], 1.0e-16);
assertEquals( -19.0 / 720.0, coeffs4[4], 1.0e-16);
double[] coeffs5 = new AdamsMoultonIntegrator(5, 0.01).getCorrectorCoeffs();
assertEquals(6, coeffs5.length);
assertEquals( 475.0 / 1440.0, coeffs5[0], 1.0e-16);
assertEquals(1427.0 / 1440.0, coeffs5[1], 1.0e-16);
assertEquals(-798.0 / 1440.0, coeffs5[2], 1.0e-16);
assertEquals( 482.0 / 1440.0, coeffs5[3], 1.0e-16);
assertEquals(-173.0 / 1440.0, coeffs5[4], 1.0e-16);
assertEquals( 27.0 / 1440.0, coeffs5[5], 1.0e-16);
double[] coeffs6 = new AdamsMoultonIntegrator(6, 0.01).getCorrectorCoeffs();
assertEquals(7, coeffs6.length);
assertEquals( 19087.0 / 60480.0, coeffs6[0], 1.0e-16);
assertEquals( 65112.0 / 60480.0, coeffs6[1], 1.0e-16);
assertEquals(-46461.0 / 60480.0, coeffs6[2], 1.0e-16);
assertEquals( 37504.0 / 60480.0, coeffs6[3], 1.0e-16);
assertEquals(-20211.0 / 60480.0, coeffs6[4], 1.0e-16);
assertEquals( 6312.0 / 60480.0, coeffs6[5], 1.0e-16);
assertEquals( -863.0 / 60480.0, coeffs6[6], 1.0e-16);
double[] coeffs7 = new AdamsMoultonIntegrator(7, 0.01).getCorrectorCoeffs();
assertEquals(8, coeffs7.length);
assertEquals( 36799.0 / 120960.0, coeffs7[0], 1.0e-16);
assertEquals( 139849.0 / 120960.0, coeffs7[1], 1.0e-16);
assertEquals(-121797.0 / 120960.0, coeffs7[2], 1.0e-16);
assertEquals( 123133.0 / 120960.0, coeffs7[3], 1.0e-16);
assertEquals( -88547.0 / 120960.0, coeffs7[4], 1.0e-16);
assertEquals( 41499.0 / 120960.0, coeffs7[5], 1.0e-16);
assertEquals( -11351.0 / 120960.0, coeffs7[6], 1.0e-16);
assertEquals( 1375.0 / 120960.0, coeffs7[7], 1.0e-16);
double[] coeffs8 = new AdamsMoultonIntegrator(8, 0.01).getCorrectorCoeffs();
assertEquals(9, coeffs8.length);
assertEquals( 1070017.0 / 3628800.0, coeffs8[0], 1.0e-16);
assertEquals( 4467094.0 / 3628800.0, coeffs8[1], 1.0e-16);
assertEquals(-4604594.0 / 3628800.0, coeffs8[2], 1.0e-16);
assertEquals( 5595358.0 / 3628800.0, coeffs8[3], 1.0e-16);
assertEquals(-5033120.0 / 3628800.0, coeffs8[4], 1.0e-16);
assertEquals( 3146338.0 / 3628800.0, coeffs8[5], 1.0e-16);
assertEquals(-1291214.0 / 3628800.0, coeffs8[6], 1.0e-16);
assertEquals( 312874.0 / 3628800.0, coeffs8[7], 1.0e-16);
assertEquals( -33953.0 / 3628800.0, coeffs8[8], 1.0e-16);
}
public void testDimensionCheck() {
try {
TestProblem1 pb = new TestProblem1();
new AdamsMoultonIntegrator(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);
if (pb instanceof TestProblem3) {
step /= 8;
}
FirstOrderIntegrator integ = new AdamsMoultonIntegrator(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 AdamsMoultonIntegrator(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() < 7.0e-12);
assertTrue(handler.getMaximalValueError() < 4.0e-11);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
assertEquals("Adams-Moulton", integ.getName());
}
public void testBigStep()
throws DerivativeException, IntegratorException {
TestProblem1 pb = new TestProblem1();
double step = (pb.getFinalTime() - pb.getInitialTime()) * 0.2;
FirstOrderIntegrator integ = new AdamsMoultonIntegrator(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.01);
assertTrue(handler.getMaximalValueError() > 0.03);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-14);
}
public static Test suite() {
return new TestSuite(AdamsMoultonIntegratorTest.class);
}
}