Numerical Integration

Unmodified
Removed
• Code

r1 r2
1 #include <iostream>
1#include <iostream>  // Not needed for integration routines
2#include <cmath>  // sqrt, exp, acos
23using namespace std;
4
5double PI = acos(-1);
36
47double x_squared(double x) {
58  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);
614}
715
816double rectangle_method(double (*func)(double), double x1, double x2, int n) {
917  /*
1018  Integrate func between x1 and x2 using the rectangle method with the midpoint
1119  approximation for n samples.
1220
1321  */
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
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));
1938  }
2039  return integral;
2140}
2241
2342int 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;
2652
2753  return 0;
2854}
2955