Merge branch 'master' of https://git-wip-us.apache.org/repos/asf/commons-math

src/changes/changes.xml

@@ -73,6 +73,9 @@ Users are encouraged to upgrade to this version as this release not
   2. A few methods in the FastMath class are in fact slower that their
   counterpart in either Math or StrictMath (cf. MATH-740 and MATH-901).
 ">
+      <action dev="luc" type="fix" issue="MATH-1160" >
+        Provide access to state derivatives in ContinuousOutputModel.
+      </action>
       <action dev="luc" type="fix" issue="MATH-1138" due-to="Hank Grabowski">
         Fixed bicubic spline interpolator, using Akima splines.
       </action>

src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java

@@ -233,12 +233,19 @@ public class ContinuousOutputModel
    * integration is over because some internal variables are set only
    * once the last step has been handled.</p>
    * <p>Setting the time outside of the integration interval is now
-   * allowed (it was not allowed up to version 5.9 of Mantissa), but
-   * should be used with care since the accuracy of the interpolator
-   * will probably be very poor far from this interval. This allowance
-   * has been added to simplify implementation of search algorithms
-   * near the interval endpoints.</p>
+   * allowed, but should be used with care since the accuracy of the
+   * interpolator will probably be very poor far from this interval.
+   * This allowance has been added to simplify implementation of search
+   * algorithms near the interval endpoints.</p>
+   * <p>Note that each time this method is called, the internal arrays
+   * returned in {@link #getInterpolatedState()}, {@link
+   * #getInterpolatedDerivatives()} and {@link #getInterpolatedSecondaryState(int)}
+   * <em>will</em> be overwritten. So if their content must be preserved
+   * across several calls, user must copy them.</p>
    * @param time time of the interpolated point
+   * @see #getInterpolatedState()
+   * @see #getInterpolatedDerivatives()
+   * @see #getInterpolatedSecondaryState(int)
    */
   public void setInterpolatedTime(final double time) {
 
@@ -330,21 +337,53 @@ public class ContinuousOutputModel
 
   /**
    * Get the state vector of the interpolated point.
+   * <p>The returned vector is a reference to a reused array, so
+   * it should not be modified and it should be copied if it needs
+   * to be preserved across several calls to the associated
+   * {@link #setInterpolatedTime(double)} method.</p>
    * @return state vector at time {@link #getInterpolatedTime}
    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
+   * @see #setInterpolatedTime(double)
+   * @see #getInterpolatedDerivatives()
    * @see #getInterpolatedSecondaryState(int)
+   * @see #getInterpolatedSecondaryDerivatives(int)
    */
   public double[] getInterpolatedState() throws MaxCountExceededException {
     return steps.get(index).getInterpolatedState();
   }
 
+  /**
+   * Get the derivatives of the state vector of the interpolated point.
+   * <p>The returned vector is a reference to a reused array, so
+   * it should not be modified and it should be copied if it needs
+   * to be preserved across several calls to the associated
+   * {@link #setInterpolatedTime(double)} method.</p>
+   * @return derivatives of the state vector at time {@link #getInterpolatedTime}
+   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
+   * @see #setInterpolatedTime(double)
+   * @see #getInterpolatedState()
+   * @see #getInterpolatedSecondaryState(int)
+   * @see #getInterpolatedSecondaryDerivatives(int)
+   * @since 3.4
+   */
+  public double[] getInterpolatedDerivatives() throws MaxCountExceededException {
+    return steps.get(index).getInterpolatedDerivatives();
+  }
+
   /** Get the interpolated secondary state corresponding to the secondary equations.
+   * <p>The returned vector is a reference to a reused array, so
+   * it should not be modified and it should be copied if it needs
+   * to be preserved across several calls to the associated
+   * {@link #setInterpolatedTime(double)} method.</p>
    * @param secondaryStateIndex index of the secondary set, as returned by {@link
    * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations(
    * org.apache.commons.math3.ode.SecondaryEquations)
    * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
    * @return interpolated secondary state at the current interpolation date
+   * @see #setInterpolatedTime(double)
    * @see #getInterpolatedState()
+   * @see #getInterpolatedDerivatives()
+   * @see #getInterpolatedSecondaryDerivatives(int)
    * @since 3.2
    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    */
@@ -353,6 +392,28 @@ public class ContinuousOutputModel
     return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex);
   }
 
