4. Chebyshev polynomials¶
As stated, Fourier series are only a good choice for periodic function. For problems with non-periodic boundary conditions, ansatz functions based on orthogonal polynomials are preferred. One popular choice are the Chebyshev polynomials (Pafnuty Lvovich Chebyshev, 1821–1894), defined on a domain x≤|1| as
The first polynomials are thus (see also Fig. 4.1)
There exists also a recursion formula

Fig. 4.1 Chebyshev polynomials Tk(x) for k=0,…,6.¶
A function u(x) is approximated via a finite series of Chebyshev polynomials as
with the ak being the Chebyshev coefficients. Note that the sum goes from 0 to N, i.e. there are N+1 coefficients. The highest order of the polynomial approximation is thus N.
A general problem of high-order interpolation of a function via high-order polynomials is the Runge phenomenon (Carl David Tolmé Runge 1856–1927), similar to the Gibbs phenomenon for Fourier series. If a function is interpolated on an equidistant grid, the error grows as 2N. However, using a non-equidistant distribution of points such that the point density is (approximately) proportional to N√1−x2−1, i.e. denser towards the domain boundaries, it can be shown that the interpolation errors decrease exponentially.
A common distribution of points in particular for Chebyshev polynomials are the Gauss–Lobatto–Chebyshev (GLC) points
The N+1 points xj correspond to the locations of the extrema of TN=±1. A more intuitive way of visualising the GLC nodes is by diving a half circle into N identical sectors and projecting their edge coordinates.

Fig. 4.2 Illustration of the Runge phenomenon: Polynomial interpolation of the function f(x)=(1+16x2)−1 on equidistant and non-equidistant (Gauss–Lobatto–Chebyshev) grids.¶
Chebyshev polynomials have a number of important properties:
Alternating even and odd functions:
Orthogonality (note the weight in the inner product):
with c0=2, ck=1 for k≥1 and the Kronecker symbol δkl.
Boundary conditions:
To evaluate the finite Chebyshev series (4.1) use the stable recursive algorithm
The values of the Chebyshev polynomials on the Gauss-Lobatto nodes are
The transformation between the physical space uN and spectral (Chebyshev) space ak is done via the so-called Chebyshev transform. Since the Chebyshev polynomials are essentially cosine functions on a transformed coordinate, there exists a fast transform based on the FFT. As usual, the linear transform can also be represented by a matrix–vector multiplication with a full matrix.
If a collocation method on the Gauss-Lobatto grid (4.2) is employed, the derivative of a discretised function uN can be written as a matrix multiplication,
The matrix D__=[Dij] is known as the (Chebyshev) derivative matrix, and represents a transformation into spectral space, spectral derivation and back-transform. An explicit form of the (N+1)×(N+1)-matrix D__=[Dij] is given by
with
Note that the computation of D__ is very sensitive to round-off errors. Therefore, it is recommended for practical computations to use specificially adapted versions of the above formula, e.g. chebdif.m
from “A Matlab Differentiation Matrix Suite” by J.A.C. Weideman and S.C. Reddy (or the corresponding Python versions as in the attached notebooks.