Rounding error in c integration program

Question

im making a program that finds the integral of a function for different values of a power n, and constant a. my program seems to be working correctly but im getting a small rounding error in my results and i cant figure out why. i know i have an error as a friend of mine is also making the same program and his results are slightly different to mine, and his a definitely the right ones as doing the integration on a calculator gives values much closer to his. below are my results and his for a=2 and n=1.

His result: 0.189070
my result: 0.189053

ive tried going through and casting just about everything i can think of but still cant work out where im getting my error from, any help in pointing out where im being an idiot would be greatly appreciated! :p

My Program:

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

#define debug 0
#define N (double)10000

double Integrand(double x, int a, int n);
double Integral(double *x, double dx, int a, int n);

int main (int argc, char* argv[])
{
    int j,a,n=0,count=0,size=(int)N;
    double dx=1/N, x[size];

    sscanf(argv[1], "%d", &a);
    for(j=0;j<N;j++) {
        x[j]=(double)(j)*dx;
    }
    for(n=1;n<=10;n++) {
        printf("n is %d integral is %lf\n",n,Integral(x,dx,a,n));
    }    
    return(EXIT_SUCCESS);
}

double Integral(double *x, double dx, int a, int n)
{
    int i;
    double result=0;

    for(i=0;i<N;i++) {
       result +=(double)((Integrand((double)x[i],a,n))*dx);
    }
    return(result);
}

double Integrand(double x, int a, int n)
{
   double result;
   result=(double)(((pow(x,(double)n))/(x+(double)a)));
   return(result);
}

Answer 1

It's not a rounding error, you just don't pick the best choice for integration points. Change the initialisation to

x[j]=(j+0.5)*dx;

so that you take the midpoint of each integration strip to calculate the value of the integrand. If you always take the left or right endpoint, you will get a systematically too large error for monotonic functions.

If you approximate the integral of a sufficiently smooth function f by a Riemann sum,

 b           n
 ∫ f(x) dx ≈ ∑ f(y_k)*(b-a)/n
 a          k=1

the choice of y_k in the interval [x_(k-1), x_k] = [a+(k-1)*(b-a)/n, a+k*(b-a)/n] influences the error and the speed of convergence. Writing

f(x) = f(y_k) + f'(y_k)*(x-y_k) + 1/2*f''(y_k)*(x-y_k)² + O((x-y_k)³)

in that interval, you find that

x_k                                   x_k                            x_k
 ∫ f(x) dx = f(y_k)*(b-a)/n + f'(y_k)* ∫ (x-y_k) dx + 1/2*f''(y_k) * ∫ (x-y_k)² dx + O(1/n^4)
x_(k-1)                               x_(k-1)                       x_(k-1)

           = f(y_k)*(b-a)/n + 1/2*f'(y_k)*(b-a)/n*((x_k-y_k)-(y_k-x_(k-1))) + O(1/n³)

and the first and largest error term with respect to the approximation f(y_k)*(b-a)/n vanishes for

y_k = (x_k + x_(k-1))/2

giving you an overall O(1/n³) error for that strip, and a total O(1/n²) error for the entire Riemann sum.

If you choose y_k = x_(k-1) (or y_k = x_k), the first error term becomes

±1/2*f'(y_k)*[(b-a)/n]²

leading to an O(1/n) total error.

Answer 2

At a Linux terminal prompt, type:

man fegetenv