[jira] [Work logged] (MATH-1558) MidpointIntegrator doesn't implement the midpoint method correctly

["ASF GitHub Bot (Jira)" <ji...@apache.org>]

     [ https://issues.apache.org/jira/browse/MATH-1558?focusedWorklogId=502877&page=com.atlassian.jira.plugin.system.issuetabpanels:worklog-tabpanel#worklog-502877 ]

ASF GitHub Bot logged work on MATH-1558:
----------------------------------------

                Author: ASF GitHub Bot
            Created on: 20/Oct/20 22:26
            Start Date: 20/Oct/20 22:26
    Worklog Time Spent: 10m 
      Work Description: sritchie opened a new pull request #161:
URL: https://github.com/apache/commons-math/pull/161


   This PR implements the fix described in https://issues.apache.org/jira/browse/MATH-1558, and fixes the implementation of the "stage" method in MidPointIntegrator.


----------------------------------------------------------------
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.

For queries about this service, please contact Infrastructure at:
users@infra.apache.org


Issue Time Tracking
-------------------

            Worklog Id:     (was: 502877)
    Remaining Estimate: 0h
            Time Spent: 10m

> MidpointIntegrator doesn't implement the midpoint method correctly
> ------------------------------------------------------------------
>
>                 Key: MATH-1558
>                 URL: https://issues.apache.org/jira/browse/MATH-1558
>             Project: Commons Math
>          Issue Type: Bug
>    Affects Versions: 3.6.1
>            Reporter: Sam Ritchie
>            Priority: Major
>              Labels: integration, midpoint
>         Attachments: midpoint.patch
>
>          Time Spent: 10m
>  Remaining Estimate: 0h
>
> 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 {{(b-a)/n}} wide and as tall as the midpoint of the interval - ie, {{((b-a)/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.
> h2. 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:
> {noformat}
> n = 1 |--------x--------|
> n = 3 |--x--|--x--|--x–-|
> {noformat}
> 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:
> {noformat}
> n = 1 |-------x-------|
> n = 2 |---x---|---x---|
> n = 4 |-x-|-x-|-x-|-x-|
> {noformat}
> The implementation of "stage" in MidpointIntegrator is combining them like this:
> {noformat}
> n = 1 |-------x-------|
> n = 2 |---x---x---x---|
> n = 4 |-x-x-x-x-x-x-x-|
> {noformat}
> 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 left-most point.
> Here are the evaluation points for a left Riemann sum for comparison:
> {noformat}
> n = 2 x-------x--------
> n = 4 x---x---x---x----
> n = 8 x-x-x-x-x-x-x-x--
> {noformat}
> Here's the code from "stage" implementing this:
> {code}
> // number of new points in this stage... 2^n points at this stage, so we
> // have 2^n-1 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);
> {code}
> h2. 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/f4-4.pdf] and notes that the savings you get still make it more efficient than NOT reusing function evaluations and implementing a correct scale-by-2 each time.
> Happy to provide more information if I can be helpful!
>  



--
This message was sent by Atlassian Jira
(v8.3.4#803005)