quadrature.cpp

Go to the documentation of this file.

/// @file quadrature.cpp
#include "analysis/quadrature.hpp"
#include "analysis/types.hpp"
#include <cmath>
#include <stdexcept>
#include <vector>
 
#ifdef NUMERICS_HAS_OMP
#include <omp.h>
#endif
 
namespace num {
 
real trapz(ScalarFn f, real a, real b, idx n, Backend backend) {
    real h   = (b - a) / n;
    real sum = 0.0;
#ifdef NUMERICS_HAS_OMP
    #pragma omp parallel for reduction(+:sum) schedule(static) if(backend == Backend::omp)
#endif
    for (idx i = 1; i < n; ++i)
        sum += f(a + i * h);
    return h * (0.5 * (f(a) + f(b)) + sum);
}
 
real simpson(ScalarFn f, real a, real b, idx n, Backend backend) {
    if (n % 2 != 0)
        throw std::invalid_argument("simpson: n must be even");
    real h   = (b - a) / n;
    real sum = f(a) + f(b);
#ifdef NUMERICS_HAS_OMP
    #pragma omp parallel for reduction(+:sum) schedule(static) if(backend == Backend::omp)
#endif
    for (idx i = 1; i < n; ++i)
        sum += f(a + i * h) * (i % 2 == 0 ? 2.0 : 4.0);
    return h / 3.0 * sum;
}
 
// Gauss-Legendre nodes and weights on [-1, 1] for p = 1..5
// Source: Abramowitz & Stegun, Table 25.4
static constexpr real GL_NODES[5][5] = {
    {0.0,                  0, 0, 0, 0},
    {-0.5773502691896257,  0.5773502691896257, 0, 0, 0},
    {-0.7745966692414834,  0.0, 0.7745966692414834, 0, 0},
    {-0.8611363115940526, -0.3399810435848563,
      0.3399810435848563,  0.8611363115940526, 0},
    {-0.9061798459386640, -0.5384693101056831,  0.0,
      0.5384693101056831,  0.9061798459386640}
};
static constexpr real GL_WEIGHTS[5][5] = {
    {2.0,                  0, 0, 0, 0},
    {1.0,                  1.0, 0, 0, 0},
    {0.5555555555555556,   0.8888888888888889, 0.5555555555555556, 0, 0},
    {0.3478548451374538,   0.6521451548625461,
     0.6521451548625461,   0.3478548451374538, 0},
    {0.2369268850561891,   0.4786286704993665, 0.5688888888888889,
     0.4786286704993665,   0.2369268850561891}
};
 
real gauss_legendre(ScalarFn f, real a, real b, idx p) {
    if (p < 1 || p > 5)
        throw std::invalid_argument("gauss_legendre: p must be 1..5");
    real mid  = 0.5 * (a + b);
    real half = 0.5 * (b - a);
    real sum  = 0.0;
    for (idx i = 0; i < p; ++i)
        sum += GL_WEIGHTS[p - 1][i] * f(mid + half * GL_NODES[p - 1][i]);
    return half * sum;
}
 
// Recursive adaptive Simpson helper
static real adaptive_helper(ScalarFn f, real a, real b,
                             real fa, real fm, real fb,
                             real whole, real tol, idx depth) {
    real mid_l = 0.5 * (a + b * 0.5 + a * 0.5);  // midpoint of [a, mid]
    real mid_r = 0.5 * (0.5 * (a + b) + b);        // midpoint of [mid, b]
    real mid   = 0.5 * (a + b);
    real fl    = f(0.5 * (a + mid));
    real fr    = f(0.5 * (mid + b));
    (void)mid_l; (void)mid_r;
 
    real left  = (b - a) / 12.0 * (fa + 4.0 * fl + fm);
    real right = (b - a) / 12.0 * (fm + 4.0 * fr + fb);
    real delta = left + right - whole;
 
    if (depth == 0 || std::abs(delta) <= 15.0 * tol)
        return left + right + delta / 15.0;
 
    return adaptive_helper(f, a,   mid, fa, fl, fm, left,  tol * 0.5, depth - 1)
         + adaptive_helper(f, mid, b,   fm, fr, fb, right, tol * 0.5, depth - 1);
}
 
real adaptive_simpson(ScalarFn f, real a, real b, real tol, idx max_depth) {
    real fa  = f(a), fb = f(b), fm = f(0.5 * (a + b));
    real est = (b - a) / 6.0 * (fa + 4.0 * fm + fb);
    return adaptive_helper(f, a, b, fa, fm, fb, est, tol, max_depth);
}
 
real romberg(ScalarFn f, real a, real b, real tol, idx max_levels) {
    std::vector<std::vector<real>> R(max_levels, std::vector<real>(max_levels, 0.0));
    R[0][0] = 0.5 * (b - a) * (f(a) + f(b));
 
    for (idx i = 1; i < max_levels; ++i) {
        idx   n   = idx(1) << i;  // 2^i panels
        real  h   = (b - a) / n;
        real  sum = 0.0;
        for (idx k = 1; k < n; k += 2)
            sum += f(a + k * h);
        R[i][0] = 0.5 * R[i - 1][0] + h * sum;
 
        // Richardson extrapolation
        real factor = 1.0;
        for (idx j = 1; j <= i; ++j) {
            factor *= 4.0;
            R[i][j] = R[i][j - 1] + (R[i][j - 1] - R[i - 1][j - 1]) / (factor - 1.0);
        }
 
        if (i > 0 && std::abs(R[i][i] - R[i - 1][i - 1]) < tol)
            return R[i][i];
    }
    return R[max_levels - 1][max_levels - 1];
}
 
} // namespace num