Fixed a problem with getInterpolatedDerivatives returning zero derivatives when

an ODE step handler is configured to not use interpolation. It now returns a
constant but non-zero value consistent with at least one point inside the step


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@919479 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2010-03-05 16:35:56 +00:00
parent 1bb6d2ad6b
commit 81e8b4f0d1
9 changed files with 46 additions and 15 deletions

View File

@ -214,7 +214,7 @@ public abstract class EmbeddedRungeKuttaIntegrator
rki.reinitialize(this, yTmp, yDotK, forward); rki.reinitialize(this, yTmp, yDotK, forward);
interpolator = rki; interpolator = rki;
} else { } else {
interpolator = new DummyStepInterpolator(yTmp, forward); interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward);
} }
interpolator.storeTime(t0); interpolator.storeTime(t0);

View File

@ -624,7 +624,7 @@ public class GraggBulirschStoerIntegrator extends AdaptiveStepsizeIntegrator {
y1, yDot1, y1, yDot1,
yMidDots, forward); yMidDots, forward);
} else { } else {
interpolator = new DummyStepInterpolator(y, forward); interpolator = new DummyStepInterpolator(y, yDot1, forward);
} }
interpolator.storeTime(t0); interpolator.storeTime(t0);

View File

@ -122,7 +122,7 @@ class GraggBulirschStoerStepInterpolator
* @param y1 reference to the integrator array holding the state at * @param y1 reference to the integrator array holding the state at
* the end of the step * the end of the step
* @param y1Dot reference to the integrator array holding the slope * @param y1Dot reference to the integrator array holding the slope
* at theend of the step * at the end of the step
* @param yMidDots reference to the integrator array holding the * @param yMidDots reference to the integrator array holding the
* derivatives at the middle point of the step * derivatives at the middle point of the step
* @param forward integration direction indicator * @param forward integration direction indicator

View File

@ -120,7 +120,7 @@ public abstract class RungeKuttaIntegrator extends AbstractIntegrator {
rki.reinitialize(this, yTmp, yDotK, forward); rki.reinitialize(this, yTmp, yDotK, forward);
interpolator = rki; interpolator = rki;
} else { } else {
interpolator = new DummyStepInterpolator(yTmp, forward); interpolator = new DummyStepInterpolator(yTmp, yDotK[stages - 1], forward);
} }
interpolator.storeTime(t0); interpolator.storeTime(t0);

View File

