Derivation of the Gaussian Quadrature #

The Gaussian quadrature method is the optimal numerical integration method for integrating functions of the form

It is ideal for cases where the integrand is composed of a specific weight function multiplied by a polynomial.
Although the positions of the points used for numerical integration are restricted, by referring to only $N$ points, it is possible to exactly integrate an integrand of the form (specific weight function) × (polynomial of degree $2N-1$).

The Gaussian quadrature method is especially effective in the field of scientific and technical computation because it produces highly accurate results with a small number of points. This article will discuss how the Gaussian quadrature method is derived.

1. Problem Setting #

Consider the numerical integration of the product of a weight function $\omega(x)$ defined over the interval $x=[a,b]$ and a function $f(x)$, which is a polynomial of degree $M$ defined over the same interval, using $N$ points.

Let’s express the problem in mathematical terms. Suppose $f(x)$ can be written using $M+1$ known coefficients $a_m$ as

$$ \begin{align}\label{e1} f(x)=\sum_{m=0}^{M}a_m x^m. \end{align} $$

The objective is to solve the integral

$$ \begin{align}\label{e2} \int_{a}^{b}\omega(x)f(x)dx. \end{align} $$

Numerically integrating, computers cannot handle continuous functions, so we attempt to compute using only discrete points of the integrand.

Thus, the problem considered in this article is to find $x_i, \omega_i, (i=1,2,\cdots,N)$ such that

$$ \begin{align} \label{e3} \int_{a}^{b}w(x)f(x)dx=\sum_{i=1}^N \omega_i f(x_i). \end{align} $$

Here, $\omega_i$ are referred to as “weights,” representing the degree of influence of the integrand evaluated at each point $x_i$ on the result of the integration.

2. Gaussian Quadrature #

2.1. Relationship Between the Number of Discrete Points and the Achievable Degree #

For example, if it were possible to reference an infinite number of values of $f(x)$, one could perform an exact integration. Similarly, if $f(x)$ is a linear function, knowing the values at the starting point ($x=a$) and the endpoint ($x=b$) would allow for an exact integration.

This implies that there is a relationship between the degree $M$ of the polynomial $f(x)$ and the number of discrete points $N$, so let’s derive that relationship first.

Substituting Eq. \eqref{e1} into Eq. \eqref{e2} yields

$$ \begin{align} \int_a^b w(x)f(x)dx&=\int_a^b w(x)\sum_{m=0}^M a_m x^m dx, \nonumber\\ &= \sum_{m=0}^N a_m g_m. \label{e4b} \end{align} $$

Here, we have defined

$$ \begin{align} g_m \equiv \int_a^b w(x) x^m dx. \end{align} $$

On the other hand, substituting Eq. \eqref{e1} into the right side of Eq. \eqref{e3}, we get

$$ \begin{align} \sum_{i=1}^N f(x_i)\omega_i &= \sum_{i=1}^N \sum_{m=0}^M a_m x_i^m\omega_i \nonumber\\ &= \sum_{m=0}^M a_m \left(\sum_{i=1}^N \omega_i x_i^m\right) \label{e5b} \end{align} $$

If Eq. \eqref{e3} is to hold exactly for any $a_m$, for each coefficient $m=0, 1, \cdots, M$,

$$ \begin{align} g_m = \sum_{i=1}^N\omega_i x_i^m,~~(m=0, 1,\cdots, M), \end{align} $$

must hold. Since this condition must be true for all $m$, $M+1$ equations can be established.
On the other hand, the unknowns on the right side include $N$ each of $x_i, \omega_i$, thus the system has $2N$ degrees of freedom.

If there exists a combination of coefficients that satisfies the condition for both sides to match, the system’s degrees of freedom $2N$ must exceed the number of equations $M+1$, thus necessitating that

$$ \begin{align} M+1\leq 2N~~\to~~M \leq 2N-1 \end{align} $$

be fulfilled. When aiming to minimize the number of evaluation points, equality in the above equation allows for the exact integration of the desired function.

Therefore, considering when equality holds, we regard

$$ \begin{align} M = 2N-1 \end{align} $$

as the most efficient scenario. This implies that, by using $N$ points, it is possible to find positions $x_i$ and weights $w_i$ that allow for the exact calculation of a polynomial of degree $2N-1$.

This establishes the relationship between the number of points $N$ used to reference the integrand and the highest degree $M$ of the polynomial that can be exactly integrated. Moving forward, we will explore the specific values of $x_i$ and $\omega_i$.

2.2. Positions of the Points #

Now that we have examined the relationship between the number of points $N$ and the degree $M$, no specific values have been determined yet.
Therefore, we will now determine the specific positions $x_i$ and weights $\omega_i$.

To calculate the integral

$$ \begin{align} \label{ed0} \int_a^b w(x)f(x)dx, \end{align} $$

we will use the values of $f(x)$ at discrete points $x=x_i$. In other words,

$$ \begin{align} \label{ed1} (x_1, f_1), ~(x_2, f_2), \cdots , (x_N, f_N) \end{align} $$

