Chebyshev polynomials on the unit circle

AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.

These notebooks are all available on Github.

This notebook illustrates the relation between Fourier modes on the unit circle and Chebyshev polynomials.

If running interactively, change the first cell to

%matplotlib notebook

so that you can rotate the 3D plots.

In [1]:
%matplotlib inline
In [2]:
from pylab import *
from mpl_toolkits.mplot3d import Axes3D
In [3]:
k = 15  # wave number.  Try 1, 2

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
theta = linspace(0, 2*pi, 1000)
x = cos(theta)
y = sin(theta)
z = cos(k*theta)
ax.plot(x,y,z,'b')
ax.plot(x,y,0*z,'k')  # plot unit circle

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

draw()

If you run this notebook you should be able to rotate the figure above. If you view it in the $x$-$z$ plane, you get the figure below.

Chebyshev polynomial

Plotting $z = \cos(k \theta)$ as a function of $x = \cos(\theta)$ gives the Chebyshev polynomial $T_k(x)$:

In [4]:
figure(figsize=(8,4))
plot(x,z)
plot([-1,1],[0,0],'k');
In [5]:
xfine = linspace(-1.1, 1.1, 1000)
Tnm = ones(xfine.shape)
Tn = xfine
for n in range(2,k+1):
    Tnp = 2*xfine*Tn - Tnm
    Tnm = Tn.copy()
    Tn = Tnp.copy()
    
figure(figsize=(8,4))
plot(xfine,Tn,'r', label='T_k(x)')
plot(x,z,'b--', label='cos(k arcos(x))')
grid(True)
legend(loc='lower right')
title('Chebyshev polynomial grows rapidly outside [-1,1]')

figure(figsize=(8,4))
plot(xfine,Tn,'r', label='T_k(x)')
plot(x,z,'b--', label='cos(k arcos(x))')
ylim(-2,2)
grid(True)
legend(loc='lower right')
title('Zooming in')
Out[5]:
Text(0.5, 1.0, 'Zooming in')

Aliasing

Because of the relation between Chebyshev polynomials and Fourier modes on the unit circle, aliasing takes place as described in Theorem 4.1 of Trefethen's Approximation Theory and Approximation Practice and illustrated below.

In [6]:
n = 7   # Use Chebyshev grid with n+1 points
         
k = 17  # Start with T_k
m = abs(mod(k+n-1,2*n) - (n-1))  # aliased to T_m

# Fine grid for plotting polynomials:
theta = linspace(0, pi, 500)
x = cos(theta)
Tk = cos(k*theta)
Tm = cos(m*theta)
         
# Grid points:
theta_j = linspace(0, pi, n+1)
xj = cos(theta_j)
Tk_j = cos(k*theta_j)

plt.figure()
plt.plot(x,Tk,'b')
plt.plot(xj, Tk_j, 'bo')
plt.plot(x,Tm,'r')

plt.title("The Chebyshev polynomial T_{%s} gets aliased to T_{%s}" % (k,m))
Out[6]:
Text(0.5, 1.0, 'The Chebyshev polynomial T_{17} gets aliased to T_{3}')

On the unit circle

Note that each Chebyshev point in the interval (with the exception of the endpoints) maps to two points on the unit circle. Below only one of these points is plotted.

In [7]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
theta = linspace(0, 2*pi, 1000)
x = cos(theta)
y = sin(theta)
Tk = cos(k*theta)
Tm = cos(m*theta)
ax.plot(x,y,Tm,'b')
ax.plot(x,y,Tk,'r')
ax.plot(x,y,0*Tm,'k')  # plot unit circle

theta_j = linspace(0, pi, n+1)
xj = cos(theta_j)
yj = sin(theta_j)
Tk_j = cos(k*theta_j)
ax.plot(xj, yj, Tk_j, 'bo')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.draw()

Fourier Aliasing

Aliasing occurs when a high frequency Fourier component is evaluated on a grid that is too coarse to represent it, and agrees with a lower frequency Fourier mode on the grid.

If we have a grid with $N$ nodes $x_j = 2\pi j/N$ (and suppose $N$ is even), then the highest wavenumber that can be represented by a Fourier mode of the form $e^{ikx}$ on this grid is $k = \pm N/2$. This is the "sawtooth" mode since $e^{i(\pm N/2)x_j} = e^{\pm i\pi j} = \left(e^{\pm i\pi}\right)^j = (-1)^j$.

If $k_1$ and $k_2$ differ by an integer multiple of $N$, then the two are indistinguishable on a grid of $N$ points, since $k_2 = k_1 + mN$ implies that $e^{ik_2x_j} = e^{ik_1x_j + imNjh} = e^{ik_1x_j} \left(e^{2\pi i}\right)^m = e^{ik_1x_j}$.

The code below illustrates this for $\sin(k_2 x)$, the imaginary part of the complex exponential (see also Figure 2.1 in Spectral Methods in Matlab).

In [8]:
N = 10
h = 2*pi/N
x = linspace(h, 2*pi, N)
x_fine = linspace(0, 2*pi, 1000)

def alias_plot(k1):
    # assumes k1 > 0
    k2 = mod(k1, N)
    if k2>N/2: 
        k2 = k2-N
    figure(figsize=(8,4))
    plot(x_fine, sin(k1*x_fine),'b')
    plot(x, sin(k1*x), 'bo')
    plot(x_fine, sin(k2*x_fine), 'r')
    xlim(0, 2*pi)
    title('sin(%s x) (blue) and sin(%s x) (red)' % (k1,k2))
In [9]:
alias_plot(19)