The Newton-Cotes family of quadrature formulas are characterized by the  n+1 nodes  (x_0,x_1,\hdots,x_n) beign equally spaced. Hence, their quadrature formulas were easy to construct. But, it also limits the accuracy of the formula. When  n is odd, the algebraic accuracy is  n . When  n is even, the algebraic accuracy is  n+1 . It can be shown that the algebraic accuracy can be up to  2n+1 , if we choose the nodes right. 

Gauss quadrature vs trapezoidal rule
Gauss quadrature vs trapezoidal rule

Consider an exmaple to integrate the following polynomial,

 f(x)=7x^3-8x^2-3x+3

The integral value of the polynomial is

\displaystyle \int_{-1}^{1}(7x^3-8x^2-3x+3)dx=\frac{2}{3}
Gauss quadrature vs trapezoidal rule

Using the trapezoidal rule, the orange ine approximates the polynomial and gives  f(-1)+f(1)=-10

The 2-point gauss quadrature rule uses the dashed-black line to approximate the polynomial. It gives the exact same result. teh green region has the same area size as the sum of the red regions.

 \displaystyle f\Bigg(-\sqrt{\frac{1}{3}}\Bigg)+f\Bigg(\sqrt{\frac{1}{3}}\Bigg)=\frac{2}{3}

In general, the gauss quadrature formula on an integral of this particular form

 \displaystyle \int_{a}^{b}\rho(x)f(x)dx\approx\sum_{k=0}^{n}w_kf(x_k)

can achieve the algebraic accuracy of atleast  n and upto  2n+1 , if we can choose the right set of nodes. The nodes  \{x_k\} are called the Gauss points and the coefficients  \{w_k\} are called Gauss coefficients. Both the nodes and the coefficients do not depend on the integrand function. the nodes  \{x_k\} are computed from the roots of a polynomial, depending on the form of the integrand or  \rho(x) .

Gauss-Legendre Quadrature Formula

The Gauss-Legendre Quadrature formula olves the integral of the following form for [late] \rho=1 [/latex].

 \displaystyle \int_{-1}^{1}\rho(x)f(x)dx\approx\sum_{k=0}^{n}w_kf(x_k)

The Gauss nodes are the roots of the Legendre polynomial  P_n(x) .

 P_n(x) is a polynomial of degree  n such that they are orthogonal. Two polynomials or functions are said to be orthogonal if and only if they satisfy the following condition,

 \displaystyle \int_{-1}^{1}P_m(x)P_n(x)dx=0,\,\, if\, n\neq m

The first 11 Legendre polynomials (for  n=0,\hdots,10 ) are shown below.

the first 11 Legendre polynomials

The first 6 of the polynomials are plotted below.

the first 6 Legendre polynomials

