[NUMERICAL ANALYSIS] Gauss-Legendre Integration

NUMERICAL INTEGRATION VIA GAUSS-LEGENDRE INTEGRATION

The Gaussian integration technique transforms the function and integrates the function f(x) to another domain by weights and nodes of an interpolating function g(x). Nodes and weights are calculated from the Legendre polynomials and a weighting function.  Basically, if we have the integral:

We can transform x to y via the relation (and of course, we transform dx into dy as well):

The integral then looks something like this:

The integral is equal to 

Which can then be solved via a summation: 

Therefore, we can say that the entire thing is equivalent to:

So, how do we determine the values for w and f? w is actually a set of tabulated values easily available in generic numerical analysis text. Also,  f is a function of y, which allows us to evaluate f given values of y, which, as seen above, are calculated from values of x, which are also found in generic numerical analysis text. 

The values w and x are called weights, and nodes, respectively. The nodes are actually roots of the resulting solutions of the Legendre differential equation between (-1, 1):

The weights are then calculated via:

P_n(x) there is the n^th-order Legendre polynomial. Say for a fifth-order Legendre polynomial, the weights and nodes are:

To generate Gauss-Legendre coefficients in MATLAB, you can paste this code:

%-----

syms x

n = 100

nodes = vpasolve(legendreP(n,x) == 0)

weights = (2 .* (1 - nodes .^ 2)) ./ (n .* legendreP(n - 1,nodes)) .^2

%-----

Anyway, that's enough deriving. So how do we use the method, exactly? Let's take the sample integral,

and use the 5-point Gauss-Legendre integration scheme to calculate the integral. The value of the integral is 844.5754.

First, we compute x_i = (b - a)/2 * y_i + (b + a)/2. 

Compute f(x_i): we then substitute the computed x_i to the function being integrated.

Then, we compute w_i * f(x_i) for each i.

Finally, to get the area, we compute the total of w_i * f(x_i) * dy, which is (b-a)/2 * Σ[w_i * f(x_i)], which is equal to 844.286.

So how do we get more accuracy here? We can use more points in calculating the integral, but I leave you with that 😂

References: 
Kreyszig, Erwin (2011). Advanced Engineering Mathematics, 10th Ed., Wiley

Chapra, S. C., & Canale, R. P. (2006). Numerical methods for engineers. Boston: McGraw-Hill Higher Education.