Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.
NIntegrate[x^2, {x, 5, 9}, 
      Method -> {"TrapezoidalRule", "Points" -> 2, "RombergQuadrature" -> False}, 
      MaxRecursion -> 1]

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

The correct trapezoid 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.

share|improve this question
    
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 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 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 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 at 17:04
    
The issue seems to be only with Points->2 but persists for any MaxRecursion .. (except 0) –  george2079 Jul 11 at 17:32

3 Answers 3

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, 
   "RombergQuadrature" -> False},
 MaxRecursion -> 10,
 PrecisionGoal -> 6]

201.333

share|improve this answer
    
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 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@
    MapThread[{{Opacity[0.3], EdgeForm[Black], FaceForm[Yellow], 
        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]]

enter image description here

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}}]

enter image description here

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

share|improve this answer
    
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 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 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..*)
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.