Details
Minor Description
org.apache.commons.math3.analysis.integration.SimpsonIntergrator
When used with minimalIterationCount == 1 the integrator computes the wrong value due to the inlining of computation of stage 1 and stage 0 of the TrapezoidIntegrator. Each stage is successive since it relies on the result of the previous stage. So stage 0 must be computed first. This inlining causes stage 1 to be computed before stage 0:
return (4 * qtrap.stage(this, 1) - qtrap.stage(this, 0)) / 3.0;
This can be fixed using:
final double s0 = qtrap.stage(this, 0); return (4 * qtrap.stage(this, 1) - s0) / 3.0;
What is not clear is why setting minimum iterations to 1 results in no iteration. This is not documented. I would expect setting it to 1 would compute the first Simpson sum and then perform 1 refinement. This would make it functionality equivalent to the other Integrator classes which compute two sums for the first iteration and allow them to be compared if minimum iterations = 1. If convergence fails then each additional iteration computes an additional sum.
Note when used with minimalIterationCount > 1 the SimpsonIntegrator wastes a stage since it computes the following stages: 0, 0, 1, 2, 3. i.e. stage 0 is computed twice. This is because the iteration is incremented after the stage is computed:
final double t = qtrap.stage(this, getIterations()); incrementCount();
This should be:
incrementCount(); final double t = qtrap.stage(this, getIterations());
On the first iteration it thus computes the trapezoid sum and compares it to zero. This would result in a bad computation if integrating a function whose trapezoid sum is zero (e.g. y=x^3 in the range -1 to 1). However since iteration only occurs for minimalIterationCount>1 no termination comparison is made on the first loop. The first termination comparison can be made at iteration=2 where the comparison will be between the Trapezoid sum and the first Simpson sum. This is a bug.
However I do not want to submit a formal patch as there is a lack of consistency across all the integrators in their doIntegrate() method with the use of incrementCount() and what the iteration number should be at the start of the while (true) loop:
- IterativeLegendreGauss integrator uses getIterations()+1 to mark the current iteration inside the loop and calls incrementCount() at the end.
- TrapezoidIntegrator calls incrementCount() outside the loop, uses getIterations() to mark the current iteration and calls incrementCount() at the end.
- The MidpointIntegrator calls incrementCount() at the start of the loop and uses getIterations() to mark the current iteration.
- The RombergIntegrator calls incrementCount() outside the loop, uses getIterations() to mark the current iteration and calls incrementCount() in the middle of the loop before termination conditions have been checked. This allows it to fail when the iterations are equal to the maximum iterations even if convergence has been achieved (see Note*).
- The SimpsonIntegrator uses getIterations() to mark the current iteration and calls incrementCount() immediately after.
Note*: This may not be discovered in a unit test since the incrementCount() uses a backing Incrementor where the Incrementor.increment() method calls Incrementor.increment(1) which ends up calling canIncrement(0) { instead of canIncrement(nTimes) } to check if the maxCountCallback should be triggered. I expect that all uses of the Incrementor actually trigger the maxCountCallback when the count has actually exceeded the maximalCount. I don't want to get into debugging that class since it also has an iterator using hasNext() with a call to canIncrement(0) and I do not know the contract that the iterator is working under.
A consistent approach would be:
- Compute the first sum before the loop
- Enter the loop and increment the iteration (so the first loop execution would be iteration 1)
- Compute the next sum
- Check termination conditions
An example for the SimpsonIntegrator is below:
protected double doIntegrate() throws TooManyEvaluationsException, MaxCountExceededException { // This is a modification from the default SimpsonIntegrator. // That only computed a single iteration if // getMinimalIterationCount() == 1. // Simpson's rule requires at least two trapezoid stages. // 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 now the first refinement of the sum. // This matches how the MidPointIntegrator works. incrementCount(); final int i = getIterations(); // 1-stage ahead of the iteration final double t = qtrap.stage(this, i + 1); final double s = (4 * t - oldt) / 3.0; 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())) { return s; } } olds = s; oldt = t; } }
Note: If this method is accepted then the SIMPSON_MAX_ITERATIONS_COUNT must be reduced by 1 to 63, since the stage method of the TrapezoidIntegrator has a maximum valid input argument of 64.
