# NIntegrate Trapezoid rule with even subdivision producing poor results

NIntegrate[x^2, {x, 5, 9},
Method -> {"TrapezoidalRule", "Points" -> 2, "RombergQuadrature" -> False},
MaxRecursion -> 1]

returns the value 190 while the correct trapezoidal estimate is 5^2/2 + 6^2 + 7^2 + 8^2 + 9^2/2 = 202.

The correct trapezoidal weights for a width 4 interval are 1/2, 1, 1, 1, 1/2. I believe that the incorrect weights of 1, 1, 1/2, 1, 1/2 are being used instead.

Similar errors are made for TrapezoidRule with any even number of points.

For example, NIntegrate[x^2, {x, 5, 9}, Method -> {"TrapezoidalRule", "Points" -> 22}] returns warning messages about slow convergence.

I reported this to Wolfram several months ago and expected to see a fix in 10.0. However, the same results appear in 10.0. So now what I'm asking is if others can verify that this really is an error in the implementation of NIntegrate or if I'm making some mistake in my analysis.

-
Hi ! Unless, I am missing something this is not a question. You can edit it to add more information, if you want to ask something, otherwise it will be closed – Sektor Jul 11 '14 at 14:53
What reply to your report did you get from WRI? Did they acknowledge your finding as a bug? – m_goldberg Jul 11 '14 at 16:50
i concur this is a bug with the weights you have surmised. Seems fairly innocuous though, Points-3,MaxRecursion->0 gives the result you want. – george2079 Jul 11 '14 at 16:53
Here's the response I received from Wolfram: "I have filed a report with our developers and attached your name. If and when the situation is corrected in a future version of Mathematica you will be notified." – Todd Will Jul 11 '14 at 17:04
The issue seems to be only with Points->2 but persists for any MaxRecursion .. (except 0) – george2079 Jul 11 '14 at 17:32

It seems people misinterpreted the settings of "TrapezoidalRule". Allow me to set the record straight.

First, to get a "pure" trapezoidal rule evaluation, one should thwart the automatic interval splitting that is usually done by NIntegrate[]. This means a setting of MaxRecursion -> 0.

Second, the "Points" setting has also not been read right. In fact, setting "Points" -> k tells the internal rule generator to produce 2 k - 1 points, thus guaranteeing an even number of panels! That this is the case can be seen by calling the internal rule generator itself:

