Gauss-Legendre Quadrature

In the method of Newton-Cotes quadrature based on $n+1$ equally spaced node points, we are free to choose the weights $c_0,\cdots,c_n$ as $n+1$ variables to achieve the highest degree of accuracy of $n$ if it is odd, or $n+1$ if $n$ is even. In comparison, in the method of Gauss-Legendre quadrature considered below, we are free to choose the $n+1$ node points, as well as the $n+1$ weights. As number of free variables in the method is doubled, its degree of accuracy is also doubled to reach $2n+1$.

Reconsider the quadrature rule based on the Lagrange interpolation of the integrand $f(x)$:

$$I[f]=\int_a^b f(x)\,dx\approx \int_a^b L_n(x)\,dx=\sum_{i=0}^n f(x_i)\int_a^b l_i(x)\,dx =\sum_{i=0}^n c_i f(x_i)$$ (155)

where $c_i\;(i=0,\cdots,n)$ are the weights independent of the integrand $f(x)$:

$$c_i=\int_a^b l_i(x)\,dx=\int_a^b \prod_{j=0,\;j\ne i}^n \frac{x-x_j}{x_i-x_j}\,dx$$ (156)

To derive the algorithm, we first make the following assumptions:

• The integrand $f(x)$ is a polynomial of degree no higher than $2n+1$,
• The integral limits are $a=-1$ and $b=1$,
• The node points $x_0,\cdots,x_n$, called Gauss points, are the $n+1$ roots of an (n+1)th polynomial $p_{n+1}(x)$ in an orthogonal polynomial family, here assumed to be the Legendre polynomial, i.e., $p_{n+1}(x_i)=0,\;(i=0,\cdots,n)$.

Dividing $f(x)$ by $p_{n+1}(x)$, we get

$$\frac{f(x)}{p_{n+1}(x)}=Q(x)+\frac{R(x)}{p_{n+1}(x)}, \;\;\;\;\;\;\; \;\;\;\;\;\;\;f(x)=p_{n+1}(x)Q(x)+R(x)$$ (157)

where the quotient $Q(x)$ is a polynomial of degree $(2n+1)-(n+1)=n$, and the remainder $R(x)$ is a polynomial of degree no higher than $n$. Evaluating $f(x)$ at any Gauss point $x_i$, we have

$$f(x_i)=p_{n+1}(x_i)Q(x_i)+R(x_i)=R(x_i)$$ (158)

Also, as $Q(x)$ is a polynomial of degree no higher than $n$, it can be expressed as a linear combination of the Legendre polynomials of degree no higher than $n$, its inner product with $p_{n+1}(x)$ is zero due to the orthogonality of the Legendre polynomials:

$$\int_{-1}^1 p_{n+1}(x)Q(x)\,dx =\int_{-1}^1 p_{n+1}(x)\left[\sum_{j=0}^n a_jp_j(x)\right]\,dx =\sum_{j=0}^n a_j\int_{-1}^1 p_{n+1}(x)p_j(x)\,dx=0$$ (159)

Now the integral can be carried out by

$$\int_{-1}^1 f(x)\,dx = \int_{-1}^1 p_{n+1}(x)Q(x)\,dx +\int_{-1}^1 R(x)\,dx$$

$$= \int_{-1}^1 R(x)\,dx=\sum_{i=0}^n c_iR(x_i)=\sum_{i=0}^n c_if(x_i)$$

We see that under the assumptions above, the integral can be calculated exactly as a linear combination of $f(x_i)$, the integrand evaluated at the $n+1$ Gauss node points $x_0,\cdots,x_n$, the roots of the Legendre polynomial $p_{n+1}(x)$.

While the weights $c_0,\cdots,c_n$ of the Gauss-Legendre quadrature can be found by Eq. (?) by integrating the Lagrange basis polynomials $l_i(x)$, they can also be found more conveniently by the method of undetermined coefficients. Specifically, byassuming the integrand to be polynomials of different orders $f(x)=x^k,\;(k=0,\cdots,n)$, we get the following $n+1$ equations:

$$\sum_{i=0}^n x_i^k c_i=\int_{-1}^1x^k\,dx =\frac{1-(-1)^{k+1}}{k+... ...s odd}\\ 2/(k+1)&\mbox{$k$\ is even}\end{array}\right. ,\;\;\;\;(k=0,\cdots,n)$$ (160)

Or, alternatively, by assuming $f(x)=p_k(x),\;(k=0,\cdots,n)$, we get

