NUMERICAL INTEGRATION: SIMPSONS 1/3 RULE







#include "stdafx.h"
#include "iostream"
#include "conio.h"
#include "math.h"
using namespace std;

class Integration{
private:
       double A; //Area under curve
       double xmin, xmax; //limits of integration
       int numberOfPoints;
       double(*f)(double x); // Function to be integratied
public:
       Integration(double(*F)(double x), double a, double b, int n)
       {
             f = F;
             xmin = a;
             xmax = b;
             numberOfPoints = n;
       }
       double Trapezoidal();
       double Simpson();
};


double Integration::Trapezoidal()
{
       double dx = (xmax - xmin) / double(numberOfPoints);
       double sum = 0.0;
       double x = xmin + dx;
       for (int i = 1; i < numberOfPoints; i++)
       {
             sum += (*f)(x);
             x += dx;
       }
       A = (*f)(xmin)+(*f)(xmax)+2.0*sum;
       A *= dx / double(2.0);
       return A;
}
double Integration::Simpson(){
       double dx = (xmax - xmin) / double(numberOfPoints);
       double sum1 = 0, sum2 = 0.0;
       double x = xmin + dx;
       for (int i = 1; i < numberOfPoints; i++)
       {
             if (i % 2 != 0)
                    sum1 += (*f)(x);
             else
                    sum2 += (*f)(x);
             x += dx;
       }
       A = (*f)(xmin)+(*f)(xmax)+4.0*sum1 + 2.0*sum2;
       A *= dx / double(3.0);
       return A;
}
double FUN(double x) // USER supplied function
{
       return x* sqrt(8.0 - x*x*x);
}
int _tmain(int argc, _TCHAR* argv[])
{
       Integration I(FUN, 0, 2, 8);
       cout << I.Trapezoidal() << endl;
       cout << I.Simpson() << endl;
       system("pause");
       return 0;
}



Share on Google Plus