commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1206723 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ode/nonstiff/ site/xdoc/ test/java/org/apache/commons/math/ode/nonstiff/
Date Sun, 27 Nov 2011 14:32:03 GMT
Author: luc
Date: Sun Nov 27 14:32:00 2011
New Revision: 1206723

URL: http://svn.apache.org/viewvc?rev=1206723&view=rev
Log:
Improved accuracy of Runge-Kutta based step interpolators near step
start.

JIRA: MATH-705

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolator.java
Sun Nov 27 14:32:00 2011
@@ -26,16 +26,25 @@ import org.apache.commons.math.ode.sampl
  * <p>This interpolator allows to compute dense output inside the last
  * step computed. The interpolation equation is consistent with the
  * integration scheme :
-
- * <pre>
- *   y(t_n + theta h) = y (t_n + h)
- *                    + (1 - theta) (h/6) [ (-4 theta^2 + 5 theta - 1) y'_1
- *                                          +(4 theta^2 - 2 theta - 2) (y'_2 + y'_3)
- *                                          -(4 theta^2 +   theta + 1) y'_4
+ * <ul>
+ *   <li>Using reference point at step start:<br>
+ *   y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub>)
+ *                    + &theta; (h/6) [  (6 - 9 &theta; + 4 &theta;<sup>2</sup>)
y'<sub>1</sub>
+ *                                     + (    6 &theta; - 4 &theta;<sup>2</sup>)
(y'<sub>2</sub> + y'<sub>3</sub>)
+ *                                     + (   -3 &theta; + 4 &theta;<sup>2</sup>)
y'<sub>4</sub>
+ *                                    ]
+ *   </li>
+ *   <li>Using reference point at step end:<br>
+ *   y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub> + h)
+ *                    + (1 - &theta;) (h/6) [ (-4 &theta;^2 + 5 &theta; - 1)
y'<sub>1</sub>
+ *                                          +(4 &theta;^2 - 2 &theta; - 2) (y'<sub>2</sub>
+ y'<sub>3</sub>)
+ *                                          -(4 &theta;^2 +   &theta; + 1) y'<sub>4</sub>
  *                                        ]
- * </pre>
+ *   </li>
+ * </ul>
+ * </p>
  *
- * where theta belongs to [0 ; 1] and where y'_1 to y'_4 are the four
+ * where &theta; belongs to [0 ; 1] and where y'<sub>1</sub> to y'<sub>4</sub>
are the four
  * evaluations of the derivatives already computed during the
  * step.</p>
  *
@@ -83,24 +92,41 @@ class ClassicalRungeKuttaStepInterpolato
     protected void computeInterpolatedStateAndDerivatives(final double theta,
                                             final double oneMinusThetaH) {
 
-        final double fourTheta      = 4 * theta;
         final double oneMinusTheta  = 1 - theta;
         final double oneMinus2Theta = 1 - 2 * theta;
-        final double s             = oneMinusThetaH / 6.0;
-        final double coeff1        = s * ((-fourTheta + 5) * theta - 1);
-        final double coeff23       = s * (( fourTheta - 2) * theta - 2);
-        final double coeff4        = s * ((-fourTheta - 1) * theta - 1);
         final double coeffDot1     = oneMinusTheta * oneMinus2Theta;
         final double coeffDot23    = 2 * theta * oneMinusTheta;
         final double coeffDot4     = -theta * oneMinus2Theta;
-        for (int i = 0; i < interpolatedState.length; ++i) {
-            final double yDot1  = yDotK[0][i];
-            final double yDot23 = yDotK[1][i] + yDotK[2][i];
-            final double yDot4  = yDotK[3][i];
-            interpolatedState[i] =
-                currentState[i] + coeff1  * yDot1 + coeff23 * yDot23 + coeff4  * yDot4;
-            interpolatedDerivatives[i] =
-                coeffDot1 * yDot1 + coeffDot23 * yDot23 + coeffDot4 * yDot4;
+        if ((previousState != null) && (theta <= 0.5)) {
+            final double fourTheta2     = 4 * theta * theta;
+            final double s             = theta * h / 6.0;
+            final double coeff1        = s * ( 6 - 9 * theta + fourTheta2);
+            final double coeff23       = s * ( 6 * theta - fourTheta2);
+            final double coeff4        = s * (-3 * theta + fourTheta2);
+            for (int i = 0; i < interpolatedState.length; ++i) {
+                final double yDot1  = yDotK[0][i];
+                final double yDot23 = yDotK[1][i] + yDotK[2][i];
+                final double yDot4  = yDotK[3][i];
+                interpolatedState[i] =
+                        previousState[i] + coeff1  * yDot1 + coeff23 * yDot23 + coeff4  *
yDot4;
+                interpolatedDerivatives[i] =
+                        coeffDot1 * yDot1 + coeffDot23 * yDot23 + coeffDot4 * yDot4;
+            }
+        } else {
+            final double fourTheta      = 4 * theta;
+            final double s             = oneMinusThetaH / 6.0;
+            final double coeff1        = s * ((-fourTheta + 5) * theta - 1);
+            final double coeff23       = s * (( fourTheta - 2) * theta - 2);
+            final double coeff4        = s * ((-fourTheta - 1) * theta - 1);
+            for (int i = 0; i < interpolatedState.length; ++i) {
+                final double yDot1  = yDotK[0][i];
+                final double yDot23 = yDotK[1][i] + yDotK[2][i];
+                final double yDot4  = yDotK[3][i];
+                interpolatedState[i] =
+                        currentState[i] + coeff1  * yDot1 + coeff23 * yDot23 + coeff4  *
yDot4;
+                interpolatedDerivatives[i] =
+                        coeffDot1 * yDot1 + coeffDot23 * yDot23 + coeffDot4 * yDot4;
+            }
         }
 
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
Sun Nov 27 14:32:00 2011
@@ -204,10 +204,18 @@ class DormandPrince54StepInterpolator
     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] - 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];
+    if ((previousState != null) && (theta <= 0.5)) {
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            interpolatedState[i] =
+                    previousState[i] + theta * h * (v1[i] + eta * (v2[i] + theta * (v3[i]
+ eta * v4[i])));
+            interpolatedDerivatives[i] = v1[i] + dot2 * v2[i] + dot3 * v3[i] + dot4 * v4[i];
+        }
+    } else {
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            interpolatedState[i] =
+                    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/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolator.java
Sun Nov 27 14:32:00 2011
@@ -25,12 +25,17 @@ import org.apache.commons.math.ode.sampl
  * <p>This interpolator computes dense output inside the last
  * step computed. The interpolation equation is consistent with the
  * integration scheme :
+ * <ul>
+ *   <li>Using reference point at step start:<br>
+ *     y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub>) + &theta;
h y'
+ *   </li>
+ *   <li>Using reference point at step end:<br>
+ *     y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub> + h) -
(1-&theta;) h y'
+ *   </li>
+ * </ul>
+ * </p>
  *
