You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2018/05/08 15:51:01 UTC

[3/8] [math] MATH-1458: Fixed iteration of the SimpsonIntegrator

MATH-1458: Fixed iteration of the SimpsonIntegrator

Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/d97cc96e
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/d97cc96e
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/d97cc96e

Branch: refs/heads/master
Commit: d97cc96e285b6cf382d96f7213e5a62355bc0aa2
Parents: 3cce9ed
Author: aherbert <a....@sussex.ac.uk>
Authored: Tue May 8 10:47:00 2018 +0100
Committer: aherbert <a....@sussex.ac.uk>
Committed: Tue May 8 10:47:00 2018 +0100

----------------------------------------------------------------------
 .../analysis/integration/SimpsonIntegrator.java | 31 +++++++++++---------
 1 file changed, 17 insertions(+), 14 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/d97cc96e/src/main/java/org/apache/commons/math4/analysis/integration/SimpsonIntegrator.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/analysis/integration/SimpsonIntegrator.java b/src/main/java/org/apache/commons/math4/analysis/integration/SimpsonIntegrator.java
index 5cfea07..24325a0 100644
--- a/src/main/java/org/apache/commons/math4/analysis/integration/SimpsonIntegrator.java
+++ b/src/main/java/org/apache/commons/math4/analysis/integration/SimpsonIntegrator.java
@@ -37,7 +37,7 @@ import org.apache.commons.math4.util.FastMath;
 public class SimpsonIntegrator extends BaseAbstractUnivariateIntegrator {
 
     /** Maximal number of iterations for Simpson. */
-    public static final int SIMPSON_MAX_ITERATIONS_COUNT = 64;
+    public static final int SIMPSON_MAX_ITERATIONS_COUNT = 63;
 
     /**
      * Build a Simpson integrator with given accuracies and iterations counts.
@@ -100,30 +100,33 @@ public class SimpsonIntegrator extends BaseAbstractUnivariateIntegrator {
     protected double doIntegrate()
         throws TooManyEvaluationsException, MaxCountExceededException {
 
+        // Simpson's rule requires at least two trapezoid stages.
+        // So we set the first sum using two trapezoid stages.
         TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
-        if (getMinimalIterationCount() == 1) {
-            return (4 * qtrap.stage(this, 1) - qtrap.stage(this, 0)) / 3.0;
-        }
 
-        // Simpson's rule requires at least two trapezoid stages.
-        double olds = 0;
-        double oldt = qtrap.stage(this, 0);
-        while (true) {
-            final double t = qtrap.stage(this, iterations.getCount());
+        final double s0 = qtrap.stage(this, 0);
+        double oldt = qtrap.stage(this, 1);
+        double olds = (4 * oldt - s0) / 3.0;
+        while (true)
+        {
+            // The first iteration is the first refinement of the sum.
             iterations.incrementCount();
+            final int i = getIterations();
+            final double t = qtrap.stage(this, i + 1); // 1-stage ahead of the iteration
             final double s = (4 * t - oldt) / 3.0;
-            if (iterations.getCount() >= getMinimalIterationCount()) {
+            if (i >= getMinimalIterationCount())
+            {
                 final double delta = FastMath.abs(s - olds);
-                final double rLimit =
-                    getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
-                if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
+                final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
+                if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy()))
+                {
                     return s;
                 }
             }
             olds = s;
             oldt = t;
         }
-
+        
     }
 
 }