quadrature.hpp

Go to the documentation of this file.

 
#pragma once
 
namespace teqp{
 
template<int N, typename T, typename Double=double>
inline auto quad(const std::function<T(Double)>& F, const Double& a, const Double& b){
    
    // Locations x in [-1,1] where the function is to be evaluated, and the corresponding weight
    static const std::map<int, std::tuple<std::vector<Double>, std::vector<Double>>> xw = {
        {
            3,
            {{ 0, sqrt(3.0/5), -sqrt(3.0/5) },
            {8.0/9.0, 5.0/9.0, 5.0/9.0}}
        },
        {
            4,
            {{sqrt(3.0/7.0-2.0/7.0*sqrt(6.0/5)), -sqrt(3.0/7.0-2.0/7.0*sqrt(6.0/5)), sqrt(3.0/7.0+2.0/7.0*sqrt(6.0/5)), -sqrt(3.0/7.0+2.0/7.0*sqrt(6.0/5))},
            {(18.0+sqrt(30.0))/36.0, (18.0+sqrt(30.0))/36.0, (18.0-sqrt(30.0))/36.0, (18.0-sqrt(30.0))/36.0}}
        },
        {
            5,
            {{0, 1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)), -1.0/3.0*sqrt(5-2*sqrt(10.0/7.0)), 1.0/3.0*sqrt(5+2*sqrt(10.0/7.0)), -1.0/3.0*sqrt(5+2*sqrt(10.0/7.0))},
            {128/225.0, (322.0+13*sqrt(70))/900.0, (322.0+13*sqrt(70))/900.0, (322.0-13*sqrt(70))/900.0, (322.0-13*sqrt(70))/900.0}}
        },
        {
            7,
            {{0.0000000000000000,0.4058451513773972,-0.4058451513773972,-0.7415311855993945,0.7415311855993945,-0.9491079123427585,0.9491079123427585},
            { 0.4179591836734694,0.3818300505051189,0.3818300505051189,0.2797053914892766,0.2797053914892766,0.1294849661688697,0.1294849661688697}}
        },
        {
            10,
            {{-0.1488743389816312,0.1488743389816312,-0.4333953941292472,0.4333953941292472,-0.6794095682990244,0.6794095682990244,-0.8650633666889845,0.8650633666889845,-0.9739065285171717,0.9739065285171717},
            { 0.2955242247147529,0.2955242247147529,0.2692667193099963,0.2692667193099963,0.2190863625159820,0.2190863625159820,0.1494513491505806,0.1494513491505806,0.0666713443086881,0.0666713443086881}}
        },
        {
            15,
            {{0.0000000000000000,-0.2011940939974345,0.2011940939974345,-0.3941513470775634,0.3941513470775634,-0.5709721726085388,0.5709721726085388,-0.7244177313601701,0.7244177313601701,-0.8482065834104272,0.8482065834104272,-0.9372733924007060,0.9372733924007060,-0.9879925180204854,0.9879925180204854},
            {0.2025782419255613,0.1984314853271116,0.1984314853271116,0.1861610000155622,0.1861610000155622,0.1662692058169939,0.1662692058169939,0.1395706779261543,0.1395706779261543,0.1071592204671719,0.1071592204671719,0.0703660474881081,0.0703660474881081,0.0307532419961173,0.0307532419961173}}
        },
        {
            30,
            {{-0.0514718425553177,0.0514718425553177,-0.1538699136085835,0.1538699136085835,-0.2546369261678899,0.2546369261678899,-0.3527047255308781,0.3527047255308781,-0.4470337695380892,0.4470337695380892,-0.5366241481420199,0.5366241481420199,-0.6205261829892429,0.6205261829892429,-0.6978504947933158,0.6978504947933158,-0.7677774321048262,0.7677774321048262,-0.8295657623827684,0.8295657623827684,-0.8825605357920527,0.8825605357920527,-0.9262000474292743,0.9262000474292743,-0.9600218649683075,0.9600218649683075,-0.9836681232797472,0.9836681232797472,-0.9968934840746495,0.9968934840746495
            },
            {0.1028526528935588,0.1028526528935588,0.1017623897484055,0.1017623897484055,0.0995934205867953,0.0995934205867953,0.0963687371746443,0.0963687371746443,0.0921225222377861,0.0921225222377861,0.0868997872010830,0.0868997872010830,0.0807558952294202,0.0807558952294202,0.0737559747377052,0.0737559747377052,0.0659742298821805,0.0659742298821805,0.0574931562176191,0.0574931562176191,0.0484026728305941,0.0484026728305941,0.0387991925696271,0.0387991925696271,0.0287847078833234,0.0287847078833234,0.0184664683110910,0.0184664683110910,0.0079681924961666,0.0079681924961666}}
        }
    };
    
    T summer = 0.0;
    // Lookup the coefficients without making a copy or re-allocating
    const std::tuple<std::vector<Double>, std::vector<Double>> &pair = xw.at(N);
    const std::vector<Double> &x = std::get<0>(pair);
    const std::vector<Double> &w = std::get<1>(pair);
    for (auto i = 0; i < N; ++i){
        Double arg = (b-a)/2.0*x[i] + (a+b)/2.0;
        summer += w[i]*F(arg);
    }
    T retval = (b-a)/2.0*summer; // Forces a flattening if T is an autodiff type
    return retval;
}
}