gauss_legendre.hpp

Go to the documentation of this file.

#pragma once
 
#include <cmath>
#include <functional>
#include <vector>
 
#include "NumCpp/Core/Constants.hpp"
#include "NumCpp/Core/Types.hpp"
#include "NumCpp/Utils/sqr.hpp"
 
namespace nc::integrate
{
    //============================================================================
    // Class Description:
    class LegendrePolynomial
    {
    public:
        //============================================================================
        // Method Description:
        explicit LegendrePolynomial(const uint32 numIterations) noexcept :
            numIterations_(numIterations),
            weight_(numIterations + 1),
            root_(numIterations + 1)
        {
            calculateWeightAndRoot();
        }
 
        //============================================================================
        // Method Description:
        [[nodiscard]] const std::vector<double>& getWeight() const noexcept
        {
            return weight_;
        }
 
        //============================================================================
        // Method Description:
        [[nodiscard]] const std::vector<double>& getRoot() const noexcept
        {
            return root_;
        }
 
    private:
        //============================================================================
        // Class Description:
        struct Result
        {
            double value{ 0. };
            double derivative{ 0. };
 
            //============================================================================
            // Method Description:
            Result(const double val, const double deriv) noexcept :
                value(val),
                derivative(deriv)
            {
            }
        };
 
        //============================================================================
        // Method Description:
        void calculateWeightAndRoot() noexcept
        {
            const auto numIterationsDouble = static_cast<double>(numIterations_);
            for (uint32 step = 0; step <= numIterations_; ++step)
            {
                double root =
                    std::cos(constants::pi * (static_cast<double>(step) - 0.25) / (numIterationsDouble + 0.5));
                Result result = calculatePolynomialValueAndDerivative(root);
 
                double newtonRaphsonRatio;
                do
                {
                    newtonRaphsonRatio = result.value / result.derivative;
                    root -= newtonRaphsonRatio;
                    result = calculatePolynomialValueAndDerivative(root);
                } while (std::fabs(newtonRaphsonRatio) > EPSILON);
 
                root_[step]   = root;
                weight_[step] = 2. / ((1. - utils::sqr(root)) * result.derivative * result.derivative);
            }
        }
 
        //============================================================================
        // Method Description:
        Result calculatePolynomialValueAndDerivative(const double x) noexcept
        {
            Result result(x, 0.);
 
            double       value_minus_1 = 1.;
            const double f             = 1. / (utils::sqr(x) - 1.);
            for (uint32 step = 2; step <= numIterations_; ++step)
            {
                const auto   stepDouble = static_cast<double>(step);
                const double value =
                    ((2. * stepDouble - 1.) * x * result.value - (stepDouble - 1.) * value_minus_1) / stepDouble;
                result.derivative = stepDouble * f * (x * value - result.value);
 
                value_minus_1 = result.value;
                result.value  = value;
            }
 
            return result;
        }
 
        //===================================Attributes==============================
        const double EPSILON{ 1e-15 };
 
        const uint32        numIterations_{};
        std::vector<double> weight_{};
        std::vector<double> root_{};
    };
 
    //============================================================================
    // Method Description:
    inline double
        gauss_legendre(const double low, const double high, const uint32 n, const std::function<double(double)>& f)
    {
        const LegendrePolynomial   legendrePolynomial(n);
        const std::vector<double>& weight = legendrePolynomial.getWeight();
        const std::vector<double>& root   = legendrePolynomial.getRoot();
 
        const double width = 0.5 * (high - low);
        const double mean  = 0.5 * (low + high);
 
        double gaussLegendre = 0.;
        for (uint32 step = 1; step <= n; ++step)
        {
            gaussLegendre += weight[step] * f(width * root[step] + mean);
        }
 
        return gaussLegendre * width;
    }
} // namespace nc::integrate