Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
SimpsonsRule.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012-2021 Euclid Science Ground Segment
3 *
4 * This library is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU Lesser General Public License as published by the Free
6 * Software Foundation; either version 3.0 of the License, or (at your option)
7 * any later version.
8 *
9 * This library is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU Lesser General Public License
15 * along with this library; if not, write to the Free Software Foundation, Inc.,
16 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17 */
18
27#include <cmath>
28
29namespace Euclid {
30namespace MathUtils {
31
32double SimpsonsRule::operator()(const Function& function, double min, double max, int order) {
33 if (order < 3) {
34 throw Elements::Exception() << "Simpson's Rule integration is define only for order bigger than 2";
35 }
36
37 int N = pow(2, order);
38 double h = (max - min) / N;
39
40 double partial_sum = 0;
41 for (int i = 3; i < N - 2; i++) {
42 partial_sum += function(min + i * h);
43 }
44
45 partial_sum += 0.375 * (function(min) + function(max));
46 partial_sum += 7. * (function(min + h) + function(max - h)) / 6.;
47 partial_sum += 23. * (function(min + 2. * h) + function(max - 2 * h)) / 24.;
48
49 return partial_sum * h;
50}
51
52double SimpsonsRule::operator()(const Function& function, double min, double max, double previous_value, int order) {
53 if (order < 4) {
54 throw Elements::Exception() << "Simpson's Rule integration with recursion is define only for order bigger than 3";
55 }
56
57 int N = pow(2, order);
58 double h = (max - min) / N;
59
60 double partial_sum = 0;
61
62 for (int j = 1; j < N / 2 - 1; j++) {
63 int i = 2 * j + 1;
64 partial_sum += function(min + i * h);
65 }
66
67 partial_sum += 7. * (function(min + h) + function(max - h)) / 6.;
68 partial_sum -= 5. * (function(min + 2. * h) + function(max - 2. * h)) / 24.;
69 partial_sum += (function(min + 4. * h) + function(max - 4. * h)) / 24.;
70 return partial_sum * h + previous_value / 2.;
71}
72
73} // namespace MathUtils
74} // end of namespace Euclid
Interface class representing a function with an arbitrary number of parameters.
Definition Function.h:104
double operator()(const Function &function, double min, double max, int order)
Integrate a function between min and max using a mesh of 2^order steps. order>=3.
std::array< std::vector< double >, N > Coordinates
Used to pass the grid coordinates to interpn. Internally will make a copy of the required values.