The weights  \{w_k\} are given as,

 \displaystyle w_k=\frac{2}{(1-x^2_k)[P'_n(x_k)]^2}

Some of the low-order quadrature rules are tabulated below over the interval [-1,1].

low-order Gauss-Legendre quadrature rules

With the Gauss nodes and Gauss weights, the Gauss-Legendre quadrature formula is given as,

 \displaystyle \int_{-1}^{1}f(x)dx\approx\sum_{k=0}^{n}w_kf(x_k)

For the functions that are polynomials and with degree upto  2n+1, the equation is exact.

The following code computes the following integral.

 \displaystyle I=\int_{-1}^{1}(x^4+x^2+x)dx=2
				
					val f: final UnivariateRealFunction f = new AbstractUnivariateRealFunction() {
override public double evaluate(double x) {
  return 4 * x * x * x + 2 * x + 1; // x^4 + x^2 + x
  }
};

// the integrators
Integrator integrator1 = new Trapezoidal(1e-8, 20); // using the trapezoidal rule
Integrator integrator2 = new Simpson(1e-8, 20); // using the Simpson rule
Integrator integrator3 = new GaussLegendreQuadrature(2);

// the limits
double a = -1, b = 1;

// the integrations
double I1 = integrator1.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the trapezoidal rule", a, b, I1));
double I2 = integrator2.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the Simpson rule", a, b, I2));
double I3 = integrator3.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the Gauss Legendre Quadrature", a, b, I3));
				
			

The output is:

				
					S_[-1,1] f(x) dx = 2.0000000000000000, using the trapezoidal rule
S_[-1,1] f(x) dx = 2.0000000000000000, using the Simpson rule
S_[-1,1] f(x) dx = 2.0000000000000000, using the Gauss Legendre Quadrature
				
			

Gauss-Laguerre Quadrature Formula

The Gauss-Laguerre Quadrature formula numerically computes the integrals of the following form,

 \displaystyle \int_{0}^{\infty}e^{-x}f(x)dx

The gauss nodes,  \{x_k\} , re the roots of the Laguerre polynomials  L_n(x) .

Below is the list of the first few Laguerre polynomials.

Laguerre polynomials

The Gauss coefficients,  \{w_k\} , are given by

 \displaystyle w_k=\frac{x_k}{(n+1)^2[L_{n+1}(x_k)]^2}

The following code evaluates this particular integral.

 \displaystyle I=\int_{0}^{\infty}e^{-2}(x^2+2x+1)dx=5
				
					final Polynomial poly = new Polynomial(1, 2, 1); // x^2 + 2x + 1
val f: final UnivariateRealFunction f = new AbstractUnivariateRealFunction() {
override public double evaluate(double x) {
  return exp(-x) * poly.evaluate(x); // e^-x * (x^2 + 2x + 1)
  }
};

// the integrators
Integrator integrator1 = new Trapezoidal(1e-8, 20); // using the trapezoidal rule
Integrator integrator2 = new Simpson(1e-8, 20); // using the Simpson rule
Integrator integrator3 = new GaussLaguerreQuadrature(2, 1e-8);

// the limits
double a = 0, b = Double.POSITIVE_INFINITY;

// the integrations
double I1 = integrator1.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the trapezoidal rule", a, b, I1));
double I2 = integrator2.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the Simpson rule", a, b, I2));
double I3 = integrator3.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the Gauss Laguerre Quadrature", a, b, I3));
				
			

The output is:

				
					S_[0,infinity] f(x) dx = NaN, using the trapezoidal rule
S_[0,infinity] f(x) dx = NaN, using the trapezoidal rule
S_[0,Infinity] f(x) dx = 5.0000000000000000, using the Gauss Laguerre Quadrature
				
			

It is interesting to note that neither the trapezoidal or the Simpson’s rule can compute this improper integral of which on of the limit is positive infinity. We can’ t divide the half real life  [0,\infty] into partitions.

Gauss-Hermite Quadrature Formula

The Gauss-Hermite Quadrature formula numerically computes the integrals in the following form,

 \displaystyle \int_{-\infty}^{\infty}e^{-x^2}f(x)dx

The Gauss nodes,  \{x_k\} , are the roots of the Hermite polynomials  H_n(x) .

Below is the list of the first few Hermite polynomials.

the first 10 Hermite polynomials
the first 6 Hermite polynomials

The Gauss coefficients,  \{w_k\} , are given by,

 \displaystyle w_k=\frac{2^{n-1}n!\sqrt{\pi}}{n^2[H_{n-1}(x_k)]^2}

The following code computes the integral.

 \displaystyle I=\int_{0}^{\infty}e^{-x^2}(x^2+2x+1)dx=\frac{3\sqrt(\pi)}{2}\approx2.658681
				
					final Polynomial poly = new Polynomial(1, 2, 1); // x^2 + 2x + 1
val f: final UnivariateRealFunction f = new AbstractUnivariateRealFunction() {
override public double evaluate(double x) {
  return exp(-(x * x)) * poly.evaluate(x); // e^(-x^2) * (x^2 + 2x + 1)
  }
};

// the integrators
Integrator integrator1 = new Trapezoidal(1e-8, 20); // using the trapezoidal rule
Integrator integrator2 = new Simpson(1e-8, 20); // using the Simpson rule
Integrator integrator3 = new GaussHermiteQuadrature(2);

// the limits
double a = Double.NEGATIVE_INFINITY, b = Double.POSITIVE_INFINITY;

// the integrations
double I1 = integrator1.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the trapezoidal rule", a, b, I1));
double I2 = integrator2.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the Simpson rule", a, b, I2));
double I3 = integrator3.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the Gauss Hermite Quadrature", a, b, I3));
				
			

The output is:

				
					S_[-Infinity,infinity] f(x) dx = NaN, using the trapezoidal rule
S_[-Infinity,infinity] f(x) dx = NaN, using the trapezoidal rule
S_[-Infinity,Infinity] f(x) dx = 2.6586807763582730, using the Gauss Hermite Quadrature
				
			

Note that neither the trapezoidal or the Simpson’s rule can compute this improper integral of which both limits are infinity. We can’ t divide the half real life  [-\infty,\infty] into partitions.

Gauss-Chebyshev Quadrature Formula

The Gauss-Chebyshev quadrature ormula numerically computes the integrals in the following form,

 \displaystyle \int_{-1}^{1}\frac{f(x)}{\sqrt{1-x^2}}dx

The Gauss nodes,  \{x_k\} are given by,

 \displaystyle x_k=cos\Big(\frac{2i-1}{2n}\pi\Big)

The Gauss coefficients,  \{w_k\} , are given by

 \displaystyle w_k=\frac{\pi}{n}

The following code computes this integral.

 \displaystyle I=\int_{-1}^{1}\frac{x^2+2x+1}{\sqrt{1-x^2}}dx=\frac{3\pi}{2}\approx4.712389

 

				
					final Polynomial poly = new Polynomial(1, 2, 1); // x^2 + 2x + 1
val f: final UnivariateRealFunction f = new AbstractUnivariateRealFunction() {
override public double evaluate(double x) {
  return poly.evaluate(x) / sqrt(1 - x* x);
  }
};

// the integrators
Integrator integrator1 = new Trapezoidal(1e-8, 20); // using the trapezoidal rule
Integrator integrator2 = new Simpson(1e-8, 20); // using the Simpson rule
Integrator integrator3 = new GaussChebyshevQuadrature(2);

// the limits
double a = -1, b = 1;

// the integrations
double I1 = integrator1.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the trapezoidal rule", a, b, I1));
double I2 = integrator2.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the Simpson rule", a, b, I2));
double I3 = integrator3.integrate(f, a, b);
System.out.println(String.format("S_[%.0f,%.0f] f(x) x = %.16f, using the Gauss Chebyshev Quadrature", a, b, I3));
				
			

The output is:

				
					S_[-1,1] f(x) dx = NaN, using the trapezoidal rule
S_[-1,1] f(x) dx = NaN, using the trapezoidal rule
S_[-1,1] f(x) dx = 4.7123889803846900, using the Gauss Chebyshev Quadrature
				
			

Note that neither trapezoidal or Simpson’s rule can compute this integral.