Trapezoid Rule C++ programming help

[jai6638]

hey there.

Posted in the wrong forum so reposting this again.
_______


Hello there..

I'm writing a trapezoidal rule program but even though it runs, it outputs the wrong value for finding the area under the curve for e^((-(x)^2)/2) from 1 to 2.

double trapRule(double limit_B,double limit_A,double number_Partitions);

int main()
{
double number_Partitions=0.0;
double limit_A=1.0;
double limit_B=2.0;
double area=0.0;

cout<<"Please enter the number of partitions"<<endl;
cin>>number_Partitions;
area=trapRule(limit_A,limit_B,number_Partitions);
cout<<"The area under the curve for the function is "<<area<<endl;

return 0;


}

double trapRule(double limit_A, double limit_B,double number_Partitions)
{
double delta_X=0.0;
double function=0.0;
double height=0.0;
double x=0.0;
double j=0.0;

delta_X=(limit_B-limit_A)/number_Partitions;
function=exp(pow(-x,2.0)/2.0);
x=(limit_A+delta_X);

for(j=x;j<2.0;j=j+delta_X)
{

height=height+exp((pow(-(j),2.0))/2.0);

}

double height_A=(0.5)*(exp(-1.0/2.0));
double height_B=(0.5)*(exp(-2.0));

double area_Function=delta_X*(height_A+height_B+height);

return area_Function;
}

I get the following results when i put the following number of partitions:

2 Partitions: .684566
3 Partitions: 1.12364
4 Partitions: 0.842733
100 Partitions: 0.993709
1000 Parititions: 1.0037

These numbers dont look right however..

Any help is much appreciated. I've tried debugging and it seems that my problem is in the for loop but not sure how to fix it.. I've tried different stuff but to no avail.

 

[QuixoticOne]

Just put your f(x) as a function in *one* place in the code.

Also (-x)^2 == x^2 just to simplify.

You'll find that your height_A and height_B formulae are incorrect, probably an oversight that was easier to make because you duplicated the function code there rather than calling a single function function.

height_A = f(a)/2 which isn't what you've written.
height_B = f(b)/2 which isn't what you've written.

That should make your answer more correct.

Also be sure you're including the proper headers that define your math functions and are linking with the appropriate libraries.

I assume that you're doing this as a learning / programming experiment as opposed to just trying to get the value of such integrals since there are better available library programs to calculate integrals for you with less work involved in coding them, of course.

Originally posted by: jai6638

function=exp(pow(-x,2.0)/2.0);
height=height+exp((pow(-(j),2.0))/2.0);
double height_A=(0.5)*(exp(-1.0/2.0));
double height_B=(0.5)*(exp(-2.0));

 

[MoMeanMugs]

QuixoticOne,
I'm pretty sure he wants -x^2 and not (-x)^2. At least that's what I gather from the form he's given. He never used "function." Also, how are his equations for Height_A and Height_B incorrect? It's not dynamic for A and B, but there's nothing wrong with the evaluation.

 

[QuixoticOne]

Originally posted by: MoMeanMugs
QuixoticOne,
I'm pretty sure he wants -x^2 and not (-x)^2. At least that's what I gather from the form he's given. Also, how are his equations for Height_A and Height_B incorrect? Last time I checked, 1/2=.5.

Well he wrote in two places evaluations of the f(x) as follows:

function=exp(pow(-x,2.0)/2.0);
height=height+exp((pow(-(j),2.0))/2.0);

pow(a,b) returns a^b where the expressions of (a) and (b) being totally separate formal arguments to the function are evaluated independently and previously to the exponentiation operation.

So pow(-x,2) = (-x) * (-x) = x^2. If he MEANT -x^2 and NOT (-x)^2 then he's certainly made a mistake in those function expressions.

In other places he wrote:

double height_A=(0.5)*(exp(-1.0/2.0));
double height_B=(0.5)*(exp(-2.0));

where the correct (AFAIK) formulae would be
height_A = f(a)/2
height_B = f(b)/2
since a = 1 and b = 2
f(a) = f(1) = exp( (a^2)/2) = exp((1^2)/2) = exp(1/2).
f(b) = f(2) = exp((b^2)/2) = exp((2^2)/2) = exp(2).

However that's totally inconsistent with:

double height_A=(0.5)*(exp(-1.0/2.0));
double height_B=(0.5)*(exp(-2.0));

