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 
    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