trapzd returning half of answer?


CWenger
07-31-2004, 01:35 PM
I must be missing something obvious (remember it's been 2 years since I finished my college calculus ;)), but the trapzd numerical integration function always seems to return half of the true answer. I am just trying simple functions, such as y=x, for testing purposes. I have retyped the entire funtion twice and I am sure it does not have any typos. Here is an example code of calling the function, which is so simple I don't think I could be screwing it up:

int main()
{
printf("x,{x,0,5}=%f\n",trapzd(yeqx,0,5,10));

return(0);
}

float yeqx(float x)
{
return(x);
}

The function returns 6.25, although it should be 25/2=12.5.

CWenger
08-02-2004, 09:38 AM
Now I'm really confused. It seems that while the trapzd function only returns half of the true answer, when you use the qtrap function, it gives the correct answer. I am new to Numerical Recipes, although I am fairly advanced with C, and I can't figure out why this would be the case.

pagastra
07-02-2005, 05:26 PM
Cwenger, I think you're right.
I just joined this forum, so I'm a year late. I think you're right, trapzd is only returning half of the value. I have tried it with different functions, and the result is still the same. I think the 0.5 should not be there.
peace up!
ps. i'm using fortran and i know next to nothing about C ;)

pagastra
07-06-2005, 08:24 PM
After tinkering with trapzd.for long enough, the 0.5 in
s = 0.5*(s+(b-a)*sum/tnm)
is correct. When used in qtrap, the iteration is performed on the degree of refinement, not simply on the number of time the program is called. So, if degree of refinement is N=3, then s is accumulated over N=1..3.
good luck!