From df6766e1e3a569163b52dbda8e7cfbc2998e3742 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Thu, 10 Jul 2008 14:06:17 +0000 Subject: [PATCH] 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 --- .../DormandPrince54StepInterpolator.java | 17 ++++----- .../DormandPrince853StepInterpolator.java | 38 ++++++++++--------- .../GraggBulirschStoerStepInterpolator.java | 8 +++- 3 files changed, 36 insertions(+), 27 deletions(-) diff --git a/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java b/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java index 40c84c5ad..00a648dc4 100644 --- a/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java +++ b/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java @@ -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]; + } } diff --git a/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java b/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java index 05b39d077..bf7d7a755 100644 --- a/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java +++ b/src/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java @@ -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]; } } diff --git a/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java b/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java index 1be91fcb0..5ff680a7d 100644 --- a/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java +++ b/src/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java @@ -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} */