2
$\begingroup$

A bit of a simple question here. Recently, I was evaluating some line integrals using a Gaussian quadrature rule on $[-1,1]$ where the abscissae $x_{i}$ are just the roots of the Legendre polynomial $P_{n}(x)$ and the weights are $$ w_{i} = \frac{-2}{(1-x_{i}^{2})[P_{n}^{'}(x_{i})]^{2}} $$

Embarrassingly, I just spent an entire day trying to debug my code and painstakingly doing line integrals out by hand for comparison before realizing that I wasn't using enough quadrature points.

From the wiki page (http://en.wikipedia.org/wiki/Gaussian_quadrature#Error_estimates), I thought a Gaussian quadrature rule of order $n$ was accurate to polynomials of order $2n-1$.

I was using a 12-point rule and my integrands were all polynomials with order of at most 3. It turns out this wasn't enough, but increasing the number of quadrature points to 20 fixed everything.

Can anybody explain or suggest some readings on what exactly is meant by the order of the quadrature rule?

$\endgroup$
  • $\begingroup$ You have the idea of the order correct. There must be something wrong with your implementation. Are you using single-precision floating point (float in C, real*4 in Fortran, or some other 32-bit representation)? $\endgroup$ – Bill Barth Apr 2 '15 at 3:02
  • $\begingroup$ Can you also describe your integrals? If you've parametrized them somehow, the order of that parametrization also affects your end integral. $\endgroup$ – Jesse Chan Apr 2 '15 at 4:40
  • $\begingroup$ @BillBarth This was double precision in C (double). I had copied the weights and abscissae from a table, which had 256 digits of accuracy. I was only using the first 16 (15 + last digit rounded), which I assumed would be double-precision accuracy. I replaced the original 12-point rule with all the digits available. Now the integration works without a problem. $\endgroup$ – Justin Dong Apr 2 '15 at 14:55
  • 2
    $\begingroup$ I don't completely understand why that made a difference. No matter, you realize that 12th order Gauss quadrature is massive overkill for integrating cubics, right? You are doing way more work than you need to. $\endgroup$ – Bill Barth Apr 2 '15 at 15:16
2
$\begingroup$

If your integrals were all with order of at most 3, then certainly there is an implementation error. Besides, all your weights should be positive. So the minus sign in your definition of weights is wrong.

Here is a python code snippet, in which I integrate with 2 quadrature points arbitrary polynoms of degree 3 within machine precision.

import numpy.polynomial.legendre as leg
import numpy.polynomial
import numpy as np

def F(x,c):
    V=numpy.polynomial.polynomial.polyvander(x,c.size)
    i=1./(np.arange(c.size)+1) # integration factor                          
    return np.dot(i*V[:,1:],c)                      

def f(x,c):
    V=numpy.polynomial.polynomial.polyvander(x,c.size-1)
    return np.dot(V,c)                          

[x,w]=leg.leggauss(2)

c=np.random.randn(4)
a1=F(1,c)-F(-1,c)
a2=np.sum(f(x,c)*w)
print a1
print a2
print a1-a2

Out[13]: array([ 0.])
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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