<jupyter_start><jupyter_text>Chebyshev polynomials on the unit circle[AMath 585, Winter Quarter 2020](http://staff.washington.edu/rjl/classes/am585w2020/) at the University of Washington. Developed by <NAME> and distributed under the [BSD license](https://github.com/rjleveque/amath585w2020/blob/master/LICENSE).  You are free to modify and use as you please, with attribution.These notebooks are all [available on Github](https://github.com/rjleveque/amath585w2020/).-------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.<jupyter_code>%matplotlib inline
from pylab import *
from mpl_toolkits.mplot3d import Axes3D
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()<jupyter_output><empty_output><jupyter_text>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 polynomialPlotting $z = \cos(k \theta)$ as a function of $x = \cos(\theta)$ gives the Chebyshev polynomial $T_k(x)$:<jupyter_code>figure(figsize=(8,4))
plot(x,z)
plot([-1,1],[0,0],'k');
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')<jupyter_output><empty_output><jupyter_text>AliasingBecause 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](https://my.siam.org/Store/Product/viewproduct/?ProductId=31254869) and illustrated below.<jupyter_code>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))<jupyter_output><empty_output><jupyter_text>On the unit circleNote 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.<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>Fourier AliasingAliasing 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](https://my.siam.org/Store/Product/viewproduct/?ProductId=1186)).<jupyter_code>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))
    
alias_plot(19)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15565847764800644, 'fraction_numerical': 0.03161498187676198, 'mean_word_length': 4.768873403019745, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 24, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '304889', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2086, 'n_tokens_neox': 1950, 'n_words': 547, 'path': 'notebooks/ChebPolyOnCircle.ipynb', 'repo_name': 'rjleveque/amath585w2020', 'size': '4972'}