is considered, where the points are known. Here, we have defined $f_i \equiv f(x_i)$ (at this point, $x_i$ are not yet determined). The function $w(x)$ is known for any point $x$ due to the orthogonality of the orthogonal polynomials used, which determines it automatically.

From here, many considerations are made in a somewhat top-down manner. We will attempt to explain as much as possible, but note that the explanation may not be complete.

First, we define new functions $F(x)$ and $g(x)$. $F(x)$ is defined as the unique polynomial of degree $N-1$ that passes through all $N$ points, specifically using Lagrange polynomials¹ as

$$ \begin{align} \label{ed2} F(x)\equiv \sum_{k=1}^{N} f_k \left(\prod_{\substack{i=1,\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\right) \end{align} $$

Furthermore, we define the function of the difference between $f(x)$ and $F(x)$ as

$$ \begin{align}\label{ed3} g(x)\equiv w(x) [f(x)-F(x)]. \end{align} $$

Let’s check the degree of each function here:

• $f(x)$: Polynomial of degree $2N-1$
• $F(x)$: Polynomial of degree $N-1$

Since the Lagrange polynomial $F(x)$ satisfies $F(x_i) = f(x_i),~~(i=1,2,\cdots,N)$ at points $x=x_i$, $g(x)$ will have the property that

$$ \begin{align} \label{ed4} g(x_i) = 0, \hspace{1em}(i=1,2,\cdots,N). \end{align} $$

Considering the integral of $g(x)$ over $x$’s domain, we get

$$ \begin{align} \label{ed5} \int_{a}^b g(x)dx = \int_{a}^b w(x)f(x)dx-\int_{a}^b w(x)F(x)dx. \end{align} $$

Since $F(x)$ is a polynomial of degree $M$, if it can be expressed discretely as in Eq. \eqref{e3}, then

$$ \begin{align} \label{ed6} \int_{a}^b g(x)dx = \int_{a}^b w(x)f(x)dx-\sum_{i=1}^N F(x_i)\omega_i. \end{align} $$

Now, since the right side $f(x)$ is a polynomial of degree $2N-1$ and $F(x)$ is expressed using $N$ points, from the previous discussion, it can be shown that there exists a combination of $(x_i, \omega_i)$ that makes Eq. \eqref{ed6} equal to zero.
If we can determine $(x_i, \omega_i)$ such that Eq. \eqref{ed6} becomes zero, then we can precisely execute the integration of $w(x)f(x)$, which was the initial goal.

Given that $g(x)$ is expressed by Eq. \eqref{ed3} and has the properties described in Eq. \eqref{ed4}, for the $N$ points $x_i$, we can rewrite it as

$$ \begin{align} g(x)&\equiv w(x)\cdot(x-x_1)(x-x_2)\cdots(x-x_N)G(x)\nonumber\\ &=w(x)G(x)\prod_{i=1}^N (x-x_i) \end{align} $$

where $G(x)$ is newly defined as a polynomial of degree $N-1$. $(x-x_1)(x-x_2)\cdots(x-x_N)$ is a polynomial of degree $N$, which fits the degree when excluding $w(x)$.

Considering that the integral of $g(x)$ is zero, the integration becomes as follows

$$ \begin{align} \label{c1} \int_{a}^b g(x)dx = \int_{a}^b w(x)\cdot (x-x_1)(x-x_2)\cdots(x-x_N)G(x) dx=0. \end{align} $$

Since $G(x)$ is a polynomial of any degree up to $N-1$, it can be expanded as a linear combination of orthogonal polynomials $p_n(x)$ of up to degree $N-1$. Thus, the $\alpha_j, \beta_j$ are exist which satisfy

$$ \begin{align}\label{c2} G(x)=w(x)\sum_{j=0}^{N-1}\alpha_j x^j=w(x) \sum_{j=0}^{N-1}\beta_j p_j(x). \end{align} $$

Furthermore, is of degree $N$, so it can be written as a linear combination of polynomials up to degree $N$.

Also, since

$$ \begin{align} (x-x_1)(x-x_2)\cdots(x-x_N) \end{align} $$

is of degree $N$, it can be written as a linear combination of polynomials up to degree $N$.

Generally, when expanding an arbitrary polynomial of degree $N$, it is necessary to write it as a linear combination of all orthogonal polynomials of degree $N$ or less. However, let’s choose points $x_1, x_2, \cdots, x_N$ such that they satisfy $p_N(x)=0$ for $N$ points which called quadrature points. By doing this, it can be described solely using the orthogonal polynomial of degree $N$, and a constant $c$ can be found from

$$ \begin{align}\label{c3} (x-x_1)(x-x_2)\cdots(x-x_N)=c p_N(x). \end{align} $$

Although the value of $c$ can be determined, it is not necessary to specify its exact value as it will be canceled out later. The positions of the points are determined from this condition.

Substituting equations \eqref{c2}, \eqref{c3} into \eqref{c1} yields

$$ \begin{align} \int_a^b g(x)dx &=\int_a^b dx w(x)\cdot (x-x_1)(x-x_2)\cdots(x-x_N)\cdot G(x) \nonumber\\ &=c\int_a^b dx w(x)p_N(x)\cdot \sum_{j=0}^{N-1}\beta_j p_j(x) \nonumber\\ &=c\sum_{j=0}^{N-1}\beta_j\int_a^b dx w(x) p_N(x) p_j(x) \nonumber\\ &=c h_N \sum_{j=0}^{N-1}\beta_j\delta_{N, j} \nonumber\\ &=0 \label{c4} \end{align} $$

where $\delta_{i, j}$ denotes the Kronecker delta. Furthermore, $h_n$ refers to the quantity that appears on the right side of the orthogonality condition for orthogonal polynomials,

$$ \begin{align}\label{qh} \int_a^b dx w(x) p_n(x) p_m(x)=h_n \delta_{n,m}. \end{align} $$

The result of Eq. \eqref{c4} implies that if you choose $N$ points such that $p_N(x)=0$, Eq. \eqref{c1} is automatically satisfied, i.e., $\int_a^b g(x)dx = 0$, meaning that an exact integration of Eq. \eqref{ed0} can be performed at these $N$ points.

Thus, by selecting the zeros of the $N$th degree orthogonal polynomial as the positions for the Gaussian quadrature points, we find that

$$ \begin{gather} \label{e4} \int_a^b g(x)dx = \int_a^b w(x)f(x)dx-\sum_{i=1}^N F(x_i)\omega_i =0 \\ p_N(x)=0 \to x=x_i ,~~(i=1, 2, \cdots, N) \end{gather} $$

can be satisfied.

From the discussion so far, we have determined the positions of the points $x_i$. What remains is to specifically determine the weights $\omega_i$.

2.3. Weights #

From the discussion so far, we have shown that there exists a set of points and weights $x_i, \omega_i$ that can exactly integrate the product of the weight function $w(x)$ and a polynomial $f(x)$ of degree $2N-1$ using $N$ points, as

$$ \begin{align} \int_a^b f(x)dx=\sum_{i=1}^N f(x_i)\omega_i. \end{align} $$

Subsequently, we found that the positions of the points $x_i$ are the zeros of the $N$th degree orthogonal polynomial, that is,

$$ \begin{align} p_{N}(x)=0\to x=x_i ,~~(i=1, 2, \cdots, N). \end{align} $$

Next, we will determine the weights $\omega_i$.

Since we have seen that Eq. \eqref{e4} holds, the integration of a polynomial $f(x)$ of degree $2N-1$ can be replaced with the integration through the zeros of the $N$th degree orthogonal polynomial. Using the fact that Eq. \eqref{ed5} is zero and Eq. \eqref{ed2},

$$ \begin{align} \int_a^bw(x)f(x)dx&= \int_a^b w(x)F(x)dx\nonumber\\ &= \int_a^b w(x)\sum_{k=1}^{N} f_k \Biggl(\prod_{\substack{i=1,\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx \nonumber\\ &= \sum_{k=1}^{N} f_k \int_a^b w(x)\Biggl(\prod_{\substack{i=1,\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx \nonumber\\ &\equiv \sum_{k=1}^{N}f_k \omega_k. \end{align} $$

Thus, the weights can be expressed as

$$ \begin{align} \omega_k=\int_{a}^b w(x)\Biggl(\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx. \label{omega0} \end{align} $$

Since the positions $x_i$ are already determined, Eq. \eqref{omega0} is uniquely defined. However, as the formula involves complex integrals and summations, making it cumbersome to use, we attempt to transform the equation for easier application.

From Eq. \eqref{c3}, we know that

$$ \begin{align}\label{u1} (x-x_1)(x-x_2)\cdots(x-x_N)=\frac{1}{k_N} p_N(x) \end{align} $$

Specifically, for the $k$th equation, it can be written as

$$ \begin{align}\label{uu1} \prod_{\substack{i=1\\ i\ne k}}^N (x-x_i)= \frac{cp_N(x)}{x-x_k} \end{align} $$

Furthermore, differentiating both sides of Eq. \eqref{u1} yields

$$ \begin{align} cp’_N(x) &=\frac{\partial}{\partial x}\Bigl[(x-x_1)(x-x_2)\cdots(x-x_N)\Bigr]\nonumber\\ &=\Bigl[(x-x_2)(x-x_3)\cdots(x-x_N)\nonumber\\ &\hspace{3em}+(x-x_1)(x-x_3)\cdots(x-x_N)\nonumber\\ &\hspace{3em}+(x-x_1)(x-x_2)\cdots(x-x_{N-1})\Bigr] \nonumber\\ &=\sum_{j=1}^N\prod_{\substack{i=1\\ i\ne j}}^{N}(x-x_i) \end{align} $$

Especially when $x=x_k$, only the term where $j=k$ remains, so the summation sign disappears, leading to

$$ \begin{align}\label{u2} cp’_N(x_k) = \prod_{\substack{i=1\\ i\ne k}}^{N} (x_k-x_i) \end{align} $$

Combining equations \eqref{uu1} and \eqref{u2}, we can proceed with the calculation as

$$ \begin{align} \prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i} &= \frac{\prod_{\substack{i=1\\ i\ne k}}^{N}{x-x_i}}{\prod_{\substack{i=1\\ i\ne k}}^{N}{x_k-x_i}}\nonumber\\ &= \frac{cp_N(x)}{x-x_k}\cdot\frac{1}{cp’_N(x_k)} \nonumber\\ &=\frac{p_N(x)}{(x-x_k)p’_N(x_k)} \end{align} $$

Therefore, the weight $\omega_k$ can be transformed into

$$ \begin{align} \omega_k&=\int_a^b w(x)\Biggl(\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx \nonumber\\ &=\int_a^b w(x)\frac{p_N(x)}{(x-x_k)p’_N(x_k)} dx \label{w1} \end{align} $$

The integration \eqref{w1} can be advanced using the properties of polynomials.

Though it is somewhat top-down, the right-hand side can be inferred to derive from the Christoffel–Darboux theorem for orthogonal polynomials (such as Legendre polynomials or Laguerre polynomials). The Christoffel–Darboux theorem for an $n$th degree orthogonal polynomial $p_n(x)$ exhibits the following property

$$ \begin{align}\label{CDt} \sum_{j=0}^N \frac{p_j(x)p_j(y)}{h_j} = \frac{k_N}{h_N k_{N+1}}\frac{p_{N+1}(x)p_{N}(y)-p_{N}(x)p_{N+1}(y)}{x-y}. \end{align} $$

Here, $k_n$ denotes the coefficient of $x^n$ in the orthogonal polynomial $p_n(x)$, meaning

$$ \begin{align} p_n(x)=k_n x^n+\cdots \end{align} $$

Furthermore, $h_n$, as defined in Eq. \eqref{qh}, is the coefficient derived from orthogonality.

Considering the case where $y$ is precisely at $x_k$ in the context of Eq. \eqref{CDt}, we obtain the relationship

$$ \begin{align} \sum_{j=0}^{N-1} \frac{p_j(x)p_j(x_k)}{h_j} &= \frac{k_{N-1}}{h_{N-1} k_{N}}\frac{p_{N}(x)p_{N-1}(x_k)-p_{N-1}(x)p_{N}(x_k)}{x-x_k} \nonumber\\ &= \frac{k_{N-1}}{h_{N-1} k_{N}}\frac{p_{N}(x)p_{N-1}(x_k)}{x-x_k} \label{PNx} \end{align} $$

Here, the positions of the points $x_k$ are chosen such that $p_{N}(x_k)=0$, utilizing this condition. Subsequently, rearranging to get the relationship

$$ \begin{align} \frac{p_{N}(x)}{x-x_k}&=\frac{k_{N}h_{N-1} }{k_{N-1}} \frac{1}{p_{N-1}(x_k)}\sum_{j=0}^{N-1} \frac{p_j(x)p_j(x_k)}{h_j} \label{w2} \end{align} $$

Substituting Eq. \eqref{w2} into the formula for weights \eqref{w1},

$$ \begin{align} \omega_k&=\int_{a}^b w(x)\frac{p_N(x)}{(x-x_k)p’_N(x_k)} dx\nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\int_{a}^b dx w(x)\sum_{j=0}^{N-1} \frac{p_j(x)p_j(x_k)}{h_j} \nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\sum_{j=0}^{N-1} \frac{p_j(x_k)}{h_j}\cdot \int_{a}^b dx w(x) p_j(x)\nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\sum_{j=0}^{N-1} \frac{p_j(x_k)}{h_j}\cdot \frac{1}{c_0}\int_{a}^b dx w(x) p_j(x)p_0(x)\nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\sum_{j=0}^{N-1} \frac{p_j(x_k)}{h_j}\cdot \frac{h_0}{c_0}\delta_{j,0}\nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\label{w3} \\ &= \left[\sum_{n=0}^{N-1} \frac{p_n^2(x_k)}{h_n}\right]^{-1}\label{w3x} \end{align} $$

Here, we used the fact that the zeroth degree of orthogonal polynomials is a constant ($p_0(x)=c_0$). Furthermore, in the last transformation of Eq. \eqref{w3x}, the relationship in the limit of $y\to x$ in the Christoffel–Darboux theorem

$$ \begin{align} \sum_{n=0}^{N-1} \frac{[p_n(x)]^2}{h_n} = \frac{k_{N-1}}{h_{N-1} k_{N}}\left[-p_N(x)p’_{N-1}(x)+p_{N-1}(x)p’_N(x)\right] \end{align} $$

and the property of the points $p_N(x_k)=0$ are used to derive the result.

When actually calculating the weights, one would use either Eq. \eqref{w3} or \eqref{w3x}, whichever is more convenient for the given situation.

3. Gauss-Lobatto Quadrature #

The Gauss-Lobatto quadrature method is designed to have points at the endpoints, as an adaptation of the Gaussian quadrature method.

Up to Eq. \eqref{c1}, the process is the same, but from here, the positions of the points are chosen such that the points $x=a$ and $x=b$ are selected.

Generally, when referring to the Gauss-Lobatto quadrature method, it is based on the Gauss-Legendre quadrature method (for the interval $x=[-1,1]$) modified to include the endpoints as points. However, this approach is not limited to Legendre polynomials and can similarly be applied to Jacobi polynomials (Jacobi polynomials)²³.

In this article, we will consider the quadrature method based on Legendre polynomials.

3.1. Point Constraints #

$$ \begin{align} \int_{a}^b g(x)dx = \int_{a}^b w(x)\cdot (x-x_1)(x-x_2)\cdots(x-x_N)G(x) dx=0 \end{align} $$

In previous derivations, we fulfilled Eq. \eqref{c1} by choosing the zeros of the polynomial $(x-x_1)\cdots(x-x_N)$. Now, we modify this approach.

That is, among

$$ \begin{align} (x-x_1)(x-x_2)\cdots(x-x_N) \end{align} $$

we consider two of the points as the endpoints, leading to

$$ \begin{align} (x-x_1)(x-x_2)\cdots(x-x_{N-3})(x-a)(x-b) \end{align} $$

This can be rewritten using the product of a polynomial $q_{N-2}(x)$ of degree $N-2$ and $(x-a)(x-b)$ as

$$ \begin{align} (x-x_1)(x-x_2)\cdots(x-x_{N-3})&(x-a)(x-b)\nonumber \\ &=c q_{N-2}(x)(x-a)(x-b)\label{k3a} \end{align} $$

Here, $c$ is an appropriate constant.

$q_{N-2}(x)$ is a polynomial of degree $N-2$. Due to the constraint that $q_{N}(x)$ must be zero at the endpoints, we lose “2” degrees of freedom for the endpoint positions. Thus, the relationship between the degree of integration $M$ that can be exactly integrated and the number of points $N$ changes to

$$ M=(2N-1) - 2 = 2N-3 $$

Although the positions of the endpoints are fixed, the degrees of freedom for the weights at the endpoints do not decrease, so only subtracting 2 is necessary.

With this decrease in degree, while we have previously treated $G(x)$ as a polynomial of any degree up to $N-1$, we now reduce its degree to any polynomial up to $N-3$. Therefore, we set

$$ \begin{align} G(x)=w(x)\sum_{j=0}^{N-3}\alpha_j x^j=w(x) \sum_{j=0}^{N-3}\beta_j q_j(x) \end{align} $$

Then, the integral can be transformed as

$$ \begin{align} \int_a^b g(x)dx \label{k2r} &=\int_a^b dx w(x)\cdot (x-x_1)(x-x_2)\cdots(x-x_{N-3})(x-a)(x-b)\cdot G(x) \\ &=-c\int_a^b dx [w(x)(x-a)(b-x)] q_{N-2}(x)\cdot \sum_{j=0}^{N-3}\beta_j q_j(x) \nonumber\\ &=-c\sum_{j=0}^{N-3}\beta_j\int_a^b dx w’(x) q_{N-2}(x) q_j(x) \\ \end{align} $$

where we have introduced a new weight function

$$ \begin{align} w’(x) = w(x)(x-a)(b-x) \end{align} $$

If $q_n(x)$ is chosen such that it satisfies the orthogonality condition

$$ \begin{align} \label{k4} \int_a^b dx w’(x) q_n(x) q_m(x)=h_n\delta_{n,m}, \end{align} $$

then

$$ \begin{align} \int_a^b g(x)dx&=c\sum_{j=0}^{N-3} h_n \beta_j \delta_{N-2,j}\\ &=0 \end{align} $$

can be achieved, allowing for the properties of the Gaussian quadrature method to be satisfied.

Now, having set Eq. \eqref{k4} in this way, can we find a function $q_{N-2}(x)$ that satisfies the orthogonality \eqref{k4}? The conditions for $q_{N-2}(x)$ must include:

• $q_{N-2}(x)$: A polynomial of degree $N-2$
• $q_{N-2}(a)<\infty$, $q_{N-2}(b)<\infty$
• Satisfies the orthogonality condition \eqref{k4}

Moving forward, to avoid complexity, we will only consider the case of the Gauss-Lobatto quadrature method based on Legendre polynomials⁴.

3.2. Positions of the quadrature points #

In the case of Legendre polynomials, we have $a=-1, b=1$.

The $N$th degree Legendre polynomial $P_N(x)$ becomes an orthogonal polynomial $P_{N-2}(x)$ with its degree reduced by two. We design the polynomial to satisfy the condition of being zero at the endpoints using the information from the reduced degree.

Normally, in Gaussian quadrature, we search for $N$ zeros of

$$ \begin{align} P_N(x)=0 \hspace{1em}\to\hspace{1em} x=x_j,\hspace{1em}(j=1,\cdots,N) \end{align} $$

to use as the quadrature points. However, in Gauss-Lobatto quadrature, we modify this to the following,

$$ \begin{align} &P_{N}(x)+ u P_{N-1}(x)+v P_{N-2}(x)=0, \nonumber \\ &\hspace{3em}\to\hspace{1em} x=-1, 1, x_j,\hspace{1em}(j=1,\cdots,N-2). \label{k6} \end{align} $$

Here, the constants $u, v$ are determined based on the condition that the equation becomes zero at the endpoints $x=\pm 1$.

After this, from Eq. \eqref{k3a}, we define the $N-2$ degree polynomial as

$$ \begin{align} \label{k7} q_{N-2}(x)=\frac{P_{N}(x)+ u P_{N-1}(x)+v P_{N-2}(x)}{(1+x)(1-x)},~~q_{N-2}(\pm 1)<\infty \end{align} $$

This becomes the polynomial considered in the Gauss-Lobatto quadrature method. Just to clarify, in Eq. \eqref{k7}, because we are dividing the roots of the numerator ($x=\pm1$) by the denominator, the result remains a polynomial, and $q_{N-2}(x)$ is also a polynomial (the operation decreases the degree, which is known as “degree reduction” ⁵).

Finally, let’s verify if $q_{N-2}(x)$ satisfies the orthogonality condition \eqref{k4}.
First, we seek the explicit form of $q_{N-2}(x)$. Given the condition that Eq. \eqref{k6} becomes zero at $x=\pm 1$,

$$ \begin{align} x=1 &:& 1 + u + v &= 0 \\ x=-1 &:& -1 + u - v &= 0 \end{align} $$

we deduce $u=0, v=-1$, utilizing the Legendre polynomial properties $P_n(1)=1, P_n(-1)=(-1)^n$. From this,

$$ \begin{align} \label{k8} q_{N-2}(x)=\frac{P_{N}(x)- P_{N-2}(x)}{(1-x^2)}=-\frac{2N-1}{N(N-1)}P’_{N-1}(x) \end{align} $$

To generalize further, let’s shift the indices and consider

$$ \begin{align} \label{k8a} q_m(x)&=\frac{P_{m+2}(x)- P_{m}(x)}{(1-x^2)} \\ &=-\frac{2m+3}{(m+2)(m+1)}P’_{m+1}(x). \end{align} $$

Then, the orthogonality becomes

$$ \begin{align} \label{k9} &\int_{-1}^{1} (1-x^2)q_n(x)q_m(x)dx \nonumber \\ &=\int_{-1}^{1} \left[P_{n+2}(x)- P_{n}(x)\right]\cdot \left[-\frac{2m+3}{(m+2)(m+1)}P^{\prime}_{m+1}(x)\right]dx \nonumber \\ &=-\frac{2m+3}{(m+2)(m+1)}\Biggl\{ \Bigl[\bigl(P_{n+2}(x)- P_{n}(x)\bigr)P_{m+1}(x)\Bigr]_{x=-1}^{x=1} \nonumber\\ &\hspace{9em} -\int_{-1}^{1}\frac{d}{dx}\Bigl[P_{n+2}(x)- P_{n}(x)\Bigr] \cdot P_{m+1}(x) dx \Biggr\} \nonumber\\ &=\frac{(2m+3)(2n+3)}{(m+2)(m+1)}\int_{-1}^{1}P_{n+1}(x) \cdot P_{m+1}(x) dx \nonumber\\ &=\frac{(2m+3)(2n+3)}{(m+2)(m+1)}\frac{2}{2n+3}\delta_{n,m} \nonumber\\ &=\frac{2(2n+3)}{(n+2)(n+1)}\delta_{n,m} \end{align} $$

where we use the property of Legendre polynomials (Appendix \eqref{ap2}).

3.3. Properties of Polynomials Used in the Gauss-Lobatto Quadrature Method #

Let’s summarize the properties of the polynomial $q_n(x)$ used in the Gauss-Lobatto quadrature method, as discussed so far:

• Relation with Legendre Polynomials

$$ \begin{gather} q_n(x)=\frac{P_{n+2}(x)- P_{n}(x)}{(1-x^2)} =-\frac{2n+3}{(n+2)(n+1)}P’_{n+1}(x) \\ \int_{-1}^{1} (1-x^2)q_n(x)q_m(x)dx = \frac{2(2n+3)}{(n+2)(n+1)}\delta_{n,m} \end{gather} $$

• Symmetry

$$ \begin{gather} q_{n}(-x)=(-1)^n q_n(x) \end{gather} $$

• Special Values

$$ \begin{gather} q_{n}(1)=-\frac{2n+3}{2} \end{gather} $$

$$ \begin{gather} \frac{d}{dx}q_{n}(x)\Bigr|_{x=1}=-\frac{n(n+3)(2n+3)}{8} \end{gather} $$

※See Appendix.

• Three-Term Recurrence Relation

$$ \begin{gather} \frac{n+3}{2n+5}q_{n+1}(x)=xq_n(x)-\frac{n}{2n+1}q_{n-1}(x) \end{gather} $$

• Differentiation

$$ \begin{gather} (1-x^2)q’_{n}(x)=-nxq_n(x)+\frac{n(2n+3)}{2n+1}q_{n-1}(x) \end{gather} $$

• Leading Coefficient

$$ \begin{gather} q_n(x)=-\frac{2n+3}{n+2} \left(\frac{(2n+2)!}{2^{n+1}((n+1)!)^2}\right)x^n+\cdots \end{gather} $$

• Explicit Form of $q_n(x)$

3.4. Weights #

To calculate the explicit form of weights, we compute for the function $c q_{N-2}(x)(x-a)(x-b)$. Therefore, the weights can be calculated as

$$ \begin{align} \omega_k=\int_{a}^b w(x)\Biggl(\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx. \end{align} $$

For convenience, we have set $x_{N-1}=a$ and $x_{N}=b$. At this stage, $w(x)=1$, but to avoid confusion, we express it as follows,

$$ \begin{align} \omega_k=\int_{a}^b \Biggl(\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx.\label{omega1} \end{align} $$

From Eq. \eqref{c3}, we can rewrite

$$ \begin{align}\label{u1x} (x-x_1)(x-x_2)\cdots(x-x_N)=c q_{N-2}(x) (x-x_{N-1})(x-x_{N}) \end{align} $$

Proceeding as before, for the $k$th term, we summarize

$$ \begin{align} \prod_{\substack{i=1\\ i\ne k}}^{N} (x-x_i)= \left\{ \begin{aligned} &\frac{cq_{N-2}(x)}{x-x_k}(x-x_{N-1})(x-x_{N}),&&(k=1,\cdots,N-2) \\ &cq_{N-2}(x)(x-x_{N}),&&(k=N-1) \\ &cq_{N-2}(x)(x-x_{N-1}),&&(k=N) \\ \end{aligned} \right.\label{uu1x} \end{align} $$

Differentiating both sides of Eq. \eqref{u1x} yields

$$ \begin{align} (d/dx \text{ of left-hand side} ) =\sum_{j=1}^N\prod_{\substack{i=1\\ i\ne j}}^{N}(x-x_i) \end{align} $$

so

$$ \begin{align} &cq’_{N-2}(x)(x-x_{N-1})(x-x_{N}) + cq_{N-2}(x)(x-x_{N}) + cq_{N-2}(x)(x-x_{N-1}) \nonumber\\ &\hspace{18em}=\sum_{j=1}^N\prod_{\substack{i=1\\ i\ne j}}^{N}(x-x_i) \end{align} $$

Especially when $x=x_k$, only the term where $j=k$ remains, so the summation disappears, resulting in

$$ \begin{align} \prod_{\substack{i=1\\ i\ne k}}^{N} (x_k-x_i) =\left\{ \begin{aligned} &cq’_{N-2}(x_k)(x_k-x_{N-1})(x_k-x_{N}), &&(k=1,\cdots,N-2) \\ &cq_{N-2}(x_k)(x_k-x_{N}),&&(k=N-1) \\ &cq_{N-2}(x_k)(x_k-x_{N-1}),&&(k=N) \end{aligned} \right.\label{u2x} \end{align} $$

ombining equations \eqref{uu1x} and \eqref{u2x}, we can calculate

$$ \begin{align} &\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i} \nonumber\\ &=\left\{ \begin{aligned} &\frac{q_{N-2}(x)}{(x-x_k)q’_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})(x-x_{N})}{(x_k-x_{N-1})(x_k-x_{N})}, &&(k=1,\cdots,N-2) \\ &\frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N})}{(x_k-x_{N})},&&(k=N-1) \\ &\frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})}{(x_k-x_{N-1})},&&(k=N) \\ \end{aligned} \right. \end{align} $$

Therefore, the weights $\omega_k$ can be calculated as follows:

• For $k=1,2,\cdots,N-2$

  $$ \begin{align} \omega_k=\int_a^b \frac{q_{N-2}(x)}{(x-x_k)q’_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})(x-x_{N})}{(x_k-x_{N-1})(x_k-x_{N})} dx \end{align} $$

• For $k=N-1$ (where $x_{N-1}=-1$)

  $$ \begin{align} \omega_k=\int_a^b \frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N})}{(x_k-x_{N})} dx \end{align} $$

