From 25aa4bd3665d8b265f03fa2b3e7ab6ee68256367 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Wed, 22 Oct 2014 17:34:29 +0200 Subject: [PATCH 1/3] Provide access to state derivatives in ContinuousOutputModel. JIRA: MATH-1160 --- src/changes/changes.xml | 3 +++ .../math3/ode/ContinuousOutputModel.java | 14 +++++++++++++ .../math3/ode/ContinuousOutputModelTest.java | 21 ++++++++++++------- 3 files changed, 31 insertions(+), 7 deletions(-) diff --git a/src/changes/changes.xml b/src/changes/changes.xml index bb5c5257e..9350d6b1a 100644 --- a/src/changes/changes.xml +++ b/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). "> + + Provide access to state derivatives in ContinuousOutputModel. + Fixed bicubic spline interpolator, using Akima splines. diff --git a/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java b/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java index 599eab900..d9f619261 100644 --- a/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java +++ b/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java @@ -332,12 +332,25 @@ public class ContinuousOutputModel * Get the state vector of the interpolated point. * @return state vector at time {@link #getInterpolatedTime} * @exception MaxCountExceededException if the number of functions evaluations is exceeded + * @see #getInterpolatedDerivatives() * @see #getInterpolatedSecondaryState(int) */ public double[] getInterpolatedState() throws MaxCountExceededException { return steps.get(index).getInterpolatedState(); } + /** + * Get the derivatives of the state vector of the interpolated point. + * @return derivatives of the state vector at time {@link #getInterpolatedTime} + * @exception MaxCountExceededException if the number of functions evaluations is exceeded + * @see #getInterpolatedState() + * @see #getInterpolatedSecondaryState(int) + * @since 3.4 + */ + public double[] getInterpolatedDerivatives() throws MaxCountExceededException { + return steps.get(index).getInterpolatedDerivatives(); + } + /** Get the interpolated secondary state corresponding to the secondary equations. * @param secondaryStateIndex index of the secondary set, as returned by {@link * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations( @@ -345,6 +358,7 @@ public class ContinuousOutputModel * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)} * @return interpolated secondary state at the current interpolation date * @see #getInterpolatedState() + * @see #getInterpolatedDerivatives() * @since 3.2 * @exception MaxCountExceededException if the number of functions evaluations is exceeded */ diff --git a/src/test/java/org/apache/commons/math3/ode/ContinuousOutputModelTest.java b/src/test/java/org/apache/commons/math3/ode/ContinuousOutputModelTest.java index 2f4053d93..3a098baf3 100644 --- a/src/test/java/org/apache/commons/math3/ode/ContinuousOutputModelTest.java +++ b/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); } From f8e6bc8ec72fd02cde8c705ca17609c76a07d3b8 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Wed, 22 Oct 2014 21:36:40 +0200 Subject: [PATCH 2/3] Improved javadoc for ContinuousOutputModel and StepInterpolator. The fact the arrays are reused is stated more clearly, as weel as the fact changing the time changes the array content. Yes, it is a bad design, we should change it in the next major version to allocate fresh new arrays :-( --- .../math3/ode/ContinuousOutputModel.java | 32 ++++++++++++++++--- .../math3/ode/sampling/StepInterpolator.java | 28 ++++++++++++---- 2 files changed, 49 insertions(+), 11 deletions(-) diff --git a/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java b/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java index d9f619261..6d324d16e 100644 --- a/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java +++ b/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.

*

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.

+ * 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.

+ *

Note that each time this method is called, the internal arrays + * returned in {@link #getInterpolatedState()}, {@link + * #getInterpolatedDerivatives()} and {@link #getInterpolatedSecondaryState(int)} + * will be overwritten. So if their content must be preserved + * across several calls, user must copy them.

* @param time time of the interpolated point + * @see #getInterpolatedState() + * @see #getInterpolatedDerivatives() + * @see #getInterpolatedSecondaryState(int) */ public void setInterpolatedTime(final double time) { @@ -330,8 +337,13 @@ public class ContinuousOutputModel /** * Get the state vector of the interpolated point. + *

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.

* @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) */ @@ -341,8 +353,13 @@ public class ContinuousOutputModel /** * Get the derivatives of the state vector of the interpolated point. + *

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.

* @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) * @since 3.4 @@ -352,11 +369,16 @@ public class ContinuousOutputModel } /** Get the interpolated secondary state corresponding to the secondary equations. + *

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.

* @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() * @since 3.2 diff --git a/src/main/java/org/apache/commons/math3/ode/sampling/StepInterpolator.java b/src/main/java/org/apache/commons/math3/ode/sampling/StepInterpolator.java index 44e77ebeb..5d27bf2c3 100644 --- a/src/main/java/org/apache/commons/math3/ode/sampling/StepInterpolator.java +++ b/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.

- *

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()}.

+ *

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.

* @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. *

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 be preserved across several calls to the associated + * {@link #setInterpolatedTime(double)} method.

* @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. *

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 be preserved across several calls to the associated + * {@link #setInterpolatedTime(double)} method.

* @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. *

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 be preserved across several calls to the associated + * {@link #setInterpolatedTime(double)} method.

* @param index index of the secondary set, as returned by {@link * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations( * org.apache.commons.math3.ode.SecondaryEquations) From 45ae5c7e42b07ff24e2385de57f70e9484a46ae3 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Wed, 22 Oct 2014 21:40:32 +0200 Subject: [PATCH 3/3] Added getInterpolatedSecondaryDerivatives to ContinuousOutputModel. This method is a close relative to getInterpolatedDerivatives, but is associated with the secondary state. JIRA: MATH-1160 --- .../math3/ode/ContinuousOutputModel.java | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java b/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java index 6d324d16e..75e497977 100644 --- a/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java +++ b/src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java @@ -346,6 +346,7 @@ public class ContinuousOutputModel * @see #setInterpolatedTime(double) * @see #getInterpolatedDerivatives() * @see #getInterpolatedSecondaryState(int) + * @see #getInterpolatedSecondaryDerivatives(int) */ public double[] getInterpolatedState() throws MaxCountExceededException { return steps.get(index).getInterpolatedState(); @@ -362,6 +363,7 @@ public class ContinuousOutputModel * @see #setInterpolatedTime(double) * @see #getInterpolatedState() * @see #getInterpolatedSecondaryState(int) + * @see #getInterpolatedSecondaryDerivatives(int) * @since 3.4 */ public double[] getInterpolatedDerivatives() throws MaxCountExceededException { @@ -381,6 +383,7 @@ public class ContinuousOutputModel * @see #setInterpolatedTime(double) * @see #getInterpolatedState() * @see #getInterpolatedDerivatives() + * @see #getInterpolatedSecondaryDerivatives(int) * @since 3.2 * @exception MaxCountExceededException if the number of functions evaluations is exceeded */ @@ -389,6 +392,28 @@ public class ContinuousOutputModel return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex); } + /** Get the interpolated secondary derivatives corresponding to the secondary equations. + *

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.

+ * @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