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/10 16:06:17 UTC
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
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} */