• For $k=N$ (where $x_{N}=1$)

  $$ \begin{align} \omega_k=\int_a^b \frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})}{(x_k-x_{N-1})} dx \end{align} $$

Given that $q_n(x)$ is an orthogonal polynomial, the Christoffel-Darboux theorem provides the following relation

$$ \begin{align} \frac{q_{N-2}(x)}{x-x_k}=\frac{k_{N-2}h_{N-3}}{k_{N-3}} \frac{1}{q_{N-3}(x_k)}\sum_{j=0}^{N-3} \frac{q_j(x)q_j(x_k)}{h_j} \label{w2x} \end{align} $$

Using this, we can calculate the weights for each case as follows:

• For $k=1,2,\cdots,N-2$

  $$ \begin{align} \omega_k&=\int_a^b \frac{q_{N-2}(x)}{(x-x_k)q’_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})(x-x_{N})}{(x_k-x_{N-1})(x_k-x_{N})} dx \nonumber\\ &=\frac{1}{(1-x_k^2)q’_{N-2}(x_k)}\int_{-1}^1 \frac{q_{N-2}(x)}{(x-x_k)}\cdot(1-x^2) dx \nonumber\\ &=\frac{k_{N-2}h_{N-3}}{k_{N-3}}\frac{1}{(1-x_k^2)q’_{N-2}(x_k)q_{N-3}(x_k)} \sum_{j=0}^{N-3} \frac{q_j(x_k)}{h_j}\int_{-1}^1 (1-x^2)q_j(x)dx \nonumber\\ &=\frac{k_{N-2}h_{N-3}}{k_{N-3}}\frac{1}{(1-x^2_k)q’_{N-2}(x_k)q_{N-3}(x_k)} \\ &=\frac{2(2N-1)(2N-3)}{N(N-1)(N-2)}\frac{1}{(1-x^2_i)q’_{N-2}(x_i)q_{N-3}(x_i)}\nonumber\\ &= \frac{2}{N(N-1)P^2_{N-1}(x_i)} \end{align} $$