- * <pre>
- *   y(t_n + theta h) = y (t_n + h) - (1-theta) h y'
- * </pre>
- *
- * where theta belongs to [0 ; 1] and where y' is the evaluation of
+ * where &theta; belongs to [0 ; 1] and where y' is the evaluation of
  * the derivatives already computed during the step.</p>
  *
  * @see EulerIntegrator
@@ -78,11 +83,17 @@ class EulerStepInterpolator
   @Override
   protected void computeInterpolatedStateAndDerivatives(final double theta,
                                           final double oneMinusThetaH) {
-
-    for (int i = 0; i < interpolatedState.length; ++i) {
-      interpolatedState[i] = currentState[i] - oneMinusThetaH * yDotK[0][i];
-    }
-    System.arraycopy(yDotK[0], 0, interpolatedDerivatives, 0, interpolatedDerivatives.length);
+      if ((previousState != null) && (theta <= 0.5)) {
+          for (int i = 0; i < interpolatedState.length; ++i) {
+              interpolatedState[i] = previousState[i] + theta * h * yDotK[0][i];
+          }
+          System.arraycopy(yDotK[0], 0, interpolatedDerivatives, 0, interpolatedDerivatives.length);
+      } else {
+          for (int i = 0; i < interpolatedState.length; ++i) {
+              interpolatedState[i] = currentState[i] - oneMinusThetaH * yDotK[0][i];
+          }
+          System.arraycopy(yDotK[0], 0, interpolatedDerivatives, 0, interpolatedDerivatives.length);
+      }
 
   }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolator.java
