commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r675578 - in /commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff: DormandPrince54StepInterpolator.java DormandPrince853StepInterpolator.java GraggBulirschStoerStepInterpolator.java
Date Thu, 10 Jul 2008 14:06:17 GMT
Author: luc
Date: Thu Jul 10 07:06:17 2008
New Revision: 675578

URL: http://svn.apache.org/viewvc?rev=675578&view=rev
Log:
prevent zero-length steps from generating NaN

Modified:
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
    commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java

Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java?rev=675578&r1=675577&r2=675578&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
Thu Jul 10 07:06:17 2008
@@ -128,10 +128,10 @@
           final double yDot4 = yDotK[4][i];
           final double yDot5 = yDotK[5][i];
           final double yDot6 = yDotK[6][i];
-          v1[i] = h * (a70 * yDot0 + a72 * yDot2 + a73 * yDot3 + a74 * yDot4 + a75 * yDot5);
-          v2[i] = h * yDot0 - v1[i];
-          v3[i] = v1[i] - v2[i] - h * yDot6;
-          v4[i] = h * (d0 * yDot0 + d2 * yDot2 + d3 * yDot3 + d4 * yDot4 + d5 * yDot5 + d6
* yDot6);
+          v1[i] = a70 * yDot0 + a72 * yDot2 + a73 * yDot3 + a74 * yDot4 + a75 * yDot5;
+          v2[i] = yDot0 - v1[i];
+          v3[i] = v1[i] - v2[i] - yDot6;
+          v4[i] = d0 * yDot0 + d2 * yDot2 + d3 * yDot3 + d4 * yDot4 + d5 * yDot5 + d6 * yDot6;
       }
 
       vectorsInitialized = true;
@@ -139,17 +139,16 @@
     }
 
     // interpolate
-    final double eta = oneMinusThetaH / h;
+    final double eta = 1 - theta;
     final double twoTheta = 2 * theta;
     final double dot2 = 1 - twoTheta;
     final double dot3 = theta * (2 - 3 * theta);
     final double dot4 = twoTheta * (1 + theta * (twoTheta - 3));
     for (int i = 0; i < interpolatedState.length; ++i) {
       interpolatedState[i] =
-          currentState[i] - eta * (v1[i] - theta * (v2[i] + theta * (v3[i] + eta * v4[i])));
-      interpolatedDerivatives[i] = 
-          (v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i]) / h;
-      }
+          currentState[i] - oneMinusThetaH * (v1[i] - theta * (v2[i] + theta * (v3[i] + eta
* v4[i])));
+      interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i];
+    }
 
   }
 

Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java?rev=675578&r1=675577&r2=675578&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
Thu Jul 10 07:06:17 2008
@@ -156,16 +156,16 @@
           final double yDot14 = yDotKLast[0][i];
           final double yDot15 = yDotKLast[1][i];
           final double yDot16 = yDotKLast[2][i];
-          v[0][i] = h * (b_01 * yDot1  + b_06 * yDot6 + b_07 * yDot7 +
-                         b_08 * yDot8  + b_09 * yDot9 + b_10 * yDot10 +
-                         b_11 * yDot11 + b_12 * yDot12);
-          v[1][i] = h * yDot1 - v[0][i];
-          v[2][i] = v[0][i] - v[1][i] - h * yDotK[12][i];
+          v[0][i] = b_01 * yDot1  + b_06 * yDot6 + b_07 * yDot7 +
+                    b_08 * yDot8  + b_09 * yDot9 + b_10 * yDot10 +
+                    b_11 * yDot11 + b_12 * yDot12;
+          v[1][i] = yDot1 - v[0][i];
+          v[2][i] = v[0][i] - v[1][i] - yDotK[12][i];
           for (int k = 0; k < d.length; ++k) {
-              v[k+3][i] = h * (d[k][0] * yDot1  + d[k][1]  * yDot6  + d[k][2]  * yDot7  +
-                               d[k][3] * yDot8  + d[k][4]  * yDot9  + d[k][5]  * yDot10 +
-                               d[k][6] * yDot11 + d[k][7]  * yDot12 + d[k][8]  * yDot13 +
-                               d[k][9] * yDot14 + d[k][10] * yDot15 + d[k][11] * yDot16);
+              v[k+3][i] = d[k][0] * yDot1  + d[k][1]  * yDot6  + d[k][2]  * yDot7  +
+                          d[k][3] * yDot8  + d[k][4]  * yDot9  + d[k][5]  * yDot10 +
+                          d[k][6] * yDot11 + d[k][7]  * yDot12 + d[k][8]  * yDot13 +
+                          d[k][9] * yDot14 + d[k][10] * yDot15 + d[k][11] * yDot16;
           }
       }
 
@@ -173,7 +173,7 @@
 
     }
 
-    final double eta      = oneMinusThetaH / h;
+    final double eta      = 1 - theta;
     final double twoTheta = 2 * theta;
     final double theta2   = theta * theta;
     final double dot1 = 1 - twoTheta;
@@ -184,13 +184,17 @@
     final double dot6 = theta2 * theta * (4 + theta * (-15 + theta * (18 - 7 * theta)));
 
     for (int i = 0; i < interpolatedState.length; ++i) {
-      interpolatedState[i] =
-          currentState[i] - eta * (v[0][i] - theta * (v[1][i] +
-                  theta * (v[2][i] + eta * (v[3][i] + theta * (v[4][i] +
-                          eta * (v[5][i] + theta * (v[6][i])))))));
-      interpolatedDerivatives[i] = 
-          (v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] + dot3 * v[3][i] +
-                     dot4 * v[4][i] + dot5 * v[5][i] + dot6 * v[6][i]) / h;
+      interpolatedState[i] = currentState[i] -
+                             oneMinusThetaH * (v[0][i] -
+                                               theta * (v[1][i] +
+                                                        theta * (v[2][i] +
+                                                                 eta * (v[3][i] +
+                                                                        theta * (v[4][i]
+
+                                                                                 eta * (v[5][i]
+
+                                                                                        theta
* (v[6][i])))))));
+      interpolatedDerivatives[i] =  v[0][i] + dot1 * v[1][i] + dot2 * v[2][i] +
+                                    dot3 * v[3][i] + dot4 * v[4][i] +
+                                    dot5 * v[5][i] + dot6 * v[6][i];
     }
 
   }

Modified: commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java?rev=675578&r1=675577&r2=675578&view=diff
==============================================================================
--- commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
(original)
+++ commons/proper/math/branches/MATH_2_0/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
Thu Jul 10 07:06:17 2008
@@ -229,7 +229,7 @@
 
 
   /** Compute the interpolation coefficients for dense output.
-   * @param mu degree of the interpolation polynom
+   * @param mu degree of the interpolation polynomial
    * @param h current step
    */
   public void computeCoefficients(final int mu, final double h) {
@@ -342,6 +342,12 @@
 
     }
 
+    if (h == 0) {
+        // in this degenerated case, the previous computation leads to NaN for derivatives
+        // we fix this by using the derivatives at midpoint
+        System.arraycopy(yMidDots[1], 0, interpolatedDerivatives, 0, dimension);
+    }
+
   }
     
   /** {@inheritDoc} */



Mime
View raw message