Phil Steitz a écrit :
Luc,
I am about to commit the code. Actually was working on that when I
saw this. Once I have the code in, submit the patch against the code
in svn.
Here is the patch. I hope it is not too large to be submitted through
this list. If anybody is upset with this, please accept my apologies.
The commit message I suggest is the following one:
fixed a problem when switching functions triggered derivatives
discontinuities
Luc
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java (working copy)
@@ -99,22 +99,6 @@
}
/** Build a rotation from the quaternion coordinates.
- * @param q0 scalar part of the quaternion
- * @param q1 first coordinate of the vectorial part of the quaternion
- * @param q2 second coordinate of the vectorial part of the quaternion
- * @param q3 third coordinate of the vectorial part of the quaternion
- * @deprecated since Mantissa 6.3, this method as been deprecated as it
- * does not properly handles non-normalized quaternions, it should be
- * replaced by [EMAIL PROTECTED] #Rotation(double, double, double, double, boolean)}
- */
- public Rotation(double q0, double q1, double q2, double q3) {
- this.q0 = q0;
- this.q1 = q1;
- this.q2 = q2;
- this.q3 = q3;
- }
-
- /** Build a rotation from the quaternion coordinates.
* <p>A rotation can be built from a <em>normalized</em> quaternion,
* i.e. a quaternion for which q<sub>0</sub><sup>2</sup> +
* q<sub>1</sub><sup>2</sup> + q<sub>2</sub><sup>2</sup> +
@@ -120,9 +104,6 @@
* q<sub>1</sub><sup>2</sup> + q<sub>2</sub><sup>2</sup> +
* q<sub>3</sub><sup>2</sup> = 1. If the quaternion is not normalized,
* the constructor can normalize it in a preprocessing step.</p>
- * <p>This method replaces the [EMAIL PROTECTED] #Rotation(double, double,
- * double, double) constructor using only 4 doubles} which was deprecated
- * as of version 6.3 of Mantissa.</p>
* @param q0 scalar part of the quaternion
* @param q1 first coordinate of the vectorial part of the quaternion
* @param q2 second coordinate of the vectorial part of the quaternion
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java (working copy)
@@ -91,7 +91,7 @@
{ "dimensions mismatch: ODE problem has dimension {0},"
+ " state vector has dimension {1}",
"incompatibilit\u00e9 de dimensions entre le probl\u00e8me ODE ({0}),"
- + " et le vecteur d'\u00e9tat ({1})" },
+ + " et le vecteur d''\u00e9tat ({1})" },
{ "too small integration interval: length = {0}",
"intervalle d''int\u00e9gration trop petit : {0}" },
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java (working copy)
@@ -880,7 +880,11 @@
interpolator.storeTime(currentT);
handler.handleStep(interpolator, lastStep);
- switchesHandler.reset(currentT, y);
+ if (switchesHandler.reset(currentT, y) && ! lastStep) {
+ // some switching function has triggered changes that
+ // invalidate the derivatives, we need to recompute them
+ firstStepAlreadyComputed = false;
+ }
int optimalIter;
if (k == 1) {
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java (working copy)
@@ -311,7 +311,11 @@
System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
}
- switchesHandler.reset(currentT, y);
+ if (switchesHandler.reset(currentT, y) && ! lastStep) {
+ // some switching function has triggered changes that
+ // invalidate the derivatives, we need to recompute them
+ equations.computeDerivatives(currentT, y, yDotK[0]);
+ }
if (! lastStep) {
// stepsize control for next step
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java (working copy)
@@ -234,7 +234,11 @@
System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
}
- switchesHandler.reset(currentT, y);
+ if (switchesHandler.reset(currentT, y) && ! lastStep) {
+ // some switching function has triggered changes that
+ // invalidate the derivatives, we need to recompute them
+ equations.computeDerivatives(currentT, y, yDotK[0]);
+ }
if (needUpdate) {
// a switching function has changed the step
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java (working copy)
@@ -22,18 +22,20 @@
* <p>A switching function allows to handle discrete events in
* integration problems. These events occur for example when the
* integration process should be stopped as some value is reached
- * (G-stop facility), or when the derivatives has
- * discontinuities. These events are traditionally defined as
- * occurring when a <code>g</code> function sign changes.</p>
+ * (G-stop facility), or when the derivatives have
+ * discontinuities, or simply when the user wants to monitor some
+ * states boundaries crossings. These events are traditionally defined
+ * as occurring when a <code>g</code> function sign changes, hence
+ * the name <em>switching functions</em>.</p>
*
* <p>Since events are only problem-dependent and are triggered by the
* independant <i>time</i> variable and the state vector, they can
- * occur at virtually any time. The integrators will take care to
- * avoid sign changes inside the steps, they will reduce the step size
- * when such an event is detected in order to put this event exactly
- * at the end of the current step. This guarantees that step
- * interpolation (which always has a one step scope) is relevant even
- * in presence of discontinuities. This is independent from the
+ * occur at virtually any time, unknown in advance. The integrators will
+ * take care to avoid sign changes inside the steps, they will reduce
+ * the step size when such an event is detected in order to put this
+ * event exactly at the end of the current step. This guarantees that
+ * step interpolation (which always has a one step scope) is relevant
+ * even in presence of discontinuities. This is independent from the
* stepsize control provided by integrators that monitor the local
* error (this feature is available on all integrators, including
* fixed step ones).</p>
@@ -52,14 +54,23 @@
*/
public static final int STOP = 0;
- /** Reset indicator.
+ /** Reset state indicator.
* <p>This value should be used as the return value of the [EMAIL PROTECTED]
* #eventOccurred eventOccurred} method when the integration should
* go on after the event ending the current step, with a new state
- * vector (which will be retrieved through the [EMAIL PROTECTED] #resetState
+ * vector (which will be retrieved thanks to the [EMAIL PROTECTED] #resetState
* resetState} method).</p>
*/
- public static final int RESET = 1;
+ public static final int RESET_STATE = 1;
+
+ /** Reset derivatives indicator.
+ * <p>This value should be used as the return value of the [EMAIL PROTECTED]
+ * #eventOccurred eventOccurred} method when the integration should
+ * go on after the event ending the current step, with a new derivatives
+ * vector (which will be retrieved thanks to the [EMAIL PROTECTED]
+ * FirstOrderDifferentialEquations#computeDerivatives} method).</p>
+ */
+ public static final int RESET_DERIVATIVES = 2;
/** Continue indicator.
* <p>This value should be used as the return value of the [EMAIL PROTECTED]
@@ -66,7 +77,7 @@
* #eventOccurred eventOccurred} method when the integration should go
* on after the event ending the current step.</p>
*/
- public static final int CONTINUE = 2;
+ public static final int CONTINUE = 3;
/** Compute the value of the switching function.
@@ -73,8 +84,8 @@
* <p>Discrete events are generated when the sign of this function
* changes, the integrator will take care to change the stepsize in
* such a way these events occur exactly at step boundaries. This
- * function must be continuous, as the integrator will need to find
- * its roots to locate the events.</p>
+ * function must be continuous (at least in its roots neighborhood),
+ * as the integrator will need to find its roots to locate the events.</p>
* @param t current value of the independant <i>time</i> variable
* @param y array containing the current value of the state vector
@@ -88,16 +99,27 @@
* ending exactly on a sign change of the function, just before the
* step handler itself is called. It allows the user to update his
* internal data to acknowledge the fact the event has been handled
- * (for example setting a flag to switch the derivatives computation
- * in case of discontinuity), and it allows to direct the integrator
- * to either stop or continue integration, possibly with a reset
- * state.</p>
+ * (for example setting a flag in the [EMAIL PROTECTED]
+ * FirstOrderDifferentialEquations differential equations} to switch
+ * the derivatives computation in case of discontinuity), or to
+ * direct the integrator to either stop or continue integration,
+ * possibly with a reset state or derivatives.</p>
- * <p>If [EMAIL PROTECTED] #STOP} is returned, the step handler will be called
- * with the <code>isLast</code> flag of the [EMAIL PROTECTED]
- * StepHandler#handleStep handleStep} method set to true. If [EMAIL PROTECTED]
- * #RESET} is returned, the [EMAIL PROTECTED] #resetState resetState} method
- * will be called once the step handler has finished its task.</p>
+ * <ul>
+ * <li>if [EMAIL PROTECTED] #STOP} is returned, the step handler will be called
+ * with the <code>isLast</code> flag of the [EMAIL PROTECTED]
+ * StepHandler#handleStep handleStep} method set to true and the
+ * integration will be stopped,</li>
+ * <li>if [EMAIL PROTECTED] #RESET_STATE} is returned, the [EMAIL PROTECTED] #resetState
+ * resetState} method will be called once the step handler has
+ * finished its task, and the integrator will also recompute the
+ * derivatives,</li>
+ * <li>if [EMAIL PROTECTED] #RESET_DERIVATIVES} is returned, the integrator
+ * will recompute the derivatives,
+ * <li>if [EMAIL PROTECTED] #CONTINUE} is returned, no specific action will
+ * be taken (apart from having called this method) and integration
+ * will continue.</li>
+ * </ul>
* @param t current value of the independant <i>time</i> variable
* @param y array containing the current value of the state vector
@@ -102,8 +124,8 @@
* @param t current value of the independant <i>time</i> variable
* @param y array containing the current value of the state vector
* @return indication of what the integrator should do next, this
- * value must be one of [EMAIL PROTECTED] #STOP}, [EMAIL PROTECTED] #RESET} or [EMAIL PROTECTED]
- * #CONTINUE}
+ * value must be one of [EMAIL PROTECTED] #STOP}, [EMAIL PROTECTED] #RESET_STATE},
+ * [EMAIL PROTECTED] #RESET_DERIVATIVES} or [EMAIL PROTECTED] #CONTINUE}
*/
public int eventOccurred(double t, double[] y);
@@ -111,11 +133,11 @@
* <p>This method is called after the step handler has returned and
* before the next step is started, but only when [EMAIL PROTECTED]
- * #eventOccurred} has itself returned the [EMAIL PROTECTED] #RESET}
+ * #eventOccurred} has itself returned the [EMAIL PROTECTED] #RESET_STATE}
* indicator. It allows the user to reset the state vector for the
* next step, without perturbing the step handler of the finishing
* step. If the [EMAIL PROTECTED] #eventOccurred} never returns the [EMAIL PROTECTED]
- * #RESET} indicator, this function will never be called, and it is
+ * #RESET_STATE} indicator, this function will never be called, and it is
* safe to leave its body empty.</p>
* @param t current value of the independant <i>time</i> variable
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java (working copy)
@@ -166,11 +166,16 @@
* beginning of the next step
* @param y array were to put the desired state vector at the beginning
* of the next step
+ * @return true if the integrator should reset the derivatives too
*/
- public void reset(double t, double[] y) {
+ public boolean reset(double t, double[] y) {
+ boolean resetDerivatives = false;
for (Iterator iter = functions.iterator(); iter.hasNext();) {
- ((SwitchState) iter.next()).reset(t, y);
+ if (((SwitchState) iter.next()).reset(t, y)) {
+ resetDerivatives = true;
+ }
}
+ return resetDerivatives;
}
/** Switching functions. */
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java (working copy)
@@ -239,15 +239,23 @@
* beginning of the next step
* @param y array were to put the desired state vector at the beginning
* of the next step
+ * @return true if the integrator should reset the derivatives too
*/
- public void reset(double t, double[] y) {
- if (pendingEvent) {
- if (nextAction == SwitchingFunction.RESET) {
- function.resetState(t, y);
- }
- pendingEvent = false;
- pendingEventTime = Double.NaN;
+ public boolean reset(double t, double[] y) {
+
+ if (! pendingEvent) {
+ return false;
+ }
+
+ if (nextAction == SwitchingFunction.RESET_STATE) {
+ function.resetState(t, y);
}
+ pendingEvent = false;
+ pendingEventTime = Double.NaN;
+
+ return (nextAction == SwitchingFunction.RESET_STATE)
+ || (nextAction == SwitchingFunction.RESET_DERIVATIVES);
+
}
/** Get the value of the g function at the specified time.
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java (working copy)
@@ -305,6 +305,17 @@
}
+ public void testUnstableDerivative()
+ throws DerivativeException, IntegratorException {
+ final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+ FirstOrderIntegrator integ =
+ new DormandPrince853Integrator(0.1, 10, 1.0e-12, 0.0);
+ integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+ double[] y = { Double.NaN };
+ integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+ assertEquals(8.0, y[0], 1.0e-12);
+ }
+
public static Test suite() {
return new TestSuite(DormandPrince853IntegratorTest.class);
}
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java (working copy)
@@ -196,6 +196,16 @@
pb.getFinalTime(), new double[pb.getDimension()]);
}
+ public void testUnstableDerivative()
+ throws DerivativeException, IntegratorException {
+ final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+ FirstOrderIntegrator integ = new GillIntegrator(0.3);
+ integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+ double[] y = { Double.NaN };
+ integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+ assertEquals(8.0, y[0], 1.0e-12);
+ }
+
public static Test suite() {
return new TestSuite(GillIntegratorTest.class);
}
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java (working copy)
@@ -256,6 +256,17 @@
pb.getFinalTime(), new double[pb.getDimension()]);
}
+ public void testUnstableDerivative()
+ throws DerivativeException, IntegratorException {
+ final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+ FirstOrderIntegrator integ =
+ new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);
+ integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+ double[] y = { Double.NaN };
+ integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+ assertEquals(8.0, y[0], 1.0e-12);
+ }
+
public static Test suite() {
return new TestSuite(GraggBulirschStoerIntegratorTest.class);
}
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/StepProblem.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/StepProblem.java (revision 0)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/StepProblem.java (revision 0)
@@ -0,0 +1,42 @@
+package org.spaceroots.mantissa.ode;
+
+
+public class StepProblem
+ implements FirstOrderDifferentialEquations, SwitchingFunction {
+
+ public StepProblem(double rateBefore, double rateAfter,
+ double switchTime) {
+ this.rateAfter = rateAfter;
+ this.switchTime = switchTime;
+ setRate(rateBefore);
+ }
+
+ public void computeDerivatives(double t, double[] y, double[] yDot) {
+ yDot[0] = rate;
+ }
+
+ public int getDimension() {
+ return 1;
+ }
+
+ public void setRate(double rate) {
+ this.rate = rate;
+ }
+
+ public int eventOccurred(double t, double[] y) {
+ setRate(rateAfter);
+ return RESET_DERIVATIVES;
+ }
+
+ public double g(double t, double[] y) {
+ return t - switchTime;
+ }
+
+ public void resetState(double t, double[] y) {
+ }
+
+ private double rate;
+ private double rateAfter;
+ private double switchTime;
+
+}
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java (revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java (working copy)
@@ -104,7 +104,7 @@
public int eventOccurred(double t, double[] y) {
// this sign change is needed because the state will be reset soon
sign = -sign;
- return SwitchingFunction.RESET;
+ return SwitchingFunction.RESET_STATE;
}
public void resetState(double t, double[] y) {
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