{λ, ω} = Most[NIntegrate`TrapezoidalRuleData[3, $MachinePrecision]] // Rationalize {{0, 1/4, 1/2, 3/4, 1}, {1/8, 1/4, 1/4, 1/4, 1/8}} and we recall that the internal generators are configured by default for the integral$\int_0^1 f(x)\,\mathrm dx$. So, to evaluate$\int_5^9 x^2\,\mathrm dx$with the$5\$-point trapezoidal rule, we can do this:

a = 5; b = 9;
(b - a) ω.(Function[x, x^2] /@ Rescale[λ, {0, 1}, {a, b}])
202

Compare:

NIntegrate[x^2, {x, 5, 9}, Method -> {"TrapezoidalRule", "Points" -> 3,
"RombergQuadrature" -> False, "SymbolicProcessing" -> 0},
MaxRecursion -> 0] // Quiet
202.

Here's another confirmatory check, using a different implementation of the trapezoidal rule:

n = 5;
tr = Interpolation[Table[{x, x^2}, {x, a, b, (b - a)/(n - 1)}],
InterpolationOrder -> 1];
Integrate[tr[x], {x, 5, 9}]
202
-

Your assumption is wrong. The correct result is

Integrate[x^2, {x, 5, 9}]

604/3 (201.333...)

With TrapezoidalRule you can only approximate this result:

NIntegrate[x^2, {x, 5, 9},
Method -> {"TrapezoidalRule", "Points" -> 3,
MaxRecursion -> 10,
PrecisionGoal -> 6]

201.333

-
Yes, I'm sorry to have said that the correct value was 202. What I meant was that the correct Trapezoid estimate using 2 subintervals with 1 level of recursion would be 202 and not 190. – Todd Will Jul 11 '14 at 15:35

I post this not as a specific answer but I think it may provide some insights. Experts and WRI would have to answer. The following will only work for positive valued functions (polygons breaking on x axis ->problems...remediable but I just post this as a quick insight). NIntegrate aims to provide best approximation within working precision. It seems there are local adaptive algorithms to refine. The out of the box answers for this example (a smooth function) is accurate.

trap[u_, n_, lim_] := Module[{int, h, p, pg, area},
int = Integrate[u, {x, lim[[1]], lim[[2]]}];
h = Range @@ ({#1, #2, (#2 - #1)/n}) & @@ lim;
p = Join @@ ({{#, 0}, {#, u /. x -> #}} & /@ h);
pg = #[[{1, 2, 4, 3}]] & /@ Partition[p, 4, 2];
area = Area /@ (Polygon /@ pg);
Show[Plot[Evaluate@u, {x, lim[[1]], lim[[2]]},
PlotStyle -> {Red, Thick},
PlotLabel ->
Grid[{{"Trapeziod", N[Total@area]}, {"Integrate",
N@int}, {"Relative error (%)",
N[100 Abs[Total@area/int - 1], 2]}}],
PlotRange -> {0, 1.1 Max[(u /. x -> #) & /@ h]}],
Graphics@
Polygon[#1]}, {Text[N[#2, 3],
RegionCentroid[Polygon@#1]]}} &, {pg, area}]]]

In the following I tabulate results for this specific example: Integrate, NIntegrate with methods specified (note max recursions->0...this seems to be first pass of trapezod area) and my "implementation" of trapezoid method just summing areas (using MMA 10 features):

Grid[Partition[
Column[#, Frame -> All] & /@
Table[{Row[{"Number of segments: ", j}], trap[x^2, j, {5, 9}],
Row[{"NIntegrate result for ", j - 1, " points: ",
Quiet@NIntegrate[x^2, {x, 5, 9},
Method -> {"TrapezoidalRule", "Points" -> j - 1,
"RombergQuadrature" -> False}, MaxRecursion -> 0]}]}, {j,
3, 10}], 4]]

Just for fun:

Manipulate[
trap[fun,
num, {5, 9}], {{fun, x^2, "function"}, {x^2, 1, x, Abs[Sin[x]], x^4,
Sin[2 x] + 2}}, {num, {2, 3, 4, 5, 6, 10, 20, 30, 50}}]

Again this is NOT meant to be a general trapezoid rule implementation but insight into NIntegrate in relation to the question.

-
really nice piece of work, but it doesn't relate to the issue at all since the BUG doesn't appear with maxrecursion->0. If you do points->2 recursion->1 you get a completely different result than points-3 recursion->0 based on the exact same sample points. – george2079 Jul 12 '14 at 11:46
@george2079 I accept your comment...I just sought to compare another way with what NIntegrate does to get some insight but I leave it to you and other numericists who know the underlying algorithms to clarify...I guess it was just some fun after a stressful week – ubpdqn Jul 12 '14 at 11:53

not an answer but a neat trick to pull out the weights that are used:

method = {"TrapezoidalRule", "Points" -> 2, "RombergQuadrature" -> False};
r = 2;

integrate once to learn the values:

xvals = Reap[i0=Quiet@NIntegrate[x^2, {x, 1, 5},
Method -> method, MaxRecursion -> r,
EvaluationMonitor :> Sow[x]]] // Last // Flatten // Union;

add a small imaginary value to one integration point at a time to learn the associated weight. ( the aim is to minimally impact the subdivision algorithm -- note this does break for high levels of recursion )

wts = (
f[x_?NumericQ] := x^2 + Boole[ x == # ] I/10^10 ;
10^10 Im@
Quiet@NIntegrate[f[x], {x, 1, 5},
Method -> method, MaxRecursion -> r] ) & /@ xvals

{1., 1., 0.5, 0.5, 0.25, 0.5, 0.25}

i0 (* 36 *)

verify these are the weights NIntegrate used originally.

wts.xvals^2 == i0 (* true *)

These weights are clearly wrong (not trapezoidal rule or anything else i can make sense of). The correct weights for trapezoidal integration are: (note the sample points are not uniformly spaced for 2 recursions )

trule= Mean /@
Partition[Differences@{xvals[[1]], Sequence @@ xvals, xvals[[-1]]}, 2, 1]

{0.5, 1., 0.75, 0.5, 0.5, 0.5, 0.25}

trule.xvals^2 (* 41.75 , much closer to exact 41.33..*)
-