

In a general Gaussian quadrature rule, an definite integral of 



f
(
x
)


{\displaystyle f(x)}

 is first approximated over the interval 



[
−
1
,
1
]


{\displaystyle [-1,1]}

 by a polynomial approximable function 



g
(
x
)


{\displaystyle g(x)}

 and a known weighting function 



W
(
x
)


{\displaystyle W(x)}

.






∫

−
1


1


f
(
x
)

d
x
=

∫

−
1


1


W
(
x
)
g
(
x
)

d
x


{\displaystyle \int _{-1}^{1}f(x)\,dx=\int _{-1}^{1}W(x)g(x)\,dx}




Those are then approximated by a sum of function values at specified points 




x

i




{\displaystyle x_{i}}

 multiplied by some weights 




w

i




{\displaystyle w_{i}}

:






∫

−
1


1


W
(
x
)
g
(
x
)

d
x
≈

∑

i
=
1


n



w

i


g
(

x

i


)


{\displaystyle \int _{-1}^{1}W(x)g(x)\,dx\approx \sum _{i=1}^{n}w_{i}g(x_{i})}




In the case of Gauss-Legendre quadrature, the weighting function 



W
(
x
)
=
1


{\displaystyle W(x)=1}

, so we can approximate an integral of 



f
(
x
)


{\displaystyle f(x)}

 with:






∫

−
1


1


f
(
x
)

d
x
≈

∑

i
=
1


n



w

i


f
(

x

i


)


{\displaystyle \int _{-1}^{1}f(x)\,dx\approx \sum _{i=1}^{n}w_{i}f(x_{i})}





For this, we first need to calculate the nodes and the weights, but after we have them, we can reuse them for numerious integral evaluations, which greatly speeds up the calculation compared to more simple numerical integration methods.



The 



n


{\displaystyle n}

 evaluation points 




x

i




{\displaystyle x_{i}}

 for a n-point rule, also called "nodes", are roots of n-th order Legendre Polynomials 




P

n


(
x
)


{\displaystyle P_{n}(x)}

. Legendre polynomials are defined by the following recursive rule:






P

0


(
x
)
=
1


{\displaystyle P_{0}(x)=1}







P

1


(
x
)
=
x


{\displaystyle P_{1}(x)=x}






n

P

n


(
x
)
=
(
2
n
−
1
)
x

P

n
−
1


(
x
)
−
(
n
−
1
)

P

n
−
2


(
x
)


{\displaystyle nP_{n}(x)=(2n-1)xP_{n-1}(x)-(n-1)P_{n-2}(x)}





There is also a recursive equation for their derivative:






P

n

′

(
x
)
=


n


x

2


−
1




(
x

P

n


(
x
)
−

P

n
−
1


(
x
)
)



