Details

Bug

Status: Open

Major

Resolution: Unresolved

3.6.1

None

None
Description
Hi all,
I've been reading through the implementation of MidpointIntegrator, and I've discovered an issue with the algorithm as implemented.
The midpoint method generates an estimate for the definite integral of f between a and b by
 subdividing (a, b) into n intervals
 estimating the area of each interval as a rectangle (ba)/n wide and as tall as the midpoint of the interval  ie, ((ba)/n) * f((a + b) / 2)
 adding up all estimates.
The MidpointIntegrator implementation here sticks to n = powers of 2 for this reason, stated in the comments:
The interval is divided equally into 2^n sections rather than an arbitrary m sections because this configuration can best utilize the already computed values.
Here is the problem: the midpoint method can't reuse values if you split an interval in half. It can only reuse previous values if you split the interval into 3; ie, use 3^n sections, not 2^n.
What's actually implemented in `MidpointIntegrator` is a left Riemann sum without the leftmost point. Or, identically, a right Riemann sum without the rightmost point: https://en.wikipedia.org/wiki/Riemann_sum
This matters because the error of a (left, right) Riemann sum is proportional to 1/n, while the error of the midpoint method is proportional to 1/n^2... a big difference.
Explanation
The idea behind the midpoint method's point reuse is that if you triple the number of integral slices, one of the midpoints with n slices will overlap with the n/3 slice evaluation:
n = 1 x n = 3 xxx–
You can incrementally update the n=1 estimate by
 sampling the points at (1/6) and (5/6) of the way across the n=1 interval
 adding them to the n=1 estimate
 scaling this sum down by 3, to scale down the widths of the rectangles
For values of n != 1, the same trick applies... just sample the 1/6, 5/6 point for every slice in the n/3 evaluation.
What happens when you try and scale from n => 2n? The old function evaluations all fall on the boundaries between the new cells, and can't be reused:
n = 1 x n = 2 xx n = 4 xxxx
The implementation of "stage" in MidpointIntegrator is combining them like this:
n = 1 x n = 2 xxx n = 4 xxxxxxx
which is actually:
 tripling the number of integration slices at each step, not doubling, and
 moving the function evaluation points out of the midpoint.
In fact, what "stage" is actually doing is calculating a left Riemann sum, just without the leftmost point.
Here are the evaluation points for a left Riemann sum for comparison:
n = 2 xx n = 4 xxxx n = 8 xxxxxxxx
Here's the code from "stage" implementing this:
// number of new points in this stage... 2^n points at this stage, so we // have 2^n1 points. final long np = 1L << (n  1); double sum = 0; // spacing between adjacent new points final double spacing = diffMaxMin / np; // the first new point}}}} double x = min + 0.5 * spacing;}}}} for (long i = 0; i < np; i++) {}}}} sum += computeObjectiveValue(x); x += spacing;}}}} } // add the new sum to previously calculated result return 0.5 * (previousStageResult + sum * spacing);
Suggested Resolution
To keep the exact same results, I think the only solution is to remove the incorrect incremental "stage"; or convert it to the algorithm that implements the correct incremental increase by 3, and then don't call it by default.
Everything does work wonderfully if you expand n by a factor of 3 each time. Press discusses this trick in Numerical Recipes, section 4.4 (p136): http://phys.uri.edu/nigh/NumRec/bookfpdf/f44.pdf and notes that the savings you get still make it more efficient than NOT reusing function evaluations and implementing a correct scaleby2 each time.
Happy to provide more information if I can be helpful!