Numerical Methods with C++ Part 2: Gauss-Legendre Integration

Welcome back to a new post at thoughts-on-coding.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.

This diagram shows the first 5 orders of the Legendre Polynomial except order 0.

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

class LegendrePolynomial {
public:
    LegendrePolynomial(double lowerBound, double upperBound, size_t numberOfIterations)
        : mLowerBound(lowerBound), mUpperBound(upperBound), mNumberOfIterations(numberOfIterations), mWeight(numberOfIterations+1), mRoot(numberOfIterations+1) {
        calculateWeightAndRoot();
    }

    const std::vector<double> & getWeight() const {
        return mWeight;
    }

    const std::vector<double> & getRoot() const {
        return mRoot;
    }

private:
    const static double EPSILON;

    struct Result {
        double value;
        double derivative;

        Result() : value(0), derivative(0) {}
        Result(double val, double deriv) : value(val), derivative(deriv) {}
    };

    void calculateWeightAndRoot() {
        for(int step = 0; step <= mNumberOfIterations; step++) {
            double root = cos(M_PI * (step-0.25)/(mNumberOfIterations+0.5));
            Result result = calculatePolynomialValueAndDerivative(root);

            double newtonRaphsonRatio;
            do {
                newtonRaphsonRatio = result.value/result.derivative;
                root -= newtonRaphsonRatio;
                result = calculatePolynomialValueAndDerivative(root);
            } while (fabs(newtonRaphsonRatio) > EPSILON);

            mRoot[step] = root;
            mWeight[step] = 2.0/((1-root*root)*result.derivative*result.derivative);
        }
    }

    Result calculatePolynomialValueAndDerivative(double x) {
        Result result(x, 0);

        double value_minus_1 = 1;
        const double f = 1/(x*x-1);
        for(int step = 2; step <= mNumberOfIterations; step++) {
            const double value = ((2*step-1)*x*result.value-(step-1)*value_minus_1)/step;
            result.derivative = step*f*(x*value - result.value);

            value_minus_1 = result.value;
            result.value = value;
        }

        return result;
    }

    const double mLowerBound;
    const double mUpperBound;
    const int mNumberOfIterations;
    std::vector<double> mWeight;
    std::vector<double> mRoot;
};

const double LegendrePolynomial::EPSILON = 1e-15;

double gaussLegendreIntegral(double a, double b, int n, const std::function<double (double)> &f) {
    const LegendrePolynomial legendrePolynomial(a, b, n);
    const std::vector<double> & weight = legendrePolynomial.getWeight();
    const std::vector<double> & root = legendrePolynomial.getRoot();

    const double width = 0.5*(b-a);
    const double mean = 0.5*(a+b);

    double gaussLegendre = 0;
    for(int step = 1; step <= n; step++) {
        gaussLegendre += weight[step]*f(width * root[step] + mean);
    }

    return gaussLegendre * width;
}

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.

Since you've made it this far, sharing this article on your favorite social media network and giving feedback would be highly appreciated 💖!

Published 25 Apr 2019