Determining Roots of Legendre Polynomials and their weights for Gaussian Quadrature – C PROGRAM

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 x 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 n , 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 n , then you can refer to the program below.
The below program calculates and returns the roots & weights for a given n .

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 [-1,1]
Then weights are calculated by integrating the Lagrange interpolation terms from -1 to 1:
\boxed{C_i=\int_{-1}^{1}Li(x)dx}
\boxed{L_i(x)=\prod^n_{j=0;j \neq i}\frac{x-x_j}{x_i-x_j}}
where x_i is the ith root of the Legendre polynomial, and n 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:

Test 1 for n=3

Test 2 for n=5

Test 3 for n=12

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.

PhD researcher at Friedrich-Schiller University Jena, Germany. I'm a physicist specializing in theoretical, computational and experimental condensed matter physics. I like to develop Physics related apps and softwares from time to time. Can code in most of the popular languages. Like to share my knowledge in Physics and applications using this Blog and a YouTube channel.



Leave a Reply

Your email address will not be published. Required fields are marked *