Integral V (método do trapézio extendido)

#include <stdio.h>
#include <math.h>

#define del 1e-5

double f(double x)
{
    return 1/x;
}

int main(void)
{
    int N = 10;
    double h, I, IOld, eps;

    // ------------------------- I_1 -------------------------

    int i;
    double soma=0.0f,a=0.1;

    h=(1.0-a)/N;

    for (i=0;i<N;i++)
    {
        double x1 = a + h*i, x2= a + h*(i+1);
        soma += (h/2)*(f(x1) + f(x2));
    }

    I = soma;

    // ------------------------- END -------------------------

    do {
        IOld = I;
        N = 2*N;
        h = h/2;

        double soma = 0;
        for(i = 1; i < N; i += 2) {
            soma += f(a + i*h);
        }

        I = 0.5*IOld + h*soma;

        eps = 1.0/3 * fabs(I - IOld);

        printf("%d %e\n", N, eps);
    } while (eps > del);

    return 0;
}

// erro = 1/12 * h^2 * [f'(a) - f'(b)]