From 8717704d9a12b8c9754f47434c033696d3c7500c Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Sat, 18 Jul 2009 18:21:36 +0000 Subject: [PATCH] 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 --- .../commons/math/ode/MultistepIntegrator.java | 32 +++++++++++++++++-- .../AdamsBashforthIntegratorTest.java | 10 +++--- .../nonstiff/AdamsMoultonIntegratorTest.java | 6 ++-- 3 files changed, 38 insertions(+), 10 deletions(-) diff --git a/src/java/org/apache/commons/math/ode/MultistepIntegrator.java b/src/java/org/apache/commons/math/ode/MultistepIntegrator.java index 5629e16b4..8212e119f 100644 --- a/src/java/org/apache/commons/math/ode/MultistepIntegrator.java +++ b/src/java/org/apache/commons/math/ode/MultistepIntegrator.java @@ -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. + *

We define scaled derivatives si(n) at step n as: + * 

+ * s1(n) = h y'n for first derivative
+ * s2(n) = h2/2 y''n for second derivative
+ * s3(n) = h3/6 y'''n for third derivative
+ * ...
+ * sk(n) = hk/k! y(k)n for kth derivative
+ *

+ *

Rather than storing several previous steps separately, this implementation uses + * the Nordsieck vector with higher degrees scaled derivatives all taken at the same + * step (yn, s1(n) and rn) where rn is defined as: + * 

+ * rn = [ s2(n), s3(n) ... sk(n) ]T
+ *
+ * (we omit the k index in the notation for clarity)

+ *

+ * Multistep integrators with Nordsieck representation are highly sensitive to + * large step changes because when the step is multiplied by a factor a, the + * kth component of the Nordsieck vector is multiplied by ak + * and the last components are the least accurate ones. The default max growth + * factor is therefore set to a quite low value: 21/order. + *

* * @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 { *

The default starter integrator is set to the {@link * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with * some defaults settings.

+ *

+ * The default max growth factor is set to a quite low value: 21/order. + *

* @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 { *

The default starter integrator is set to the {@link * DormandPrince853Integrator Dormand-Prince 8(5,3)} integrator with * some defaults settings.

+ *

+ * The default max growth factor is set to a quite low value: 21/order. + *

* @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)); } diff --git a/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java b/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java index 6842b3a5d..8c4edac6b 100644 --- a/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java +++ b/src/test/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java @@ -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); } } diff --git a/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java b/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java index d0f6d27df..0355800ac 100644 --- a/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java +++ b/src/test/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegratorTest.java @@ -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); } }