Alexandria
2.16
Please provide a description of the project.
MathUtils
src
lib
numericalIntegration
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>
26
#include "
MathUtils/numericalIntegration/SimpsonsRule.h
"
27
#include "
ElementsKernel/Exception.h
"
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
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
N
constexpr std::size_t N
Euclid
Definition:
InstOrRefHolder.h:29
Generated by
1.8.18