riemann sums in c++

[cygh]

void float_Exp (int maxiter)
{
float sections = 6.0/maxiter, total=0.0, maxitercopy=maxiter;
while (maxiter>0)
{
float x = (maxiter*sections)-3.0;
total = total + exp(-pow(x,2.0));
maxiter--;
}
cout << "exp = " << total << " for "<< maxitercopy << " iterations \n" ;
}
void float_T2 (int maxiter)
{
float sections = 6.0/maxiter, total=0.0, maxitercopy=maxiter;
while (maxiter>0)
{
float x = (maxiter*sections)-3.0;
total = total + (pow (x,4.0)/2.0)- pow (x,2.0) + 1.0;
maxiter--;
}
cout << "float t2 = " << total << " for "<< maxitercopy << " iterations \n" ;
}

i am passing in 100, 1000 then 10,000 for each function. the first is supposed to use the c++ exp() method to get me the right result, the second needs to use the first two terms of the Taylor series of the function e ^ -1(x)^2. both methods are supposed to take my function and compute the Riemann sum of them. im getting wrong results, the number of iterations is drastically changing my answer. and i cant find my error. Was hoping i could get some help with this please.

[D80DSL3A]

Just to check that I am on the same page here.
A Riemann sum is used to approximate a definite integral, right?
If so, how does the integration interval enter your calcs here?
Is sections supposed to be the length of the sub-intervals?
And a Taylor series is used to approximate a function.
Taylor: f(x) = f(x0) + f'(x0)*(x-x0) +(1/2)* f''(x0)*(x-x0)^2 + ...
How is that related to evaluating an integral?

[cygh]

basically we got a function, e^(-1(x)^2) for -3 to 3. we have to evaluate this function using exp built in method as well as using taylor series with 2,6,and 8 terms. once we have one of the taylor series we use Riemann sum on it to approximate the area under the curve. yes sections s length of interval. also the range, and first couple of terms of the taylor series are all hard coded.

void float_Exp (int maxiter)
{
     float sections = 6.0/maxiter, total=0.0, maxitercopy=maxiter;
     while (maxiter>0)
     {
          float x = -3.0+maxiter*sections;
          total = total + exp(-pow(x,2.0));
          maxiter--;
     }
     cout <<  "float exp = " << total << " for "<< maxitercopy << " iterations \n" ;
}
void float_T2 (int maxiter)
{
     float sections = 6.0/maxiter, total=0.0, maxitercopy=maxiter;
     while (maxiter>0)
     {
          float x = -3.0+maxiter*sections;
          total = total + (pow (x,4.0)/2.0)- pow (x,2.0) + 1.0;
          maxiter--;
     }
     cout <<  "float t2 = " << total << " for "<< maxitercopy << " iterations \n" ;
}
void float_T6 (int maxiter)

the code is slightly different now but i still have the same issues, using 100 intervals of riemann gives me 1/10th of the value of using 1000 intervals same for 1000 and 10000. also I am not getting any difference in output between my code using floats and my code using doubles even though seeing the difference is supposed to be the point of my lab.

thx for the help btw

[ne555]

Riemann approximates the function with rectangles. The area under the curve is then the area under the rectangles.

Area of a rectangle = base * height

You are just calculating the height

[D80DSL3A]

@ cygh: You're welcome. Was ne555's hint sufficient?
I think it concerns a missing factor of sections on lines 7 and 18. You need integrand*deltaX there.

[cygh]

yeah i added i *sections on lines 7/18 it fixed the issue i had at the time with diffrent numbers of iterations giving me different answers. the rest of the code is still completely messed up though getting ridiculous answers, but thats another story. thanks guys!

[D80DSL3A]

In case you come back again...
I believe that the Taylor series you are using is incorrect.
The series for e^u = 1 + u + u*u/2 + ...
For your problem u = -(X^2)/2 though, so the 1st 3 terms should be:
CORRECTION: u = -(x^2)
T2(x) = 1.0 - x*x + x*x*x*x/2.0;// EDIT

Also, since I had wanted to add a Riemann sum function to my collection of functions anyways, I thought I'd share the result.

This function takes a function pointer as an argument so that it doesn't have to be modified when you want it to use a different integrand. a and b are the lower and upper integration limits (respectively) and n is the # of sub-intervals.

double RiemannSum(const double& a, const double& b, double (*pIntegrand)(double), const int& n )
{
	double sum = 0.0;	
	double subInt = (b-a)/static_cast<double>(n);
	double x = a + subInt/2.0;// evaluate integrand at middle of sub-interval

	for(int k=0; k<n; k++)
	{
		sum += pIntegrand(x)*subInt;
		x += subInt;
	}

	return sum;
}

You could then write whatever function for the integrand you like (as long as it takes one double as argument and returns a double result).
Two examples could be:

double eMinusXsq(double x)
{
    return exp(-x*x);// EDIT
}

double T2(double x)
{
    double u = -x*x;// EDIT
    return (1.0 + u + u*u/2.0);
}

Examples of calls to RiemannSum() are:

    // 10 terms
    cout << "Using exact = " << RiemannSum( -3.0, 3.0, eMinusXsq, 10 ) << endl;
    cout << "Using T2    = " << RiemannSum( -3.0, 3.0, T2, 10 ) << endl;

[cygh]

"For your problem u = -(X^2)/2 though"
i dont think it should be divided by 2 unless i am missing something else.

[D80DSL3A]

You're right! Not sure where I got that from. Omit the division by 2.0 in the 2 functions then (line 3 in eXsqOver2 and line 8 in T2 - but keep it in line 9). This makes your series correct then.

[cygh]

i took my taylor series expansion from wolframalpha here
http://www.wolframalpha.com/input/?i=taylor+e^(-1*(x^2))

not really sure where the difference between the two comes from but there is a difference.

[D80DSL3A]

The series you had was fine. I had been looking at another function while working on my program. Sorry for the confusion.
The only other difference is that I am not using the pow() as pow(x,2.0) is the same as x*x. I have modified the 2 functions I gave to reflect this.

[tition]

Umm, when you do Riemanian sums, don't you have to multiply by the size of the interval?

like,

total = total + ((pow (x,4.0)/2.0)- pow (x,2.0) + 1.0)*sections;

oops just saw you guys already found that issue...

[Edit:] Btw, as far as I remember, the integral of -infinity to +infinity of e^{-x^2} was something like \sqrt (\pi) (not 100% sure). Now, at the end points of your interval, e^{-x^2} is like e^{-9}, which is of the order of 1/3^9, damn small.

So, you should be getting as a final answer something close to sqrt(pi) ~ 1.77 .