Skip to content

Gaussian Quardrature

Qiang Zhu edited this page Oct 5, 2017 · 1 revision

Gaussian Quadrature

Gaussian Quadrature is a numerical integration method to approximate the definite integral of a function. The Gaussian Quadrature guarantees a precise calculation of integration of all polynomial up to degree of (2n-1) given a weighted sum of function values at n points within the domain of integration [-1,1].

An integral function can be expressed with the following general approximation:

First Equation

where is the corresponding weight for . In the case of the trapezoidal method the weights are (h/2, h/2) and (h/3, 4h/3, h/3) in the case of Simpson’s method. In addition, the points are separated equidistantly. Gauss-Legendre integration is introduced in this section. However, a basic introduction of Legendre polynomial is required.

Legendre Polynomial

Let’s consider an ordinary differential equation of this form:

second equation

where l∈. The equation (2) can be solved by the standard power series method. The equation is called Legendre’s differential equation.

The solution of Legendre’s differential equation is called Legendre polynomials:

In fact, these solutions are orthogonal polynomials which have the following property:

Gauss-Legendre Integration

Since we had introduced some property of Legendre polynomials (equation 3), we can talk about Gauss-Legendre integrations. As a starting point for our discussion, we will assume that the function f(x) can be approximated by a polynomial of order 2n-1, f(x)=P_(2n-1) (x) because the Gaussian Quadrature guarantees a precise calculation of integration of all polynomial up to degree of (2n-1). As described in equation 1, we can take the integral by substituting a and b to -1 and 1:

In other words, knowing the n abscissas and the corresponding weight values we can compute the integral. In the first step, let’s decompose (x) into two terms:

where are Legendre polynomials of order n. and are polynomials of order (n-1). Let’s integrate both sides of equation 5:

Sixth equation

Based on equation 3, we know that the first term of the right-hand side is zero because (x) is orthogonal to (x). Apparently, Gauss know that (x) has n zeros in [-1,1], which we will label . For these points, we have the relation:

Expressing (x) with Legendre polynomials:

Eighth equation

Thus, we can express equation 7 in the form of equation 8:

Ninth equation

Expressing equation 9 into matrix form:

Tenth equation

In order to solve for the , we need to invert the matrix :

Eleventh equation

Finally, we can get back to our integral equation (4):

Twelveth equation

Due to this Legendre polynomials are only defined in [-1,1], we can change our integration limits from [a,b] [-1,1] by the following substitution:

Thirteenth equation

This transformation leads to:

Fourteenth equation

Python code

The code of Gaussian Quadrature method is written in python script:

from numpy import ones,copy,cos,tan,pi,linspace

def gaussxw(N):

    # Initial approximation to roots of the Legendre polynomial
    a = linspace(3,4*N-1,N)/(4*N+2)
    x = cos(pi*a+1/(8*N*N*tan(a)))

    # Find roots using Newton's method
    epsilon = 1e-15
    delta = 1.0
    while delta>epsilon:
        p0 = ones(N,float)
        p1 = copy(x)
        for k in range(1,N):
            p0,p1 = p1,((2*k+1)*x*p1-k*p0)/(k+1)
        dp = (N+1)*(p0-x*p1)/(1-x*x)
        dx = p1/dp
        x -= dx
        delta = max(abs(dx))

    # Calculate the weights
    w = 2*(N+1)*(N+1)/(N*N*(1-x*x)*dp*dp)

    return x,w

def f(x):
    return x**5 - 2*x + 1

N = 3
a,b = 0.0, 2.0

# Calculate the sample points and weights, then map them
# to the required integration domain
x,w = gaussxw(N)
xp = 0.5*(b-a)*x + 0.5*(b+a)
wp = 0.5*(b-a)*w

# Perform the integration
s = 0.0
for k in range(N):
    s += wp[k]*f(xp[k])
print(s)
8.66666666667

Acknowledgement

The original wiki was created by Howard Yanxon, a graduate student in my class.