Gauss-Legendre quadrature

Gauss-Legendre quadrature is set of coordinates \(x_n\) and weights \(w_n\), n=1..N, characterized in that, that formula

\(F=\displaystyle \sum_{n=1}^N w_n f(x_n)\)

approximates the integral

\(\Phi= \displaystyle \int_{-1}^1 f(x) \mathrm d x\)

in such a way, that for any polinomail \(f\) of order up to \(2N-1\), the approximation is exact, id est, \(F=\Phi\).

Also, the same term Gauss-Legendre quadrature may be used for the modification of the formulas above for arbitrary interval of integration.

Example in C++

One example adopted from the Numerical recipes in C is shown below:

#include <stdio.h>
#include <math.h>
#include<stdlib.h>
#define DB double
#define EPS 1.0e-18
// Be careful: numeration begins with UNITY !!!!
void gauleg(long double x1, long double x2, long double x[], long double w[], int n)
{ int m,j,i; long double z1,z,xm,xl,pp,p3,p2,p1;
 m=(n+1)/2;
 xm=0.5*(x2+x1);
 xl=0.5*(x2-x1);
 for(i=1;i<=m;i++)
 { z=cos(M_PI*(i-0.25)/(n+0.5));
   do{ p1=1.0;
        p2=0.0;
        for(j=1;j<=n;j++)
        { p3=p2; p2=p1;
          p1=((2.*j-1.)*z*p2-(j-1.)*p3)/j;
        }
        pp=n*(z*p1-p2)/(z*z-1.0);
        z1=z;
        z=z1-p1/pp;
      } while (fabs(z-z1) > EPS);
    x[i]=xm-xl*z;
    x[n+1-i]=xm+xl*z;
    w[i]=2.0*xl/((1.0-z*z)*pp*pp);
    w[n+1-i]=w[i];
  }
}
#undef EPS
//#include "gauleg.cin"
#define NPOINT 8
#define X1 -1.0
#define X2 1.0
int main(void){ int i; long double xx=0.0; long double *x,*w;
x=(long double *)malloc((size_t)(NPOINT+1)*sizeof(long double));
w=(long double *)malloc((size_t)(NPOINT+1)*sizeof(long double));
gauleg(X1,X2,x,w,NPOINT);
printf("\n%2s %10s %12s\n","#","x[i]","w[i]");
for (i=1;i<=NPOINT;i++) { printf("%4d ",i);
                                printf("%18.16Lf ", x[i]);
                                printf("%18.16Lf\n",w[i]);
                        }
FILE *o; o=fopen("GLxw8.inc","w");
fprintf(o,"int NPO=%3d;\n",NPOINT);
fprintf(o,"DB GLx[%3d]={\n",NPOINT);
for (i=1;i<=NPOINT;i++)
{ if(i>1) fprintf(o,","); else fprintf(o," ");
        fprintf(o,"%24.21Lf\n",x[i]);}
fprintf(o,"};\n");
fprintf(o,"DB GLw[%3d]={\n",NPOINT);
for (i=1;i<=NPOINT;i++)
{ if(i>1) fprintf(o,","); else fprintf(o," ");
        fprintf(o,"%24.22Lf\n",w[i]);}
fprintf(o,"};\n");
fclose(o);
//small test:
        for (i=1;i<=NPOINT;i++)
                xx += (w[i]*cos(x[i]));
        printf("from GAULEG: %32.30Lf\n",xx);
        printf("should be : %32.30f\n", 2*sin(1.));
        printf("diff : % 9.3Le\n", xx-2*sin(1.));
        printf("failure in last diguits? Why are not they different?\n");
       free((char*)(w));
       free((char*)(x));
        return 0;
}

More advanced service to evaluate the coordinates and weights is suggested by John Burghardt

Keywords

John Burghardt Numerical recipes in C Gauss quadrature