commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject 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
Date Mon, 14 Jul 2008 15:08:52 GMT
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);
             }
           }
 



Mime
View raw message