Alexandria
2.25.0
SDC-CH common library for the Euclid project
MathUtils
src
lib
numericalIntegration
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
25
#include "
MathUtils/numericalIntegration/SimpsonsRule.h
"
26
#include "
ElementsKernel/Exception.h
"
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
Exception.h
N
constexpr std::size_t N
SimpsonsRule.h
Elements::Exception
Euclid::MathUtils::NAryFunction
Interface class representing a function with an arbitrary number of parameters.
Definition:
Function.h:104
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
Euclid
Definition:
index_sequence.h:27
Generated by
1.9.1