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 2008/07/14 17:08:52 UTC
svn commit: r676615 - in
/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff:
AdaptiveStepsizeIntegrator.java EmbeddedRungeKuttaIntegrator.java
GraggBulirschStoerIntegrator.java
Author: luc
Date: Mon Jul 14 08:08:51 2008
New Revision: 676615
URL: http://svn.apache.org/viewvc?rev=676615&view=rev
Log:
fixed step size handling in borderline cases.
When an even occurred at step start, the step size dropped to zero
which put integration in an infinite loop
Modified:
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java?rev=676615&r1=676614&r2=676615&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java Mon Jul 14 08:08:51 2008
@@ -249,19 +249,20 @@
/** 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 @@
if (filteredH > maxStep) {
filteredH = maxStep;
- } else if (h < -maxStep) {
+ } else if (filteredH < -maxStep) {
filteredH = -maxStep;
}
Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java?rev=676615&r1=676614&r2=676615&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java Mon Jul 14 08:08:51 2008
@@ -274,7 +274,7 @@
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 @@
}
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 @@
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);
}
}
Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java?rev=676615&r1=676614&r2=676615&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java (original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java Mon Jul 14 08:08:51 2008
@@ -638,7 +638,7 @@
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 @@
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 @@
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 @@
} 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);
}
}