
Type: Bug

Status: Open

Priority: Major

Resolution: Unresolved

Affects Version/s: 3.6.1

Fix Version/s: None

Labels:
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!