So either way he means -x^2 or (-x)^2 he's using one form in a couple of places and the OTHER in a couple of OTHER places, hence my suggestion that his main problem was inconsistent multiple definitions of f(x) and incorrect definitions in half of them.

 

[MoMeanMugs]

I've traced through your program again, and I'm pretty sure your error lies in "height=height+exp((pow(-(j),2.0))/2.0); " like was stated in your Off Topic thread. The integral from 1 to 2 for that function is approximately .341. In order for that to be true with only two intervals, height needs to be approximately .31. However, your height is about .998.

 

[MoMeanMugs]

He's multiplying by .5 instead of doing it in the trapezoid formula. The trapezoid formula is:
area=delta_X/2*(Height_A+Height_B+2*Height)
He just distributed the .5. As for the height function, you and I both agree that it's incorrect. I was just looking at your quote for the eqate of "function," which he never uses.

 

[MoMeanMugs]

OP, are you studying to be an engineer? Although I don't recall much, if anything, from the numerical methods class I took many eons ago, it was really interesting. FWIW, I think TI calculators use either Simpson's 1/3 Rule, or possibly Simpson's 2/3 Rule. My preference is Newton, when I need it. Although, computers do everything for me.

 

[QuixoticOne]

Originally posted by: MoMeanMugs
He's multiplying by .5 instead of doing it in the trapezoid formula. The trapezoid formula is:
area=delta_X/2*(Height_A+Height_B+2*Height)
He just distributed the .5. As for the height function, you and I both agree that it's incorrect. I was just looking at your quote for the eqate of "function," which he never uses.

Agreed. Actually I shouldn't have implied one set was 'correct' and the other 'incorrect' since there's really no way of telling which of the two forms are correct, only that they're inconsistent implmentations of f(x) due to one involving -(x^2) and the other involving (-x)^2.

"Numerical Recipes in C" which is a book that can be read freely online or gotten hardcopy is handy for learning about basic integration techniques and other numerical methods.

http://www.nr.com/
http://www.fizyka.umk.pl/nrbook/bookcpdf.html

Now if I could just solve my bug with calculating derivative points around singularities due to a domain mapping relevant to a 3-dimensional integration....

 

[MoMeanMugs]

You can't put in some continuity test?

 

[QuixoticOne]

Originally posted by: MoMeanMugs
You can't put in some continuity test?

Well I know WHERE the singularities are, though I'm not quite sure what value to attribute to the functions at the points of singularity so that the precision of the calculation (basically a polynomial fitting to data points) that comes out of the sampled data points of the functions is suitable.

Maybe I just need to find a series expansion to calculate that works around the signularities instead of trying to evaluate the values of the singular functions themselves. I think there's somwhere I should be clever / cheating when I'm being literal which doesn't give very usefully precise results.

 

[MoMeanMugs]

I leave that up to you. Can you try approaching it from the left and right once it's detected using small intervals? Personally, I hate programming for anything other than playing around with it.

 

[jai6638]

Originally posted by: QuixoticOne

Originally posted by: MoMeanMugs
QuixoticOne,
I'm pretty sure he wants -x^2 and not (-x)^2. At least that's what I gather from the form he's given. Also, how are his equations for Height_A and Height_B incorrect? Last time I checked, 1/2=.5.

..........

So either way he means -x^2 or (-x)^2 he's using one form in a couple of places and the OTHER in a couple of OTHER places, hence my suggestion that his main problem was inconsistent multiple definitions of f(x) and incorrect definitions in half of them.

I mean -x^2 and not (-x)^(2).

Actually, I realized that i wasnt using the variable function so as per your advice, I put that in... Here is what the code looks like now:

delta_X=(limit_B-limit_A)/number_Partitions;
function=exp(pow(j,2.0)/-2.0);
x=(limit_A+delta_X);

for(j=x;j<2.0;j=j+delta_X)
{

height=height+function;


}




double height_A=(0.5)*(exp(-1.0/2.0));
double height_B=(0.5)*(exp(-2.0));

The height_A and height_B equatoins you have written are assuming that its (-x)^(2) which it isnt. Hence, e will be raised to the power of the negative of 0.5 and the 2 as opposed to just positive 0.5 and 2.

Is there some alternative logic that you guys can suggest with which I could do away with this for loop which is not working? the height it posts, as No mean Mugs said, is wrong

No mean Mugs, I am an engineering major and am doing this for a programming class.. god I hate programming. lol

 

[jai6638]

Ok nevermind. got my roommate to debug it too and it seems that the j=0.0 was screwing everything up. I took that out and voila, it worked!