Metody numeryczne 01.05.18

NewtonCotesQuadrature.cpp

#include <iostream>

using namespace std;

double compositeTrapezoidRule(double * pointsX, double *pointsY, int N){
	double h = (pointsX[N - 1] - pointsX[0]) / (N - 1);
	double area = 0.0;
	area += pointsY[0] + pointsY[N - 1];
	for (int i = 1; i < N - 1; i ++){
		area += 2 * pointsY[i];
	}
	return ((area*h) / 2.0);
}

double compositeSimpsonRule(double * pointsX, double *pointsY, int N){
	if (N > 2){
		double h = (pointsX[N - 1] - pointsX[0]) / (N - 1);
		double area = 0.0;
		double area2 = 0.0;
		if (N % 2 == 0){	//Simpson 1/3 doesn't work when the number of points is even so we use Simpson 3/8 at the end of range
			area2 += pointsY[N - 4]+ 3 * pointsY[N - 3]+ 3 *pointsY[N - 2]+ pointsY[N - 1];
			N -= 3;	//we used last 4 points so we no longer need them in our calculations EDIT: i have no idea why you have to remove last 3 points and not 4 but 3 works.
			cout << "even";
		}
		area += pointsY[0]+ pointsY[N - 1];
		for (int i = 1; i < N - 1; i += 2){
			area += 4 * pointsY[i];
		}
		for (int i = 2; i < N - 2; i += 2){
			area += 2 * pointsY[i];
		}
		return ((area*h) / 3)+((area2*h*3)/8);
	}
	return 0;
}

double cube(double x){
	return x*x*x;
}


int main(){

	//when odd
	const int N = 9;	//number of points degree is N-1
	double  givenPointsX[N] = {0,0.5,1,1.5,2,2.5,3,3.5,4};
	double  givenPointsY[N] = { 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4 };
	cout << "Simpson: " <<compositeSimpsonRule(givenPointsX, givenPointsY, N)<<endl;
	cout << "Trapezoid: " << compositeTrapezoidRule(givenPointsX, givenPointsY, N) << endl;


	//when even
	const int N2 = 12;	//number of points degree is N-1
	double  givenPointsX2[N2] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11 };
	double  givenPointsY2[N2] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11 };
	cout << "Simpson: " <<compositeSimpsonRule(givenPointsX2, givenPointsY2, N2) << endl;
	cout << "Trapezoid: " << compositeTrapezoidRule(givenPointsX2, givenPointsY2, N2) << endl;

	//when only 4 points (to test 3/8 simpson)
	const int N3 = 4;	//number of points degree is N-1
	double  givenPointsX3[N3] = { 0, 1, 2, 3};
	double  givenPointsY3[N3] = { 0, 1, 2, 3};
	cout << "Simpson: " << compositeSimpsonRule(givenPointsX3, givenPointsY3, N3) << endl;
	cout << "Trapezoid: " << compositeTrapezoidRule(givenPointsX3, givenPointsY3, N3) << endl;

	//when f(x)= x^3
	const int N4 = 10;	//number of points degree is N-1
	double  givenPointsX4[N4];
	double  givenPointsY4[N4];
	for (int i = 0; i < N4; i++){
		givenPointsX4[i]=i;
		givenPointsY4[i] = cube(i);
	}
	cout << "Simpson: " << compositeSimpsonRule(givenPointsX4, givenPointsY4, N4) << endl;
	cout << "Trapezoid: " << compositeTrapezoidRule(givenPointsX4, givenPointsY4, N4) << endl;



	system("pause");
}