• For $k=N-1$ (where $x_{N-1}=-1$)

  $$ \begin{align} \omega_k&=\int_a^b \frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N})}{(x_k-x_{N})} dx \nonumber\\ &=\frac{1}{-2q_{N-2}(-1)}\int_{-1}^1 q_{N-2}(x)\cdot (x-1) dx \nonumber\\ &=\frac{1}{2 \cdot(-c_{N-2}P’_{N-1}(-1))}(-c_{N-2})\int_{-1}^1 (1-x)P’_{N-1}(x) dx \nonumber\\ &=\frac{1}{2 \cdot (-1)^N \frac{N(N-1)}{2}}\left(\Bigl[(1-x)P_{N-1}(x)\Bigr]^1_{-1} +\int_{-1}^1 P_{N-1}(x) dx \right) \nonumber\\ &=\frac{2}{N(N-1)}\frac{1}{2 \cdot (-1)^N }\bigl(-2P_{N-1}(-1) \bigr) \\ &=\frac{2}{N(N-1)} \end{align} $$

• For $k=N$ (where $x_{N}=1$) $$ \begin{align} \omega_k=\frac{2}{N(N-1)} \end{align} $$

For $k=N$, although direct computation is possible, it’s not necessary due to symmetry.

During these calculations, $c_n$ is set such that $q_n(x) = -c_n P’_{n+1}(x)$ to satisfy the relation.

