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

prevent zero-length steps from generating NaN 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/branches/MATH_2_0@675578 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-07-10 14:06:17 +00:00
parent e7a11765d4
commit df6766e1e3
3 changed files with 36 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
View File

@ -128,10 +128,10 @@ class DormandPrince54StepInterpolator
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 @@ class DormandPrince54StepInterpolator
}
// 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];
}
}

38
src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
View File

@ -156,16 +156,16 @@ class DormandPrince853StepInterpolator
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 @@ class DormandPrince853StepInterpolator
}
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 @@ class DormandPrince853StepInterpolator
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];
}
}

8
src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
View File

@ -229,7 +229,7 @@ class GraggBulirschStoerStepInterpolator
/** 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 @@ class GraggBulirschStoerStepInterpolator
}
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} */