From d97cc96e285b6cf382d96f7213e5a62355bc0aa2 Mon Sep 17 00:00:00 2001 From: aherbert Date: Tue, 8 May 2018 10:47:00 +0100 Subject: [PATCH] MATH-1458: Fixed iteration of the SimpsonIntegrator --- .../integration/SimpsonIntegrator.java | 33 ++++++++++--------- 1 file changed, 18 insertions(+), 15 deletions(-) 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 5cfea0701..24325a097 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 { - 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()); + // So we set the first sum using two trapezoid stages. + TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); + + 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; } - + } }