Trapezoidal Riemann Sum in C?

Question

I've been trying to do Riemann Sums to approximate integrals in C. In my code below, I'm trying to approximate by both the trapezoidal way and the rectangular way (The trapezoidal way should be better, obviously). I tried making an algorithm for this on paper, and I got the following: NOTE: N is the number of rectangles (or trapezoids) and dx is calculated using a, b and N ( dx = (b-a)/N ). f(x) = x^2 Rectangular Method: <!-- begin snippet: js hide: false console: true babel: false --> <!-- language: lang-html --> <a href="http://www.codecogs.com/eqnedit.php?latex=\\int_a^b&space;x^2&space;dx&space;\\approx&space;\\sum_{i=1}^N&space;f(a&space;&plus;&space;(i-1)dx)dx" target="_blank"><img src="http://latex.codecogs.com/png.latex?\\int_a^b&space;x^2&space;dx&space;\\approx&space;\\sum_{i=1}^N&space;f(a&space;&plus;&space;(i-1)dx)dx" title="\\int_a^b x^2 dx \\approx \\sum_{i=1}^N f(a + (i-1)dx)dx" /></a> <!-- end snippet --> Trapezoidal Method: <!-- begin snippet: js hide: false console: true babel: false --> <!-- language: lang-html --> <a href="http://www.codecogs.com/eqnedit.php?latex=\\int_a^b&space;x^2&space;dx&space;\\approx&space;\\sum_{i=1}^N&space;[f(a&space;&plus;&space;(i-1)dx)&space;&plus;&space;f(a&space;&plus;&space;i\\cdot&space;dx)]dx" target="_blank"><img src="http://latex.codecogs.com/png.latex?\\int_a^b&space;x^2&space;dx&space;\\approx&space;\\sum_{i=1}^N&space;[f(a&space;&plus;&space;(i-1)dx)&space;&plus;&space;f(a&space;&plus;&space;i\\cdot&space;dx)]dx" title="\\int_a^b x^2 dx \\approx \\sum_{i=1}^N [f(a + (i-1)dx) + f(a + i\\cdot dx)]dx" /></a> <!-- end snippet --> Code (In the following code, f(x)=x^2 and F(x) is it's antiderivative (x^3/3): int main() {   int no_of_rects;   double a, b;   printf("Number of subdivisions = ");   scanf("%d", &no_of_rects);   printf("a = ");   scanf("%lf", &a);   printf("b = ");   scanf("%lf", &b);   double dx = (b-a)/no_of_rects;   double rectangular_riemann_sum = 0;   int i;   for (i=1;i<=no_of_rects;i++) {    rectangular_riemann_sum += (f(a + (i-1)*dx)*dx);   }   double trapezoidal_riemann_sum = 0;   int j;   for (j=1;j<=no_of_rects;j++) {    trapezoidal_riemann_sum += (1/2)*(dx)*(f(a + (j-1)*dx) + f(a + j*dx));    printf("trapezoidal_riemann_sum: %lf\\n", trapezoidal_riemann_sum);   }   double exact_integral = F(b) - F(a);   double rect_error = exact_integral - rectangular_riemann_sum;   double trap_error = exact_integral - trapezoidal_riemann_sum;   printf("\\n\\nExact Integral: %lf", exact_integral);   printf("\\nRectangular Riemann Sum: %lf", rectangular_riemann_sum);   printf("\\nTrapezoidal Riemann Sum: %lf", trapezoidal_riemann_sum);   printf("\\n\\nRectangular Error: %lf", rect_error);   printf("\\nTrapezoidal Error: %lf\\n", trap_error);   return 0;  } Where:  double f(double x) {   return x*x;  }  double F(double x) {   return x*x*x/3;  } I have included the math and stdio header files. What is happening is that the rectangular riemann sum is okay, but the trapezoidal riemann sum is always 0 for some reason. What is the problem? Is it something in my formulas? Or my code? (I am a newbie in C by the way) Thanks in advance.

Answer 1

THe bug is in this statement containing the one half:

trapezoidal_riemann_sum += (1/2)*(dx)*(f(a + (j-1)*dx) + f(a + j*dx)); 

integer 1 divided by integer 2 gives a quotient of zero, which is why they're getting blown away.

change it to: 1.0f/2

that will force a typecase to float instead of doing integer division.

Don't worry, I've gotten burned on that very same bug myself. It is very easy to do.....

I got a better approximation with error -0.45............

I also split up the function calls as y1 and y2;

Here's the listing:

==================================================

//Code (In the following code, f(x)=x^2 and F(x) is it's antiderivative (x^3/3)

#include <stdio.h>

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

double F(double x) { return x*x*x/3; }

int main()

{

int no_of_rects;

double a, b;

printf("Number of subdivisions = ");

scanf("%d", &no_of_rects);

printf("a = "); scanf("%lf", &a);

printf("b = "); scanf("%lf", &b);

double dx = (b-a)/no_of_rects;

printf(" dx = %lf \n",dx);

double rectangular_riemann_sum = 0;

int i;

for (i=1;i<=no_of_rects;i++)

{

rectangular_riemann_sum += (f(a + (i-1)*dx)*dx);

}

double trapezoidal_riemann_sum = 0;

int j;

for (j=1;j<=no_of_rects;j++)

{

double y1 = f(a + (j-1)*dx);

double y2 = f(a + j*dx);

printf(" y1 = %12.6lf and y2 = %12.6lf \n",y1,y2);

trapezoidal_riemann_sum += (1.0f/2)*(dx)*( y1+y2 ); //HERE"S THE FIX !!!!

printf("trapezoidal_riemann_sum: %12.6lf\n", trapezoidal_riemann_sum);

}

double exact_integral = F(b) - F(a);

double rect_error = exact_integral - rectangular_riemann_sum;

double trap_error = exact_integral - trapezoidal_riemann_sum;

printf("\n\n Exact Integral: %lf", exact_integral);

printf("\n Rectangular Riemann Sum: %lf", rectangular_riemann_sum);

printf("\nTrapezoidal Riemann Sum: %lf", trapezoidal_riemann_sum);

printf("\n \n Rectangular Error: %lf", rect_error);

printf("\nTrapezoidal Error: %lf\n", trap_error);

return 0;

}