Sun Nov 27 14:32:00 2011
@@ -27,15 +27,24 @@ import org.apache.commons.math.util.Fast
  * <p>This interpolator allows to compute dense output inside the last
  * step computed. The interpolation equation is consistent with the
  * integration scheme :
- *
- * <pre>
- *   y(t_n + theta h) = y (t_n + h)
- *                    - (1 - theta) (h/6) [ (1 - theta) (1 - 4 theta) y'_1
- *                                        + (1 - theta) (1 + 2 theta) ((2-q) y'_2 + (2+q)
y'_3)
- *                                        + (1 + theta + 4 theta^2) y'_4
- *                                        ]
- * </pre>
- * where theta belongs to [0 ; 1], q = sqrt(2) and where y'_1 to y'_4
+ * <ul>
+ *   <li>Using reference point at step start:<br>
+ *   y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub>)
+ *                    + &theta; (h/6) [ (6 - 9 &theta; + 4 &theta;<sup>2</sup>)
y'<sub>1</sub>
+ *                                    + (    6 &theta; - 4 &theta;<sup>2</sup>)
((1-1/&radic;2) y'<sub>2</sub> + (1+1/&radic;2)) y'<sub>3</sub>)
+ *                                    + (  - 3 &theta; + 4 &theta;<sup>2</sup>)
y'<sub>4</sub>
+ *                                    ]
+ *   </li>
+ *   <li>Using reference point at step start:<br>
+ *   y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub> + h)
+ *                    - (1 - &theta;) (h/6) [ (1 - 5 &theta; + 4 &theta;<sup>2</sup>)
y'<sub>1</sub>
+ *                                          + (2 + 2 &theta; - 4 &theta;<sup>2</sup>)
((1-1/&radic;2) y'<sub>2</sub> + (1+1/&radic;2)) y'<sub>3</sub>)
+ *                                          + (1 +   &theta; + 4 &theta;<sup>2</sup>)
y'<sub>4</sub>
+ *                                          ]
+ *   </li>
+ * </ul>
+ * </p>
+ * where &theta; belongs to [0 ; 1] and where y'<sub>1</sub> to y'<sub>4</sub>
  * are the four evaluations of the derivatives already computed during
  * the step.</p>
  *