3.5. Summary #

The Gauss-Lobatto quadrature method allows for the exact evaluation of the integral of a polynomial $f(x)$ of degree $2N-3$ using $N (\geq 2)$ discrete function values.

$$ \begin{align} \int_{-1}^{1}f(x)dx=\sum_{i=1}^N \omega_i f(x_i) \end{align} $$

• Quadrature Points

For the polynomial of degree $N-2$ defined in the interval $x=[-1,1]$

$$ \begin{align} q_{N-2}(x)=\frac{P_{N}(x)+ u P_{N-1}(x)+v P_{N-2}(x)}{(1+x)(1-x)}=0,~~q_{N-2}(\pm 1)<\infty \end{align} $$

the total of $N$ points include $N-2$ points where $q_{N-2}(x)=0$ and the two points at $x=\pm 1$. When numbered in ascending order of $x$,

$$ \begin{gather} x_1=-1, x_N=1, \\ q_{N-2}(x)=0\to x=x_{i},~~(i=2,3,\cdots, N-1) \end{gather} $$

• Weights

$$ \begin{align} \omega_1 &= \omega_N = \frac{2}{N(N-1)}, \\ \omega_i &= \frac{2}{N(N-1)P^2_{N-1}(x_i)},~~(i=2,3,\cdots, N-1) \end{align} $$

