Introduction to pseudospectral methods
Waves on a grid
- Periodic function \(u(x)\) with period \(L\)
- Domain with \(N\) points indexed \(j=0\) to \(j=N-1\)
- Sine and cosine waves on the grid.
- How many will fit?
\[\cos(2\pi nx/L)\]
- Smallest wave that fits has a period twice the grid spacing: \(\cos{\pi Nx/L}\)
- \(n_\text{max} = N/2\)
Sine and cosine combinations
\[u(x) = a_0 + \sum_{n=1}^{N/2}[a_n\cos(2\pi nx/L)) + b_n\sin(2\pi nx/L)]\]
Complex form using Euler identities \[e^{i\theta} = \cos(\theta) + i\sin(\theta)\] \[\cos(\theta) = \frac{1}{2}(e^{i\theta} + e^{-i\theta})\] \[\sin(\theta) = \frac{1}{2i}(e^{i\theta} - e^{-i\theta})\] Insert into the above equation \(\Rightarrow\)
\[u(x) = \sum_{n=-N/2}^{N/2}c_ne^{2\pi inx/L}\]
\[ c_n = \begin{cases} a_0 & \text{for } n=0 \\ a_n/2 + b_ni/2 & \text{for } n<0, \\ a_n/2 - b_ni/2 & \text{for } n>0 \end{cases} \]
- \(e^{2\pi inx/L}\) is point on the unit circle in the complex plane
- Write this as \(e^{\theta_xn i}\), where \(\theta_x=2\pi x/L\), and is the angle measured from the real line
- Corresponding \(\pm n\) will give points on the unit circle with conjugate symmetry
- For example, for \(x=\Delta x\), we have \(e^{2\pi in/N}\), where \(\Delta x/L=1/N\)
- Points on the complex unit circle for \(n=-N/2\ldots N/2\) are shown for \(N=6\)
- This figure represents the basis functions \(e^{2\pi inx/L}\) for \(x=\Delta x\)
- In general, \(x_j = j\Delta x = jL/N\), giving \(e^{2\pi inj/N}\) for grid points \(j\).
- For other \(j\) values (grid points), the angle stepped from one value of \(n\) to the next is proportional to \(j\).
- This is illustrated in the next for each \(x_j=j\Delta x\)
- The red curve shows the angular size of the step between each \(n\) value
- The locations of the \(n\) indexes are labeled
Odd and even \(N\)
- We have considered even \(N\).
- For odd \(N\):
- smallest period is again \(2\Delta x=2L/N\)
- but the max number \(n\) of these waves that fit fully on the grid is \(\lfloor{N/2}\rfloor\)
- For \(N=5\) we have \(n=\) -2, -1, 0, 1, 2
IDFT and DFT
IDFT: at \(x_j\), our equation is the inverse discrete Fourier Transform,
\[\color{blue}{u_j = \sum_{n=-N/2+1}^{N/2}c_ne^{2\pi inj/N}}\]
DFT: the discrete Fourier transform is
\[\color{blue}{c_n = \frac{1}{N}\sum_{j=0}^{N-1}u_je^{-2\pi inj/N}}\]
Evaluate using fast transforms (fft, ifft), denoted \(\mathbf{c} = \mathcal{F}(\mathbf{u})\), and \(\mathbf{u} = \mathcal{F}^{-1}(\mathbf{c})\).
\[\sum_{n=-N/2+1}^{N/2}c_ne^{2\pi inj/N} = \sum_{m=0}^{N-1}\hat{c}_m e^{2\pi imj/N}\]
Points shown for sums over \(n\) and \(m\) on the left and right for \(j=1\)
- The same points appear on the circle for both index sets.
- For the sums to be equal, the same complex coefficients need to multiply the same points on the unit circle.
- Hence, \(c_{-2}=\hat{c}_4\), \(c_{-1}=\hat{c}_5\), \(c_{0}=\hat{c}_0\), \(c_{1}=\hat{c}_1\), \(c_{2}=\hat{c}_2\), \(c_{3}=\hat{c}_3\)
- For other values of \(j\), the points on the unit circle are the same as those shown in Figure 1
Compare spectral and finite difference
- For an error of 4\(\times 10^6\), the spectral method requires 13 times fewer points (40 versus 512) than the FD central method
- Spectral methods give exponential convergence; FD methods give power-law convergence.
Chebyshev polynomials
The first seven Chebyshev polynomials of the first kind are
\[\begin{align} &T_0(x) = 1 \\ &T_1(x) = x \\ &T_2(x) = 2x^2-1 \\ &T_3(x) = 4x^3-3x \\ &T_4(x) = 8x^4-8x^2+1 \\ &T_5(x) = 16x^5-20x^3+5x \\ &T_6(x) = 32x^6-48x^4+18x^2-1 \\ \end{align}\]
Trigonometric relationship
- \(T_n(x) = T_n(\cos\theta) = \cos(n\theta)\)
- as \(\theta\in[0,\pi]\) around the half circle
- \(x=\cos(0)\) is projected to the axis
- The value of \(T_n\) on the unit circle, and projected to the axis is \(\cos(n\theta)\)
- We don’t see the variation because it is out-of-plane
- tip the plane…
Chebyshev roots and extrema
Roots \(N\) roots of \(T_N(x)\) are at \[x_j = \cos\theta_j,\] \[\theta_j = \frac{\pi}{N}(j+1/2)\] with \(j=0,\,\ldots,\,N-1\)
- uniform \(\theta_j\in\left[\frac{\pi}{2N}, \pi-\frac{\pi}{2N}\right]\)
- roots projected to \(x_j=\cos\theta_j\)
- roots ordered from high to low
Extrema \(N+1\) extrema of \(T_N(x)\), (min/max values) at \[x_j = \cos(\theta_j),\; \theta_j = \frac{\pi j}{N},\; j=0,\,\ldots,\,N\]
- includes \(x=1\) and \(x=-1\).
- When using pseudospectral methods, we choose a grid
- The roots and extrema grids are the most common
- The extrema grid is convenient for treating boundary conditions since the boundary points are included
def make_plot(N, aspect=(8,8), do_unif=True):
#--------- roots grid, option 1, direct calculation (same as option 2 below)
xcr = cheb.chebpts1(N)[::-1] # roots grid
ucr = urunge(xcr)
T = np.zeros((N,N))
for n in range(N):
for j in range(N):
T[n,j] = np.cos(np.pi*n/N*(j+0.5))
cr = np.linalg.inv(T.T) @ ucr
#--------- roots grid, option 2, us numpy.polynomial.chebyshev module
cr = cheb.chebinterpolate(urunge, N-1) # TT = T^T = cheb.chebvander(xcr, N-1); not needed...
#--------- extrema grid
xce = cheb.chebpts2(N)[::-1] # extrema grid
uce = urunge(xce)
for n in range(N):
for j in range(N):
T[n,j] = np.cos(np.pi*n*j/(N-1))
ce = np.linalg.inv(T.T) @ uce
#--------- uniform grid, polyfit
xu = np.linspace(-1,1,N)
p = np.polyfit(xu, urunge(xu), N-1)
#--------- plot results
plt.rcParams.update({'font.size': 16})
fig, ax = plt.subplots(1,1, figsize=aspect)
ax.plot(xx,urunge(xx), color='black', label='Exact', lw=2)
ax.plot(xx, chebyterp(cr,xx), '-', color='blue', label='Cheby roots')
ax.plot(xcr, ucr, 'o', ms=7, fillstyle='full', color='blue')
ax.plot(xx, chebyterp(ce,xx), '-', color='green', label='Cheby extrema')
ax.plot(xce, uce, '^', ms=7, fillstyle='full', color='green')
if do_unif:
ax.plot(xx, np.polyval(p,xx), '-', color='red', label='Uniform')
ax.plot(xu, urunge(xu), 's', ms=7, fillstyle='none', color='red')
ax.legend(frameon=False)
ax.set_xlabel('x')
ax.set_ylabel('u(x)');
#-------------
make_plot(9)
\(N=19\) without the uniform grid
Example, Runge function with DCT\(_\text{II}\)
Example, Runge function with DCT\(_\text{I}\)
Example, Runge derivative
- derivative of \(u_v(x) = 0.1e^{-2x}\sin(6x) - 8x\) shown for variety
