Approximating Integrals with the Gaussian Quadrature (w/ example)

 

This is more like a memory dump so I will have a backup in case I’ll ever need it again. And if someone else finds this information useful, the better it is.

Intro

The Gaussian Quadrature is a method used to approximate the value of a given integral by choosing a set of points (x₁, x₂, x₃, … x_n) that will maximize the accuracy. Basically the integral can be approximated using coefficients and known values of our function. It’s a pretty neat method since it doesn’t require many points and it works for a set of integrals - not only for one.

∫1−1f(x)w(x)dx=∑i=1nAif(xi)

On the example above, we have an integral that we want to approximate and on the right side of the equal sign is the Gaussian Quadrature.

f(x) is our function (it’s not required to know how it looks like), w(x) is a weight function.

As you noticed we’ll be working within the interval [-1;1] - and we’ll consider w(x) = 1 (Legendre).
- But…what if the integral has a different interval?; use this little trick:

∫baf(x)dx=b−a2∫1−1f(b−a2x+b+a2)dx

Our purpose is to find a set of coefficients (A₁, A₂, … A_n) and a set of points (x₁, x₂, … x_n) that will make this quadrature as exact as possible, for a set of integrals.

To keep everything simple we’ll take n=2 so our equation will look like this:

∫1−1f(x)w(x)dx=A1f(x1)+A2f(x2)

From now on I’ll be working with this example, however it’s the same method when n has a different value.

How to find the ‘x’ values

In order to calculate the coefficients (A₁, A₂) we need to know the values of x₁ and x₂. If these are not given, we will have to calculate them using some theory implying orthogonal polynomials.

Now here's how you do this:

Create a polynomial (usually noted π), that has n=2 roots: x₁ and x₂.
It will look like this:

π=∏i=1n=2x−xi⇒π=(x−x1)(x−x2)

Here you need to know this little part of theory related to orthogonal polynomials:

1. polynomials P and Q are orthogonal if their inner product is 0 (<p, q>= 0)
2. a polynomial of degree n is orthogonal to any other polynomial of degree lower than n

These properties are required in order to produce a system of equations which will provide us the values for x₁ and x₂.

We define that inner product that I was talking about earlier:

<P,Q>=∫1−1P(x)Q(x)w(x)dx

where P and Q are 2 polynomials and w(x) is the weight function.

Ok enough with the theory, moving back to our example: we had that π polynomial.
Now we pick the first n=2 terms from the monomial basis (1, x, x², x³…).

Notice that those 2 “polynomials” are of degree 0 and 1 - both lower than n=2(the degree of our π polynomial) => these are orthogonal to π. Knowing this we can create the following system of equations:

<π,1>=∫1−1π(x)⋅1⋅1dx=0

<π,x>=∫1−1π(x)⋅x⋅1dx=0

Note that in there I already did the substitution w(x)=1.
By solving the integrals…

23+2x1x2=0

−23(x1+x2)=0

So we get these values:

x1=−3–√3

x2=3–√3

Finding the coefficients

This is the easy part - if you managed to get here then you’re almost done.
To find the values of the coefficients, we’ll use the same Gauss Quadrature, but this time on these functions:

f₁(x) = 1
f₂(x) = x

where 1 and x are the monomials we chose before.

If we rewrite the Gaussian Quadrature for these 2 functions we get:

∫1−11⋅1dx=A1f1(x1)+A2f1(x2)⇒A1+A2=2

∫1−1x⋅1dx=A1f2(x1)+A2f2(x2)⇒A1x1+A2x2=0

(substituted w(x) = 1)

So A₁ = A₂ = 1, hence the final equation:

∫1−1f(x)d=f(−3–√3)+f(3–√3)

Precision

The degree of precision for this example is 2*n-1 = 3. This means if the integrand is a polynomial of degree 3 (or lower) the result will be exact. If not, the result will be an approximation.

Note that this also works for non-polynomial functions; if we integrate cos(x) from -1 to 1

• original result: 2*sin(1) ~ 1.6829
• result approximated with the Gaussian Quadrature: cos(-sqrt(3)/3) + cos(sqrt(3)/3) = 1.6758

Not the best result, but will do - usually adding more points helps gaining precision.