Welcome back to a new post at thoughts-on-cpp.com. In this post, I would like to discuss the Gauss integration algorithm, more precisely the Gauss-Legendre integration algorithm. The Gauss-Legendre integration is the most known form of the Gauss integrations. Others are

• Gauss-Tschebyschow
• Gauss-Hermite
• Gauss-Laguerre
• Gauss-Lobatto
• Gauss-Kronrod

The idea of the Gauss integration algorithm is to approximate, similar to the Simpson Rule, the function f(x) by $f(x)=w(x) \cdot \phi (x)$

While w(x) is a weighting function, $\phi(x)$ is a polynomial function (Legendre-Polynomials) with defined nodes $x_i$ which can be exactly integrated. A general form for a range of a-b looks like the following. $\int_{a}^{b} f(x) \mathrm{d} x \approx \frac{b-a}{2} \sum_{i=1}^{n} f\left(\frac{b-a}{2} x_{i}+\frac{a+b}{2}\right) w_{i}$

The Legendre-Polynomials are defined by the general formula and its derivative $P_n(x)=\sum_{k=0}^{n/2}=(-1)^k\frac{(2n-2k)!}{(n-k)!(n-2k)!k!2^n}x^{n-2k}$ $P_{n}^{\prime}(x)=\frac{n}{x^{2}-1}\left(x P_{n}(x)-P_{n-1}(x)\right)$

The following image is showing the 3rd until the 7th Legendre Polynomials, the 1st and 2nd polynomials are just 1 and x and therefore not necessary to show.

Let’s have a closer look at the source code:

The integral is done by the gaussLegendreIntegral (line 69) function which is initializing the LegendrePolynomial class and afterward solving the integral (line 77 – 80). Something very interesting to note: We need to calculate the Legendre-Polynomials only once and can use them for any function of order n in the range a-b. The Gauss-Legendre integration is therefore extremely fast for all subsequent integrations.

The method calculatePolynomialValueAndDerivative is calculating the value (line 50) at a certain node $x_i$ and its derivative (line 51). Both results are used at method calculateWeightAndRoot to calculate the the node $x_i$ by the Newton-Raphson method (line 33 – 37). $x_{n+1}=x_{n}-\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)}$

The weight w(x)  will be calculated (line 40) by $w_{i}=\frac{2}{\left(1-x_{i}^{2}\right)\left[P_{n}^{\prime}\left(x_{i}\right)\right]^{2}}$

As we can see in the screen capture below, the resulting approximation of $I=\int_0^{\pi/2} \frac{5}{e^\pi-2}\exp(2x)\cos(x)dx=1.0$

is very accurate. We end up with an error of only $-3.8151e^{-6}$. Gauss-Legendre integration works very good for integrating smooth functions and result in higher accuracy with the same number of nodes compared to Newton-Cotes Integration. A drawback of Gauss-Legendre integration might be the performance in case of dynamic integration where the number of nodes are changing.

Did you like the post?

Processing…
Success! You're on the list.

## 2 thoughts on “Numerical Methods with C++ Part 2: Gauss-Legendre Integration”

1. Weili Cheng says:

On line 37, the variable ratio is not defined, which is likely defined as newtonRaphsonRatio on line 32.

Like

1. bmahr says:

Yes your right, it’s the newtonRaphsonRatio. Thanks for finding

Like

This site uses Akismet to reduce spam. Learn how your comment data is processed.