Author: luc
Date: Sat Jul 18 18:21:36 2009
New Revision: 795407
URL: http://svn.apache.org/viewvc?rev=795407&view=rev
Log:
Changed the default max growth factor for multistep methods using Nordsieck
representation.
The previous value (10.0) was far too big and lead to numerical instability at
high orders
because the last component of the Nordsieck vector, which has a low accuracy,
could be
multiplied by 10^k which was ... huge.
These integrators are at least usable now!
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
Modified:
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java?rev=795407&r1=795406&r2=795407&view=diff
==============================================================================
---
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java
(original)
+++
commons/proper/math/trunk/src/java/org/apache/commons/math/ode/MultistepIntegrator.java
Sat Jul 18 18:21:36 2009
@@ -28,6 +28,28 @@
/**
* This class is the base class for multistep integrators for Ordinary
* Differential Equations.
+ * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
+ * <pre>
+ * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
+ * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
+ * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
+ * ...
+ * s<sub>k</sub>(n) = h<sup>k</sup>/k! y(k)<sub>n</sub> for k<sup>th</sup>
derivative
+ * </pre></p>
+ * <p>Rather than storing several previous steps separately, this
implementation uses
+ * the Nordsieck vector with higher degrees scaled derivatives all taken at
the same
+ * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where
r<sub>n</sub> is defined as:
+ * <pre>
+ * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n)
]<sup>T</sup>
+ * </pre>
+ * (we omit the k index in the notation for clarity)</p>
+ * <p>
+ * Multistep integrators with Nordsieck representation are highly sensitive to
+ * large step changes because when the step is multiplied by a factor a, the
+ * k<sup>th</sup> component of the Nordsieck vector is multiplied by
a<sup>k</sup>
+ * and the last components are the least accurate ones. The default max growth
+ * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
+ * </p>
*
* @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
* @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
@@ -67,6 +89,9 @@
* <p>The default starter integrator is set to the {...@link
* DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
* some defaults settings.</p>
+ * <p>
+ * The default max growth factor is set to a quite low value:
2<sup>1/order</sup>.
+ * </p>
* @param name name of the method
* @param nSteps number of steps of the multistep method
* (excluding the one being computed)
@@ -102,7 +127,7 @@
// set the default values of the algorithm control parameters
setSafety(0.9);
setMinReduction(0.2);
- setMaxGrowth(10.0);
+ setMaxGrowth(Math.pow(2.0, -exp));
}
@@ -111,6 +136,9 @@
* <p>The default starter integrator is set to the {...@link
* DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with
* some defaults settings.</p>
+ * <p>
+ * The default max growth factor is set to a quite low value:
2<sup>1/order</sup>.
+ * </p>
* @param name name of the method
* @param nSteps number of steps of the multistep method
* (excluding the one being computed)
@@ -138,7 +166,7 @@
// set the default values of the algorithm control parameters
setSafety(0.9);
setMinReduction(0.2);
- setMaxGrowth(10.0);
+ setMaxGrowth(Math.pow(2.0, -exp));
}
Modified:
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java?rev=795407&r1=795406&r2=795407&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
(original)
+++
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
Sat Jul 18 18:21:36 2009
@@ -66,7 +66,7 @@
throws DerivativeException, IntegratorException {
int previousCalls = Integer.MAX_VALUE;
- for (int i = -12; i < -2; ++i) {
+ for (int i = -12; i < -5; ++i) {
TestProblem1 pb = new TestProblem1();
double minStep = 0;
double maxStep = pb.getFinalTime() - pb.getInitialTime();
@@ -82,11 +82,11 @@
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- // the 33 and 45 factors are only valid for this test
+ // the 31 and 36 factors are only valid for this test
// and has been obtained from trial and error
// there is no general relation between local and global errors
- assertTrue(handler.getMaximalValueError() > (33.0 *
scalAbsoluteTolerance));
- assertTrue(handler.getMaximalValueError() < (45.0 *
scalAbsoluteTolerance));
+ assertTrue(handler.getMaximalValueError() > (31.0 *
scalAbsoluteTolerance));
+ assertTrue(handler.getMaximalValueError() < (36.0 *
scalAbsoluteTolerance));
assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
int calls = pb.getCalls();
@@ -147,7 +147,7 @@
if (nSteps < 4) {
assertTrue(integ.getEvaluations() > 160);
} else {
- assertTrue(integ.getEvaluations() < 70);
+ assertTrue(integ.getEvaluations() < 80);
}
}
Modified:
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java?rev=795407&r1=795406&r2=795407&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
(original)
+++
commons/proper/math/trunk/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
Sat Jul 18 18:21:36 2009
@@ -82,10 +82,10 @@
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
- // the 0.4 and 3.0 factors are only valid for this test
+ // the 0.15 and 3.0 factors are only valid for this test
// and has been obtained from trial and error
// there is no general relation between local and global errors
- assertTrue(handler.getMaximalValueError() > (0.4 *
scalAbsoluteTolerance));
+ assertTrue(handler.getMaximalValueError() > (0.15 *
scalAbsoluteTolerance));
assertTrue(handler.getMaximalValueError() < (3.0 *
scalAbsoluteTolerance));
assertEquals(0, handler.getMaximalTimeError(), 1.0e-16);
@@ -147,7 +147,7 @@
if (nSteps < 4) {
assertTrue(integ.getEvaluations() > 150);
} else {
- assertTrue(integ.getEvaluations() < 90);
+ assertTrue(integ.getEvaluations() < 100);
}
}
