Author: luc
Date: Wed Mar 24 22:11:51 2010
New Revision: 927202
URL: http://svn.apache.org/viewvc?rev=927202&view=rev
Log:
Fixed an error in events handling in ODE solvers. In some rare cases, events
occurring close to a step start were handled without truncating the step,
making them appear as is they occurred close to the step end
JIRA: MATH-358
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.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/DormandPrince853IntegratorTest.java
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=927202&r1=927201&r2=927202&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
Wed Mar 24 22:11:51 2010
@@ -271,8 +271,16 @@ public class AdamsBashforthIntegrator ex
if (manager.evaluateStep(interpolatorTmp)) {
final double dt = manager.getEventTime() - stepStart;
if (Math.abs(dt) <= Math.ulp(stepStart)) {
- // rejecting the step would lead to a too small
next step, we accept it
- loop = false;
+ // we cannot simply truncate the step, reject the
current computation
+ // and let the loop compute another state with the
truncated step.
+ // it is so small (much probably exactly 0 due to
limited accuracy)
+ // that the code above would fail handling it.
+ // So we set up an artificial 0 size step by
copying states
+ interpolator.storeTime(stepStart);
+ System.arraycopy(y, 0, yTmp, 0, y0.length);
+ hNew = 0;
+ stepSize = 0;
+ loop = false;
} else {
// reject the step to match exactly the next
switch time
hNew = dt;
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=927202&r1=927201&r2=927202&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
Wed Mar 24 22:11:51 2010
@@ -289,8 +289,16 @@ public class AdamsMoultonIntegrator exte
if (manager.evaluateStep(interpolatorTmp)) {
final double dt = manager.getEventTime() - stepStart;
if (Math.abs(dt) <= Math.ulp(stepStart)) {
- // rejecting the step would lead to a too small
next step, we accept it
- loop = false;
+ // we cannot simply truncate the step, reject the
current computation
+ // and let the loop compute another state with the
truncated step.
+ // it is so small (much probably exactly 0 due to
limited accuracy)
+ // that the code above would fail handling it.
+ // So we set up an artificial 0 size step by
copying states
+ interpolator.storeTime(stepStart);
+ System.arraycopy(y, 0, yTmp, 0, y0.length);
+ hNew = 0;
+ stepSize = 0;
+ loop = false;
} else {
// reject the step to match exactly the next
switch time
hNew = dt;
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java?rev=927202&r1=927201&r2=927202&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
Wed Mar 24 22:11:51 2010
@@ -292,8 +292,16 @@ public abstract class EmbeddedRungeKutta
if (manager.evaluateStep(interpolator)) {
final double dt = manager.getEventTime() - stepStart;
if (Math.abs(dt) <= Math.ulp(stepStart)) {
- // rejecting the step would lead to a too small next step,
we accept it
- loop = false;
+ // we cannot simply truncate the step, reject the current
computation
+ // and let the loop compute another state with the truncated
step.
+ // it is so small (much probably exactly 0 due to limited
accuracy)
+ // that the code above would fail handling it.
+ // So we set up an artificial 0 size step by copying states
+ interpolator.storeTime(stepStart);
+ System.arraycopy(y, 0, yTmp, 0, y0.length);
+ hNew = 0;
+ stepSize = 0;
+ loop = false;
} else {
// reject the step to match exactly the next switch time
hNew = dt;
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java?rev=927202&r1=927201&r2=927202&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
Wed Mar 24 22:11:51 2010
@@ -172,8 +172,15 @@ public abstract class RungeKuttaIntegrat
if (manager.evaluateStep(interpolator)) {
final double dt = manager.getEventTime() - stepStart;
if (Math.abs(dt) <= Math.ulp(stepStart)) {
- // rejecting the step would lead to a too small next step, we
accept it
- loop = false;
+ // we cannot simply truncate the step, reject the current
computation
+ // and let the loop compute another state with the truncated
step.
+ // it is so small (much probably exactly 0 due to limited
accuracy)
+ // that the code above would fail handling it.
+ // So we set up an artificial 0 size step by copying states
+ interpolator.storeTime(stepStart);
+ System.arraycopy(y, 0, yTmp, 0, y0.length);
+ stepSize = 0;
+ loop = false;
} else {
// reject the step to match exactly the next switch time
stepSize = dt;
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=927202&r1=927201&r2=927202&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Wed Mar 24 22:11:51 2010
@@ -90,6 +90,11 @@ The <action> type attribute can be add,u
used to solve Boundary Value Problems (BVP). There are no
implementations (yet)
of BVP solvers in the library.
</action>
+ <action dev="luc" type="fix" issue="MATH-358" >
+ Fixed an error in events handling in ODE solvers. In some rare cases,
events
+ occurring close to a step start were handled without truncating the
step, making
+ them appear as is they occurred close to the step end
+ </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
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=927202&r1=927201&r2=927202&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
Wed Mar 24 22:11:51 2010
@@ -41,6 +41,59 @@ public class ClassicalRungeKuttaIntegrat
super(name);
}
+ public void testMissedEndEvent() throws IntegratorException,
DerivativeException {
+ final double t0 = 1878250320.0000029;
+ final double tEvent = 1878250379.9999986;
+ final double[] k = { 1.0e-4, 1.0e-5, 1.0e-6 };
+ FirstOrderDifferentialEquations ode = new
FirstOrderDifferentialEquations() {
+
+ public int getDimension() {
+ return k.length;
+ }
+
+ public void computeDerivatives(double t, double[] y, double[] yDot) {
+ for (int i = 0; i < y.length; ++i) {
+ yDot[i] = k[i] * y[i];
+ }
+ }
+ };
+
+ ClassicalRungeKuttaIntegrator integrator = new
ClassicalRungeKuttaIntegrator(60.0);
+
+ double[] y0 = new double[k.length];
+ for (int i = 0; i < y0.length; ++i) {
+ y0[i] = i + 1;
+ }
+ double[] y = new double[k.length];
+
+ double finalT = integrator.integrate(ode, t0, y0, tEvent, y);
+ Assert.assertEquals(tEvent, finalT, 5.0e-6);
+ for (int i = 0; i < y.length; ++i) {
+ Assert.assertEquals(y0[i] * Math.exp(k[i] * (finalT - t0)), y[i],
1.0e-9);
+ }
+
+ integrator.addEventHandler(new EventHandler() {
+
+ public void resetState(double t, double[] y) {
+ }
+
+ public double g(double t, double[] y) {
+ return t - tEvent;
+ }
+
+ public int eventOccurred(double t, double[] y, boolean increasing) {
+ Assert.assertEquals(tEvent, t, 5.0e-6);
+ return CONTINUE;
+ }
+ }, Double.POSITIVE_INFINITY, 1.0e-20, 100);
+ finalT = integrator.integrate(ode, t0, y0, tEvent + 120, y);
+ Assert.assertEquals(tEvent + 120, finalT, 5.0e-6);
+ for (int i = 0; i < y.length; ++i) {
+ Assert.assertEquals(y0[i] * Math.exp(k[i] * (finalT - t0)), y[i],
1.0e-9);
+ }
+
+ }
+
public void testSanityChecks() {
try {
TestProblem1 pb = new TestProblem1();
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java?rev=927202&r1=927201&r2=927202&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
Wed Mar 24 22:11:51 2010
@@ -18,6 +18,7 @@
package org.apache.commons.math.ode.nonstiff;
import org.apache.commons.math.ode.DerivativeException;
+import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
import org.apache.commons.math.ode.TestProblem1;
@@ -40,6 +41,62 @@ public class DormandPrince853IntegratorT
super(name);
}
+ public void testMissedEndEvent() throws IntegratorException,
DerivativeException {
+ final double t0 = 1878250320.0000029;
+ final double tEvent = 1878250379.9999986;
+ final double[] k = { 1.0e-4, 1.0e-5, 1.0e-6 };
+ FirstOrderDifferentialEquations ode = new
FirstOrderDifferentialEquations() {
+
+ public int getDimension() {
+ return k.length;
+ }
+
+ public void computeDerivatives(double t, double[] y, double[] yDot) {
+ for (int i = 0; i < y.length; ++i) {
+ yDot[i] = k[i] * y[i];
+ }
+ }
+ };
+
+ DormandPrince853Integrator integrator = new
DormandPrince853Integrator(0.0, 100.0,
+
1.0e-10, 1.0e-10);
+
+ double[] y0 = new double[k.length];
+ for (int i = 0; i < y0.length; ++i) {
+ y0[i] = i + 1;
+ }
+ double[] y = new double[k.length];
+
+ integrator.setInitialStepSize(60.0);
+ double finalT = integrator.integrate(ode, t0, y0, tEvent, y);
+ Assert.assertEquals(tEvent, finalT, 5.0e-6);
+ for (int i = 0; i < y.length; ++i) {
+ Assert.assertEquals(y0[i] * Math.exp(k[i] * (finalT - t0)), y[i],
1.0e-9);
+ }
+
+ integrator.setInitialStepSize(60.0);
+ integrator.addEventHandler(new EventHandler() {
+
+ public void resetState(double t, double[] y) {
+ }
+
+ public double g(double t, double[] y) {
+ return t - tEvent;
+ }
+
+ public int eventOccurred(double t, double[] y, boolean increasing) {
+ Assert.assertEquals(tEvent, t, 5.0e-6);
+ return CONTINUE;
+ }
+ }, Double.POSITIVE_INFINITY, 1.0e-20, 100);
+ finalT = integrator.integrate(ode, t0, y0, tEvent + 120, y);
+ Assert.assertEquals(tEvent + 120, finalT, 5.0e-6);
+ for (int i = 0; i < y.length; ++i) {
+ Assert.assertEquals(y0[i] * Math.exp(k[i] * (finalT - t0)), y[i],
1.0e-9);
+ }
+
+ }
+
public void testDimensionCheck() {
try {
TestProblem1 pb = new TestProblem1();
