Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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!

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@795407 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-07-18 18:21:36 +00:00
parent 56a4d632c5
commit 8717704d9a
3 changed files with 38 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
src/java/org/apache/commons/math/ode/MultistepIntegrator.java
View File

@ -28,6 +28,28 @@ import org.apache.commons.math.ode.sampling.StepInterpolator;
/**
* 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 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
* <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 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
// 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 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
* <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 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
// set the default values of the algorithm control parameters
setSafety(0.9);
setMinReduction(0.2);
setMaxGrowth(10.0);
setMaxGrowth(Math.pow(2.0, -exp));
}

10
src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
View File

@ -66,7 +66,7 @@ public class AdamsBashforthIntegratorTest {
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 @@ public class AdamsBashforthIntegratorTest {
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 @@ public class AdamsBashforthIntegratorTest {
if (nSteps < 4) {
assertTrue(integ.getEvaluations() > 160);
} else {
assertTrue(integ.getEvaluations() < 70);
assertTrue(integ.getEvaluations() < 80);
}
}

6
src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java
View File

@ -82,10 +82,10 @@ public class AdamsMoultonIntegratorTest {
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 @@ public class AdamsMoultonIntegratorTest {
if (nSteps < 4) {
assertTrue(integ.getEvaluations() > 150);
} else {
assertTrue(integ.getEvaluations() < 90);
assertTrue(integ.getEvaluations() < 100);
}
}