$$\sum_{i=0}^n p_k(x_i)c_i=\int_{-1}^1 p_k(x)\,dx =\left\{\begin{array}{cc}2 & k=0\\ 0 & k>0\end{array}\right. \;\;\;\;(k=0,\cdots,n)$$ (161)

Solving either of these systems of $n+1$ equations we get the $n+1$ weights $c_0,\cdots,c_n$.

We can also assume the integrand to be a special function $f(x)=p_{n+1}(x)p'_{n+1}(x)/(x-x_i)$ with $x_i$ to be one of the roots of $p_{n+1}(x)$, i.e., $p_{n+1}(x_i)=0$. The integral of this function can be found by the Gauss-Legendre quadrature rule:

$$\int_{-1}^1 f(x)\,dx = \int_{-1}^1\frac{p_{n+1}(x)p'_{n+1}(x)}{x-x_i}\,dx =\sum_{j=0}^nc_j \frac{p_{n+1}(x_j)p'_{n+1}(x_j)}{x_j-x_i}$$

As $x_j\;(j=0,\cdots,n)$ are the roots of $p_{n+1}(x)$, all terms in the summation are zero except the ith one, an indeterminate form $0/0$ which can be evaluated by L'Hopital's rule:

$$\frac{p_{n+1}(x_i)p'_{n+1}(x_i)}{x_i-x_i} =\lim\limits_{x\rightar... ..._{x\rightarrow x_i}\frac{[p_{n+1}(x)p'_{n+1}(x)]'}{(x-x_i)'} =[p'_{n+1}(x_i)]^2$$ (162)

Now the equation above becomes

$$\int_{-1}^1 \frac{p_{n+1}(x)p'_{n+1}(x)}{x-x_i}\,dx=c_i[p'_{n+1}(x_i)]^2$$ (163)

On the other hand, the integral can be carried out using integration by parts with $u(x)=p_{n+1}(x)/(x-x_i)$ and $dv(x)=p'_{n+1}(x)dx$, we get $u(x)\,v(x)=p^2_{n+1}(x)/(x-x_i)$, $v(x)\,du(x)=p_{n+1}u'(x)dx$ and

$$\int_{-1}^1 \frac{p_{n+1}(x)p'_{n+1}(x)}{x-x_i}\,dx = \frac{p_{n+1}^2(x)}{x-x_i}\bigg\vert _{-1}^1-\int_{-1}^1 p_{n+1}(x)u'(x)\,dx$$

$$= \frac{p_{n+1}^2(1)}{1-x_i}-\frac{p_{n+1}^2(-1)}{-1-x_i} =\frac{1}{1-x_i}+\frac{1}{1+x_i}=\frac{2}{1-x_i^2}$$

Here the integral of the second term on the right-hand side is an inner product of $u'(x)$, a polynomial of degree lower than $n$, and $p_{n+1}(x)$, which is zero due to the orthogonality of the Legendre polynomials, and we have also used the fact that $p^2_{n+1}(\pm 1)=1$. Equating the two expressions for the integral, we get

$$c_i=\frac{2}{(1-x_i^2)[p'_{n+1}(x_i)]^2},\;\;\;\;\;(i=0,\cdots,n)$$ (164)

Finally we show that the integral limits $[-1,\;1]$ can be generalized to $[a,\;b]$ by the following linear mapping:

$$x=\frac{b-a}{2}u+\frac{a+b}{2} =\left\{\begin{array}{ll}a & u=-1\\ b & u=1\end{array}\right.$$ (165)

with $dx=(b-a)du/2$. Now we get the Gauss-Legendre quadrature rule:

$$I[f]=\int_a^b f(x)\,dx = \int_{-1}^1 f\left(\frac{b-a}{2}u+\frac{a+b}{2}\right)\frac{b-a}{2}\,du$$

$$\approx \frac{b-a}{2}\;\sum_{i=0}^{n-1} c_if\left(\frac{b-a}{2}x_i+\frac{b+a}{2}\right)$$

which is exact if the integrand is a polynomial of degree no higher than $2n+1$, twice as high as the the degree of accuracy of the Newton-Cotes quadrature based on $n+1$ equally spaced node points.

Example: When $n=2$, the Gauss nodes are the roots of $p_{n+1}(x)=p_3(x)=(5x^3-3x)/2$: $x_0=-\sqrt{3/5},\;x_1=0,\; x_2=\sqrt{3/5}$, and the weights can be found Eq. (?):

$$c_0 = \int_{-1}^1\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}\,dx =\frac{5}{6}\int_{-1}^1 x(x-\sqrt{3/5})\,dx=\frac{5}{9}$$

$$c_1 = \int_{-1}^1\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}\,dx =-\frac{5}{3}\int_{-1}^1 (x^2-3/5)\,dx=\frac{8}{9}$$

$$c_2 = \int_{-1}^1\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}\,dx =\frac{5}{6}\int_{-1}^1 x(x+\sqrt{3/5})\,dx=\frac{5}{9}$$

Alternatively and more conveniently, we can also obtain these weights by Eq. (??). We first find

$$p'_{n+1}(x)=p'_3(x)=\frac{d}{dx}\left[\frac{1}{2}(5x^3-3x)\right] =\frac{1}{2}(15x^2-3)$$ (166)

and then evaluate it at the Gauss nodes to get $p'_3(x_0)=p'_3(x_2)=3$, and $p'_3(x_1)=-3/2$. Now we can get the weights:

$$c_0 = \frac{2}{(1-x_0^2)[p'_3(x_0)]^2}=\frac{2}{2/5\times 9}=\frac{5}{9}$$

$$c_0 = \frac{2}{(1-x_1^2)[p'_3(x_1)]^2}=\frac{2}{9/4}=\frac{8}{9}$$

$$c_0 = \frac{2}{(1-x_2^2)[p'_3(x_2)]^2}=\frac{2}{2/5\times 9}=\frac{5}{9}$$

There weights can also be obtained by solving either of the following equivalent equation systems:

$$\left\{\begin{array}{l}x_0^0c_0+x_1^0c_1+x_2^0c_2=2\\ x_0^1c_0+x... ...)\sqrt{3/5}=0\\ x_0^2c_0+x_1^2c_1+x_2^2c_2=(c_0+c_2)3/5=2/3 \end{array}\right.$$ (167)

and

$$\left\{\begin{array}{l} p_0(x_0)c_0+p_0(x_1)c_1+p_0(x_2)c_2=c_0+c... ... p_2(x_0)c_0+p_2(x_1)c_1+p_2(x_2)c_2=4c_0/5-c_1+4c_2/5=0\\ \end{array}\right.$$ (168)

where $p_0(x)=1,\;p_1(x)=x,\;p_2(x)=(3x^2-1)/2$.

Now consider integrating a simple polynomial $f(x)=x^5$ from $a=0$ to $b=1$. First, to linearly map the integral limits to $[-1,\,1]$, we find $(b-a)/2=(b+a)/2=1/2$ and get

$$\int_0^1 x^5\,dx =\int_{-1}^1\left(\frac{u+1}{2}\right)^5\,\frac{1}{2}\,du =\frac{1}{6}$$ (169)

Of course, we are interested in calculating the integral by the Gauss-Legendre quadrature by Eq. (???)

$$\int_0^1 x^5\,dx = \frac{1}{2}\left[c_0f\left(\frac{x_0+1}{2}\right) +c_1f\left(\frac{x_1+1}{2}\right) +c_2f\left(\frac{x_2+1}{2}\right)\right]$$

$$= \frac{1}{2}\left[ \frac{5}{9}\left(\frac{-\sqrt{3/5}+1}{2}\right)... ...c{5}{9}\left(\frac{\sqrt{3/5}+1}{2}\right)^5\right] =\frac{1}{6}\approx 1.66667$$

The result is exact as the degree of the integrand $f(x)=x^5$ is no higher than $2n+1=5$. In comparison, the result by the Newton-Cotes quadrature based on three equally spaced points $x_0=0$, $x_1=0.5$ and $x_2=1$ with $h=0.5$ is $h(x_0^5+4x_1^5+x_2^5)/3=0.1875$. If the integrand is $f(x)=x^6$ of degree higher than $2n+1=5$, the result by the Gauss-Legendre quadrature is no longer exact:

$$\frac{1}{2}\left[ \frac{5}{9}\left(\frac{-\sqrt{3/5}+1}{2}\right)... ...1}{2}\right)^6 +\frac{5}{9}\left(\frac{\sqrt{3/5}+1}{2}\right)^6\right] =0.1425$$

$$\ne \int_0^1 x^6\,dx=\frac{1}{7}=0.1429$$

But this is still more accurate than the result by the Newton-Cotes quadrature: $h(x_0^6+4x_1^6+x_2^6)/3=0.1771$.

Here is the Matlab function to generate the Gauss points and the coefficients for the Gauss-Legendre quadrature:

function [x,w]=GaussLegendre(n,method)  % n: total number of points
                    % x and w are the Gauss nodes and weights
    P=zeros(n+1);   % coefficients of n+1 Legendre polynomials
    b=zeros(n,1);   % vector on right side for undetermined coefficients
    w=zeros(n,1);   % weights for the quadrature
    P([1,2],1)=1;   % coefficients of p_0(x)=1 and p_1(x)=x in first 2 rows
    for k=1:n-1     % recursively generat remaining p_2(x) to p_{n+1}(x)
        P(k+2,1:k+2)=((2*k+1)*[P(k+1,1:k+1) 0]-k*[0 0 P(k,1:k)])/(k+1);
    end
    % (i+1)th row of P contains coefficients of p_i(x) (i=0,...,n)
    x=sort(roots(P(n+1,1:n+1)))  % find all n+1 roots of p_{n+1}(x) (Gauss points)

    switch method
    % Method 1
    dp=P(n+1,:).*[n:-1:0];       % coefficients of p'_{n+1}(x)
    for i=1:n
        t=x(i).^(n-1:-1:0);      % p'_n(x_i)
        t=dp(1:n)*t';
        w(i)=2/(1-x(i)^2)/t/t;   % weights for the quadrature
    end

    % method of undetermined coefficients: f(x)=x^k
    A=zeros(n);
    for k=1:n
        A(k,:)=x.^(k-1);              % matrix on left side
        b(k)=(1-(-1)^k)/k;            % vector on right side
    end

    % Method of undetermined coefficients: f(x)=p_k(x)
    A=zeros(n);
    for k=1:n    % evaluate p_0(x)...p_{n-1}(x) at roots of p_n(x)
        A(k,:)=polyval(P(k,1:k),x)';  % matrix on left side
    end
    b=[2;zeros(n-1,1)];               % vector on right side
    w=inv(A)*b   % find all weights 
end