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.
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} \)
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 6 of the polynomials are plotted below.
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].
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.
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 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.