I am trying to make Mathematica do the following: Split up the interval [0,1] into $n$ equal intervals. Then on each interval apply the Gaussian Quadrature for 2 points.

I tried to use the method from mky question here:Using Composite Newton-Cotes integration rules in Mathematica

But the results I am getting are way too accurate.

Here is the code:

n = 1;                                                                                                            
NIntegrate[Sin[x]/x, Evaluate@Flatten@{x, Subdivide[0., 1., n]}, 
 Method -> {"GaussBerntsenEspelidRule", "Points" -> 2}, 
 MaxRecursion -> 0]

So what this should be doing is just do the Gaussian Quadrature with 2 points on $[0,1]$ but I am getting the answer with at least $10^{-5}$ accuracy which should not be happening. What did I do wrong?


1 Answer 1


You're just lucky (as far as getting low error):

{abs, wts, err} = 
 NIntegrate`GaussBerntsenEspelidRuleData[2, MachinePrecision]
  {{0.0469101, 0.230765, 0.5, 0.769235, 0.95309},
   {0.118463, 0.239314, 0.284444, 0.239314, 0.118463},
   {0.155257, -0.439701, 0.568889, -0.439701, 0.155257}}
(Sin[x]/x /. x -> abs).wts
(Sin[x]/x /. x -> abs).err
  0.946083       <-- integral estimate
  0.0000639286   <-- estimated error bound
(Sin[x]/x /. x -> abs).wts - Integrate[Sin[x]/x, {x, 0, 1}]
  3.31957*10^-14  <-- actual error (less than the bound)

The above code reproduces the NIntegrate result:

(Sin[x]/x /. x -> abs).wts -
 NIntegrate[Sin[x]/x, Evaluate@Flatten@{x, Subdivide[0., 1., nn]}, 
  Method -> {"GaussBerntsenEspelidRule", "Points" -> 2}, 
  MaxRecursion -> 0]

Why are we lucky in this case? The error is equal to the integral of the difference of the function and the interpolating polynomial through abscissae abs, which roughly has the same area above and below the x axis:

 InterpolatingPolynomial[Transpose@{abs, (Sin[x]/x /. x -> abs)}, x] -
  Sin[x]/x // Evaluate,
 {x, 0, 1}]
  • 1
    $\begingroup$ Does this condition happen more often when the y-values are so small? Good catch. +1 $\endgroup$
    – BBirdsell
    Commented Nov 14, 2020 at 2:30
  • $\begingroup$ @BBirdsell No, I don't think so. The smaller the values, the smaller the expected error. You can (theoretically) construct functions that violate error estimates and ones for which the positive and negative errors in the function approximation cancel out in the integration, as they almost do in the graph. In some sense, I think the OP chose an example in which the error cancels out better than one might expect. (I also wonder if the OP realized that "Points -> n results in 2n+1 abscissae and not just n.) Here's a bad function for the OP's NIntegrate: 1/(10^4 (x - 36/100)^2 + 1). $\endgroup$
    – Michael E2
    Commented Nov 14, 2020 at 5:39
  • $\begingroup$ @MichaelE2 Thank you for your answer. So Mathematica is doing what I want it to do, but in this case we get lucky. That is interesting. Thank you very much for explaining it to me. $\endgroup$
    – 2132123
    Commented Nov 14, 2020 at 16:03
  • $\begingroup$ @MichaelE2 I just noticed your comment and you are correct in regards to my confusion with Points option. However, I could not specify Points less than $2$ since it would give me an error. How do I force it to do Gaussian Quadrature rule with 2 points on each subinterval, not $2*2+1=5$ points? $\endgroup$
    – 2132123
    Commented Nov 14, 2020 at 16:41
  • $\begingroup$ @2132123 NIntegrate`GaussRuleData[2, MachinePrecision], but it has no error estimate weights. $\endgroup$
    – Michael E2
    Commented Nov 14, 2020 at 18:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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