For $N=2$, with $x_1=-1, x_2=1$, the weights are $\omega_1=1, \omega_2=1$, making it equivalent to the trapezoidal rule (exact for linear functions)

4. Appendix #

4.1. Legendre Polynomials Relations #

The relationship between Legendre polynomials is expressed as follows

Legendre polynomials -wikipedia

$$ \begin{align}\label{ap2} (2n+1)P_n(x)=\frac{d}{dx}\left(P_{n+1}(x)-P_{n-1}(x)\right), \end{align} $$

$$ \begin{align}\label{ap3} nP’_{n+1}(x)-(2n+1)xP’_{n}+(n+1)P’_{n-1}(x)=0. \end{align} $$

The differentiation of Legendre polynomials is given by⁶

$$ \begin{align}\label{ap5} \frac{d^m}{dx^m}P_n(x) = 1 \cdot 3 \cdots (2m-1)C^{(m+\frac{1}{2})}_{n-m}(x),~~(m\leq n) \end{align} $$

where $C^{(\alpha)}_{n}(x)$ is the Gegenbauer polynomial (or called Ultraspherical polynomial).

The special values and symmetry at $x=1$ are⁷

$$ \begin{align}\label{ap6} C^{(\alpha)}_{n}(1)&=\binom{n+2\alpha-1}{n}\\ C^{(\alpha)}_{n}(-x)&=(-1)^n C^{(\alpha)}_{n}(x) \end{align} $$

Consequently,

$$ \begin{align}\label{ap7} \frac{d^m}{dx^m}P_n(x)\Bigr|_{x=1} = 1 \cdot 3 \cdots (2m-1)\cdot \binom{n+m}{n-m},~~(m\leq n). \end{align} $$

5. References #