Search code examples
cintegralcalculus

Trapezoidal Riemann Sum in C

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:

<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>

Trapezoidal Method: 

<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>

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.

Solution

  • In this statement:

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

    1/2 == zero, so the whole statement is zero. Change at least the numerator, or the denominator to the form of a double to get a double value back. i.e. 1/2.0 or 1.0/2 or 1.0/2.0 will all work.