[
https://issues.apache.org/jira/browse/MATH-1458?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=16462262#comment-16462262
]
Alex D Herbert commented on MATH-1458:
--------------------------------------
I've branched the git repo and will create tests that the current
SimpsonItegrator fails. I'll then fix it and submit a pull request for review.
> Simpson Integrator computes incorrect value at minimum iterations=1 and
> wastes an iteration
> -------------------------------------------------------------------------------------------
>
> Key: MATH-1458
> URL: https://issues.apache.org/jira/browse/MATH-1458
> Project: Commons Math
> Issue Type: Bug
> Affects Versions: 3.6.1
> Environment: openjdk version "1.8.0_162"
>
> OpenJDK Runtime Environment (build 1.8.0_162-8u162-b12-0ubuntu0.16.04.2-b12)
>
> OpenJDK 64-Bit Server VM (build 25.162-b12, mixed mode)
> Reporter: Alex D Herbert
> Priority: Minor
> Labels: documentation, easyfix, newbie, patch
>
> 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:
> {code:java}
> return (4 * qtrap.stage(this, 1) - qtrap.stage(this, 0)) / 3.0;
> {code}
> This can be fixed using:
> {code:java}
> final double s0 = qtrap.stage(this, 0);
> return (4 * qtrap.stage(this, 1) - s0) / 3.0;{code}
> 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:
> {code:java}
> final double t = qtrap.stage(this, getIterations());
> incrementCount();
> {code}
> This should be:
> {code:java}
> incrementCount();
> final double t = qtrap.stage(this, getIterations());
> {code}
> 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:
> {code:java}
> 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;
> }
> }
> {code}
> 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.
>
--
This message was sent by Atlassian JIRA
(v7.6.3#76005)
