You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2009/06/09 21:33:20 UTC
svn commit: r783103 - in /commons/proper/math/trunk/src:
java/org/apache/commons/math/ode/ java/org/apache/commons/math/ode/nonstiff/
java/org/apache/commons/math/ode/sampling/ site/xdoc/ site/xdoc/userguide/
test/org/apache/commons/math/ode/nonstiff/
Author: luc
Date: Tue Jun 9 19:33:19 2009
New Revision: 783103
URL: http://svn.apache.org/viewvc?rev=783103&view=rev
Log:
added support for max number of evaluations to ODE integrators
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AbstractIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/ODEIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
commons/proper/math/trunk/src/site/xdoc/changes.xml
commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince54IntegratorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AbstractIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AbstractIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AbstractIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/AbstractIntegrator.java Tue Jun 9 19:33:19 2009
@@ -21,6 +21,7 @@
import java.util.Collection;
import java.util.Collections;
+import org.apache.commons.math.MaxEvaluationsExceededException;
import org.apache.commons.math.ode.events.CombinedEventsManager;
import org.apache.commons.math.ode.events.EventHandler;
import org.apache.commons.math.ode.events.EventState;
@@ -39,6 +40,15 @@
/** Name of the method. */
private final String name;
+ /** Maximal number of evaluations allowed. */
+ private int maxEvaluations;
+
+ /** Number of evaluations already performed. */
+ private int evaluations;
+
+ /** Differential equations to integrate. */
+ private transient FirstOrderDifferentialEquations equations;
+
/** Step handler. */
protected Collection<StepHandler> stepHandlers;
@@ -60,6 +70,8 @@
stepStart = Double.NaN;
stepSize = Double.NaN;
eventsHandlersManager = new CombinedEventsManager();
+ setMaxEvaluations(-1);
+ resetEvaluations();
}
/** {@inheritDoc} */
@@ -123,6 +135,49 @@
return stepSize;
}
+ /** {@inheritDoc} */
+ public void setMaxEvaluations(int maxEvaluations) {
+ this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations;
+ }
+
+ /** {@inheritDoc} */
+ public int getMaxEvaluations() {
+ return maxEvaluations;
+ }
+
+ /** {@inheritDoc} */
+ public int getEvaluations() {
+ return evaluations;
+ }
+
+ /** Reset the number of evaluations to zero.
+ */
+ protected void resetEvaluations() {
+ evaluations = 0;
+ }
+
+ /** Set the differential equations.
+ * @see #computeDerivatives(double, double[], double[])
+ */
+ protected void setEquations(final FirstOrderDifferentialEquations equations) {
+ this.equations = equations;
+ }
+
+ /** Compute the derivatives and check the number of evaluations.
+ * @param t current value of the independent <I>time</I> variable
+ * @param y array containing the current value of the state vector
+ * @param yDot placeholder array where to put the time derivative of the state vector
+ * @throws DerivativeException this exception is propagated to the caller if the
+ * underlying user function triggers one
+ */
+ public void computeDerivatives(final double t, final double[] y, final double[] yDot)
+ throws DerivativeException {
+ if (++evaluations > maxEvaluations) {
+ throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations));
+ }
+ equations.computeDerivatives(t, y, yDot);
+ }
+
/** Perform some sanity checks on the integration parameters.
* @param equations differential equations set
* @param t0 start time
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java Tue Jun 9 19:33:19 2009
@@ -126,7 +126,6 @@
*/
protected double start(final int n, final double h,
final CombinedEventsManager manager,
- final FirstOrderDifferentialEquations equations,
final double t0, final double[] y)
throws DerivativeException, IntegratorException {
@@ -156,6 +155,8 @@
// integrate over the first few steps, ensuring no intermediate reset occurs
double t = t0;
double stopTime = Double.NaN;
+ FirstOrderDifferentialEquations equations =
+ new CountingDifferentialEquations(y.length);
do {
resetTime = Double.NaN;
final double dt = (n - 0.9999) * h;
@@ -341,4 +342,32 @@
}
+ /** Wrapper for differential equations, ensuring start evaluations are counted. */
+ private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {
+
+ /** Serializable uid. */
+ private static final long serialVersionUID = -6329212616396607764L;
+
+ /** Dimension of the problem. */
+ private final int dimension;
+
+ /** Simple constructor.
+ * @param dimension dimension of the problem
+ */
+ public CountingDifferentialEquations(final int dimension) {
+ this.dimension = dimension;
+ }
+
+ /** {@inheritDoc} */
+ public void computeDerivatives(double t, double[] y, double[] dot)
+ throws DerivativeException {
+ MultistepIntegrator.this.computeDerivatives(t, y, dot);
+ }
+
+ /** {@inheritDoc} */
+ public int getDimension() {
+ return dimension;
+ }
+ }
+
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/ODEIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/ODEIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/ODEIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/ODEIntegrator.java Tue Jun 9 19:33:19 2009
@@ -114,4 +114,28 @@
*/
double getCurrentSignedStepsize();
+ /** Set the maximal number of differential equations function evaluations.
+ * <p>The purpose of this method is to avoid infinite loops which can occur
+ * for example when stringent error constraints are set or when lots of
+ * discrete events are triggered, thus leading to many rejected steps.</p>
+ * @param maxEvaluations maximal number of function evaluations (negative
+ * values are silently converted to maximal integer value, thus representing
+ * almost unlimited evaluations)
+ */
+ void setMaxEvaluations(int maxEvaluations);
+
+ /** Get the maximal number of functions evaluations.
+ * @return maximal number of functions evaluations
+ */
+ int getMaxEvaluations();
+
+ /** Get the number of evaluations of the differential equations function.
+ * <p>
+ * The number of evaluations corresponds to the last call to the
+ * <code>integrate</code> method. It is 0 if the method has not been called yet.
+ * </p>
+ * @return number of evaluations of the differential equations function
+ */
+ int getEvaluations();
+
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Tue Jun 9 19:33:19 2009
@@ -209,6 +209,8 @@
final int n = y0.length;
sanityChecks(equations, t0, y0, t, y);
+ setEquations(equations);
+ resetEvaluations();
final boolean forward = (t > t0);
// initialize working arrays
@@ -229,8 +231,7 @@
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);
+ double stopTime = start(previousF.length, stepSize, manager, stepStart, y);
if (Double.isNaN(previousT[0])) {
return stopTime;
}
@@ -264,7 +265,7 @@
// update the Nordsieck vector
final double[] f0 = previousF[0];
previousT[0] = nextStep;
- equations.computeDerivatives(nextStep, y, f0);
+ computeDerivatives(nextStep, y, f0);
nordsieck = coefficients.msUpdate.multiply(nordsieck);
final double[] end = new double[y0.length];
for (int j = 0; j < y0.length; ++j) {
@@ -284,8 +285,7 @@
// 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);
+ stopTime = start(previousF.length, stepSize, manager, stepStart, y);
if (Double.isNaN(previousT[0])) {
return stopTime;
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Tue Jun 9 19:33:19 2009
@@ -223,6 +223,8 @@
final int n = y0.length;
sanityChecks(equations, t0, y0, t, y);
+ setEquations(equations);
+ resetEvaluations();
final boolean forward = (t > t0);
// initialize working arrays
@@ -246,8 +248,7 @@
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);
+ double stopTime = start(previousF.length, stepSize, manager, stepStart, y);
if (Double.isNaN(previousT[0])) {
return stopTime;
}
@@ -279,7 +280,7 @@
// evaluate a first estimate of the derivative (first E in the PECE sequence)
final double[] f0 = previousF[0];
previousT[0] = stepEnd;
- equations.computeDerivatives(stepEnd, yTmp, f0);
+ computeDerivatives(stepEnd, yTmp, f0);
// update Nordsieck vector
final RealMatrix nordsieckTmp = coefficients.msUpdate.multiply(nordsieck);
@@ -293,7 +294,7 @@
nordsieckTmp.walkInOptimizedOrder(new Corrector(y, predictedScaled, yTmp));
// evaluate a final estimate of the derivative (second E in the PECE sequence)
- equations.computeDerivatives(stepEnd, yTmp, f0);
+ computeDerivatives(stepEnd, yTmp, f0);
// update Nordsieck vector
final double[] correctedScaled = new double[y0.length];
@@ -338,8 +339,7 @@
// 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);
+ stopTime = start(previousF.length, stepSize, manager, stepStart, y);
if (Double.isNaN(previousT[0])) {
return stopTime;
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java Tue Jun 9 19:33:19 2009
@@ -215,7 +215,7 @@
for (int j = 0; j < y0.length; ++j) {
y1[j] = y0[j] + h * yDot0[j];
}
- equations.computeDerivatives(t0 + h, y1, yDot1);
+ computeDerivatives(t0 + h, y1, yDot1);
// estimate the second derivative of the solution
double yDDotOnScale = 0;
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java Tue Jun 9 19:33:19 2009
@@ -17,8 +17,8 @@
package org.apache.commons.math.ode.nonstiff;
+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.sampling.StepInterpolator;
/**
@@ -90,9 +90,9 @@
/** {@inheritDoc} */
@Override
- public void reinitialize(final FirstOrderDifferentialEquations equations,
+ public void reinitialize(final AbstractIntegrator integrator,
final double[] y, final double[][] yDotK, final boolean forward) {
- super.reinitialize(equations, y, yDotK, forward);
+ super.reinitialize(integrator, y, yDotK, forward);
v1 = null;
v2 = null;
v3 = null;
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java Tue Jun 9 19:33:19 2009
@@ -22,8 +22,8 @@
import java.io.ObjectOutput;
import org.apache.commons.math.MathRuntimeException;
+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.sampling.StepInterpolator;
/**
@@ -101,10 +101,10 @@
/** {@inheritDoc} */
@Override
- public void reinitialize(final FirstOrderDifferentialEquations equations,
+ public void reinitialize(final AbstractIntegrator integrator,
final double[] y, final double[][] yDotK, final boolean forward) {
- super.reinitialize(equations, y, yDotK, forward);
+ super.reinitialize(integrator, y, yDotK, forward);
final int dimension = currentState.length;
@@ -224,7 +224,7 @@
k14_11 * yDotK[10][j] + k14_12 * yDotK[11][j] + k14_13 * yDotK[12][j];
yTmp[j] = currentState[j] + h * s;
}
- equations.computeDerivatives(previousTime + c14 * h, yTmp, yDotKLast[0]);
+ integrator.computeDerivatives(previousTime + c14 * h, yTmp, yDotKLast[0]);
// k15
for (int j = 0; j < currentState.length; ++j) {
@@ -234,7 +234,7 @@
k15_14 * yDotKLast[0][j];
yTmp[j] = currentState[j] + h * s;
}
- equations.computeDerivatives(previousTime + c15 * h, yTmp, yDotKLast[1]);
+ integrator.computeDerivatives(previousTime + c15 * h, yTmp, yDotKLast[1]);
// k16
for (int j = 0; j < currentState.length; ++j) {
@@ -244,7 +244,7 @@
k16_14 * yDotKLast[0][j] + k16_15 * yDotKLast[1][j];
yTmp[j] = currentState[j] + h * s;
}
- equations.computeDerivatives(previousTime + c16 * h, yTmp, yDotKLast[2]);
+ integrator.computeDerivatives(previousTime + c16 * h, yTmp, yDotKLast[2]);
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java Tue Jun 9 19:33:19 2009
@@ -171,6 +171,8 @@
throws DerivativeException, IntegratorException {
sanityChecks(equations, t0, y0, t, y);
+ setEquations(equations);
+ resetEvaluations();
final boolean forward = (t > t0);
// create some internal working arrays
@@ -188,7 +190,7 @@
AbstractStepInterpolator interpolator;
if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
- rki.reinitialize(equations, yTmp, yDotK, forward);
+ rki.reinitialize(this, yTmp, yDotK, forward);
interpolator = rki;
} else {
interpolator = new DummyStepInterpolator(yTmp, forward);
@@ -215,7 +217,7 @@
if (firstTime || !fsal) {
// first stage
- equations.computeDerivatives(stepStart, y, yDotK[0]);
+ computeDerivatives(stepStart, y, yDotK[0]);
}
if (firstTime) {
@@ -246,7 +248,7 @@
yTmp[j] = y[j] + stepSize * sum;
}
- equations.computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
+ computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
}
@@ -304,7 +306,7 @@
if (manager.reset(stepStart, y) && ! lastStep) {
// some event handler has triggered changes that
// invalidate the derivatives, we need to recompute them
- equations.computeDerivatives(stepStart, y, yDotK[0]);
+ computeDerivatives(stepStart, y, yDotK[0]);
}
if (! lastStep) {
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java Tue Jun 9 19:33:19 2009
@@ -395,7 +395,6 @@
/** Perform integration over one step using substeps of a modified
* midpoint method.
- * @param equations differential equations to integrate
* @param t0 initial time
* @param y0 initial value of the state vector at t0
* @param step global step
@@ -411,8 +410,7 @@
* @throws DerivativeException this exception is propagated to the caller if the
* underlying user function triggers one
*/
- private boolean tryStep(final FirstOrderDifferentialEquations equations,
- final double t0, final double[] y0, final double step, final int k,
+ private boolean tryStep(final double t0, final double[] y0, final double step, final int k,
final double[] scale, final double[][] f,
final double[] yMiddle, final double[] yEnd,
final double[] yTmp)
@@ -428,7 +426,7 @@
yTmp[i] = y0[i];
yEnd[i] = y0[i] + subStep * f[0][i];
}
- equations.computeDerivatives(t, yEnd, f[1]);
+ computeDerivatives(t, yEnd, f[1]);
// other substeps
for (int j = 1; j < n; ++j) {
@@ -445,7 +443,7 @@
yTmp[i] = middle;
}
- equations.computeDerivatives(t, yEnd, f[j+1]);
+ computeDerivatives(t, yEnd, f[j+1]);
// stability check
if (performTest && (j <= maxChecks) && (k < maxIter)) {
@@ -508,6 +506,8 @@
throws DerivativeException, IntegratorException {
sanityChecks(equations, t0, y0, t, y);
+ setEquations(equations);
+ resetEvaluations();
final boolean forward = (t > t0);
// create some internal working arrays
@@ -599,7 +599,7 @@
// first evaluation, at the beginning of the step
if (! firstStepAlreadyComputed) {
- equations.computeDerivatives(stepStart, y, yDot0);
+ computeDerivatives(stepStart, y, yDot0);
}
if (firstTime) {
@@ -635,7 +635,7 @@
++k;
// modified midpoint integration with the current substep
- if ( ! tryStep(equations, stepStart, y, stepSize, k, scale, fk[k],
+ if ( ! tryStep(stepStart, y, stepSize, k, scale, fk[k],
(k == 0) ? yMidDots[0] : diagonal[k-1],
(k == 0) ? y1 : y1Diag[k-1],
yTmp)) {
@@ -773,7 +773,7 @@
}
// derivative at end of step
- equations.computeDerivatives(stepStart + stepSize, y1, yDot1);
+ computeDerivatives(stepStart + stepSize, y1, yDot1);
final int mu = 2 * k - mudif + 3;
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java Tue Jun 9 19:33:19 2009
@@ -86,6 +86,8 @@
throws DerivativeException, IntegratorException {
sanityChecks(equations, t0, y0, t, y);
+ setEquations(equations);
+ resetEvaluations();
final boolean forward = (t > t0);
// create some internal working arrays
@@ -103,7 +105,7 @@
AbstractStepInterpolator interpolator;
if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
- rki.reinitialize(equations, yTmp, yDotK, forward);
+ rki.reinitialize(this, yTmp, yDotK, forward);
interpolator = rki;
} else {
interpolator = new DummyStepInterpolator(yTmp, forward);
@@ -127,7 +129,7 @@
for (boolean loop = true; loop;) {
// first stage
- equations.computeDerivatives(stepStart, y, yDotK[0]);
+ computeDerivatives(stepStart, y, yDotK[0]);
// next stages
for (int k = 1; k < stages; ++k) {
@@ -140,7 +142,7 @@
yTmp[j] = y[j] + stepSize * sum;
}
- equations.computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
+ computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
}
@@ -179,7 +181,7 @@
if (manager.reset(stepStart, y) && ! lastStep) {
// some events handler has triggered changes that
// invalidate the derivatives, we need to recompute them
- equations.computeDerivatives(stepStart, y, yDotK[0]);
+ computeDerivatives(stepStart, y, yDotK[0]);
}
// make sure step size is set to default before next step
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java Tue Jun 9 19:33:19 2009
@@ -21,7 +21,7 @@
import java.io.ObjectInput;
import java.io.ObjectOutput;
-import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.AbstractIntegrator;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
/** This class represents an interpolator over the last step during an
@@ -49,8 +49,8 @@
*/
protected RungeKuttaStepInterpolator() {
super();
- yDotK = null;
- equations = null;
+ yDotK = null;
+ integrator = null;
}
/** Copy constructor.
@@ -90,7 +90,7 @@
// we cannot keep any reference to the equations in the copy
// the interpolator should have been finalized before
- equations = null;
+ integrator = null;
}
@@ -108,18 +108,18 @@
* {@link AbstractStepInterpolator#getInterpolatedState
* getInterpolatedState} method (for an interpolator which needs a
* finalization) or if it clones the step interpolator.</p>
- * @param equations set of differential equations being integrated
+ * @param integrator integrator being used
* @param y reference to the integrator array holding the state at
* the end of the step
* @param yDotK reference to the integrator array holding all the
* intermediate slopes
* @param forward integration direction indicator
*/
- public void reinitialize(final FirstOrderDifferentialEquations equations,
+ public void reinitialize(final AbstractIntegrator integrator,
final double[] y, final double[][] yDotK, final boolean forward) {
reinitialize(y, forward);
this.yDotK = yDotK;
- this.equations = equations;
+ this.integrator = integrator;
}
/** {@inheritDoc} */
@@ -163,7 +163,7 @@
}
}
- equations = null;
+ integrator = null;
if (currentState != null) {
// we can now set the interpolated time and state
@@ -177,7 +177,7 @@
/** Slopes at the intermediate points */
protected double[][] yDotK;
- /** Reference to the differential equations being integrated. */
- protected FirstOrderDifferentialEquations equations;
+ /** Reference to the integrator. */
+ protected AbstractIntegrator integrator;
}
Modified: commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java (original)
+++ commons/proper/math/trunk/src/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java Tue Jun 9 19:33:19 2009
@@ -17,9 +17,9 @@
package org.apache.commons.math.ode.sampling;
+import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;
-import java.io.IOException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.ode.DerivativeException;
@@ -337,7 +337,6 @@
* @throws DerivativeException this exception is propagated to the
* caller if the underlying user function triggers one
-
*/
public final void finalizeStep()
throws DerivativeException {
Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Tue Jun 9 19:33:19 2009
@@ -309,6 +309,12 @@
and a clearStepHandlers method has been added.
</action>
<action dev="luc" type="add">
+ All ODE integrators now support setting a maximal number of evaluations of differential
+ equations function. If this number is exceeded, an exception will be thrown during
+ integration. This can be used to prevent infinite loops if for example error control or
+ discrete events create a really large number of extremely small steps.
+ </action>
+ <action dev="luc" type="add">
All step interpolators for ODE integrators now provide interpolation for
both the state and its time derivatives. The interpolated derivatives are
the exact derivatives of the interpolated state, thus preserving consistency.
Modified: commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/ode.xml Tue Jun 9 19:33:19 2009
@@ -54,6 +54,15 @@
equations gracefully and get accurate dense output even close to the discontinuity.
</p>
<p>
+ All integrators support setting a maximal number of evaluations of differential
+ equations function. If this number is exceeded, an exception will be thrown during
+ integration. This can be used to prevent infinite loops if for example error control or
+ discrete events create a really large number of extremely small steps. By default, the
+ maximal number of evaluation is set to <code>Integer.MAX_VALUE</code> (i.e. 2<sup>31</sup>-1
+ or 2147483647). It is recommended to set this maximal number to a value suited to the ODE
+ problem, integration range, and step size or error control settings.
+ </p>
+ <p>
The user should describe his problem in his own classes which should implement the
<a href="../apidocs/org/apache/commons/math/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
interface. Then he should pass it to the integrator he prefers among all the classes that implement
Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java Tue Jun 9 19:33:19 2009
@@ -99,6 +99,26 @@
assertTrue(handler.getMaximalValueError() < 9.0e-9);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
assertEquals("Adams-Bashforth", integ.getName());
+ assertTrue(integ.getEvaluations() > 1000);
+ assertEquals(Integer.MAX_VALUE, integ.getMaxEvaluations());
+
+ }
+
+ @Test(expected = DerivativeException.class)
+ public void exceedMaxEvaluations() throws DerivativeException, IntegratorException {
+
+ TestProblem1 pb = new TestProblem1();
+ double range = pb.getFinalTime() - pb.getInitialTime();
+ double step = range * 0.001;
+
+ AdamsBashforthIntegrator integ = new AdamsBashforthIntegrator(3, step);
+ integ.setStarterIntegrator(new DormandPrince853Integrator(0, range, 1.0e-12, 1.0e-12));
+ TestProblemHandler handler = new TestProblemHandler(pb, integ);
+ integ.addStepHandler(handler);
+ integ.setMaxEvaluations(1000);
+ integ.integrate(pb,
+ pb.getInitialTime(), pb.getInitialState(),
+ pb.getFinalTime(), new double[pb.getDimension()]);
}
@@ -172,8 +192,8 @@
ByteArrayOutputStream bos = new ByteArrayOutputStream();
ObjectOutputStream oos = new ObjectOutputStream(bos);
oos.writeObject(new AdamsBashforthIntegrator(8, step));
- assertTrue(bos.size() > 2800);
- assertTrue(bos.size() < 2900);
+ assertTrue(bos.size() > 2900);
+ assertTrue(bos.size() < 3000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);
Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java Tue Jun 9 19:33:19 2009
@@ -176,8 +176,8 @@
ByteArrayOutputStream bos = new ByteArrayOutputStream();
ObjectOutputStream oos = new ObjectOutputStream(bos);
oos.writeObject(new AdamsMoultonIntegrator(8, step));
- assertTrue(bos.size() > 2800);
- assertTrue(bos.size() < 2900);
+ assertTrue(bos.size() > 2900);
+ assertTrue(bos.size() < 3000);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);
Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince54IntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince54IntegratorTest.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince54IntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince54IntegratorTest.java Tue Jun 9 19:33:19 2009
@@ -189,6 +189,7 @@
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
int calls = pb.getCalls();
+ assertEquals(integ.getEvaluations(), calls);
assertTrue(calls <= previousCalls);
previousCalls = calls;
@@ -245,6 +246,7 @@
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
+ assertEquals(integ.getEvaluations(), pb.getCalls());
assertTrue(pb.getCalls() < 2800);
}
Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java Tue Jun 9 19:33:19 2009
@@ -117,6 +117,7 @@
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
int calls = pb.getCalls();
+ assertEquals(integ.getEvaluations(), calls);
assertTrue(calls <= previousCalls);
previousCalls = calls;
@@ -196,6 +197,7 @@
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
+ assertEquals(integ.getEvaluations(), pb.getCalls());
assertTrue(pb.getCalls() < 3300);
}
@@ -237,12 +239,14 @@
pb1.getInitialTime(), pb1.getInitialState(),
pb1.getFinalTime(), new double[pb1.getDimension()]);
int callsWithoutDenseOutput = pb1.getCalls();
+ assertEquals(integ.getEvaluations(), callsWithoutDenseOutput);
integ.addStepHandler(new InterpolatingStepHandler());
integ.integrate(pb2,
pb2.getInitialTime(), pb2.getInitialState(),
pb2.getFinalTime(), new double[pb2.getDimension()]);
int callsWithDenseOutput = pb2.getCalls();
+ assertEquals(integ.getEvaluations(), callsWithDenseOutput);
assertTrue(callsWithDenseOutput > callsWithoutDenseOutput);
Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java Tue Jun 9 19:33:19 2009
@@ -28,7 +28,6 @@
import org.apache.commons.math.ode.ContinuousOutputModel;
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.sampling.StepHandler;
import org.apache.commons.math.ode.sampling.StepInterpolatorTestUtils;
@@ -42,7 +41,7 @@
double[] y = { 0.0, 1.0, -2.0 };
double[][] yDot = { { 1.0, 2.0, -2.0 } };
EulerStepInterpolator interpolator = new EulerStepInterpolator();
- interpolator.reinitialize(new DummyEquations(), y, yDot, true);
+ interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
interpolator.storeTime(0);
interpolator.shift();
interpolator.storeTime(1);
@@ -64,7 +63,7 @@
double[] y = y0.clone();
double[][] yDot = { new double[y0.length] };
EulerStepInterpolator interpolator = new EulerStepInterpolator();
- interpolator.reinitialize(new DummyEquations(), y, yDot, true);
+ interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
interpolator.storeTime(t0);
double dt = 1.0;
@@ -98,7 +97,7 @@
double[] y = { 1.0, 3.0, -4.0 };
double[][] yDot = { { 1.0, 2.0, -2.0 } };
EulerStepInterpolator interpolator = new EulerStepInterpolator();
- interpolator.reinitialize(new DummyEquations(), y, yDot, true);
+ interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true);
interpolator.storeTime(0);
interpolator.shift();
interpolator.storeTime(1);
@@ -172,14 +171,14 @@
}
- private static class DummyEquations
- implements FirstOrderDifferentialEquations {
- private static final long serialVersionUID = 291437140744677100L;
- public int getDimension() {
- return 0;
- }
- public void computeDerivatives(double t, double[] y, double[] yDot) {
- }
+ private static class DummyIntegrator extends RungeKuttaIntegrator {
+
+ private static final long serialVersionUID = -6936405965711773334L;
+
+ protected DummyIntegrator(RungeKuttaStepInterpolator prototype) {
+ super("dummy", new double[0], new double[0][0], new double[0], prototype, Double.NaN);
+ }
+
}
}
Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java Tue Jun 9 19:33:19 2009
@@ -141,6 +141,7 @@
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
int calls = pb.getCalls();
+ assertEquals(integ.getEvaluations(), calls);
assertTrue(calls <= previousCalls);
previousCalls = calls;
@@ -236,6 +237,7 @@
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
+ assertEquals(integ.getEvaluations(), pb.getCalls());
assertTrue(pb.getCalls() < 2150);
}
Modified: commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java?rev=783103&r1=783102&r2=783103&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java Tue Jun 9 19:33:19 2009
@@ -130,6 +130,7 @@
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
int calls = pb.getCalls();
+ assertEquals(integ.getEvaluations(), calls);
assertTrue(calls <= previousCalls);
previousCalls = calls;