@@ -48,10 +57,10 @@ class GillStepInterpolator
   extends RungeKuttaStepInterpolator {
 
     /** First Gill coefficient. */
-    private static final double TWO_MINUS_SQRT_2 = 2 - FastMath.sqrt(2.0);
+    private static final double ONE_MINUS_INV_SQRT_2 = 1 - FastMath.sqrt(0.5);
 
     /** Second Gill coefficient. */
-    private static final double TWO_PLUS_SQRT_2 = 2 + FastMath.sqrt(2.0);
+    private static final double ONE_PLUS_INV_SQRT_2 = 1 + FastMath.sqrt(0.5);
 
     /** Serializable version identifier. */
     private static final long serialVersionUID = 20111120L;
@@ -91,32 +100,49 @@ class GillStepInterpolator
   protected void computeInterpolatedStateAndDerivatives(final double theta,
                                           final double oneMinusThetaH) {
 
-    final double twoTheta  = 2 * theta;
-    final double fourTheta = 4 * theta;
-    final double s         = oneMinusThetaH / 6.0;
-    final double oMt       = 1 - theta;
-    final double soMt      = s * oMt;
-    final double c23       = soMt * (1 + twoTheta);
-    final double coeff1    = soMt * (1 - fourTheta);
-    final double coeff2    = c23  * TWO_MINUS_SQRT_2;
-    final double coeff3    = c23  * TWO_PLUS_SQRT_2;
-    final double coeff4    = s * (1 + theta * (1 + fourTheta));
-    final double coeffDot1 = theta * (twoTheta - 3) + 1;
-    final double cDot23    = theta * oMt;
-    final double coeffDot2 = cDot23  * TWO_MINUS_SQRT_2;
-    final double coeffDot3 = cDot23  * TWO_PLUS_SQRT_2;
-    final double coeffDot4 = theta * (twoTheta - 1);
-
-    for (int i = 0; i < interpolatedState.length; ++i) {
-        final double yDot1 = yDotK[0][i];
-        final double yDot2 = yDotK[1][i];
-        final double yDot3 = yDotK[2][i];
-        final double yDot4 = yDotK[3][i];
-        interpolatedState[i] =
-            currentState[i] - coeff1 * yDot1 - coeff2 * yDot2 - coeff3 * yDot3 - coeff4 *
yDot4;
-        interpolatedDerivatives[i] =
-            coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4 * yDot4;
-     }
+    final double twoTheta   = 2 * theta;
+    final double fourTheta2 = twoTheta * twoTheta;
+    final double coeffDot1  = theta * (twoTheta - 3) + 1;
+    final double cDot23     = twoTheta * (1 - theta);
+    final double coeffDot2  = cDot23  * ONE_MINUS_INV_SQRT_2;
+    final double coeffDot3  = cDot23  * ONE_PLUS_INV_SQRT_2;
+    final double coeffDot4  = theta * (twoTheta - 1);
+
+    if ((previousState != null) && (theta <= 0.5)) {
+        final double s         = theta * h / 6.0;
+        final double c23       = s * (6 * theta - fourTheta2);
+        final double coeff1    = s * (6 - 9 * theta + fourTheta2);
+        final double coeff2    = c23  * ONE_MINUS_INV_SQRT_2;
+        final double coeff3    = c23  * ONE_PLUS_INV_SQRT_2;
+        final double coeff4    = s * (-3 * theta + fourTheta2);
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            final double yDot1 = yDotK[0][i];
+            final double yDot2 = yDotK[1][i];
+            final double yDot3 = yDotK[2][i];
+            final double yDot4 = yDotK[3][i];
+            interpolatedState[i] =
+                    previousState[i] + coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3 +
coeff4 * yDot4;
+            interpolatedDerivatives[i] =
+                    coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4
* yDot4;
+        }
+    } else {
+        final double s      = oneMinusThetaH / 6.0;
+        final double c23    = s * (2 + twoTheta - fourTheta2);
+        final double coeff1 = s * (1 - 5 * theta + fourTheta2);
+        final double coeff2 = c23  * ONE_MINUS_INV_SQRT_2;
+        final double coeff3 = c23  * ONE_PLUS_INV_SQRT_2;
+        final double coeff4 = s * (1 + theta + fourTheta2);
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            final double yDot1 = yDotK[0][i];
+            final double yDot2 = yDotK[1][i];
+            final double yDot3 = yDotK[2][i];
+            final double yDot4 = yDotK[3][i];
+            interpolatedState[i] =
+                    currentState[i] - coeff1 * yDot1 - coeff2 * yDot2 - coeff3 * yDot3 -
coeff4 * yDot4;
+            interpolatedDerivatives[i] =
+                    coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4
* yDot4;
+        }
+    }
 
   }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolator.java