{\displaystyle P_{n}'(x)={\frac {n}{x^{2}-1}}\left(xP_{n}(x)-P_{n-1}(x)\right)}




The roots of those polynomials are in general not analytically solvable, so they have to be approximated numerically, for example by Newton-Raphson iteration:






x

n
+
1


=

x

n


−



f
(

x

n


)



f
′

(

x

n


)





{\displaystyle x_{n+1}=x_{n}-{\frac {f(x_{n})}{f'(x_{n})}}}




The first guess 




x

0




{\displaystyle x_{0}}

 for the 



i


{\displaystyle i}

-th root of a 



n


{\displaystyle n}

-order polynomial 




P

n




{\displaystyle P_{n}}

 can be given by






x

0


=
cos
⁡

(
π




i
−


1
4




n
+


1
2





)



{\displaystyle x_{0}=\cos \left(\pi \,{\frac {i-{\frac {1}{4}}}{n+{\frac {1}{2}}}}\right)}




After we get the nodes 




x

i




{\displaystyle x_{i}}

, we compute the appropriate weights by:






w

i


=


2


(
1
−

x

i


2


)

[

P

n

′

(

x

i


)

]

2







{\displaystyle w_{i}={\frac {2}{\left(1-x_{i}^{2}\right)[P'_{n}(x_{i})]^{2}}}}




After we have the nodes and the weights for a n-point quadrature rule, we can approximate an integral over any interval 



[
a
,
b
]


{\displaystyle [a,b]}

 by






∫

a


b


f
(
x
)

d
x
≈



b
−
a

2



∑

i
=
1


n



w

i


f

(



b
−
a

2



x

i


+



a
+
b

2


)



{\displaystyle \int _{a}^{b}f(x)\,dx\approx {\frac {b-a}{2}}\sum _{i=1}^{n}w_{i}f\left({\frac {b-a}{2}}x_{i}+{\frac {a+b}{2}}\right)}




Task description

Similar to the task Numerical Integration, the task here is to calculate the definite integral of a function 



f
(
x
)


{\displaystyle f(x)}

, but by applying an n-point Gauss-Legendre quadrature rule, as described here, for example. The input values should be an function f to integrate, the bounds of the integration interval a and b, and the number of gaussian evaluation points n. An reference implementation in Common Lisp is provided for comparison.

To demonstrate the calculation, compute the weights and nodes for an 5-point quadrature rule and then use them to compute:

         
  
    
      
        
          ∫
          
            −
            3
          
          
            3
          
        
        exp
        ⁡
        (
        x
        )
        
        d
        x
        ≈
        
          ∑
          
            i
            =
            1
          
          
            5
          
        
        
          w
          
            i
          
        
        
        exp
        ⁡
        (
        
          x
          
            i
          
        
        )
        ≈
        20.036
      
    
    {\displaystyle \int _{-3}^{3}\exp(x)\,dx\approx \sum _{i=1}^{5}w_{i}\;\exp(x_{i})\approx 20.036}
  




#include <iostream>
#include <iomanip>
#include <cmath>
 
namespace Rosetta {
 
    /*! Implementation of Gauss-Legendre quadrature
    *  http://en.wikipedia.org/wiki/Gaussian_quadrature
    *  http://rosettacode.org/wiki/Numerical_integration/Gauss-Legendre_Quadrature
    * 
    */
    template <int N>
    class GaussLegendreQuadrature {
    public:
        enum {eDEGREE = N};
 
        /*! Compute the integral of a functor
        *
        *   @param a    lower limit of integration
        *   @param b    upper limit of integration
        *   @param f    the function to integrate
        *   @param err  callback in case of problems
        */
        template <typename Function>
        double integrate(double a, double b, Function f) {
            double p = (b - a) / 2;
            double q = (b + a) / 2;
            const LegendrePolynomial& legpoly = s_LegendrePolynomial;
 
            double sum = 0;
            for (int i = 1; i <= eDEGREE; ++i) {
                sum += legpoly.weight(i) * f(p * legpoly.root(i) + q);
            }
 
            return p * sum;
        }
 
        /*! Print out roots and weights for information
        */
        void print_roots_and_weights(std::ostream& out) const {
            const LegendrePolynomial& legpoly = s_LegendrePolynomial;
            out << "Roots:  ";
            for (int i = 0; i <= eDEGREE; ++i) {
                out << ' ' << legpoly.root(i);
            }
            out << std::endl;
            out << "Weights:";
            for (int i = 0; i <= eDEGREE; ++i) {
                out << ' ' << legpoly.weight(i);
            }
            out << std::endl;
        }
    private:
        /*! Implementation of the Legendre polynomials that form
        *   the basis of this quadrature
        */
        class LegendrePolynomial {
        public:
            LegendrePolynomial () {
                // Solve roots and weights
                for (int i = 0; i <= eDEGREE; ++i) {
                    double dr = 1;
 
                    // Find zero
                    Evaluation eval(cos(M_PI * (i - 0.25) / (eDEGREE + 0.5)));
                    do {
                        dr = eval.v() / eval.d();
                        eval.evaluate(eval.x() - dr);
                    } while (fabs (dr) > 2e-16);
 
                    this->_r[i] = eval.x();
                    this->_w[i] = 2 / ((1 - eval.x() * eval.x()) * eval.d() * eval.d());
                }
            }
 
            double root(int i) const { return this->_r[i]; }
            double weight(int i) const { return this->_w[i]; }
        private:
            double _r[eDEGREE + 1];
            double _w[eDEGREE + 1];
 
            /*! Evaluate the value *and* derivative of the
            *   Legendre polynomial
            */
            class Evaluation {
            public:
                explicit Evaluation (double x) : _x(x), _v(1), _d(0) {
                    this->evaluate(x);
                }
 
                void evaluate(double x) {
                    this->_x = x;
 
                    double vsub1 = x;
                    double vsub2 = 1;
                    double f     = 1 / (x * x - 1);
 
                    for (int i = 2; i <= eDEGREE; ++i) {
                        this->_v = ((2 * i - 1) * x * vsub1 - (i - 1) * vsub2) / i;
                        this->_d = i * f * (x * this->_v - vsub1);
 
                        vsub2 = vsub1;
                        vsub1 = this->_v;
                    }
                }
 
                double v() const { return this->_v; }
                double d() const { return this->_d; }
                double x() const { return this->_x; }
 
            private:
                double _x;
                double _v;
                double _d;
            };
        };
 
        /*! Pre-compute the weights and abscissae of the Legendre polynomials
        */
        static LegendrePolynomial s_LegendrePolynomial;
    };
 
    template <int N>
    typename GaussLegendreQuadrature<N>::LegendrePolynomial GaussLegendreQuadrature<N>::s_LegendrePolynomial;
}
 
// This to avoid issues with exp being a templated function
double RosettaExp(double x) {
    return exp(x);
}
 
int main() {
    Rosetta::GaussLegendreQuadrature<5> gl5;
 
    std::cout << std::setprecision(10);
 
    gl5.print_roots_and_weights(std::cout);
    std::cout << "Integrating Exp(X) over [-3, 3]: " << gl5.integrate(-3., 3., RosettaExp) << std::endl;
    std::cout << "Actual value:                    " << RosettaExp(3) - RosettaExp(-3) << std::endl; 
}