@ -40,8 +40,11 @@ import org.apache.commons.math.ode.DerivativeException;
public class DummyStepInterpolator public class DummyStepInterpolator
extends AbstractStepInterpolator { extends AbstractStepInterpolator {
/** Serializable version identifier */ /** Serializable version identifier. */
private static final long serialVersionUID = 1708010296707839488L; private static final long serialVersionUID = 1708010296707839488L;
/** Current derivative. */
private double[] currentDerivative;
/** Simple constructor. /** Simple constructor.
* This constructor builds an instance that is not usable yet, the * This constructor builds an instance that is not usable yet, the
@ -55,15 +58,19 @@ public class DummyStepInterpolator
*/ */
public DummyStepInterpolator() { public DummyStepInterpolator() {
super(); super();
currentDerivative = null;
} }
/** Simple constructor. /** Simple constructor.
* @param y reference to the integrator array holding the state at * @param y reference to the integrator array holding the state at
* the end of the step * the end of the step
* @param yDot reference to the integrator array holding the state
* derivative at some arbitrary point within the step
* @param forward integration direction indicator * @param forward integration direction indicator
*/ */
public DummyStepInterpolator(final double[] y, final boolean forward) { public DummyStepInterpolator(final double[] y, final double[] yDot, final boolean forward) {
super(y, forward); super(y, forward);
currentDerivative = yDot;
} }
/** Copy constructor. /** Copy constructor.
@ -73,6 +80,7 @@ public class DummyStepInterpolator
*/ */
public DummyStepInterpolator(final DummyStepInterpolator interpolator) { public DummyStepInterpolator(final DummyStepInterpolator interpolator) {
super(interpolator); super(interpolator);
currentDerivative = interpolator.currentDerivative.clone();
} }
/** Really copy the finalized instance. /** Really copy the finalized instance.
@ -96,7 +104,8 @@ public class DummyStepInterpolator
@Override @Override
protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH)
throws DerivativeException { throws DerivativeException {
System.arraycopy(currentState, 0, interpolatedState, 0, currentState.length); System.arraycopy(currentState, 0, interpolatedState, 0, currentState.length);
System.arraycopy(currentDerivative, 0, interpolatedDerivatives, 0, currentDerivative.length);
} }
/** Write the instance to an output channel. /** Write the instance to an output channel.
@ -106,8 +115,16 @@ public class DummyStepInterpolator
@Override @Override
public void writeExternal(final ObjectOutput out) public void writeExternal(final ObjectOutput out)
throws IOException { throws IOException {
// save the state of the base class
// save the state of the base class
writeBaseExternal(out); writeBaseExternal(out);
if (currentDerivative != null) {
for (int i = 0; i < currentDerivative.length; ++i) {
out.writeDouble(currentDerivative[i]);
}
}
} }
/** Read the instance from an input channel. /** Read the instance from an input channel.
@ -121,6 +138,15 @@ public class DummyStepInterpolator
// read the base class // read the base class
final double t = readBaseExternal(in); final double t = readBaseExternal(in);
if (currentState == null) {
currentDerivative = null;
} else {
currentDerivative = new double[currentState.length];
for (int i = 0; i < currentDerivative.length; ++i) {
currentDerivative[i] = in.readDouble();
}
}
// we can now set the interpolated time and state // we can now set the interpolated time and state
setInterpolatedTime(t); setInterpolatedTime(t);

View File

@ -51,6 +51,11 @@ The <action> type attribute can be add,update,fix,remove.
used to solve Boundary Value Problems (BVP). There are no implementations (yet) used to solve Boundary Value Problems (BVP). There are no implementations (yet)
of BVP solvers in the library. of BVP solvers in the library.
</action> </action>
<action dev="luc" type="fix" >
Fixed a problem with getInterpolatedDerivatives returning zero derivatives when
an ODE step handler is configured to not use interpolation. It now returns a
constant but non-zero value consistent with at least one point inside the step
</action>
<action dev="luc" type="fix" issue="MATH-344" > <action dev="luc" type="fix" issue="MATH-344" >
Fixed wrong return values when enpoints are roots in Brent solver with Fixed wrong return values when enpoints are roots in Brent solver with
a user provided initial guess a user provided initial guess

View File

@ -166,7 +166,7 @@ public class ContinuousOutputModelTest
} }
private StepInterpolator buildInterpolator(double t0, double[] y0, double t1) { private StepInterpolator buildInterpolator(double t0, double[] y0, double t1) {
DummyStepInterpolator interpolator = new DummyStepInterpolator(y0, t1 >= t0); DummyStepInterpolator interpolator = new DummyStepInterpolator(y0, new double[y0.length], t1 >= t0);
interpolator.storeTime(t0); interpolator.storeTime(t0);
interpolator.shift(); interpolator.shift();
interpolator.storeTime(t1); interpolator.storeTime(t1);

View File

@ -52,7 +52,7 @@ public class EventStateTest {
double t0 = r1 - 0.5 * gap; double t0 = r1 - 0.5 * gap;
es.reinitializeBegin(t0, new double[0]); es.reinitializeBegin(t0, new double[0]);
AbstractStepInterpolator interpolator = AbstractStepInterpolator interpolator =
new DummyStepInterpolator(new double[0], true); new DummyStepInterpolator(new double[0], new double[0], true);
interpolator.storeTime(t0); interpolator.storeTime(t0);
interpolator.shift(); interpolator.shift();

View File

@ -38,7 +38,7 @@ public class DummyStepInterpolatorTest {
public void testNoReset() throws DerivativeException { public void testNoReset() throws DerivativeException {
double[] y = { 0.0, 1.0, -2.0 }; double[] y = { 0.0, 1.0, -2.0 };
DummyStepInterpolator interpolator = new DummyStepInterpolator(y, true); DummyStepInterpolator interpolator = new DummyStepInterpolator(y, new double[y.length], true);
interpolator.storeTime(0); interpolator.storeTime(0);
interpolator.shift(); interpolator.shift();
interpolator.storeTime(1); interpolator.storeTime(1);
@ -55,7 +55,7 @@ public class DummyStepInterpolatorTest {
throws DerivativeException { throws DerivativeException {
double[] y = { 1.0, 3.0, -4.0 }; double[] y = { 1.0, 3.0, -4.0 };
DummyStepInterpolator interpolator = new DummyStepInterpolator(y, true); DummyStepInterpolator interpolator = new DummyStepInterpolator(y, new double[y.length], true);
interpolator.storeTime(0); interpolator.storeTime(0);
interpolator.shift(); interpolator.shift();
interpolator.storeTime(1); interpolator.storeTime(1);
@ -79,7 +79,7 @@ public class DummyStepInterpolatorTest {
throws DerivativeException, IOException, ClassNotFoundException { throws DerivativeException, IOException, ClassNotFoundException {
double[] y = { 0.0, 1.0, -2.0 }; double[] y = { 0.0, 1.0, -2.0 };
DummyStepInterpolator interpolator = new DummyStepInterpolator(y, true); DummyStepInterpolator interpolator = new DummyStepInterpolator(y, new double[y.length], true);
interpolator.storeTime(0); interpolator.storeTime(0);
interpolator.shift(); interpolator.shift();
interpolator.storeTime(1); interpolator.storeTime(1);
@ -132,7 +132,7 @@ public class DummyStepInterpolatorTest {
public BadStepInterpolator() { public BadStepInterpolator() {
} }
public BadStepInterpolator(double[] y, boolean forward) { public BadStepInterpolator(double[] y, boolean forward) {
super(y, forward); super(y, new double[y.length], forward);
} }
@Override @Override
protected void doFinalize() protected void doFinalize()