Meet people who work on similar things as you – get help if you need it
Join Siafoo Now
or
Learn More
Numerical Integration
Revision 1 vs. Revision 2
Legend:
 Unmodified
 Added
 Removed

Code
r1 r2 1 #include <iostream> 1 #include <iostream> // Not needed for integration routines 2 #include <cmath> // sqrt, exp, acos 2 3 using namespace std; 4 5 double PI = acos(1); 3 6 4 7 double x_squared(double x) { 5 8 return x * x; 9 } 10 11 double normal_dist(double x) { 12 double m = 0, s = 1; 13 return exp((x  m) * (m  x) / 2 / s / s) / s / sqrt(2 * PI); 6 14 } 7 15 8 16 double rectangle_method(double (*func)(double), double x1, double x2, int n) { 9 17 /* 10 18 Integrate func between x1 and x2 using the rectangle method with the midpoint 11 19 approximation for n samples. 12 20 13 21 */ 14 double step = (x2  x1) / n; 15 double x = x1 + step / 2, integral = 0; 16 for (int i = 0; i < n; i++) { 17 integral += (*func)(x) * step; 18 x += step; 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 29 double 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)); 19 38 } 20 39 return integral; 21 40 } 22 41 23 42 int main() { 24 double integral = rectangle_method(x_squared, 1, 10, 1000); 25 cout << integral << endl; 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; 26 52 27 53 return 0; 28 54 } 29 55