+  /** Get the interpolated secondary derivatives corresponding to the secondary equations.
+   * <p>The returned vector is a reference to a reused array, so
+   * it should not be modified and it should be copied if it needs
+   * to be preserved across several calls to the associated
+   * {@link #setInterpolatedTime(double)} method.</p>
+   * @param secondaryStateIndex index of the secondary set, as returned by {@link
+   * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations(
+   * org.apache.commons.math3.ode.SecondaryEquations)
+   * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
+   * @return interpolated secondary derivatives at the current interpolation date
+   * @see #setInterpolatedTime(double)
+   * @see #getInterpolatedState()
+   * @see #getInterpolatedDerivatives()
+   * @see #getInterpolatedSecondaryState(int)
+   * @since 3.4
+   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
+   */
+  public double[] getInterpolatedSecondaryDerivatives(final int secondaryStateIndex)
+    throws MaxCountExceededException {
+    return steps.get(index).getInterpolatedSecondaryDerivatives(secondaryStateIndex);
+  }
+
   /** Compare a step interval and a double.
    * @param time point to locate
    * @param interval step interval

src/main/java/org/apache/commons/math3/ode/sampling/StepInterpolator.java

@@ -74,10 +74,17 @@ public interface StepInterpolator extends Externalizable {
    * probably be very poor far from this step. This allowance has been
    * added to simplify implementation of search algorithms near the
    * step endpoints.</p>
-   * <p>Setting the time changes the instance internal state. If a
-   * specific state must be preserved, a copy of the instance must be
-   * created using {@link #copy()}.</p>
+   * <p>Setting the time changes the instance internal state. This includes
+   * the internal arrays returned in {@link #getInterpolatedState()},
+   * {@link #getInterpolatedDerivatives()}, {@link
+   * #getInterpolatedSecondaryState(int)} and {@link
+   * #getInterpolatedSecondaryDerivatives(int)}. So if their content must be preserved
+   * across several calls, user must copy them.</p>
    * @param time time of the interpolated point
+   * @see #getInterpolatedState()
+   * @see #getInterpolatedDerivatives()
+   * @see #getInterpolatedSecondaryState(int)
+   * @see #getInterpolatedSecondaryDerivatives(int)
    */
   void setInterpolatedTime(double time);
 
@@ -85,9 +92,13 @@ public interface StepInterpolator extends Externalizable {
    * Get the state vector of the interpolated point.
    * <p>The returned vector is a reference to a reused array, so
    * it should not be modified and it should be copied if it needs
-   * to be preserved across several calls.</p>
+   * to be preserved across several calls to the associated
+   * {@link #setInterpolatedTime(double)} method.</p>
    * @return state vector at time {@link #getInterpolatedTime}
    * @see #getInterpolatedDerivatives()
+   * @see #getInterpolatedSecondaryState(int)
+   * @see #getInterpolatedSecondaryDerivatives(int)
+   * @see #setInterpolatedTime(double)
    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    */
   double[] getInterpolatedState() throws MaxCountExceededException;
@@ -96,9 +107,13 @@ public interface StepInterpolator extends Externalizable {
    * Get the derivatives of the state vector of the interpolated point.
    * <p>The returned vector is a reference to a reused array, so
    * it should not be modified and it should be copied if it needs
-   * to be preserved across several calls.</p>
+   * to be preserved across several calls to the associated
+   * {@link #setInterpolatedTime(double)} method.</p>
    * @return derivatives of the state vector at time {@link #getInterpolatedTime}
    * @see #getInterpolatedState()
+   * @see #getInterpolatedSecondaryState(int)
+   * @see #getInterpolatedSecondaryDerivatives(int)
+   * @see #setInterpolatedTime(double)
    * @since 2.0
    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    */
@@ -107,7 +122,8 @@ public interface StepInterpolator extends Externalizable {
   /** Get the interpolated secondary state corresponding to the secondary equations.
    * <p>The returned vector is a reference to a reused array, so
    * it should not be modified and it should be copied if it needs
-   * to be preserved across several calls.</p>
+   * to be preserved across several calls to the associated
+   * {@link #setInterpolatedTime(double)} method.</p>
    * @param index index of the secondary set, as returned by {@link
    * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations(
    * org.apache.commons.math3.ode.SecondaryEquations)

src/test/java/org/apache/commons/math3/ode/ContinuousOutputModelTest.java

@@ -63,22 +63,29 @@ public class ContinuousOutputModelTest {
                     pb.getFinalTime(), new double[pb.getDimension()]);
 
     Random random = new Random(347588535632l);
-    double maxError = 0.0;
+    double maxError    = 0.0;
+    double maxErrorDot = 0.0;
     for (int i = 0; i < 1000; ++i) {
       double r = random.nextDouble();
       double time = r * pb.getInitialTime() + (1.0 - r) * pb.getFinalTime();
       cm.setInterpolatedTime(time);
-      double[] interpolatedY = cm.getInterpolatedState ();
-      double[] theoreticalY  = pb.computeTheoreticalState(time);
+      double[] interpolatedY    = cm.getInterpolatedState();
+      double[] interpolatedYDot = cm.getInterpolatedDerivatives();
+      double[] theoreticalY     = pb.computeTheoreticalState(time);
+      double[] theoreticalYDot  = new double[pb.getDimension()];
+      pb.doComputeDerivatives(time, theoreticalY, theoreticalYDot);
       double dx = interpolatedY[0] - theoreticalY[0];
       double dy = interpolatedY[1] - theoreticalY[1];
       double error = dx * dx + dy * dy;
-      if (error > maxError) {
-        maxError = error;
-      }
+      maxError = FastMath.max(maxError, error);
+      double dxDot = interpolatedYDot[0] - theoreticalYDot[0];
+      double dyDot = interpolatedYDot[1] - theoreticalYDot[1];
+      double errorDot = dxDot * dxDot + dyDot * dyDot;
+      maxErrorDot = FastMath.max(maxErrorDot, errorDot);
     }
 
-    Assert.assertTrue(maxError < 1.0e-9);
+    Assert.assertEquals(0.0, maxError,    1.0e-9);
+    Assert.assertEquals(0.0, maxErrorDot, 4.0e-7);
 
   }