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

MATH-1458: Fixed iteration of the SimpsonIntegrator

This commit is contained in:
aherbert 2018-05-08 10:47:00 +01:00
parent 3cce9ed6c3
commit d97cc96e28
1 changed files with 18 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
src/main/java/org/apache/commons/math4/analysis/integration/SimpsonIntegrator.java
View File

@ -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;
}
}
}