License Public Domain
Lines 73
Keywords
Num_integration (3)
Permissions
Owner: jsimones
Group Owner: SnortSnort
Viewable by Everyone
Editable by All Siafoo Users

Numerical Integration 2 Atom Feed 0

# 's
 1/*
2 Rectangle method and Simpson's rule integrators for x and y data.
3
4*/
5
6#include <iostream> // Not needed for integration routines
7#include <cmath> // sqrt, exp, acos
8#include <vector>
9using namespace std;
10
11double PI = acos(-1);
12
13double x_squared(double x) {
14 return x * x;
15}
16
17double normal_dist(double x) {
18 double m = 0, s = 1;
19 return exp((x - m) * (m - x) / 2 / s / s) / s / sqrt(2 * PI);
20}
21
22double rectangle_method(vector<double>& x_list, vector<double>& y_list) {
23 /*
24 Integrate a function given vectors of x and y data using the rectangle method
25 with the midpoint approximation.
26
27 */
28 double integral = y_list[0] * (x_list[1] - x_list[0]) / 2;
29 for (int i = 1; i < x_list.size() - 1; i++) {
30 integral += y_list[i] * (x_list[i+1] - x_list[i-1]) / 2;
31 }
32 int i = x_list.size() - 1;
33 integral += y_list[i] * (x_list[i] - x_list[i - 1]) / 2;
34 return integral;
35}
36
37double simpsons_rule(vector<double>& x_list, vector<double>& y_list) {
38 /*
39 Integrate a function given vectors of x and y data using Simpson's rule.
40 Assumes that x values are equally spaced! The looping in this implementation
41 requires that x_list and y_list contain an odd number of elements, so if the
42 number of elements is even, then the very last bin will be integrated using
43 the trapezoid rule.
44
45 */
46 double integral = 0;
47 for (int i = 1; i < x_list.size() - 1; i += 2) {
48 integral += (x_list[i+1] - x_list[i-1]) / 6 *
49 (y_list[i-1] + 4 * y_list[i] + y_list[i+1]);
50 }
51 if (x_list.size() % 2 == 0) {
52 int i = x_list.size() - 1;
53 integral += (y_list[i-1] + y_list[i]) * (x_list[i] - x_list[i-1]) / 2;
54 }
55 return integral;
56}
57
58int main() {
59 // Generate x,y data for n samples:
60 int n = 1000;
61
62 double x1a = 1, x1b = 10, x2a = -10, x2b = 10;
63 double step1 = (x1b - x1a) / (n - 1), step2 = (x2b - x2a) / (n - 1);
64 vector<double> x_list1, x_list2, y_list1, y_list2;
65 for (int i = 0; i < n; i++) {
66 // y_list1 = x_squared of x_list1
67 x_list1.push_back(x1a);
68 y_list1.push_back(x_squared(x1a));
69 x1a += step1;
70
71 // y_list2 = normal_dist of x_list2
72 x_list2.push_back(x2a);
73 y_list2.push_back(normal_dist(x2a));
74 x2a += step2;
75 }
76
77 // Test using x_squared data:
78 double integral1 = rectangle_method(x_list1, y_list1);
79 double integral2 = simpsons_rule(x_list1, y_list1);
80 cout << integral1 << " " << integral2 << endl;
81
82 // Test using normal_dist data:
83 integral1 = rectangle_method(x_list2, y_list2);
84 integral2 = simpsons_rule(x_list2, y_list2);
85 cout << integral1 << " " << integral2 << endl;
86
87 return 0;
88}