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

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

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_0@676615 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-07-14 15:08:51 +00:00
parent 87d350f022
commit e1306edb51
3 changed files with 18 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
src/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
View File

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

9
src/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
View File

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

16
src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
View File

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