1

I would like to write a program to calculate abscissas and weights of Gauss-Legendre Quadrature. I found the following source

http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature.

I don't understand why the first guess $x_0$ for the $i$-th root of a $n$-order polynomial $P_n$ is defined as

$x_0=\cos\left(\pi\frac{i-\frac{1}{4}}{n+\frac{1}{2}}\right)$

Could someone explain this to me?

LiN
  • 11
  • Have a look at this paper http://www.ams.org/journals/tran/1936-039-01/S0002-9947-1936-1501831-2/S0002-9947-1936-1501831-2.pdf – uranix May 12 '15 at 15:14

1 Answers1

1

H. Bruns proved in 1881 that the arccosine of the $i$th root of $P_n$ lies in the interval $$ \left[ \frac{2i-1}{2n+1}\pi, \ \frac{2i}{2n+1}\pi\right] $$ (When working with Legendre polynomials, it is often convenient to compose them with cosine, i.e., look at $P_n(\cos\theta)$ instead of $P_n(x)$.)

So it's natural to take the cosine of the midpoint of this interval as the first guess for the root. One source for the above is Szegő's paper mentioned by uranix; another one is Szegő's book Orthogonal Polynomials, Theorem 6.21.2.

That said, I think that using a root search procedure for each root individually is not particularly efficient, as this does not take advantage of the special structure of the polynomials $P_n$. The recurrence property of $P_n$ leads to a family of tridiagonal matrices for which $P_n$ is (a multiple of) the characteristic polynomial. It's easier to find the eigenvalues of a tridiagonal matrix than it is to find the roots of a generic polynomial: see the paper Calculation of Gauss quadrature rules by G. H. Golub and J. H. Welsch.

Spoiler: the paper Is Gauss Quadrature Better than Clenshaw–Curtis? by Lloyd N. Trefethen has the following 7-line implementation of Gaussian quadrature on $[-1,1]$ in Matlab. 

function I = gauss(f,n)               % (n+1)-pt Gauss quadrature of f
beta = .5./sqrt(1-(2*(1:n)).ˆ(-2));   % 3-term recurrence coeffs
T = diag(beta,1) + diag(beta,-1);     % Jacobi matrix
[V,D] = eig(T);                       % eigenvalue decomposition
x = diag(D); [x,i] = sort(x);         % nodes (= Legendre points)
w = 2*V(1,i).ˆ2;                      % weights
I = w*feval(f,x);                     % the integral
  • but $x_0=\cos\left(\pi\frac{i-\frac{1}{4}}{n+\frac{1}{2}}\right)$ is not the middle point of that interval. So why we take the value of $\cos()$? – LiN May 15 '15 at 17:43
  • Corrected. In the book Szegő consistently writes $x=\cos\theta$, and the bounds I quoted are for $\theta$. –  May 15 '15 at 17:59
  • I have problem with understand one thing. As we know the abscissas of $n$point Gaussian quadrature in the interval (a,b) are precisely the roots of the orthogonal polynomial $p_N(x)$ for the same interval. Let $N=2$. $P_2(x)=\frac{1}{2}(3x^2-1)$ Roots of $P_2(x)$ are: $0.5773502692, -0.5773502692$. We got what we had expected.

    If we substitute $x=\cos(\theta)$ we have $P_2(cos\theta)=\frac{1}{2}(3(\cos(\theta))^2-1)$ Roots: $0.9553166181, 2.186276036$.

    It is a completely different result (?!). So why we use $x=\cos\theta$?

     – LiN May 17 '15 at 13:58
  • $\cos(0.9553166181) = 0.57735...$. It's called substitution. Like when you are trying to solve $x^4-3x^2+2=0$, you substitute $x^2=u$ and solve for $u$. Here they substitute $x=\cos\theta$ and solve for $\theta$. To see why this helps, you have to read that book on orthogonal polynomials. –  May 17 '15 at 14:51