Computational Physics - Department of Physics

5.3 Gaussian Quadrature 121

L 3 (x) = ( 5 x^3 − 3 x)/ 2 ,

and
L 4 (x) = ( 35 x^4 − 30 x^2 + 3 )/ 8.

The following simple function implements the above recursion relation of Eq. (5.11). for com-
puting Legendre polynomials of orderN.

// This function computes the Legendre polynomial of degreeN

doubleLegendre(intn,doublex)
{
doubler, s, t;
intm;
r = 0; s = 1.;
// Use recursion relation to generate p1 and p2
for(m=0; m < n; m++ )
{
t = r; r = s;
s = (2m+1)xr - mt;
s /= (m+1);
}// end of do loop
return s;
} // end of function Legendre

The variablesrepresentsLj+ 1 (x), whilerholdsLj(x)andtthe valueLj− 1 (x).

5.3.2 Integration points and weights with orthogonal polynomials

To understand how the weights and the mesh points are generated, we define first a polyno-
mial of degree 2 N− 1 (since we have 2 Nvariables at hand, the mesh points and weights forN
points). This polynomial can be represented through polynomial division by

P 2 N− 1 (x) =LN(x)PN− 1 (x)+QN− 1 (x),

wherePN− 1 (x)andQN− 1 (x)are some polynomials of degreeN− 1 or less. The functionLN(x)
is a Legendre polynomial of orderN.
Recall that we wanted to approximate an arbitrary functionf(x)with a polynomialP 2 N− 1 in
order to evaluate ∫ 1

− 1

f(x)dx≈

∫ 1

− 1

P 2 N− 1 (x)dx.

We can use Eq. (5.14) to rewrite the above integral as
∫ 1
− 1
P 2 N− 1 (x)dx=

∫ 1

− 1

(LN(x)PN− 1 (x)+QN− 1 (x))dx=

∫ 1

− 1

QN− 1 (x)dx,

due to the orthogonality properties of the Legendre polynomials. We see that it suffices to
evaluate the integral over
∫ 1
− 1 QN−^1 (x)dxin order to evaluate

∫ 1
− 1 P^2 N−^1 (x)dx. In addition, at the
pointsxkwhereLNis zero, we have

P 2 N− 1 (xk) =QN− 1 (xk) k= 0 , 1 ,...,N− 1 ,

and we see that through theseNpoints we can fully defineQN− 1 (x)and thereby the integral.
Note that we have chosen to let the numbering of the points runfrom 0 toN− 1. The reason
for this choice is that we wish to have the same numbering as the order of a polynomial of