-
Notifications
You must be signed in to change notification settings - Fork 0
/
IntegralSolver.hpp
97 lines (80 loc) · 2.38 KB
/
IntegralSolver.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#include <utility>
#include <cmath>
#include <vector>
#include <functional>
#include <fstream>
#include <iostream>
namespace IntegralSolver
{
double Trapezoid(const std::vector<double> &fvec,
const double &dx)
{
const std::size_t N = fvec.size();
const std::size_t N1 = N - 1;
double ret = 0.5 * fvec[0];
for (std::size_t i = 1; i < N1; ++i)
{
ret += fvec[i];
}
ret += 0.5 * fvec[N1];
ret *= dx;
return ret;
}
double Simpson(const std::vector<double> &fvec,
const double &dx)
{
const std::size_t N = fvec.size();
const std::size_t N1 = N - 1;
if (N % 2 == 0)
{
std::cerr << "# of quadrature points must be an odd number!" << std::endl;
std::exit(EXIT_FAILURE);
}
double ret = fvec[0] + 4.0 * fvec[1];
for (std::size_t i = 2; i < N1; i += 2)
{
ret += 2.0 * fvec[i] + 4.0 * fvec[i + 1];
}
ret += fvec[N1];
ret *= dx / 3.0;
return ret;
}
double Simpson38(const std::vector<double> &fvec,
const double &dx)
{
const std::size_t N = fvec.size();
const std::size_t N1 = N - 1;
if (N % 3 != 1)
{
std::cerr << "# of quadrature points must be a multiple of 3 + 1!" << std::endl;
std::exit(EXIT_FAILURE);
}
double ret = fvec[0] + 3.0 * (fvec[1] + fvec[2]);
for (std::size_t i = 3; i < N1; i += 3)
{
ret += 2.0 * fvec[i] + 3.0 * (fvec[i + 1] + fvec[i + 2]);
}
ret += fvec[N1];
ret *= 3.0 * dx / 8.0;
return ret;
}
double Boole(const std::vector<double> &fvec,
const double &dx)
{
const std::size_t N = fvec.size();
const std::size_t N1 = N - 1;
if (N % 4 != 1)
{
std::cerr << "# of quadrature points must be a multiple of 4 + 1!" << std::endl;
std::exit(EXIT_FAILURE);
}
double ret = 7.0 * fvec[0] + 32.0 * (fvec[1] + fvec[3]) + 12.0 * fvec[2];
for (std::size_t i = 4; i < N1; i += 4)
{
ret += 14.0 * fvec[i] + 32.0 * (fvec[i + 1] + fvec[i + 3]) + 12.0 * fvec[i + 2];
}
ret += 7.0 * fvec[N1];
ret *= dx / 22.5;
return ret;
}
}