Sun Nov 27 14:32:00 2011
@@ -71,29 +71,48 @@ class HighamHall54StepInterpolator
   protected void computeInterpolatedStateAndDerivatives(final double theta,
                                           final double oneMinusThetaH) {
 
-    final double theta2 = theta * theta;
-
-    final double b0 = h * (-1.0/12.0 + theta * (1.0 + theta * (-15.0/4.0 + theta * (16.0/3.0
+ theta * -5.0/2.0))));
-    final double b2 = h * (-27.0/32.0 + theta2 * (459.0/32.0 + theta * (-243.0/8.0 + theta
* 135.0/8.0)));
-    final double b3 = h * (4.0/3.0 + theta2 * (-22.0 + theta * (152.0/3.0  + theta * -30.0)));
-    final double b4 = h * (-125.0/96.0 + theta2 * (375.0/32.0 + theta * (-625.0/24.0 + theta
* 125.0/8.0)));
-    final double b5 = h * (-5.0/48.0 + theta2 * (-5.0/16.0 + theta * 5.0/12.0));
     final double bDot0 = 1 + theta * (-15.0/2.0 + theta * (16.0 - 10.0 * theta));
     final double bDot2 = theta * (459.0/16.0 + theta * (-729.0/8.0 + 135.0/2.0 * theta));
     final double bDot3 = theta * (-44.0 + theta * (152.0 - 120.0 * theta));
     final double bDot4 = theta * (375.0/16.0 + theta * (-625.0/8.0 + 125.0/2.0 * theta));
     final double bDot5 = theta * 5.0/8.0 * (2 * theta - 1);
 
-    for (int i = 0; i < interpolatedState.length; ++i) {
-        final double yDot0 = yDotK[0][i];
-        final double yDot2 = yDotK[2][i];
-        final double yDot3 = yDotK[3][i];
-        final double yDot4 = yDotK[4][i];
-        final double yDot5 = yDotK[5][i];
-        interpolatedState[i] =
-            currentState[i] + b0 * yDot0 + b2 * yDot2 + b3 * yDot3 + b4 * yDot4 + b5 * yDot5;
-        interpolatedDerivatives[i] =
-            bDot0 * yDot0 + bDot2 * yDot2 + bDot3 * yDot3 + bDot4 * yDot4 + bDot5 * yDot5;
+    if ((previousState != null) && (theta <= 0.5)) {
+        final double hTheta = h * theta;
+        final double b0 = hTheta * (1.0 + theta * (-15.0/4.0  + theta * (16.0/3.0 - 5.0/2.0
* theta)));
+        final double b2 = hTheta * (      theta * (459.0/32.0 + theta * (-243.0/8.0 + theta
* 135.0/8.0)));
+        final double b3 = hTheta * (      theta * (-22.0      + theta * (152.0/3.0  + theta
* -30.0)));
+        final double b4 = hTheta * (      theta * (375.0/32.0 + theta * (-625.0/24.0 + theta
* 125.0/8.0)));
+        final double b5 = hTheta * (      theta * (-5.0/16.0  + theta *  5.0/12.0));
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            final double yDot0 = yDotK[0][i];
+            final double yDot2 = yDotK[2][i];
+            final double yDot3 = yDotK[3][i];
+            final double yDot4 = yDotK[4][i];
+            final double yDot5 = yDotK[5][i];
+            interpolatedState[i] =
+                    previousState[i] + b0 * yDot0 + b2 * yDot2 + b3 * yDot3 + b4 * yDot4
+ b5 * yDot5;
+            interpolatedDerivatives[i] =
+                    bDot0 * yDot0 + bDot2 * yDot2 + bDot3 * yDot3 + bDot4 * yDot4 + bDot5
* yDot5;
+        }
+    } else {
+        final double theta2 = theta * theta;
+        final double b0 = h * (-1.0/12.0 + theta * (1.0 + theta * (-15.0/4.0 + theta * (16.0/3.0
+ theta * -5.0/2.0))));
+        final double b2 = h * (-27.0/32.0 + theta2 * (459.0/32.0 + theta * (-243.0/8.0 +
theta * 135.0/8.0)));
+        final double b3 = h * (4.0/3.0 + theta2 * (-22.0 + theta * (152.0/3.0  + theta *
-30.0)));
+        final double b4 = h * (-125.0/96.0 + theta2 * (375.0/32.0 + theta * (-625.0/24.0
+ theta * 125.0/8.0)));
+        final double b5 = h * (-5.0/48.0 + theta2 * (-5.0/16.0 + theta * 5.0/12.0));
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            final double yDot0 = yDotK[0][i];
+            final double yDot2 = yDotK[2][i];
+            final double yDot3 = yDotK[3][i];
+            final double yDot4 = yDotK[4][i];
+            final double yDot5 = yDotK[5][i];
+            interpolatedState[i] =
+                    currentState[i] + b0 * yDot0 + b2 * yDot2 + b3 * yDot3 + b4 * yDot4 +
b5 * yDot5;
+            interpolatedDerivatives[i] =
+                    bDot0 * yDot0 + bDot2 * yDot2 + bDot3 * yDot3 + bDot4 * yDot4 + bDot5
* yDot5;
+        }
     }
 
   }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolator.java
