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 [1]
