Gaussian Legendre quadrature/ Gauss-Legendre quadrature is a numerical technique used to calculate the definite integral of a function. This is done by evaluating the function at some specific values of given by the roots of the Legendre polynomials, and then multiplying that by the weight of that root.
The weight calculation is a little complicated involving an integration step.
This means that we would need to use an already existing numerical integration technique to be able to calculate the weights, which are then again going to be used for numerical integration. This might seem stupid/weird. But usually, what is done is that, the weights, and roots are calculated once and then stored for future use.
On the internet you could find these weights and roots for a large value of , say 100. You could then just use those values to perform an integration using Gauss-Quadrature.
However, if you wanted to calculate the weights and roots for a higher , then you can refer to the program below.
The below program calculates and returns the roots & weights for a given .
The program uses several concepts that I have discussed and wrote about in the last few posts, like Simpson’s 1/3rd integration, Lagrange interpolation, Bisection method, Recursion relations, etc.
The program uses recursion relation to calculate the value of the nth order Legendre polynomial.
Then finds the root using bisection method within the interval
Then weights are calculated by integrating the Lagrange interpolation terms from -1 to 1:
where is the ith root of the Legendre polynomial, and is the total number of roots.
CODE:
/************************************************************************** ******Find out the roots & weights of Gauss-Legendre Quadrature for given n ***************************************************************************/ #include<stdio.h> #include<math.h> /*Legendre Polynomial P0(x)*/ double P0(double x){ return 1; } /*Legendre Polynomial P1(x)*/ double P1(double x){ return x; } /*Nth Legendre Polynomial Pn(x)*/ double Pn(int n, double x){ if(n==0){ return P0(x); }else if(n==1){ return P1(x); }else{ //Use the recurrence relation return (double )((2*n-1)*x*Pn(n-1,x)-(n-1)*Pn(n-2,x))/n; } } /*Lagrange terms*/ double Li(int n, double x[n+1], int i, double X){ int j; double prod=1; for(j=0;j<=n;j++){ if (j!=i){ prod=prod*(X-x[j])/(x[i]-x[j]); } } return prod; } /*Function definition to perform integration by Simpson's 1/3rd Rule */ double Ci(int i, int n, double x[n], double a, double b, int N){ double h,integral,X,sum=0; int j,k; h=(b-a)/N; for(j=1;j<N;j++){ X=a+j*h; if(j%2==0){ sum=sum+2*Li(n-1,x,i,X); } else{ sum=sum+4*Li(n-1,x,i,X);; } } double Fa=Li(n-1,x,i,a);; double Fb=Li(n-1,x,i,b); integral=(h/3.0)*(Fa+Fb+sum); return integral; } /*Function definition for bisection procedure[Returns the root if found or 999 for failure]*/ double bisection(int n,double f(int n,double x),double a, double b, double eps, int maxSteps){ double c; if(f(n,a)*f(n,b)<=0){ int iter=1; /*Bisection Method begins that tabulates the various values at each iteration*/ do{ c=(a+b)/2; if(f(n,a)*f(n,c)>0){ a=c; } else if(f(n,a)*f(n,c)<0){ b=c; } else if(f(n,c)==0){ return c; } iter++; }while(fabs(a-b)>=eps&&iter<=maxSteps); return c; } else{ return 999; } } main(){ int i=0; int n; // order/terms printf("Enter the value of n(data-points):\n"); scanf("%d",&n); //Array to store the roots of Legendre polynomials double xi[n]; //window(Step-size) for bisection method double h=0.01; //dummy variable for bisection method double x; //dummy variable where the root is returned after bisection routine double root; printf("\n\nThe roots (xi's) are:\n_____________________________________________________\nAccuracy: 10^(-15)\n\n"); for(x=-1.0;x<=1.0;x=x+h){ //set the accuracy to approx. 10^-15 but there is also a limit on maxSteps. (Modify these acc. to your needs) root=bisection(n,Pn,x,x+h,0.0000000000000001,1000000); if(root!=999){ xi[i]=root; printf("x[%d] = %17.16lf\n",i+1,root); i++; } } printf("_____________________________________________________\n"); printf("\n\nThe weights (ci's) are:\n_____________________________________________________\nAccuracy: 10^(-13)\n\n"); for(i=0;i<n;i++){ //(Modify the number of sub-intervals according to your needs) printf("c[%d] = %17.16lf\n",i+1,Ci(i,n,xi,-1,1,1000000)); } }
OUTPUT:
Android Apps:
I’ve also created a few Android apps that perform various numerical methods calculations and can come in handy to those taking a course on Numerical Methods.
Download: https://play.google.com/store/apps/details?id=com.bragitoff.numericalmethods
Download: https://play.google.com/store/apps/details?id=com.bragitoff.matrixcalculator
Download: https://play.google.com/store/apps/details?id=com.bragitoff.lagrangeinterpolatingpolynomial
Download: https://play.google.com/store/apps/details?id=com.bragitoff.polynomialroots
References:
https://pomax.github.io/bezierinfo/legendre-gauss.html
http://mathworld.wolfram.com/Legendre-GaussQuadrature.html
http://keisan.casio.com/exec/system/1329114617
Well, that’s it!
I hope you guys enjoyed this post.
If you have any questions/doubts leave them in the comments section down below.
Ph.D. researcher at Friedrich-Schiller University Jena, Germany. I’m a physicist specializing in computational material science. I write efficient codes for simulating light-matter interactions at atomic scales. I like to develop Physics, DFT, and Machine Learning related apps and software from time to time. Can code in most of the popular languages. I like to share my knowledge in Physics and applications using this Blog and a YouTube channel.