Alexandria  2.18
Please provide a description of the project.
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 
29 namespace Euclid {
30 namespace MathUtils {
31 
32 double 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 
52 double 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
SimpsonsRule.h
Euclid::MathUtils::Function
Interface class representing a function.
Definition: Function.h:46
Euclid::MathUtils::SimpsonsRule::operator()
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.
Definition: SimpsonsRule.cpp:32
Exception.h
Elements::Exception
Euclid
Definition: InstOrRefHolder.h:29