NumCpp  2.1.0
A C++ implementation of the Python Numpy library
gauss_legendre.hpp
Go to the documentation of this file.
1 #pragma once
34 
36 #include "NumCpp/Core/Types.hpp"
37 #include "NumCpp/Utils/sqr.hpp"
38 
39 #include <cmath>
40 #include <functional>
41 #include <vector>
42 
43 namespace nc
44 {
45  namespace integrate
46  {
47  //============================================================================
48  // Class Description:
52  {
53  public:
54  //============================================================================
55  // Method Description:
60  explicit LegendrePolynomial(const uint32 numIterations) noexcept :
61  numIterations_(numIterations),
62  weight_(numIterations + 1),
63  root_(numIterations + 1)
64  {
65  calculateWeightAndRoot();
66  }
67 
68  //============================================================================
69  // Method Description:
74  const std::vector<double>& getWeight() const noexcept
75  {
76  return weight_;
77  }
78 
79  //============================================================================
80  // Method Description:
85  const std::vector<double>& getRoot() const noexcept
86  {
87  return root_;
88  }
89 
90  private:
91  //============================================================================
92  // Class Description:
95  struct Result
96  {
97  double value{ 0.0 };
98  double derivative{ 0.0 };
99 
100  //============================================================================
101  // Method Description:
107  Result(const double val, const double deriv) noexcept :
108  value(val),
109  derivative(deriv)
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 = std::cos(constants::pi * (static_cast<double>(step) - 0.25) / (numIterationsDouble + 0.5));
123  Result result = calculatePolynomialValueAndDerivative(root);
124 
125  double newtonRaphsonRatio;
126  do
127  {
128  newtonRaphsonRatio = result.value / result.derivative;
129  root -= newtonRaphsonRatio;
130  result = calculatePolynomialValueAndDerivative(root);
131  } while (std::fabs(newtonRaphsonRatio) > EPSILON);
132 
133  root_[step] = root;
134  weight_[step] = 2.0 / ((1.0 - utils::sqr(root)) * result.derivative * result.derivative);
135  }
136  }
137 
138  //============================================================================
139  // Method Description:
145  Result calculatePolynomialValueAndDerivative(const double x) noexcept
146  {
147  Result result(x, 0.0);
148 
149  double value_minus_1 = 1.0;
150  const double f = 1.0 / (utils::sqr(x) - 1.0);
151  for (uint32 step = 2; step <= numIterations_; ++step)
152  {
153  const auto stepDouble = static_cast<double>(step);
154  const double value = ((2.0 * stepDouble - 1.0) * x * result.value - (stepDouble - 1.0) * 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 gauss_legendre(const double low, const double high, const uint32 n,
184  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.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 integrate
202 } // namespace nc
nc::integrate::gauss_legendre
double gauss_legendre(const double low, const double high, const uint32 n, const std::function< double(double)> &f)
Definition: gauss_legendre.hpp:183
nc::integrate::LegendrePolynomial::getWeight
const std::vector< double > & getWeight() const noexcept
Definition: gauss_legendre.hpp:74
Constants.hpp
nc::mean
NdArray< double > mean(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: mean.hpp:54
nc::integrate::LegendrePolynomial::LegendrePolynomial
LegendrePolynomial(const uint32 numIterations) noexcept
Definition: gauss_legendre.hpp:60
nc::cos
auto cos(dtype inValue) noexcept
Definition: cos.hpp:52
nc::constants::e
constexpr double e
eulers number
Definition: Constants.hpp:42
nc::uint32
std::uint32_t uint32
Definition: Types.hpp:41
nc::constants::pi
constexpr double pi
Pi.
Definition: Constants.hpp:44
nc
Definition: Coordinate.hpp:45
sqr.hpp
nc::utils::sqr
constexpr dtype sqr(dtype inValue) noexcept
Definition: sqr.hpp:45
nc::integrate::LegendrePolynomial
Definition: gauss_legendre.hpp:51
Types.hpp
nc::random::f
dtype f(dtype inDofN, dtype inDofD)
Definition: f.hpp:58
nc::integrate::LegendrePolynomial::getRoot
const std::vector< double > & getRoot() const noexcept
Definition: gauss_legendre.hpp:85