diff --git a/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java b/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java index 71cdce710..336d89e22 100644 --- a/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java +++ b/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java @@ -249,19 +249,20 @@ public abstract class AdaptiveStepsizeIntegrator /** Filter the integration step. * @param h signed step + * @param forward forward integration indicator * @param acceptSmall if true, steps smaller than the minimal value * are silently increased up to this value, if false such small * steps generate an exception * @return a bounded integration step (h if no bound is reach, or a bounded value) * @exception IntegratorException if the step is too small and acceptSmall is false */ - protected double filterStep(final double h, final boolean acceptSmall) + protected double filterStep(final double h, final boolean forward, final boolean acceptSmall) throws IntegratorException { double filteredH = h; if (Math.abs(h) < minStep) { if (acceptSmall) { - filteredH = (filteredH < 0) ? -minStep : minStep; + filteredH = forward ? minStep : -minStep; } else { throw new IntegratorException("minimal step size ({0}) reached," + " integration needs {1}", @@ -274,7 +275,7 @@ public abstract class AdaptiveStepsizeIntegrator if (filteredH > maxStep) { filteredH = maxStep; - } else if (h < -maxStep) { + } else if (filteredH < -maxStep) { filteredH = -maxStep; } diff --git a/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java b/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java index db5372063..8e9686edc 100644 --- a/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java +++ b/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java @@ -274,7 +274,7 @@ public abstract class EmbeddedRungeKuttaIntegrator final double factor = Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp))); - hNew = filterStep(stepSize * factor, false); + hNew = filterStep(stepSize * factor, forward, false); } } @@ -304,6 +304,11 @@ public abstract class EmbeddedRungeKuttaIntegrator } if (! lastStep) { + // in some rare cases we may get here with stepSize = 0, for example + // when an event occurs at integration start, reducing the first step + // to zero; we have to reset the step to some safe non zero value + stepSize = filterStep(stepSize, forward, true); + // stepsize control for next step final double factor = Math.min(maxGrowth, Math.max(minReduction, @@ -311,7 +316,7 @@ public abstract class EmbeddedRungeKuttaIntegrator final double scaledH = stepSize * factor; final double nextT = stepStart + scaledH; final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t); - hNew = filterStep(scaledH, nextIsLast); + hNew = filterStep(scaledH, forward, nextIsLast); } } diff --git a/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java b/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java index 73e89b3ac..ceb387002 100644 --- a/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java +++ b/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java @@ -638,7 +638,7 @@ public class GraggBulirschStoerIntegrator yTmp)) { // the stability check failed, we reduce the global step - hNew = Math.abs(filterStep(stepSize * stabilityReduction, false)); + hNew = Math.abs(filterStep(stepSize * stabilityReduction, forward, false)); reject = true; loop = false; @@ -662,7 +662,7 @@ public class GraggBulirschStoerIntegrator if ((error > 1.0e15) || ((k > 1) && (error > maxError))) { // error is too big, we reduce the global step - hNew = Math.abs(filterStep(stepSize * stabilityReduction, false)); + hNew = Math.abs(filterStep(stepSize * stabilityReduction, forward, false)); reject = true; loop = false; } else { @@ -674,7 +674,7 @@ public class GraggBulirschStoerIntegrator double fac = stepControl2 / Math.pow(error / stepControl1, exp); final double pow = Math.pow(stepControl3, exp); fac = Math.max(pow / stepControl4, Math.min(1 / pow, fac)); - optimalStep[k] = Math.abs(filterStep(stepSize * fac, true)); + optimalStep[k] = Math.abs(filterStep(stepSize * fac, forward, true)); costPerTimeUnit[k] = costPerStep[k] / optimalStep[k]; // check convergence @@ -903,13 +903,11 @@ public class GraggBulirschStoerIntegrator } else { if ((k < targetIter) && (costPerTimeUnit[k] < orderControl2 * costPerTimeUnit[k-1])) { - hNew = filterStep(optimalStep[k] * - costPerStep[optimalIter+1] / costPerStep[k], - false); + hNew = filterStep(optimalStep[k] * costPerStep[optimalIter+1] / costPerStep[k], + forward, false); } else { - hNew = filterStep(optimalStep[k] * - costPerStep[optimalIter] / costPerStep[k], - false); + hNew = filterStep(optimalStep[k] * costPerStep[optimalIter] / costPerStep[k], + forward, false); } }