Simpson's method code has a redundancy and I can't find it!

[Aristophanes]

Hello,
I just joined, so if I infringe upon any of the forum rules, then I apologize. Anyway, I have a code for approximating integrals via the Simpson's Method. The problem with the code is that somewhere within the mere 52 lines, I have divided my entire sum by 2. Having seen this in the end result, I found I have to multiply the whole thing by 2 to get the right answer. I know this is redundant, and I fear that, depending where the devision is, it might make my integration approx. less accurate. If it isn't too much work, can someone find where I went wrong in my code? By the way, this is not the usual Simpson's approx. It is called the Composite Simpson's Rule; the algorithm can be found on Wikipedia: http://en.wikipedia.org/wiki/Simpson_method

Also, as a side note, how do I turn this source code into a class for an arbitrary function (not just, say, for cos(x)?)

/*
 * Simprule.cpp
 *
 *  Created on: Mar 14, 2009
 *      Author: keith burghardt
 */
# include <iostream>
using namespace std;
#include <cmath>
#include<iomanip>
double begin,end,area,h;

long double f(long double x){
	return(//State function for integration
	cos(x)
	);
}
long double simps(double,double,long double,double);

int main (void){

cout<<"State beginning and end number for integration: ";
cin>>begin>>end;
long int n(1e6*int(fabs(end-begin)));//make sure n is even (must be able to divide by 2, see below)
long double h(fabs(end-begin)/n);//the fabs() are useful in making sure that, in the simpsons rule,
								//I need not have 2 cases for the for loops
area=simps(begin,end,h,n);
cout<<"AREA of F(x)= "<<setprecision(19)<<area<<endl<<"Area really is "<<sin(end);//I use this to check my answer
	return 0;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//composite simp rule
long double simps (double begin, double end,long double h,double n){//using simpson's rule to approx. an area
	double b(begin),e(end);
	if(begin>end){//needed in order to make sure that I don't have to make 2 for loops
					//one for b>e, and one for e>b.
	b=end;
	e=begin;
}
	long double sum(f(b)),sum2(0),sum4(0);
for(int i=1;i<=n/2-1;i+=2){//adding up all even points
	sum2=sum2+f(b+h*(2*i));
}
for(int k=1;k<=n/2;k+=2){//adding up all odd points
	sum4=sum4+f(b+h*(2*k-1));
}
sum=(2*h/3)*(sum+2*sum2+4*sum4+f(e));//!!*can't find where I divided answer*!!<----------PROBLEM!!!
if(begin>end){//needed so that if I integrate from a larger to smaller number,
				sum=-sum;//the sum will be negitive
				}
return(sum);
}

[MattEvans]

In your for loops, you are adding 2 to i and k. You probably want to add 1 to i and k, since your trying to get 2*i and 2*k-1, i.e the odds and evens... so it looks like you are skipping half of the points by adding 2 instead of 1. I can't see where (or if) you're performing the redundant multiplication in the code that you posted, but when I made that change, the result needed to be divided by 2 to be correct (it was correct before for ranges like, 0-1, 0-10 etc, so if that's not the problem, what ranges ARE you trying?).

Second question, pass a function pointer as an argument.. the signature of your function is long double () ( long double ), so change your simps function to be:

long double simps (double begin, double end,long double h,double n, [b]long double ( *f ) ( long double )[/b] )

and that's actually all you need to do.. you can pass any function in place of cos, now, providing that the signature is exactly the same.

Btw, the signature for forward delcaring simps will change to:

long double simps(double,double,long double,double,long double(*)(long double));

Actually, I may aswell just go ahead and post the modified version, all changes emboldened:

/*
 * Simprule.cpp
 *
 *  Created on: Mar 14, 2009
 *      Author: keith burghardt
 */
# include <iostream>
using namespace std;
#include <cmath>
#include<iomanip>
double begin,end,area,h;

long double [b]coswrapper[/b](long double x){
	return(//State function for integration
	cos(x)
	);
}

long double simps(double,double,long double,double[b],long double(*)(long double)[/b]);

int main (void){

cout<<"State beginning and end number for integration: ";
cin>>begin>>end;
long int n(1e6*int(fabs(end-begin)));//make sure n is even (must be able to divide by 2, see below)
long double h(fabs(end-begin)/n);//the fabs() are useful in making sure that, in the simpsons rule,
								//I need not have 2 cases for the for loops
area=simps(begin,end,h,n[b],&coswrapper[/b]);
cout<<"AREA of F1(x)= "<<setprecision(19)<<area<<endl<<"Area really is "<<sin(end);//I use this to check my answer
	return 0;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//composite simp rule
long double simps (double begin, double end,long double h,double n[b],long double (*f) (long double)[/b]){//using simpson's rule to approx. an area
	double b(begin),e(end);
	if(begin>end){//needed in order to make sure that I don't have to make 2 for loops
					//one for b>e, and one for e>b.
	b=end;
	e=begin;
}
	long double sum(f(b)),sum2(0),sum4(0);
for(int i=1;i<=n/2-1;[b]i+=1[/b]){//adding up all even points
	sum2=sum2+f(b+h*(2*i));
}
for(int k=1;k<=n/2;[b]k+=1[/b]){//adding up all odd points
	sum4=sum4+f(b+h*(2*k-1));
}
sum=(2*h/3)*(sum+2*sum2+4*sum4+f(e));//!!*can't find where I divided answer*!!<----------PROBLEM!!!
if(begin>end){//needed so that if I integrate from a larger to smaller number,
				sum=-sum;//the sum will be negitive
				}
return(sum[b]/2[/b]);
}

[Aristophanes]

Thanks for the help, the code with your fixes can now easily give me 12 or so places of accuracy, and I don't have to double the answer!