NumCpp  2.12.1
A Templatized Header Only C++ Implementation of the Python NumPy Library
gauss_legendre.hpp
Go to the documentation of this file.
1
32#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
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: Core/Constants.hpp:39
constexpr double e
eulers number
Definition: Core/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
NdArray< double > mean(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: mean.hpp:52
auto cos(dtype inValue) noexcept
Definition: cos.hpp:49
std::uint32_t uint32
Definition: Types.hpp:40