Sun Nov 27 14:32:00 2011
@@ -26,12 +26,17 @@ import org.apache.commons.math.ode.sampl
  * <p>This interpolator computes dense output inside the last
  * step computed. The interpolation equation is consistent with the
  * integration scheme :
+ * <ul>
+ *   <li>Using reference point at step start:<br>
+ *   y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub>) + &theta;
h [(1 - &theta;) y'<sub>1</sub> + &theta; y'<sub>2</sub>]
+ *   </li>
+ *   <li>Using reference point at step end:<br>
+ *   y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub> + h) + (1-&theta;)
h [&theta; y'<sub>1</sub> - (1+&theta;) y'<sub>2</sub>]
+ *   </li>
+ * </ul>
+ * </p>
  *
- * <pre>
- *   y(t_n + theta h) = y (t_n + h) + (1-theta) h [theta y'_1 - (1+theta) y'_2]
- * </pre>
- *
- * where theta belongs to [0 ; 1] and where y'_1 and y'_2 are the two
+ * where &theta; belongs to [0 ; 1] and where y'<sub>1</sub> and y'<sub>2</sub>
are the two
  * evaluations of the derivatives already computed during the
  * step.</p>
  *
@@ -81,16 +86,27 @@ class MidpointStepInterpolator
   protected void computeInterpolatedStateAndDerivatives(final double theta,
                                           final double oneMinusThetaH) {
 
-    final double coeff1    = oneMinusThetaH * theta;
-    final double coeff2    = oneMinusThetaH * (1.0 + theta);
     final double coeffDot2 = 2 * theta;
     final double coeffDot1 = 1 - coeffDot2;
 
-    for (int i = 0; i < interpolatedState.length; ++i) {
-      final double yDot1 = yDotK[0][i];
-      final double yDot2 = yDotK[1][i];
-      interpolatedState[i] = currentState[i] + coeff1 * yDot1 - coeff2 * yDot2;
-      interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2;
+    if ((previousState != null) && (theta <= 0.5)) {
+        final double coeff1    = theta * oneMinusThetaH;
+        final double coeff2    = theta * theta * h;
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            final double yDot1 = yDotK[0][i];
+            final double yDot2 = yDotK[1][i];
+            interpolatedState[i] = previousState[i] + coeff1 * yDot1 + coeff2 * yDot2;
+            interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2;
+        }
+    } else {
+        final double coeff1    = oneMinusThetaH * theta;
+        final double coeff2    = oneMinusThetaH * (1.0 + theta);
+        for (int i = 0; i < interpolatedState.length; ++i) {
+            final double yDot1 = yDotK[0][i];
+            final double yDot2 = yDotK[1][i];
+            interpolatedState[i] = currentState[i] + coeff1 * yDot1 - coeff2 * yDot2;
+            interpolatedDerivatives[i] = coeffDot1 * yDot1 + coeffDot2 * yDot2;
+        }
     }
 
   }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolator.java
Sun Nov 27 14:32:00 2011
@@ -26,17 +26,27 @@ import org.apache.commons.math.ode.sampl
  * <p>This interpolator allows to compute dense output inside the last
  * step computed. The interpolation equation is consistent with the
  * integration scheme :
+ * <ul>
+ *   <li>Using reference point at step start:<br>
+ *     y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub>)
+ *                      + &theta; (h/8) [ (8 - 15 &theta; +  8 &theta;<sup>2</sup>)
y'<sub>1</sub>
+ *                                     +  3 * (15 &theta; - 12 &theta;<sup>2</sup>)
y'<sub>2</sub>
+ *                                     +        3 &theta;                           y'<sub>3</sub>
+ *                                     +      (-3 &theta; +  4 &theta;<sup>2</sup>)
y'<sub>4</sub>
+ *                                    ]
+ *   </li>
+ *   <li>Using reference point at step end:<br>
+ *     y(t<sub>n</sub> + &theta; h) = y (t<sub>n</sub> + h)
+ *                      - (1 - &theta;) (h/8) [(1 - 7 &theta; + 8 &theta;<sup>2</sup>)
y'<sub>1</sub>
+ *                                         + 3 (1 +   &theta; - 4 &theta;<sup>2</sup>)
y'<sub>2</sub>
+ *                                         + 3 (1 +   &theta;)                      
  y'<sub>3</sub>
+ *                                         +   (1 +   &theta; + 4 &theta;<sup>2</sup>)
y'<sub>4</sub>
+ *                                          ]
+ *   </li>
+ * </ul>
+ * </p>
  *
