NumCpp  2.10.1
A Templatized Header Only C++ Implementation of the Python NumPy Library
gauss_legendre.hpp
Go to the documentation of this file.
1 #pragma once
33 
34 #include <cmath>
35 #include <functional>
36 #include <vector>
37 
39 #include "NumCpp/Core/Types.hpp"
40 #include "NumCpp/Utils/sqr.hpp"
41 
42 namespace nc::integrate
43 {
44  //============================================================================
45  // Class Description:
49  {
50  public:
51  //============================================================================
52  // Method Description:
57  explicit LegendrePolynomial(const uint32 numIterations) noexcept :
58  numIterations_(numIterations),
59  weight_(numIterations + 1),
60  root_(numIterations + 1)
61  {
62  calculateWeightAndRoot();
63  }
64 
65  //============================================================================
66  // Method Description:
71  [[nodiscard]] const std::vector<double>& getWeight() const noexcept
72  {
73  return weight_;
74  }
75 
76  //============================================================================
77  // Method Description:
82  [[nodiscard]] const std::vector<double>& getRoot() const noexcept
83  {
84  return root_;
85  }
86 
87  private:
88  //============================================================================
89  // Class Description:
92  struct Result
93  {
94  double value{ 0. };
95  double derivative{ 0. };
96 
97  //============================================================================
98  // Method Description:
104  Result(const double val, const double deriv) noexcept :
105  value(val),
106  derivative(deriv)
107  {
108  }
109  };
110 
111  //============================================================================
112  // Method Description:
115  void calculateWeightAndRoot() noexcept
116  {
117  const auto numIterationsDouble = static_cast<double>(numIterations_);
118  for (uint32 step = 0; step <= numIterations_; ++step)
119  {
120  double root =
121  std::cos(constants::pi * (static_cast<double>(step) - 0.25) / (numIterationsDouble + 0.5));
122  Result result = calculatePolynomialValueAndDerivative(root);
123 
124  double newtonRaphsonRatio;
125  do
126  {
127  newtonRaphsonRatio = result.value / result.derivative;
128  root -= newtonRaphsonRatio;
129  result = calculatePolynomialValueAndDerivative(root);
130  } while (std::fabs(newtonRaphsonRatio) > EPSILON);
131 
132  root_[step] = root;
133  weight_[step] = 2. / ((1. - utils::sqr(root)) * result.derivative * result.derivative);
134  }
135  }
136 
137  //============================================================================
138  // Method Description:
144  Result calculatePolynomialValueAndDerivative(const double x) noexcept
145  {
146  Result result(x, 0.);
147 
148  double value_minus_1 = 1.;
149  const double f = 1. / (utils::sqr(x) - 1.);
150  for (uint32 step = 2; step <= numIterations_; ++step)
151  {
152  const auto stepDouble = static_cast<double>(step);
153  const double value =
154  ((2. * stepDouble - 1.) * x * result.value - (stepDouble - 1.) * value_minus_1) / stepDouble;
155  result.derivative = stepDouble * f * (x * value - result.value);
156 
157  value_minus_1 = result.value;
158  result.value = value;
159  }
160 
161  return result;
162  }
163 
164  //===================================Attributes==============================
165  const double EPSILON{ 1e-15 };
166 
167  const uint32 numIterations_{};
168  std::vector<double> weight_{};
169  std::vector<double> root_{};
170  };
171 
172  //============================================================================
173  // Method Description:
183  inline double
184  gauss_legendre(const double low, const double high, const uint32 n, const std::function<double(double)>& f)
185  {
186  const LegendrePolynomial legendrePolynomial(n);
187  const std::vector<double>& weight = legendrePolynomial.getWeight();
188  const std::vector<double>& root = legendrePolynomial.getRoot();
189 
190  const double width = 0.5 * (high - low);
191  const double mean = 0.5 * (low + high);
192 
193  double gaussLegendre = 0.;
194  for (uint32 step = 1; step <= n; ++step)
195  {
196  gaussLegendre += weight[step] * f(width * root[step] + mean);
197  }
198 
199  return gaussLegendre * width;
200  }
201 } // namespace nc::integrate
Definition: gauss_legendre.hpp:49
LegendrePolynomial(const uint32 numIterations) noexcept
Definition: gauss_legendre.hpp:57
const std::vector< double > & getRoot() const noexcept
Definition: gauss_legendre.hpp:82
const std::vector< double > & getWeight() const noexcept
Definition: gauss_legendre.hpp:71
constexpr double pi
Pi.
Definition: Constants.hpp:39
constexpr double e
eulers number
Definition: Constants.hpp:37
Definition: gauss_legendre.hpp:43
double gauss_legendre(const double low, const double high, const uint32 n, const std::function< double(double)> &f)
Definition: gauss_legendre.hpp:184
dtype f(GeneratorType &generator, dtype inDofN, dtype inDofD)
Definition: f.hpp:56
constexpr dtype sqr(dtype inValue) noexcept
Definition: sqr.hpp:42
auto cos(dtype inValue) noexcept
Definition: cos.hpp:49
NdArray< double > mean(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: mean.hpp:52
std::uint32_t uint32
Definition: Types.hpp:40