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) f(x)=w(x) \cdot \phi (x)](https://s0.wp.com/latex.php?latex=f%28x%29%3Dw%28x%29+%5Ccdot+%5Cphi+%28x%29++&bg=%23ffffff&fg=%23383838&s=2)
While w(x) is a weighting function,
is a polynomial function (Legendre-Polynomials) with defined nodes
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} \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}](https://s0.wp.com/latex.php?latex=%5Cint_%7Ba%7D%5E%7Bb%7D+f%28x%29+%5Cmathrm%7Bd%7D+x+%5Capprox+%5Cfrac%7Bb-a%7D%7B2%7D+%5Csum_%7Bi%3D1%7D%5E%7Bn%7D+f%5Cleft%28%5Cfrac%7Bb-a%7D%7B2%7D+x_%7Bi%7D%2B%5Cfrac%7Ba%2Bb%7D%7B2%7D%5Cright%29+w_%7Bi%7D+&bg=%23ffffff&fg=%23383838&s=2)
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(x)=\sum_{k=0}^{n/2}=(-1)^k\frac{(2n-2k)!}{(n-k)!(n-2k)!k!2^n}x^{n-2k}](https://s0.wp.com/latex.php?latex=P_n%28x%29%3D%5Csum_%7Bk%3D0%7D%5E%7Bn%2F2%7D%3D%28-1%29%5Ek%5Cfrac%7B%282n-2k%29%21%7D%7B%28n-k%29%21%28n-2k%29%21k%212%5En%7Dx%5E%7Bn-2k%7D++&bg=%23ffffff&fg=%23383838&s=2)
![P_{n}^{\prime}(x)=\frac{n}{x^{2}-1}\left(x P_{n}(x)-P_{n-1}(x)\right) P_{n}^{\prime}(x)=\frac{n}{x^{2}-1}\left(x P_{n}(x)-P_{n-1}(x)\right)](https://s0.wp.com/latex.php?latex=P_%7Bn%7D%5E%7B%5Cprime%7D%28x%29%3D%5Cfrac%7Bn%7D%7Bx%5E%7B2%7D-1%7D%5Cleft%28x+P_%7Bn%7D%28x%29-P_%7Bn-1%7D%28x%29%5Cright%29++&bg=%23ffffff&fg=%23383838&s=2)
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.
![First 5 Legendre Polynomials](https://thoughtsoncpp.files.wordpress.com/2019/04/legendrepolynoms.png?w=800)
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
and its derivative (line 51). Both results are used at method calculateWeightAndRoot
to calculate the the node
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)} x_{n+1}=x_{n}-\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)}](https://s0.wp.com/latex.php?latex=x_%7Bn%2B1%7D%3Dx_%7Bn%7D-%5Cfrac%7Bf%5Cleft%28x_%7Bn%7D%5Cright%29%7D%7Bf%5E%7B%5Cprime%7D%5Cleft%28x_%7Bn%7D%5Cright%29%7D++&bg=%23ffffff&fg=%23383838&s=2)
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}} w_{i}=\frac{2}{\left(1-x_{i}^{2}\right)\left[P_{n}^{\prime}\left(x_{i}\right)\right]^{2}}](https://s0.wp.com/latex.php?latex=w_%7Bi%7D%3D%5Cfrac%7B2%7D%7B%5Cleft%281-x_%7Bi%7D%5E%7B2%7D%5Cright%29%5Cleft%5BP_%7Bn%7D%5E%7B%5Cprime%7D%5Cleft%28x_%7Bi%7D%5Cright%29%5Cright%5D%5E%7B2%7D%7D++&bg=%23ffffff&fg=%23383838&s=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 I=\int_0^{\pi/2} \frac{5}{e^\pi-2}\exp(2x)\cos(x)dx=1.0](https://s0.wp.com/latex.php?latex=I%3D%5Cint_0%5E%7B%5Cpi%2F2%7D+%5Cfrac%7B5%7D%7Be%5E%5Cpi-2%7D%5Cexp%282x%29%5Ccos%28x%29dx%3D1.0++&bg=%23ffffff&fg=%23383838&s=2)
is very accurate. We end up with an error of only
. 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.
![Screen capture of program execution](https://thoughtsoncpp.files.wordpress.com/2019/04/numericalintegration.gif?w=800)
Did you like the post?
What are your thoughts?
Feel free to comment and share this post.
Passionate C++ developer, mechanical engineer, Head of Software Development at KISSsoft AG and host of https://thoughs-on-coding.com