[
https://issues.apache.org/jira/browse/MATH-1558?focusedWorklogId=502880&page=com.atlassian.jira.plugin.system.issuetabpanels:worklog-tabpanel#worklog-502880
]
ASF GitHub Bot logged work on MATH-1558:
----------------------------------------
Author: ASF GitHub Bot
Created on: 20/Oct/20 22:38
Start Date: 20/Oct/20 22:38
Worklog Time Spent: 10m
Work Description: sritchie commented on a change in pull request #161:
URL: https://github.com/apache/commons-math/pull/161#discussion_r508881160
##########
File path:
src/test/java/org/apache/commons/math4/analysis/integration/MidPointIntegratorTest.java
##########
@@ -48,8 +67,9 @@ public void testLowAccuracy() {
double expected = -3697001.0 / 48.0;
double tolerance = FastMath.abs(expected *
integrator.getRelativeAccuracy());
double result = integrator.integrate(Integer.MAX_VALUE, f, min, max);
- Assert.assertTrue(integrator.getEvaluations() < Integer.MAX_VALUE / 2);
+ Assert.assertTrue(integrator.getEvaluations() < Integer.MAX_VALUE / 3);
Review comment:
This test might be redundant now that we know exactly how many
evaluations we should have, as a function of iterations.
----------------------------------------------------------------
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:
[email protected]
Issue Time Tracking
-------------------
Worklog Id: (was: 502880)
Time Spent: 0.5h (was: 20m)
> 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: 0.5h
> 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)
