GSI_Toolbox 1.0.0
A toolbox for Gas-Surface Interaction simulations
Loading...
Searching...
No Matches
Hermite_Tools.h
1
10template <typename T>
11class Hermite {
12
13 private:
14 Matrix<T> roots;
15 Matrix<T> weights;
16 unsigned int max_order = 100;
17
18 public:
19
25 Hermite();
26
30 ~Hermite();
31
39 Matrix<T> evaluate(Matrix<T> coefficients, Matrix<T> values);
40
48 T evaluate(Matrix<T> coefficients, T value);
49
57
65 Matrix<T> evaluate_deriv(Matrix<T> coefficients, Matrix<T> values);
66
74 T evaluate_deriv(Matrix<T> coefficients, T value);
75
82 Matrix<T> import_data(const char* file);
83
90 Matrix<T> get_roots(unsigned int order);
91
98 Matrix<T> get_weights(unsigned int order);
99
107 Matrix<T> compute_coefficients(unsigned int order, T (*func)(T));
108};
109
115template <typename T>
117 this->roots = import_data("code_files/roots.txt");
118 this->weights = import_data("code_files/weights.txt");
119 this->weights /= exp((this->roots ^ 2.0) * (-1.0));
120}
121
125template <typename T>
127
135template <typename T>
137
138 unsigned int order = coefficients.get_rows();
139 Matrix<T> Hn_0(values), Hn_1(values), Hn(values), result(values);
140
141 // Initialize the first two Hermite polynomials
142 Hn_0 = 1.0;
143 Hn_1 = values * 2.0;
144 result = Hn_0 * coefficients(0, 0) + Hn_1 * coefficients(1, 0);
145
146 // Compute higher order Hermite polynomials iteratively
147 for (auto i = 2; i < order; i++) {
148 Hn = Hn_1 * values * 2.0 - Hn_0 * 2.0 * (i - 1.0);
149 Hn_0 = Hn_1;
150 Hn_1 = Hn;
151 result += (Hn * coefficients(i, 0));
152 }
153
154 return result;
155}
156
164template <typename T>
165T Hermite<T>::evaluate(Matrix<T> coefficients, T value) {
166
167 T Hn_0, Hn_1, Hn, result;
168 unsigned int order = coefficients.get_rows();
169 Hn_0 = 1.0;
170 Hn_1 = value * 2.0;
171 result = Hn_0 * coefficients(0, 0) + Hn_1 * coefficients(1, 0);
172
173 // Compute higher order Hermite polynomials iteratively
174 for (auto i = 2; i < order; i++) {
175 Hn = Hn_1 * value * 2.0 - Hn_0 * 2.0 * (i - 1.0);
176 Hn_0 = Hn_1;
177 Hn_1 = Hn;
178 result += Hn * coefficients(i, 0);
179 }
180
181 return result;
182}
183
190template <typename T>
192
193 unsigned int order = coefficients.get_rows();
194 Matrix<T> new_coefficients(order - 1, 1);
195
196 // Compute the derivative of the polynomial coefficients
197 for (auto i = 1; i < order; i++) {
198 new_coefficients(i - 1, 0) = 2.0 * i * coefficients(i, 0);
199 }
200
201 return new_coefficients;
202}
203
214template <typename T>
216 Matrix<T> deriv_coefficients = coefficients_deriv(coefficients);
217 return evaluate(deriv_coefficients, values);
218}
219
230template <typename T>
231T Hermite<T>::evaluate_deriv(Matrix<T> coefficients, T value) {
232 Matrix<T> deriv_coefficients = coefficients_deriv(coefficients);
233 return evaluate(deriv_coefficients, value);
234}
235
245template <typename T>
247 unsigned int rows, cols;
248 T buffer;
249 std::ifstream myfile;
250 myfile.open(file);
251 myfile >> rows >> cols;
252
253 Matrix<T> data(rows, cols);
254
255 // Read data into the matrix
256 for (auto i = 0; i < rows; i++) {
257 for (auto j = 0; j < cols; j++) {
258 myfile >> buffer;
259 data(i, j) = buffer;
260 }
261 }
262 myfile.close();
263 return data;
264}
265
275template <typename T>
276Matrix<T> Hermite<T>::get_roots(unsigned int order) {
277 Matrix<T> roots(order, 1);
278 if (order >= max_order) {
279 std::cerr << "Maximum order exceeded!\n";
280 return roots;
281 } else {
282 for (auto i = 0; i < order; i++) {
283 roots(i, 0) = this->roots(order - 1, i);
284 }
285 return roots;
286 }
287}
288
298template <typename T>
300 Matrix<T> weights(order, 1);
301 if (order >= max_order) {
302 std::cerr << "Maximum order exceeded!\n";
303 return weights;
304 } else {
305 for (auto i = 0; i < order; i++) {
306 weights(i, 0) = this->weights(order - 1, i);
307 }
308 return weights;
309 }
310}
311
323template <typename T>
324Matrix<T> Hermite<T>::compute_coefficients(unsigned int order, T (*func)(T)) {
325 Matrix<T> weights = get_weights(order);
326 Matrix<T> roots = get_roots(order);
327 Matrix<T> func_values(order, 1);
328 Matrix<T> coefficients(order, 1);
329
330 T factorial = 1.0;
331
332 // Evaluate the function at each root
333 for (auto i = 0; i < order; i++) {
334 func_values(i, 0) = *func(roots(i, 0));
335 }
336
337 // Compute the coefficients
338 for (auto i = 0; i < order; i++) {
339 coefficients(i, 0) = weights.tr().dot(func_values)(0, 0) / std::pow(2, i) / factorial / std::sqrt(M_PI);
340 factorial *= i;
341 }
342
343 return coefficients;
344}
Class for computing Hermite polynomials and related operations.
Definition Hermite_Tools.h:11
Matrix< T > evaluate(Matrix< T > coefficients, Matrix< T > values)
Evaluates the Hermite polynomial at specified values.
Definition Hermite_Tools.h:136
Matrix< T > import_data(const char *file)
Imports data for roots or weights from a file.
Definition Hermite_Tools.h:246
Matrix< T > compute_coefficients(unsigned int order, T(*func)(T))
Computes the coefficients for the Hermite polynomial using a given function.
Definition Hermite_Tools.h:324
Matrix< T > get_roots(unsigned int order)
Retrieves the roots of the Hermite polynomial of specified order.
Definition Hermite_Tools.h:276
Matrix< T > coefficients_deriv(Matrix< T > coefficients)
Computes the derivative of the Hermite polynomial coefficients.
Definition Hermite_Tools.h:191
~Hermite()
Destructor for Hermite class.
Definition Hermite_Tools.h:126
Hermite()
Default constructor for Hermite class.
Definition Hermite_Tools.h:116
Matrix< T > get_weights(unsigned int order)
Retrieves the weights of the Hermite polynomial of specified order.
Definition Hermite_Tools.h:299
Matrix< T > evaluate_deriv(Matrix< T > coefficients, Matrix< T > values)
Evaluates the derivative of the Hermite polynomial at specified values.
Definition Hermite_Tools.h:215
This class implements a matrix object used for linear algebra and vectorized operations.
Definition Matrix.h:25
Matrix dot(const Matrix &other)
Calculates the dot product of two matrices.
Definition Matrix.h:604
Matrix tr()
Transposes the matrix.
Definition Matrix.h:705