- * <pre>
- *   y(t_n + theta h) = y (t_n + h)
- *                    - (1 - theta) (h/8) [ (1 - 7 theta + 8 theta^2) y'_1
- *                                      + 3 (1 +   theta - 4 theta^2) y'_2
- *                                      + 3 (1 +   theta)             y'_3
- *                                      +   (1 +   theta + 4 theta^2) y'_4
- *                                        ]
- * </pre>
- *
- * where theta belongs to [0 ; 1] and where y'_1 to y'_4 are the four
+ * where &theta; belongs to [0 ; 1] and where y'<sub>1</sub> to y'<sub>4</sub>
are the four
  * evaluations of the derivatives already computed during the
  * step.</p>
  *
@@ -86,27 +96,47 @@ class ThreeEighthesStepInterpolator
   protected void computeInterpolatedStateAndDerivatives(final double theta,
                                           final double oneMinusThetaH) {
 
-      final double fourTheta2 = 4 * theta * theta;
-      final double s          = oneMinusThetaH / 8.0;
-      final double coeff1     = s * (1 - 7 * theta + 2 * fourTheta2);
-      final double coeff2     = 3 * s * (1 + theta - fourTheta2);
-      final double coeff3     = 3 * s * (1 + theta);
-      final double coeff4     = s * (1 + theta + fourTheta2);
       final double coeffDot3  = 0.75 * theta;
       final double coeffDot1  = coeffDot3 * (4 * theta - 5) + 1;
       final double coeffDot2  = coeffDot3 * (5 - 6 * theta);
       final double coeffDot4  = coeffDot3 * (2 * theta - 1);
 
-      for (int i = 0; i < interpolatedState.length; ++i) {
-          final double yDot1 = yDotK[0][i];
-          final double yDot2 = yDotK[1][i];
-          final double yDot3 = yDotK[2][i];
-          final double yDot4 = yDotK[3][i];
-          interpolatedState[i] =
-              currentState[i] - coeff1 * yDot1 - coeff2 * yDot2 - coeff3 * yDot3 - coeff4
* yDot4;
-          interpolatedDerivatives[i] =
-              coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4 * yDot4;
+      if ((previousState != null) && (theta <= 0.5)) {
+          final double s          = theta * h / 8.0;
+          final double fourTheta2 = 4 * theta * theta;
+          final double coeff1     = s * (8 - 15 * theta + 2 * fourTheta2);
+          final double coeff2     = 3 * s * (5 * theta - fourTheta2);
+          final double coeff3     = 3 * s * theta;
+          final double coeff4     = s * (-3 * theta + fourTheta2);
+          for (int i = 0; i < interpolatedState.length; ++i) {
+              final double yDot1 = yDotK[0][i];
+              final double yDot2 = yDotK[1][i];
+              final double yDot3 = yDotK[2][i];
+              final double yDot4 = yDotK[3][i];
+              interpolatedState[i] =
+                      previousState[i] + coeff1 * yDot1 + coeff2 * yDot2 + coeff3 * yDot3
+ coeff4 * yDot4;
+              interpolatedDerivatives[i] =
+                      coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4
* yDot4;
+
+          }
+      } else {
+          final double s          = oneMinusThetaH / 8.0;
+          final double fourTheta2 = 4 * theta * theta;
+          final double coeff1     = s * (1 - 7 * theta + 2 * fourTheta2);
+          final double coeff2     = 3 * s * (1 + theta - fourTheta2);
+          final double coeff3     = 3 * s * (1 + theta);
+          final double coeff4     = s * (1 + theta + fourTheta2);
+          for (int i = 0; i < interpolatedState.length; ++i) {
+              final double yDot1 = yDotK[0][i];
+              final double yDot2 = yDotK[1][i];
+              final double yDot3 = yDotK[2][i];
+              final double yDot4 = yDotK[3][i];
+              interpolatedState[i] =
+                      currentState[i] - coeff1 * yDot1 - coeff2 * yDot2 - coeff3 * yDot3
- coeff4 * yDot4;
+              interpolatedDerivatives[i] =
+                      coeffDot1 * yDot1 + coeffDot2 * yDot2 + coeffDot3 * yDot3 + coeffDot4
* yDot4;
 
+          }
       }
 
   }

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Sun Nov 27 14:32:00 2011
@@ -52,6 +52,9 @@ The <action> type attribute can be add,u
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-705">
+        Improved accuracy of Runge-Kutta based step interpolators near step start.
+      </action>
       <action dev="psteitz" type="fix" issue="MATH-691">
         Fixed errors in SummaryStatistics addValue causing variance, mean, or
         geometric mean statistics not to be updated if they have been overriden

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java
Sun Nov 27 14:32:00 2011
@@ -150,7 +150,7 @@ public class ClassicalRungeKuttaIntegrat
 
         double error = handler.getMaximalValueError();
         if (i > 4) {
-          Assert.assertTrue(error < FastMath.abs(previousValueError));
+          Assert.assertTrue(error < 1.01 * FastMath.abs(previousValueError));
         }
         previousValueError = error;
 

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java
Sun Nov 27 14:32:00 2011
@@ -72,19 +72,19 @@ public class EulerStepInterpolatorTest {
     interpolator.storeTime(t0);
 
     double dt = 1.0;
+    interpolator.shift();
     y[0] =  1.0;
     y[1] =  3.0;
     y[2] = -4.0;
     yDot[0][0] = (y[0] - y0[0]) / dt;
     yDot[0][1] = (y[1] - y0[1]) / dt;
     yDot[0][2] = (y[2] - y0[2]) / dt;
-    interpolator.shift();
     interpolator.storeTime(t0 + dt);
 
     interpolator.setInterpolatedTime(interpolator.getPreviousTime());
     double[] result = interpolator.getInterpolatedState();
     for (int i = 0; i < result.length; ++i) {
-      Assert.assertTrue(FastMath.abs(result[i] - y0[i]) < 1.0e-10);
+        Assert.assertTrue(FastMath.abs(result[i] - y0[i]) < 1.0e-10);
     }
 
     interpolator.setInterpolatedTime(interpolator.getCurrentTime());
@@ -98,7 +98,7 @@ public class EulerStepInterpolatorTest {
   @Test
   public void interpolationInside() {
 
-    double[]   y    =   { 1.0, 3.0, -4.0 };
+    double[]   y    =   { 0.0, 1.0, -2.0 };
     double[][] yDot = { { 1.0, 2.0, -2.0 } };
     EulerStepInterpolator interpolator = new EulerStepInterpolator();
     interpolator.reinitialize(new DummyIntegrator(interpolator), y, yDot, true,
@@ -106,6 +106,9 @@ public class EulerStepInterpolatorTest {
                               new EquationsMapper[0]);
     interpolator.storeTime(0);
     interpolator.shift();
+    y[0] =  1.0;
+    y[1] =  3.0;
+    y[2] = -4.0;
     interpolator.storeTime(1);
 
     interpolator.setInterpolatedTime(0.1);

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java
Sun Nov 27 14:32:00 2011
@@ -75,7 +75,7 @@ public class GillIntegratorTest {
 
         double valueError = handler.getMaximalValueError();
         if (i > 5) {
-          Assert.assertTrue(valueError < FastMath.abs(previousValueError));
+          Assert.assertTrue(valueError < 1.01 * FastMath.abs(previousValueError));
         }
         previousValueError = valueError;
 

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java?rev=1206723&r1=1206722&r2=1206723&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java
Sun Nov 27 14:32:00 2011
@@ -75,7 +75,7 @@ public class ThreeEighthesIntegratorTest
 
         double error = handler.getMaximalValueError();
         if (i > 4) {
-          Assert.assertTrue(error < FastMath.abs(previousValueError));
+          Assert.assertTrue(error < 1.01 * FastMath.abs(previousValueError));
         }
         previousValueError = error;
 



Mime
View raw message