License Public Domain
Lines 44
Keywords
Num_integration (3)
Permissions
Owner: jsimones
Group Owner: SnortSnort
Viewable by Everyone
Editable by All Siafoo Users
Hide
Easily highlight source code for your blog with our Syntax Highlighter. Join Siafoo Now or Learn More

Numerical Integration Atom Feed 0

# 's
 1#include <iostream>  // Not needed for integration routines
2#include <cmath> // sqrt, exp, acos
3using namespace std;
4
5double PI = acos(-1);
6
7double x_squared(double x) {
8 return x * x;
9}
10
11double normal_dist(double x) {
12 double m = 0, s = 1;
13 return exp((x - m) * (m - x) / 2 / s / s) / s / sqrt(2 * PI);
14}
15
16double rectangle_method(double (*func)(double), double x1, double x2, int n) {
17 /*
18 Integrate func between x1 and x2 using the rectangle method with the midpoint
19 approximation for n samples.
20
21 */
22 double step = (x2 - x1) / n, integral = 0;
23 for (x1 += step / 2; x1 < x2; x1 += step) {
24 integral += (*func)(x1) * step;
25 }
26 return integral;
27}
28
29double simpsons_rule(double (*func)(double), double x1, double x2, int n) {
30 /*
31 Integrate func between x1 and x2 using simpson's rule for n samples.
32
33 */
34 double step = (x2 - x1) / n, integral = 0;
35 for (x1; x1 < x2; x1 += step) {
36 integral += step / 6 * ((*func)(x1) + 4 * (*func)(x1 + step / 2) +
37 (*func)(x1 + step));
38 }
39 return integral;
40}
41
42int main() {
43 // Test using x_squared:
44 double integral1 = rectangle_method(x_squared, 1, 10, 1000);
45 double integral2 = simpsons_rule(x_squared, 1, 10, 1000);
46 cout << integral1 << " " << integral2 << endl;
47
48 // Test using normal_dist:
49 integral1 = rectangle_method(normal_dist, -10, 10, 1000);
50 integral2 = simpsons_rule(normal_dist, -10, 10, 1000);
51 cout << integral1 << " " << integral2 << endl;
52
53 return 0;
54}