NumCpp  2.8.0
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
43 {
44  namespace integrate
45  {
46  //============================================================================
47  // Class Description:
51  {
52  public:
53  //============================================================================
54  // Method Description:
59  explicit LegendrePolynomial(const uint32 numIterations) noexcept :
60  numIterations_(numIterations),
61  weight_(numIterations + 1),
62  root_(numIterations + 1)
63  {
64  calculateWeightAndRoot();
65  }
66 
67  //============================================================================
68  // Method Description:
73  const std::vector<double>& getWeight() const noexcept
74  {
75  return weight_;
76  }
77 
78  //============================================================================
79  // Method Description:
84  const std::vector<double>& getRoot() const noexcept
85  {
86  return root_;
87  }
88 
89  private:
90  //============================================================================
91  // Class Description:
94  struct Result
95  {
96  double value{ 0.0 };
97  double derivative{ 0.0 };
98 
99  //============================================================================
100  // Method Description:
106  Result(const double val, const double deriv) noexcept :
107  value(val),
108  derivative(deriv)
109  {
110  }
111  };
112 
113  //============================================================================
114  // Method Description:
117  void calculateWeightAndRoot() noexcept
118  {
119  const auto numIterationsDouble = static_cast<double>(numIterations_);
120  for (uint32 step = 0; step <= numIterations_; ++step)
121  {
122  double root =
123  std::cos(constants::pi * (static_cast<double>(step) - 0.25) / (numIterationsDouble + 0.5));
124  Result result = calculatePolynomialValueAndDerivative(root);
125 
126  double newtonRaphsonRatio;
127  do
128  {
129  newtonRaphsonRatio = result.value / result.derivative;
130  root -= newtonRaphsonRatio;
131  result = calculatePolynomialValueAndDerivative(root);
132  } while (std::fabs(newtonRaphsonRatio) > EPSILON);
133 
134  root_[step] = root;
135  weight_[step] = 2.0 / ((1.0 - utils::sqr(root)) * result.derivative * result.derivative);
136  }
137  }
138 
139  //============================================================================
140  // Method Description:
146  Result calculatePolynomialValueAndDerivative(const double x) noexcept
147  {
148  Result result(x, 0.0);
149 
150  double value_minus_1 = 1.0;
151  const double f = 1.0 / (utils::sqr(x) - 1.0);
152  for (uint32 step = 2; step <= numIterations_; ++step)
153  {
154  const auto stepDouble = static_cast<double>(step);
155  const double value =
156  ((2.0 * stepDouble - 1.0) * x * result.value - (stepDouble - 1.0) * value_minus_1) / stepDouble;
157  result.derivative = stepDouble * f * (x * value - result.value);
158 
159  value_minus_1 = result.value;
160  result.value = value;
161  }
162 
163  return result;
164  }
165 
166  //===================================Attributes==============================
167  const double EPSILON{ 1e-15 };
168 
169  const uint32 numIterations_;
170  std::vector<double> weight_;
171  std::vector<double> root_;
172  };
173 
174  //============================================================================
175  // Method Description:
185  inline double
186  gauss_legendre(const double low, const double high, const uint32 n, const std::function<double(double)>& f)
187  {
188  const LegendrePolynomial legendrePolynomial(n);
189  const std::vector<double>& weight = legendrePolynomial.getWeight();
190  const std::vector<double>& root = legendrePolynomial.getRoot();
191 
192  const double width = 0.5 * (high - low);
193  const double mean = 0.5 * (low + high);
194 
195  double gaussLegendre = 0.0;
196  for (uint32 step = 1; step <= n; ++step)
197  {
198  gaussLegendre += weight[step] * f(width * root[step] + mean);
199  }
200 
201  return gaussLegendre * width;
202  }
203  } // namespace integrate
204 } // namespace nc
Definition: gauss_legendre.hpp:51
LegendrePolynomial(const uint32 numIterations) noexcept
Definition: gauss_legendre.hpp:59
const std::vector< double > & getRoot() const noexcept
Definition: gauss_legendre.hpp:84
const std::vector< double > & getWeight() const noexcept
Definition: gauss_legendre.hpp:73
constexpr double pi
Pi.
Definition: Constants.hpp:43
constexpr double e
eulers number
Definition: Constants.hpp:41
double gauss_legendre(const double low, const double high, const uint32 n, const std::function< double(double)> &f)
Definition: gauss_legendre.hpp:186
dtype f(GeneratorType &generator, dtype inDofN, dtype inDofD)
Definition: f.hpp:58
constexpr dtype sqr(dtype inValue) noexcept
Definition: sqr.hpp:44
Definition: Coordinate.hpp:45
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