Alexandria  2.14.1
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
SimpsonsRule.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012-2020 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 
25 #include <cmath>
28 
29 namespace Euclid {
30 namespace MathUtils {
31 
32 double SimpsonsRule::operator()(const Function& function, double min,
33  double max, int order) {
34  if (order < 3) {
35  throw Elements::Exception()
36  << "Simpson's Rule integration is define only for order bigger than 2";
37  }
38 
39  int N = pow(2, order);
40  double h = (max - min) / N;
41 
42  double partial_sum = 0;
43  for (int i = 3; i < N - 2; i++) {
44  partial_sum += function(min + i * h);
45  }
46 
47  partial_sum += 0.375 * (function(min) + function(max));
48  partial_sum += 7. * (function(min + h) + function(max - h)) / 6.;
49  partial_sum += 23. * (function(min + 2. * h) + function(max - 2 * h)) / 24.;
50 
51  return partial_sum * h;
52 
53 }
54 
55 double SimpsonsRule::operator()(const Function& function, double min,
56  double max, double previous_value, int order) {
57  if (order < 4) {
58  throw Elements::Exception()
59  << "Simpson's Rule integration with recursion is define only for order bigger than 3";
60  }
61 
62  int N = pow(2, order);
63  double h = (max - min) / N;
64 
65  double partial_sum = 0;
66 
67  for (int j = 1; j < N / 2 - 1; j++) {
68  int i = 2 * j + 1;
69  partial_sum += function(min + i * h);
70  }
71 
72  partial_sum += 7. * (function(min + h) + function(max - h)) / 6.;
73  partial_sum -= 5. * (function(min + 2. * h) + function(max - 2. * h)) / 24.;
74  partial_sum += (function(min + 4. * h) + function(max - 4. * h)) / 24.;
75  return partial_sum * h + previous_value / 2.;
76 }
77 
78 } // End of MathUtils
79 } // end of namespace Euclid
Interface class representing a function.
Definition: Function.h:46
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&gt;=3. ...