quadrature.hpp

Go to the documentation of this file.

// Copyright (C) 2021-2022 Petter Nilsson. MIT License.
 
#pragma once
 
#include <array>
#include <numbers>
#include <utility>
 
#include "basis.hpp"
#include "static_matrix.hpp"
 
SMOOTH_BEGIN_NAMESPACE
 
// \cond
namespace detail {
 
template<std::size_t K>
constexpr auto pow_s(const auto x) noexcept
{
  const auto f = [&x]<auto... Is>(std::index_sequence<Is...>) {
    return ((static_cast<void>(Is), x) * ... * 1.);
  };  // NOLINT
  return f(std::make_index_sequence<K>{});
}
 
template<std::size_t K>
constexpr auto factorial_s() noexcept
{
  const auto f = []<auto... Is>(std::index_sequence<Is...>) { return ((Is + 1) * ... * 1.); };  // NOLINT
  return f(std::make_index_sequence<K>{});
}
 
constexpr auto cos_s(const auto x) noexcept
{
  const auto f = [&x]<auto... Is>(std::index_sequence<Is...>) {
    return ((pow_s<Is>(-1.) * pow_s<2 * Is>(x) / factorial_s<2 * Is>()) + ...);
  };  // NOLINT
  return f(std::make_index_sequence<8>{});
}
 
}  // namespace detail
// \endcond
 
template<std::size_t K>
constexpr std::array<double, K> cgr_nodes()
{
  std::array<double, K> x;
  for (auto i = 0u; i < K; ++i) { x[i] = -detail::cos_s(2. * std::numbers::pi_v<double> * i / (2 * K - 1)); }
  return x;
}
 
template<std::size_t K, std::size_t I = 8>
  requires(K <= 40)
constexpr std::pair<std::array<double, K>, std::array<double, K>> lgr_nodes()
{
  // initial guess
  auto xs = cgr_nodes<K>();
 
  // initialize weights
  std::array<double, K> ws;
  ws[0] = 2. / (K * K);
 
  // two rightmost column of legendre basis matrix
  constexpr auto B = polynomial_basis<PolynomialBasis::Legendre, K>().template block<K + 1, 2>(0, K - 1);
 
  for (auto i = 1u; i < K; ++i) {
    StaticMatrix<double, 1, 2> U_B;
    for (auto iter = 0u; iter < I; ++iter) {
      U_B               = monomial_derivative<K>(xs[i], 0) * B;
      const auto dU_B   = monomial_derivative<K>(xs[i], 1) * B;
      const double fxi  = U_B[0][0] + U_B[0][1];
      const double dfxi = dU_B[0][0] + dU_B[0][1];
      xs[i] -= fxi / dfxi;
    }
    ws[i] = (1. - xs[i]) / (K * K * U_B[0][0] * U_B[0][0]);
  }
 
  return {xs, ws};
}
 
SMOOTH_END_NAMESPACE