An introduction to pseudospectral methods
This document presents an introduction to pseudospectral methods for solution of partial differential equations. The focus is on
- periodic domains using Fourier series,
- nonperiodic domains using Chebyshev series,
- one and two spatial dimensions,
- boundary conditions including periodic, Dirichlet, Nuemann, and Robin.
Suppose we have a PDE for an unsteady, one-dimensional conservation law with some diffusive flux \(q=-\Gamma\partial u/\partial x\), and a source term \(S\): \[ \frac{\partial u}{\partial t} = -\frac{\partial}{\partial x}\left(\Gamma\frac{\partial u}{\partial x}\right) + S. \tag{1}\]
Here, \(u = u(t,x)\), \(S\) is a source term, and \(\Gamma\) is a diffusivity. If \(\Gamma\) is a function of \(u\), or if \(S\) is a nonlinear function of \(u\), then the problem is nonlinear.
We approximate the function \(u(t,x)\) as \[ u(t,x) = \sum_{n=0}^{N-1} c_n(t)\phi_n(x). \tag{2}\] That is, \(u\) is approximated as a linear combination of basis functions \(\phi_n(x)\), which depend only on \(x\), and \(c_n\) are coefficients that depend only on \(t\). (For simplicity, we use the same notation for the solution \(u(t,x)\) of the PDE, and its \(N\)-term series approximation.)
For periodic spatial domains with period \(L\), the basis functions \(\phi_n\) are taken as \[ \phi_n(x) = e^{2\pi inx/L}. \] We can refer to these \(\phi_n\) as Fourier modes.
For nonperiodic spatial domains with \(-1\le x\le 1\), we will use \[ \phi_n(x) = T_n(x), \] where \(T_n(x)\) is a Chebyshev polynomial of first kind.
Other basis functions can be used for \(\phi_n\), and the form of Equation 2 is used in other methods, such as Finite Element and Galerkin methods. The latter is a spectral method, as opposed to the pseudospectral method discussed here. Importantly, for spectral and pseudospectral methods, the basis functions \(\phi_n\) are global (they span the domain), rather than local (nonzero around a given location), as in Finite Element methods. This global character of spectral and pseudospectral methods contributes to the high accuracy and convergence rates of the methods. The Fourier modes and Chebyshev polynomials are discussed in more detail below.
In the pseudospectral method, the coefficients \(c_n\) are evaluated by requiring that \(u(t,x)\) satisfy the PDE at \(N\) discrete grid points called collocation points (sometimes also called called interpolation points). The pseudospectral method is also sometimes called the collocation method. We denote the collocation points as \(x_j\) for \(j=0,\,\ldots N-1\), and \(u_j\equiv u(x_j)\). Similarly, \(\Gamma_j\) and \(S_j\) are \(\Gamma\) and \(S\) evaluated at \(x_j\). Then \[ u_j = \sum_{n=0}^{N-1} c_n(t)\phi_n(x_j). \tag{3}\]
In the following sections, Equation 3 will be related to the inverse discrete Fourier transform (IDFT) or the inverse discrete cosine transform (IDCT). These then allow evaluation of spatial derivatives at the collocation (grid) points, so that a PDE, such as that in Equation 1 can be solved using, e.g., the method of lines (MOL), where the PDE is converted to a system of ODEs at each point, \[ \frac{\partial u_j}{\partial t} = -\frac{\partial}{\partial x}\left(\Gamma_j\frac{\partial u_j}{\partial x}\right)_j + S_j. \tag{4}\] Here, \(\partial u_j/\partial x\) denotes the partial derivative of \(u(t,x)\) evaluated at \(x_j\); and similarly for the \(()_j\) term. Equation 4 is then solved as a system of ODEs using some ODE solver, such as the fourth-order Runge-Kutta method. Importantly, the pseudospectral method solves the PDE in physical space, with derivatives evaluated as noted, and nonlinearity poses no special difficultly. This will be demonstrated in examples below.
Fourier pseudospectral methods
Here, we relate Equation 3 with \(\phi_n(x)=e^{2\pi inx/L}\) to the inverse discrete Fourier transform (IDFT), and show examples of computing derivatives and solving PDEs.
Discrete Fourier transform
Consider a periodic function \(u(x)\) with period \(L\), and a domain with \(N\) points indexed \(j=0\) to \(j=N-1\).
We can put sine and cosine waves on this grid. How many waves will fit? That is, what is the highest wavenumber that will fit on the grid? Consider \(\cos(2\pi nx/L)\). This function has period \(L/n\). When \(n=1\), we have a single complete wave on the domain. When \(n\) is larger, we have more waves on the grid with smaller wavelengths. The smallest wave that will fit on the grid has a period of twice the grid spacing, or \(2\Delta x=2L/N\). But this period is \(L/n_\text{max}\), so \(n_\text{max}=N/2\). This gives \(\cos(2\pi n_\text{max}x/L) = \cos(\pi N x/L)\), shown below.
For a sine wave, we have \(\sin(2\pi n_\text{max}x/L) = \sin(\pi N x/L)\).
Plotting all the cosine and sine waves gives
A function \(u(x)\) can be written as a linear combination of the sine and cosine basis functions (one for each \(n\)),
\[u(x) = a_0 + \sum_{n=1}^{N/2}[a_n\cos(2\pi nx/L)) + b_n\sin(2\pi nx/L)]. \tag{5}\]
This equation is converted to complex form using the 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 these into Equation 5 to give
\[u(x) = a_0 + \sum_{n=1}^{N/2}c_ne^{2\pi inx/L} + \sum_{n=1}^{N/2}k_ne^{-2\pi inx/L}, \tag{6}\]
where \[c_n = \frac{a_n}{2} + \frac{b_n}{2i} = \frac{a_n}{2} - \frac{b_ni}{2},\] \[k_n = \frac{a_n}{2} - \frac{b_n}{2i} = \frac{a_n}{2} + \frac{b_ni}{2}.\]
In the second summation in Equation 6, if we replace \(n\) with \(-n\), then the first and second sums are over different ranges of \(n\). This allows us to use \(c_n\) in both sums, that is, use \(c_n\) instead of \(k_n\), notationally, in the second sum. Recognizing that, for \(n=0\), \(a_0=c_ne^{2\pi inx/L}\), we can combine all three terms into one sum from \(n=-N/2\) to \(n=N/2\):
\[u(x) = \sum_{n=-N/2}^{N/2}c_ne^{2\pi inx/L}. \tag{7}\]
Here, \[ 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} \]
Note that \(e^{2\pi inx/L}\) is point on the unit circle in the complex plane. We can write this as \(e^{\theta_xn i}\), where \(\theta_x=2\pi x/L\), and is the angle measured from the real line. Corresponding positive and negative values of \(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 below for \(N=6\).
Note, this figure represents the basis functions \(e^{2\pi inx/L}\) for \(x=\Delta x\). For other \(x\) values, the angle stepped from one value of \(n\) to the next is proportional to the size of \(x\). This is illustrated below 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.
For real \(u(x)\), the \(c_n\) will have conjugate symmetry, \(c_n=c_{-n}^*\), so that the imaginary parts of the sum cancel. For example, consider two terms of the summation for some given \(n\) and \(-n.\) If \(c_n=c_{n,r} + ic_{n,i}\) and \(c_{-n} = c_{-n,r}+ic_{-n,i}\), then
\[\begin{align} c_ne^{2\pi inx/L} + c_{-n}e^{-2\pi inx/L} = & (c_{n,r}+ic_{n,i})(\cos(2\pi nx/L) + i\sin(2\pi nx/L)) + \\ & (c_{-n,r}+ic_{-n,i})(\cos(-2\pi nx/L) + i\sin(-2\pi nx/L)). \end{align}\]
The imaginary part of this expression is
\[\begin{align} &ic_{n,r}\sin(2\pi nx/L) + ic_{n,i}\cos(2\pi nx/L) + \\ &ic_{-n,r}\sin(-2\pi nx/L) + ic_{-n,i}\cos(-2\pi nx/L) \\ \end{align}\]
Since \(\sin(-\phi)=-\sin(\phi)\) and \(\cos(-\phi) = \cos(\phi)\), the imaginary part cancels when \(c_{n,r}=c_{-n,r}\) and \(c_{n,i}=-c_{-n,i}\), that is, when \(c_n=c_{-n}^*\). Note that points \(n=0\), \(n=-N/2\), and \(n=N/2\) have no imaginary component.
Note that while the \(c_n\) have conjugate symmetry, and there is a coefficient \(c_n\) associated with each basis function \(e^{2\pi inx/L}\), and the \(c_n\) are complex numbers, those points are not at the same locations as the points in the images above, nor do they have the same complex angle. That is, don’t confuse the points on the circle associated with the Fourier modes with the complex coefficients multiplying those Fourier modes.
The sum in Equation 7 is from \(n=-N/2\) to \(n=N/2\). At these two bounding values, \(e^{2\pi inx/L}\) becomes \(e^{\pm\pi iNx/L}.\) At \(x=x_j=j\Delta x\), with \(\Delta x/L = 1/N\), we have \(e^{\pm\pi ij} = \pm 1\) for all integer \(j.\) (In Figure 1, \(n=3\) and \(n=-3\) are at the same location for all \(j.\)) Hence, including both \(n=-N/2\) and \(n=N/2\) in the sum is redundant, and we remove \(n=-N/2,\) giving
\[u(x) = \sum_{n=-N/2+1}^{N/2}c_ne^{2\pi inx/L}. \tag{8}\]
Above, we have considered even \(N\). For odd \(N\), the smallest period is again \(2\Delta x=2L/N\), but the maximum number \(n\) of these waves that will fit fully on the grid is \(\lfloor{N/2}\rfloor\). For \(N=5\), we have \(n=\) -2, -1, 0, 1, 2.
Equation 8 is the inverse discrete Fourier transform (IDFT), rewritten here and evaluated at \(x=x_j,\)
\[\color{blue}u_j = \sum_{n=-N/2+1}^{N/2}c_ne^{2\pi inj/N}. \tag{9}\]
The corresponding discrete Fourier transform (DFT) is
\[\color{blue}c_n = \frac{1}{N}\sum_{j=0}^{N-1}u_je^{-2\pi inj/N}. \tag{10}\]
These are evaluated using fast Fourier transform (FFT) and its inverse (IFFT), denoted \(\mathcal{F}\) and \(\mathcal{F}^{-1},\) respectively. That is, \(\mathbf{c}=\mathcal{F}(\mathbf{u})\), and \(\mathbf{u}=\mathcal{F}^{-1}(\mathbf{c})\). Here, the bold symbols denote arrays. So, \(\mathcal{F}(\mathbf{u})\) is the discrete Fourier transform of the \(u\) array at each grid point \(j\), and \(\mathbf{c}\) is the array of coefficients with elements \(c_n\).
We can show that \(u_j = \mathcal{F}^{-1}(\mathcal{F}(u_j))\) by inserting Equation 10 into Equation 9. Before doing this, we change index \(j\) in Equation 10 to \(J\) so that it doesn’t conflict with the \(j\) in Equation 9:
\[u_j = \sum_{n=-N/2+1}^{N/2}\left[ \frac{1}{N}\sum_{J=0}^{N-1}u_Je^{-2\pi inJ/N} \right]e^{2\pi inj/N}.\]
Change the order of the sums and arrange as
\[u_j = \sum_{J=0}^{N-1}u_J\frac{1}{N}\left[ \sum_{n=-N/2+1}^{N/2}e^{2\pi in(j-J)/N}\right].\]
When \(J=j\), \(e^{2\pi in(j-k)/N}=1\) and the sum in brackets is \(N\). When \(J\ne j\), the sum is 0. Hence, \(u_j=u_j\) is recovered.
Variations
The DFT/IDFT forms given here are not unique, and several variations are common.
- The normalizing factor \(1/N\) sometimes appears on the IDFT instead of the DFT. Sometimes it is split uniformly so that a \(1/\sqrt{N}\) factor appears in both the IDFT and the DFT.
- The sign in the exponent is sometimes switched between the IDFT and the DFT.
- In the IDFT, the sum from \(n=-N/2+1\) to \(N/2\) is sometimes given from \(n=-N/2\) to \(N/2-1\).
- Finally, in the IDFT, the range \(n=-N/2+1\) to \(N/2\) is sometimes cycled to \(n=0\) to \(N-1\).
This last point can be shown as follows. Show that
\[\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}.\]
Consider \(N=6\). The sums over \(n\) and \(m\) are, respectively,
\[ c_{-2}e^{\theta(-2)} + c_{-1}e^{\theta(-1)} + c_{0}e^{\theta(0)} + c_{1}e^{\theta(1)} + c_{2}e^{\theta(2)} + c_{3}e^{\theta(3)},\]
\[ \hat{c}_{0}e^{\theta(0)} + \hat{c}_{1}e^{\theta(1)} + \hat{c}_{2}e^{\theta(2)} + \hat{c}_{3}e^{\theta(3)} + \hat{c}_{4}e^{\theta(4)} + \hat{c}_{5}e^{\theta(5)},\]
where \(\theta = 2\pi ij/N\). As before, the exponential terms are points in the complex plane on the unit circle. Those points are shown for first and second sums (over \(n\), and \(m\)), on the left and right below for \(j=1\).
Note that 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, we have \(c_{-2}=\hat{c}_4\), \(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
Derivatives
In applying the pseudospectral method to partial differential equations, we need to evaluate spatial derivatives as points \(x_j\), such as \(\partial u_j/\partial x\). The method is subject to truncation error associated with the finite number of grid points and spectral modes used, but derivatives are evaluated exactly in spectral space. We have
\[u(x) = \sum_{n=-N/2+1}^{N/2}c_ne^{2\pi inx/L},\] \[\frac{\partial u}{\partial x} = \sum_{n=-N/2+1}^{N/2}(2\pi in/L)c_ne^{2\pi inx/L},\] \[\frac{\partial^2 u}{\partial x^2} = \sum_{n=-N/2+1}^{N/2}-(2\pi n/L)^2c_ne^{2\pi inx/L}.\]
Evaluated at \(x=x_j\), the derivatives can be written as
\[\frac{\partial \mathbf{u}}{\partial x} = \mathcal{F}^{-1}[(2\pi i\mathbf{n}/L)\odot\mathcal{F}(\mathbf{u})],\] \[\frac{\partial^2 \mathbf{u}}{\partial x^2} = \mathcal{F}^{-1}[(2\pi i\mathbf{n}/L)^2\odot\mathcal{F}(\mathbf{u})].\]
Here, \(\mathbf{n}\) is a vector with elements \(-N/2+1,\,\ldots,\,N/2\). The \(\odot\) symbol indicates element-wise multiplication of two arrays. \(\partial\mathbf{u}/\partial x\) is the array of first derivatives at each grid point \(x_j\).
In a PDE with the divergence of a diffusive flux \(q=-\Gamma(u)du/dx\), we have
\[\frac{\partial \mathbf{q}}{\partial x} = \mathcal{F}^{-1}[(2\pi i\mathbf{n}/L)\odot\mathcal{F}(\mathbf{q})],\] \[\mathbf{q} = -\Gamma(\mathbf{u})\odot\mathcal{F}^{-1}[(2\pi i\mathbf{n}/L)\odot\mathcal{F}(\mathbf{u})].\]
Example, viscous Burgers equation
The viscous Burgers equation is
\[\frac{\partial u}{\partial t} = -u\frac{\partial u}{\partial x} + \nu\frac{\partial^2u}{\partial x^2}.\]
Here, \(\nu\) is viscosity. This equation is nonlinear.
- Use a periodic initial condition of \(u_0(x) = \cos(2\pi x/L)+2\).
- \(\nu=0.1\)
- \(L=10\)
- Solve to \(t=10\).
Exact solution
The exact solution is given by
\[u(x,t) = -2\nu\frac{\partial}{\partial x}\ln\left\{\frac{1}{\sqrt{4\pi\nu t}}\int_{-\infty}^\infty\exp\left[-\frac{(x-x^\prime)^2}{4\nu t}-\frac{1}{2\nu}\int_0^{x^\prime}u_0(x^{\prime\prime})dx^{\prime\prime}\right]dx^\prime\right\}.\]
This equation is evaluated here using Sympy.
import sympy as sp
v, t, L, x, xp = sp.symbols('v, t, L, x, xp')
ex = L/(4*sp.pi*v)*sp.sin(2*sp.pi*xp/L) + xp/v
ex = -2*v*sp.diff(sp.ln(1/sp.sqrt(4*sp.pi*v*t)*
sp.Integral(sp.exp(-(x-xp)**2/4/v/t - ex),
(xp,-sp.oo, sp.oo))), x)
ex = ex.subs([(v,0.1),(t,10.0),(L,10.0)]) # make sure these are the same as below
uex = sp.lambdify(x, ex) # then u_exact(x) = uex(x)Spectral solution
\[\frac{d\mathbf{u}}{dt} = -\mathbf{u}\odot\mathcal{F}^{-1}\left(\frac{2\pi i\mathbf{n}}{L}\odot\mathcal{F}(\mathbf{u})\right) + \nu \mathcal{F}^{-1}\left(\left(\frac{2\pi i\mathbf{n}}{L}\right)^2\odot\mathcal{F}(\mathbf{u})\right).\]
Here, the square of \(\mathbf{n}\) is done element-wise.
from scipy.integrate import odeint
from scipy.fft import fft, ifft
def spectral(nx):
L = 10.0
v = 0.1 # use same value as the exact solution above
tend = 10 # use same value as the exact solution above
#---------- solution domain, initial condition
dx = L/nx # not L/(nx-1)
x = np.linspace(0.0, L-dx, nx)
u0 = np.cos(2*np.pi*x/L) + 2
#---------- solve the problem
def rates(u, t):
N = len(u)
n = np.arange(N); n[int(N/2)+1:]-= N # n[int(N/2):]-=N
return -u*ifft(2*np.pi*1j*n/L*fft(u)).real - \
v*ifft((2*np.pi*n/L)**2*fft(u)).real
t = np.linspace(0,tend,11)
u = odeint(rates, u0, t)
return x, u, u0
#---------- (Now, call spectral(32) to solve, then plot result, not shown)
Finite difference solver
Discretize the spatial derivatives using second order central differences on the diffusion term, and either upwind or central differences on the advective term.
def FD(nx, upwind=False):
L = 10.0
v = 0.1 # use same value as the exact solution above
tend = 10 # use same value as the exact solution above
#---------- solution domain, initial condition
dx = L/nx # not L/(nx-1)
x = np.linspace(0.0, L-dx, nx)
u0 = np.cos(2*np.pi*x/L) + 2
#---------- solve the problem
i = np.arange(nx)
ip = i+1; ip[-1] = 0
im = i-1; im[0] = nx-1
def rates(u, t):
if upwind: # upwind on the advective term, 1st order
return -u*(u - u[im])/dx + v/dx/dx*(u[im] - 2*u + u[ip])
else: # central difference, 2nd order
return -u*(u[ip] - u[im])/2/dx + v/dx/dx*(u[im] - 2*u + u[ip])
t = np.linspace(0,tend,11)
u = odeint(rates, u0, t)
return x, u, u0Compare spectral and finite difference
Compare the methods on grids with varying number of points.
- The spectral method bottoms out at \(\le\) 80 points, where the average error is \(0.4\times 10^6\) times less than the second order finite difference method.
- For an average relative error of around 4\(\times 10^6\), the spectral method requires around 13 times fewer grid points (40 versus 512).
- Spectral methods give exponential convergence, while finite difference methods give power-law convergence.
2D transforms
In one dimension, the basis function of the Fourier series is \(\phi_n(x)=e^{2\pi inx/L}\). In two dimensions we have \(\phi_{n,m}(x,y) = e^{2\pi i(nx/L_x+my/L_y)} = e^{2\pi inx/L_x}e^{2\pi imy/L_y}\).
The IDFT and DFT are given by
\[\text{IDFT:}\phantom{xxx}u_{j,k} = \sum_{n=-N_n/2+1}^{N_n/2}\left[ \sum_{m=-N_m/2+1}^{N_m/2} c_{n,m}e^{2\pi imy_k/L_y}\right] e^{2\pi inx_j/L_x},\]
\[ \text{DFT:}\phantom{Ixxx}c_{n,m} = \frac{1}{N_x}\sum_{j=0}^{N_x-1}\underbrace{\left[ \frac{1}{N_y} \sum_{k=0}^{N_y-1} u_{j,k}e^{-2\pi imy_k/L_y} \right]}_{A_j} e^{-2\pi inx_j/L_x}. \tag{11}\]
Note, \(N_n=N_x\), and \(N_m=N_y\). The term in brackets and \(A_j\) are defined to facilitate the presentation and computation. We recognize the term in brackets in the DFT as a one-dimensional DFT in the y-direction, operating on the row \(u_{j,k}\) for given \(j\). Each \(A_j\) is an array with elements corresponding to \(m=-N_m/2+1,\ldots,N_m/2\). Hence, \(A\) is a matrix with elements \(A_{j,m}\) (rows \(j\) and columns \(m\)). The outer sum is also a one-dimensional DFT operating on columns of \(A_{j,m}\) for each \(m\).
In Python, \(c_{n,m}\) is computed as:
for j in range(Nx):
A[j,:] = fft(u[j,:])
for m in range(Nm):
c[:,m] = fft(a([:,m]))We don’t need an intermediate matrix \(A\), rather, \(c\) can be computed . The code for the transforms is given below and compared with the two-dimensional transform available directly in ths scipy.fft module. Note that the loops are the same between the DFT and the IDFT, except for notation, and trading dft and idft, as well as \(c\) and \(u\).
import numpy as np
from scipy.fft import fft, ifft, fft2, ifft2
#---------- 2D DFT
def fft2D(u):
c = np.empty_like(u, dtype=np.complex128)
for j in range(Nx):
c[j,:] = fft(u[j,:])
for m in range(Nm):
c[:,m] = fft(c[:,m])
return c
#---------- 2D IDFT
def ifft2D(c):
u = np.empty_like(c, dtype=np.complex128)
for n in range(Nn):
u[n,:] = ifft(c[n,:])
for k in range(Ny):
u[:,k] = ifft(u[:,k])
return u
#---------- (These give identical results to fft2 and ifft2)Derivatives
The partial derivative with respect to \(x\) evaluated at grid point \(j\), \(k\), is denoted \(\partial u_{j,k}/\partial x\), and written as
\[ \frac{\partial u_{j,k}}{\partial x} = \sum_{n=-N_n/2+1}^{N_n/2}\frac{2\pi in}{L_x}\left[ \sum_{m=-N_m/2+1}^{N_m/2} c_{n,m}e^{2\pi imy_k/L_y}\right] e^{2\pi inx_j/L_x}. \tag{12}\]
For higher-order derivatives, we take corresponding powers of \((2\pi in/L_x)\).
In principle, if we compute the 2D array \(c_{n,m}\) from a 2D DFT, then we can use it in Equation 12 to compute the 2D derivative array using loops over 1D IDFT.
Instead, Equation 12 can be simplified. We insert Equation 11 for \(c_{n,m}\) after changing the indices in that equation from \(j\) to \(J\) and \(k\) to \(K\) so they don’t conflict with \(j\) and \(k\) in Equation 12. This gives
\[ \frac{\partial u_{j,k}}{\partial x} = \frac{1}{N_xN_y}\sum_n\frac{2\pi i n}{L_x}e^{2\pi inx_j/L_x}\sum_J\sum_K u_{J,K}e^{-2\pi inx_J/L_x}\sum_me^{2\pi i\frac{m}{N_y}(k-K)}, \]
where we have reordered the summations and terms, and used \(y_k/L_y=k/N_y\) (and similarly for \(K\)). The sum on the far right is \(N_y\) when \(k=K\) and zero otherwise. This cancels with the leading \(1/N_y\) factor and the sum over \(K\) becomes a single term for \(k\). Rearranging gives
\[ \frac{\partial u_{j,k}}{\partial x} = \sum_{n=-N/2+1}^{N_x/2}\frac{2\pi i n}{L_x}\left[\frac{1}{N_x}\sum_{J=0}^{N_x-1}u_{J,k}e^{-2\pi inx_J/L_x}\right]e^{2\pi inx_j/L_x}. \]
The term in brackets is the DFT of a line of \(u\) in the \(x\)-direction for a given \(y\)-direction index \(k\). The outer sum is an IDFT of the term multiplying the far right exponential factor. The derivative is then
\[\frac{\partial \mathbf{u}_k}{\partial x} = \mathcal{F}^{-1}[(2\pi i\mathbf{n}/L_x)\odot\mathcal{F}(\mathbf{u}_k)].\]
Here, \(\mathbf{u}_k\) is the array of \(u\) values at points \(j\) for given \(k\). The 2D partial derivative with respect to \(x\) at a given \(y\) is evaluated the same as for 1D, and similarly for the derivative with respect to \(y\). Higher derivatives have higher powers of \(2\pi i\mathbf{n}/L_x\), where powers of \(\mathbf{n}\) are computed element-wise
Example, 2D heat equation
Solve the following PDE with a constant diffusivity \(\Gamma=1\):
\[\frac{\partial u}{\partial t} = \Gamma\frac{\partial^2u}{\partial x^2} +\Gamma\frac{\partial^2u}{\partial x^2}.\]
Use periodic boundary conditions and centered rectangular initial condition. The code below solves the problem using the pseudospectral method. A finite difference solution is also coded. Just uncomment the last line.
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, ifft, fft2
from IPython.display import display, clear_output
import time
#--------- setup
Lx = 2.0*np.pi
Ly = 4.0*np.pi
Nx = 20
Ny = 40
G = 1.0
trun = 1.0
ntimes = 31
dx = Lx/Nx
dy = Ly/Ny
x = np.linspace(0.0, Lx-dx, Nx)
y = np.linspace(0.0, Ly-dy, Ny)
X,Y = np.meshgrid(x,y)
X = X.T; Y = Y.T
IJ = np.ix_ # for convenient 2D indexing
u0 = np.zeros((Nx,Ny))
u0[IJ(np.arange(5,15), np.arange(10,30))] = 1.0
#--------- define PDE right-hand-side
def rates(uu, t):
u = np.reshape(uu, (Nx,Ny)) # convert 1D array to 2D for eval
n = np.arange(Nx); n[int(Nx/2)+1:]-= Nx
m = np.arange(Ny); m[int(Ny/2)+1:]-= Ny
d2udx2 = np.empty((Nx, Ny))
d2udy2 = np.empty((Nx, Ny))
for k in range(Ny): d2udx2[:,k] = ifft( (2*np.pi*1j*n/Lx)**2 * fft(u[:,k])).real
for j in range(Nx): d2udy2[j,:] = ifft( (2*np.pi*1j*m/Ly)**2 * fft(u[j,:])).real
dudt = G*(d2udx2 + d2udy2)
return np.reshape(dudt, shape=Nx*Ny)
def rates_FD(uu, t):
u = np.reshape(uu, (Nx,Ny)) # convert 1D array to 2D for eval
i = np.arange(0,Nx); ip=i+1; ip[-1]=0; im=i-1; im[0]=Nx-1
j = np.arange(0,Ny); jp=j+1; jp[-1]=0; jm=j-1; jm[0]=Ny-1;
dudt = G*( (u[IJ(im,j)] - 2.0*u[IJ(i,j)] + u[IJ(ip,j)])/dx**2 +
(u[IJ(i,jm)] - 2.0*u[IJ(i,j)] + u[IJ(i,jp)])/dy**2 )
return np.reshape(dudt, shape=Nx*Ny)
#--------- solve the ODE system (method of lines)
uu0 = np.reshape(u0, shape=Nx*Ny) # convert 2D array to 1D for solve
times = np.linspace(0,trun, ntimes)
uu = odeint(rates, uu0, times) # solve spectral
#uu = odeint(rates_FD, uu0, times) # solve finite difference (2nd order central)
Chebyshev pseudospectral methods
Here, we relate Equation 3 with \(\phi_n(x)=T_n(x)\) to the inverse discrete cosine transform (IDCT), and show examples of computing derivatives and solving PDEs.
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}\]
These are plotted below.
In general, a recurrence relation defines \(T_n\) given \(T_0\) and \(T_1\). That is,
\[T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x).\]
The domain of interest is \(-1\le x \le 1\) and the \(T_n\) are bounded between -1 and 1. This is due to the relation to the cosine function. The Chebyshev polynomials are defined as
\[T_n(\cos\theta) = \cos(n\theta).\]
This is seen using Euler’s identity,
\[(e^{i\theta})^n = e^{in\theta} = (\cos\theta + i\sin\theta)^n = \cos(n\theta) + i\sin(n\theta),\] which we can rearrange to \[\cos(n\theta) = e^{in\theta} - i\sin(n\theta) = (\cos\theta + i\sin\theta)^n -i\sin(n\theta).\]
Since \(\cos(n\theta)\) is real, the imaginary parts occuring from use of Euler’s identity and the above rearrangement cancel. We define
\[x = \cos\theta.\]
It is clear that \(T_0(x) = T_0(\cos\theta) = \cos(0\cdot\theta) = 1.\) Also, \(T_1(x) = T_1(\cos\theta) = \cos{\theta} = x.\)
The second Chebyshev polynomial is
\[\begin{align} T_2(x) &= T_2(\cos\theta) = \cos(2\theta) \\ &= (\cos\theta + i\sin\theta)^2 - i\sin(2\theta) \\ &= \cos^2\theta - \sin^2\theta + \cancel{[2i\cos\theta\sin\theta - i\sin(2\theta)]} \\ &= \cos^2\theta - (1-\cos^2\theta) \\ &= 2\cos^2\theta - 1 \\ &= 2x^2 - 1. \end{align}\]
Here, we are using the identity \(\cos^2\theta + \sin^2\theta = 1\), and the imaginary term in brackets is zero. Note that \(T_2 = 2xT_1 - T_0\). If we write \(y=\sin\theta\), then we can see that
\[T_n = (x+iy)^n -\cancel{i\sin(n\theta)}.\]
Evaluating this for higher \(n\), shows a pattern, such as
\[\begin{align} T_6 = \binom{6}{0}x^6y^0i^0 + &\cancel{\binom{6}{1}x^5y^1i^1} + \binom{6}{2}x^4y^2i^2 + \cancel{\binom{6}{3}x^3y^3i^3} + \\ &\binom{6}{4}x^2y^4i^4 + \cancel{\binom{6}{5}x^1y^5i^5} + \binom{6}{6}x^0y^6i^6. \end{align}\]
This pattern gives
\[ T_n = \sum_{k=0}^{n/2}\binom{n}{2k}x^{n-2k}(x^2-1)^k,\]
where \(y^2 = 1-x^2\) is used, that is, \(\sin^2\theta = 1-\cos^2\theta\).
The periodic, trigonometric character of the Chebyshev polynomials, and the relation between the angle \(\theta\) and its projection to \(x\), given by \(x=\cos\theta\) is shown in the figure below. \(\theta\) increases from 0 to \(\pi\), traversing the half circle once, while the angle \(n\theta\) moves faster, in proportion to \(n\). That is, for \(n=2\), \(n\theta\) varies from 0 to \(2\pi\) as \(\theta\) varies from 0 to \(\pi\). For \(n=3\), \(n\theta\) varies from 0 to 3$ as \(\theta\) varies from 0 to \(\pi\). \(T_n\) is \(\cos{n\theta}\), so \(T_n\) varies from 1 to -1 or -1 to 1 \(n\) times as \(\theta\) varies. The math is obvious but the point is to visualize the behavior. The figure below shows \(\cos(n\theta)\), which are the value of \(T_n\), plotted versus \(\theta\) for \(0\le\theta\le \pi\) on the left, and \(\cos(n\theta)\) plotted versus \(x=\cos\theta\) on the right. The right plot is just the plot of \(T_n(x)\) versus \(x\).
As angle \(\theta\) varies uniformly on the half-circle, the value \(T_n = \cos{n\theta}\) oscillates between 1 and -1 \(n\) times. We can think of this as a cosine wave wrapped around a half-cylinder. \(T_n\) is the vertical value. \(x=\cos\theta\) is the projection of the radial point on the circumference of the circle at angle \(\theta\) to the x-axis bisecting the circle. \(T_n(x)\) is the projection of the \(\cos(n\theta)\) function that wraps around the half-cylinder onto the vertical plane bisecting the cylinder. This is illustrated in the figure below. The left plots are angled so as to ssee the wrapped cosine functions. The left plots are angled to highlight the projection to the green vertical bisecting plane.
Chebyshev roots and extrema
The roots of the \(N^{th}\) order Chebyshev polynomial of the first kind are given by
\[x_j = \cos\theta_j,\] \[\theta_j = \frac{\pi}{N}(j+1/2),\]
with \(j=0,\,\ldots,\,N-1\). That is, \(\theta_j\) is uniformly distributed between and including \(\pi/2N\), and \(\pi-\pi/2N\). The roots are then just the projection \(x_j=\cos\theta_j\).
For example, for \(N=6\), we have \(T_6(x_j) = 32x^6 - 48x^4 + 18x^2-1 = 0\) as shown below. Note that these roots are ordered from high to low. Replace \(x_j=\cos\theta_j\) with \(x_j=-\cos\theta_j\) to order them from low to high.
For the \(N^{th}\) order Chebyshev polynomial of the first kind, the locations of the peak and minimum values, that is, the extrema, are at
\[x_j = \cos(\theta_j),\] \[\theta_j = \frac{\pi j}{N},\]
with \(j=0,\,\ldots,\,N\). These include \(x=1\) and \(x=-1\).
Chebyshev interpolation
Recall that we are appproximating our function \(u(t,x)\) as
\[ u(t,x) = \sum_{n=0}^{N-1} c_n(t)T_n(x). \]
We drop the time-dependence in the following.
For interpolation, that is, evaluating the function \(u(t,x)\) using the Chebyshev series approximation, the \(x_j\) are taken as either the roots grid or the extrema grid (noting that N is one larger for the extrema grid than the roots grid). Using the extrema grid is convenient for solving PDEs since, in that case, there are grid points on the boundaries at \(x=1\) and \(x=-1\).
Denote \(T_n(x_j)\) as \(T_{n,j}\). Then we then have
\[u_j = \sum_{n=0}^{N-1} c_nT_{n,j}.\]
This defines a matrix multiplication, \(u=T^Tc\), (that is, \(T\) transpose), with \(T_{n,j}\) components on row \(n\) and column \(j\):
\[ \left( \begin{matrix} u_0 \\ u_1 \\ u_2 \\ u_3 \\ u_4 \\ \end{matrix} \right) = \left[\begin{matrix} T_0(x_0) & T_1(x_0) & T_2(x_0) & T_3(x_0) \\ T_0(x_1) & T_1(x_1) & T_2(x_1) & T_3(x_1) \\ T_0(x_2) & T_1(x_2) & T_2(x_2) & T_3(x_2) \\ T_0(x_3) & T_1(x_3) & T_2(x_3) & T_3(x_3) \\ T_0(x_4) & T_1(x_4) & T_2(x_4) & T_3(x_4) \end{matrix}\right] \left( \begin{matrix} c_0 \\ c_1 \\ c_2 \\ c_3 \\ c_4 \end{matrix} \right). \]
The \(c_n\) can be found by inversion.
Example, Chebyshev interpolation
Boilerplate libraries and functions
For the code and examples below, use the following boilerplate libraries and functions
As an example of interpolation, consider the famous Runge function
\[u(x) = \frac{1}{1+25x^2}.\]
Use Chebyshev interpolation for \(N=9\) points using the roots grid, the extrema grid, and a uniform grid.
import numpy as np
import matplotlib.pyplot as plt
import numpy.polynomial.chebyshev as cheb
from scipy.special import eval_chebyt
from scipy.fft import fft, ifft
from scipy.fft import dct, idct, dctn, idctn
from scipy.integrate import odeint
import matplotlib.cm as cm
def urunge(x): # the Runge function for testing
return 1.0/(1.0+25*x*x)
def chebyterp(c,x):
return cheb.chebvander(x,len(c)-1) @ c # = sum_{n=0}^{N-1} c_nT_n(x)
xx = np.linspace(-1,1,1000) # lots of points for plotting
uu = urunge(xx)\[ \newcommand{\joinR}{\hspace{-.1em}} \newcommand{\RI}{\text{I}} \newcommand{\RII}{\text{I}\hspace{-.1em}\text{I}} \]
def make_plot(N, 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.plot(xx,urunge(xx), color='black', label='Exact', lw=2)
plt.plot(xx, chebyterp(cr,xx), '-', color='blue', label='Cheby roots')
plt.plot(xcr, ucr, 'o', ms=7, fillstyle='full', color='blue')
plt.plot(xx, chebyterp(ce,xx), '-', color='green', label='Cheby extrema')
plt.plot(xce, uce, '^', ms=7, fillstyle='full', color='green')
if do_unif:
plt.plot(xx, np.polyval(p,xx), '-', color='red', label='Uniform')
plt.plot(xu, urunge(xu), 's', ms=7, fillstyle='none', color='red')
plt.legend(frameon=False)
plt.xlabel('x')
plt.ylabel('u(x)');
#-------------
make_plot(9)
Note that the roots grid and extrema grid give similar results and similar errors. Interpolation using a uniform grid results in the Runge phenomena: large errors near the boundaries. These errors increase with increasing \(N\). Conversely, the Chebyshev interpolation points result in fast convergence with increasing \(N\). For example, for \(N=19\), we have the plot below. (The uniform grid gives a Runge peak near \(u(x)=30\), not shown.)
Discrete cosine transform, roots grid
The spectral coefficients \(c_n\) can be computed using the discrete cosine transform (DCT). Most fast Fourier transform librarier provide implementations of the fast DCT and its inverse (IDCT). There are several of discrete cosine transforms. The Type-II transform corresponds to Chebyshev interpolation on the roots grid, and is given by
\[\begin{align} &\text{DCT}_{\RII}:\phantom{Ixxx}a_n = 2\sum_{j=0}^{N-1} u_j\cos\left(\frac{\pi n}{N}(j+1/2)\right),\phantom{xx}n=0\ldots N-1, \\ &\text{IDCT}_{\RII}:\phantom{xxx}u_j = \frac{1}{N}\left[\frac{c_0}{2} + \sum_{n=1}^{N-1}a_n\cos\left(\frac{\pi n}{N}(j+1/2)\right)\right],\phantom{xx}j=0\ldots N-1. \end{align}\]
The subscript \(\RII\) denotes the Type-II transform. This is the form used in the scipy.fft.idct function by default. On Wikipedia, the factor of 2 on the DCT is moved to the IDCT. The leading \(c_0\) term in the IDCT can be incorporated into the sum and the IDCT written as
\[\text{IDCT}_{\RII}:\phantom{xxx}u_j = \frac{1}{N}\sum_{n=0}^{N-1}p_{\RII,n}a_n\cos\left(\frac{\pi n}{N}(j+1/2)\right), \tag{13}\]
where \[\phantom{xx} p_{\RII,0}=\frac{1}{2},\,p_{\RII,n\ne 0}=1.\]
The Chebyshev interpolation on the roots grid, denoted with subscript \(r\), is given by
\[u_j = \sum_{n=0}^{N-1}c_nT_n(x_{j,r}) = \sum_{n=0}^{N-1}c_n\cos\left(\frac{\pi n}{N}(j+1/2)\right),\]
Comparing this to Equation 13 shows that the Chebyshev interpolation is just the IDCT with \(c_n = a_np_{\RII,n}/N.\)
Denote the DCT and IDCT as \(\mathcal{C}\) and \(\mathcal{C}^{-1}\), respectively.
Summary:
\[u_j = \sum_{n=0}^{N-1}c_nT_n(x_{j,r}),\] \[c_n = \frac{a_np_{2,n}}{N}, \] \[\mathbf{a} = \mathcal{C}_{\RII}(\mathbf{u}), \] \[\mathbf{u} = \mathcal{C}_{\RII}^{-1}(\mathbf{a}).\]
Example, Runge function
Interpolate the Runge function using \(N=17\) points using a DCT.
N = 17
x = cheb.chebpts1(N) # roots grid
u = urunge(x)
#---------- the crux:
a = dct(u)
p = np.ones(N); p[0]=0.5
c = a*p/N
uu_cheby = chebyterp(c,xx) # interpolate the function everywhere
#----------- plot results (code not shown)
Domain mapping
Chebyshev interpolation is done on a domain \(-1\le x\le 1\). If we have a function on some other domain \(\hat{x}\) with \(a\le\hat{x}\le b\), we can do a simple linear mapping between the two domains. That is,
\[x = 2\left(\frac{\hat{x}-a}{b-a}\right) - 1,\] \[\hat{x} = (x+1)\left(\frac{b-a}{2}\right) + 1.\]
Discrete cosine transform, extrema grid
The Type-I DCT is given by
\[\begin{align} &\text{DCT}_\RI:\phantom{Ixxx}a_n = 2\sum_{j=0}^{N-1}u_j\cos\left(\frac{\pi nj}{N-1}\right),\phantom{xx}n=0\ldots N-1, \\ &\text{IDCT}_\RI:\phantom{xxx}u_j = \frac{1}{N-1}\sum_{n=0}^{N-1}a_n p_{\RI,n}\cos\left(\frac{\pi nj}{N-1}\right), \phantom{xx}j=0\ldots N-1, \end{align}\]
where \(p_{\RI,n}=1/2\) for \(n=0\) and \(n=N-1\), and \(p_{\RI,n}=1\) otherwise. The Chebyshev interpolation on the extrema grid, denoted by subscript \(e\), is given by
\[u_j = \sum_{n=0}^{N-1}c_nT_n(x_{j,e}) = \sum_{n=0}^{N-1}c_n\cos\left(\frac{\pi nj}{N-1}\right).\]
Comparing to the IDCT\(_1\) shows that \(c_n=a_np_{\RI,n}/(N-1).\)
Summary:
\[u_j = \sum_{n=0}^{N-1}c_nT_n(x_{j,e}),\] \[c_n = \frac{a_np_{\RI,n}}{N-1}, \] \[\mathbf{a} = \mathcal{C}_\RI(\mathbf{u}), \] \[\mathbf{u} = \mathcal{C}_\RI(\mathbf{a}).\]
Example, Runge function with DCT
Interpolate the Runge function using \(N=17\) points using a DCT on the extrema grid. The key differences with the example on the roots grid is using the type=1 argument in the dct function, and modifying the \(p\) and \(c\) arrays.
N = 17
x = cheb.chebpts2(N) # extrema grid
u = urunge(x)
#---------- the crux:
a = dct(u, type=1)
p = np.ones(N); p[[0,N-1]]=0.5
c = a*p/(N-1)
uu_cheby = chebyterp(c,xx) # interpolate the function everywhere
#---------- plot results (code not shown)
Derivatives
The derivative \(\partial u(x)/\partial x\) can be expressed as a Chebyshev series:
\[\frac{\partial u}{\partial x} = \sum_{n=0}^{N-1}c_n^\prime T_n(x),\]
where \(c_n^\prime\) are coefficients of the Chebyshev series of the derivative. These are computed using a recurrence relation from high to low \(n\); see Boyd “Chebyshev and Fourier Spectral Methods,” 2\(^{nd}\) edition, page 298.
\[(c^\prime_N=0),\] \[c^\prime_{N-1}=0,\] \[c^\prime_{n} = \gamma_n(c^\prime_{n+2} + 2(n+1)c_{n+1}),\;\; n=N-2,\ldots 0,\]
where \(\gamma_0=0.5\) and \(\gamma_{n\ne 1}=1\). In numpy.polynomial.chebyshev, the chebder function returns the size \(N-1\) array \(c^\prime\); the last coefficient (for \(n=N-1\)) is not included since it is zero, so it is convenient to append that value. Note that chebder operates on \(c_n\) and is independent of the \(x\) grid used.
N = 5
c = np.array([0.1, 0.2, -3, 4, 5])
cp = np.zeros(N)
for n in range(N-2,-1,-1):
if n < N-2: cp[n] = cp[n+2]
cp[n] += 2*(n+1)*c[n+1]
cp[0] *= 0.5
#-------- compare direct calculation to numpy.polynomial.chebyshev.chebder
print(cp)
print(np.append(cheb.chebder(c), 0.0))[12.2 28. 24. 40. 0. ]
[12.2 28. 24. 40. 0. ]Example, Runge derivative
Interpolate the derivative of the Runge function for \(N=17\) points. The Runge function is particularly challenging. For variety, the function \(u_2(x) = 0.1e^{-2x}\sin(6x) + 6x\) is also shown.
N = 17
x = cheb.chebpts1(N)[::-1] # roots grid
u = urunge(x)
p = np.ones(N); p[0] = 0.5
c = p/N*dct(u)
cp = np.append(cheb.chebder(c), 0.0)
dudx_cheby = chebyterp(cp, xx) # cheby interp derivative everywhere
#---------- plot results (code not shown)
Derivatives on the grid
Denote \(c^\prime = \mathcal{D}(c)\). Then we can write
\[\frac{\partial \mathbf{u}}{\partial x} = \mathcal{C}_{\RII}^{-1}\left[\frac{N}{\mathbf{p_{\RII}}}\odot\mathcal{D}\left(\frac{\mathbf{p_{\RII}}}{N}\odot\mathcal{C}_{\RII}(\mathbf{u})\right) \right],\]
on the roots grid, and
\[\frac{\partial \mathbf{u}}{\partial x} = \mathcal{C}_\RI^{-1}\left[\frac{N-1}{\mathbf{p_\RI}}\odot\mathcal{D}\left(\frac{\mathbf{p_\RI}}{N-1}\odot\mathcal{C}_\RI(\mathbf{u})\right) \right],\]
on the extrema grid.
Here, \(\mathbf{p_{\RII}}\) and \(\mathbf{p_\RI}\) are vectors with components \(p_{\RII,n}\) or \(p_{\RI,n}\). Using fast cosine transforms, these evaluations should have computational cost \(\mathcal{O}(N\log N)\). Note that \(N/\mathbf{p_{\RII}}\) and \((N-1)/\mathbf{p_\RI}\), are computed elementwise.
Second derivative
A second derivative would be
\[\frac{\partial^2 \mathbf{u}}{\partial x^2} = \mathcal{C}_\RI^{-1}\left[\frac{N-1}{\mathbf{p_\RI}}\odot\mathcal{D}^2\left(\frac{\mathbf{p_{\RI}}}{N-1}\odot\mathcal{C}_\RI(\mathbf{u})\right) \right],\]
on the extrema grid, and similarly on the roots grid with subscript 2 replacing subscript 1 and \(N\) replacing \(N-1\). In this equation, \(\mathcal{D}^2\) indicates two successive applications of the spectral derivative. Note that the \(N-1\) factors cancel, giving
\[\frac{\partial^2 \mathbf{u}}{\partial x^2} = \mathcal{C}_\RI^{-1}\left[\frac{1}{\mathbf{p_\RI}}\odot\mathcal{D}^2\left(\mathbf{p_\RI}\odot\mathcal{C}_\RI(\mathbf{u})\right) \right].\]
Example, 1D heat equation
Solve the 1D heat equation for \(u(t,x)\) with Dirichlet boundary conditions and a zero initial condition. Let \(S=1.\)
\[\frac{\partial u}{\partial t} = \alpha\frac{\partial^2u}{\partial x^2} + S,\] \[u(t,-1)=0,\] \[u(t,1)=1,\] \[u(0,x)= 0.\]
The pseudospectral solution is given by
\[\frac{\partial \mathbf{u}}{\partial t} = \alpha\mathcal{C}_\RI^{-1}\left[\frac{1}{\mathbf{p_\RI}}\odot\mathcal{D}^2\left(\mathbf{p_\RI}\odot\mathcal{C}_\RI(\mathbf{u})\right) \right] + S,\]
which is solved in time using the method of lines.
To include the boundary conditions we solve the problem on the extrema grid, which includes \(x=-1\) and \(x=1\). The array \(\mathbf{u}\) will have the apppropriate boundary values, and the rates at the boundaries, \(\partial u_0/\partial t\), \(\partial u_{N-1}/\partial t\) are set to zero.
def chebyspec(N, trun):
α = 1.0
S = 1.0
uinit = np.zeros(N)
uinit[0] = 1.0 # Dirichlet BC at x=1
uinit[-1] = 0.0 # Dirichlet BC at x=-1
#---------------------
def rates(u,t):
p = np.ones(N); p[0]=p[-1]=0.5
#-------- d2u/dx2
a = dct(u, type=1)
c = a*p
cpp = np.append(cheb.chebder(c, 2),[0.0,0.0])
app = cpp / p
upp = idct(app, type=1)
#-------- compute the du/dt
dudt = α*upp + S
dudt[0] = dudt[-1] = 0 # don't change BC (make du/dt = 0)
return dudt
#---------------------
times = np.linspace(0,trun,9)
u = odeint(rates, uinit, times)
return times, u
#---------------------
N = 20
trun = 2.0
times, u = chebyspec(N, trun)
x = cheb.chebpts2(N)[::-1]
#--------------------- plot results (code not shown)
Neumann BCs, example
Modify the boundary condition at \(x=1\) to be the Neumann condition \(u^\prime(t,x=1)=1\). Use a linear initial condition that satisfies the boundary conditions.
Rewrite the PDE as \[\frac{\partial u}{\partial t} = \alpha\frac{\partial}{\partial x}\left(\frac{\partial u}{\partial x}\right) + S.\]
The inner derivative is given by
\[\frac{\partial \mathbf{u}}{\partial x} \equiv\mathbf{u}^\prime = \mathcal{C}_\RI^{-1}\left[\frac{1}{\mathbf{p_\RI}}\odot\mathcal{D}\left(\mathbf{p_\RI}\odot\mathcal{C}_\RI(\mathbf{u})\right) \right].\]
A Neumann condition \(\partial u/\partial x=\beta\) for some given \(\beta\) is imposed on the first or last element of the \(u^\prime\) array once it is computed.
We then write the pseudospectral solution as
\[\frac{\partial \mathbf{u}}{\partial t} = \alpha\mathcal{C}_\RI^{-1}\left[\frac{1}{\mathbf{p_\RI}}\odot\mathcal{D}\left(\mathbf{p_\RI}\odot\mathcal{C}_\RI(\mathbf{u}^\prime)\right) \right] + S.\]
def chebyspec_Neumann(N):
α = 1.0
trun = 3 * 2.0/α
S = 1.0
x = cheb.chebpts2(N)[::-1]
uinit = np.zeros(N) # initial condition that satisfies u(-1)=0, u'(1)=1
for i in range(N):
uinit[i] = 2.0 + (x[i] - x[0])
#---------------------
def rates(u,t):
p = np.ones(N); p[0]=p[-1]=0.5
#-------- compute uprime
a = dct(u, type=1)
c = a*p
cp = np.append(cheb.chebder(c), 0.0)
ap = cp / p
up = idct(ap, type=1)
up[0] = 1.0 # impose Neumann BC at x=1
#-------- d2u/dx2
ap = dct(up, type=1)
cp = ap*p
cpp = np.append(cheb.chebder(cp), 0.0)
app = cpp/p
upp = idct(app, type=1)
#-------- compute the du/dt rates
dudt = α*upp + S
dudt[-1] = 0 # don't change the Dirichlet BC (make du/dt = 0)
return dudt
#---------------------
times = np.linspace(0,trun,9)
u = odeint(rates, uinit, times)
return times, u
#--------------------- solve the problem
N = 20
times, u = chebyspec_Neumann(N)
x = cheb.chebpts2(N)[::-1]
#--------------------- plot results (code not shown)
Robin BCs
A Robin condition is done the same way as the Neumann condition. For Neumann, we impose \(\partial u/\partial x=\beta\) at the boundary, for some given \(\beta\) that defines the Neumann boundary condition. For Robin, since we know \(u\) at the boundary (at the given time), we can impose \(\partial u/\partial x=\gamma u+\beta\) at the boundary, where \(\gamma\) and \(\beta\) are some constants for the Robin condition.
2D transforms
The Type-II cosine transforms (for the roots grid) in two dimensions are given by
\[\begin{align} &\text{IDCT}_{\RII}:\phantom{xxx}u_{j,k} = \frac{1}{N_x}\sum_{n=0}^{N_x-1}p_{\RII,n}\left[ \frac{1}{N_y}\sum_{m=0}^{N_y-1}p_{\RII,m}a_{n,m}\cos\left(\frac{\pi m}{N_x}(k+1/2)\right) \right]\cos\left(\frac{\pi n}{N_y}(j+1/2)\right), \\ &\text{DCT}_{\RII}:\phantom{Ixxx}a_{n,m} = 2\sum_{j=0}^{N_x-1}\underbrace{\left[ 2\sum_{k=0}^{N_y-1}u_{j,k}\cos\left(\frac{\pi m}{N_x}(k+1/2)\right) \right]}_{A_j}\cos\left(\frac{\pi n}{N_y}(j+1/2)\right). \\ \end{align}\]
This presentation of the transforms makes it clear that the 2D tranform is just a nested 1D transform. For the 2D DCT, the term in brackets is a 1D DCT over grid points \(k\) for given \(j\). This results in \(A_j\), which is an array with elements corresponding to \(m=0,\,\ldots,\,N_y-1\), for each \(j\). That is, \(A\) is a matrix with elements \(A_{j,m}\). The outer sum is also a 13 DCT over grid points \(j\) for each \(m\).
The Type-I cosine transforms (for the extrema grid) in two dimensions are given by
\[\begin{align} &\text{IDCT}_{\RI}:\phantom{xxx}u_{j,k} = \frac{1}{N_x-1}\sum_{n=0}^{N_x-1}p_{\RI,n}\left[ \frac{1}{N_y-1}\sum_{m=0}^{N_y-1}p_{\RI,m}a_{n,m}\cos\left(\frac{\pi mk}{N_x-1}\right) \right]\cos\left(\frac{\pi nj}{N_y-1}\right), \\ &\text{DCT}_{\RI}:\phantom{Ixxx}a_{n,m} = 2\sum_{j=0}^{N_x-1}\underbrace{\left[ 2\sum_{k=0}^{N_y-1}u_{j,k}\cos\left(\frac{\pi mk}{N_x-1}\right) \right]}_{A_j}\cos\left(\frac{\pi nj}{N_y-1}\right). \\ \end{align}\]
type = 1
#---------- 2D DCT, type I or II
def dct2D(u):
a = np.empty_like(u)
Nx,Ny = np.shape(u)
for j in range(Nx):
a[j,:] = dct(u[j,:],type=type)
for m in range(Ny):
a[:,m] = dct(a[:,m],type=type)
return a
#---------- 2D IDCT, type II
def idct2D(a):
u = np.empty_like(a)
Nx,Ny = np.shape(a)
for n in range(Nx):
u[n,:] = idct(a[n,:],type=type)
for k in range(Ny):
u[:,k] = idct(u[:,k],type=type)
return u
#---------- (These give identical results to dctn and idctn)
a = np.random.rand(4,3)
np.allclose(dct2D(a), dctn(a, type=type))TrueDerivatives
On a two-dimensional grid, partial derivatives in a direction, say \(x\), are computed along a line of \(x\) at a given \(y\) exactly as for one-dimensional cases. On the extrema grid, we have
\[\frac{\partial \mathbf{u}_k}{\partial x} = \mathcal{C}_\RI^{-1}\left[\frac{1}{\mathbf{p_\RI}}\odot\mathcal{D}\left(\mathbf{p_\RI}\odot\mathcal{C}_\RI(\mathbf{u}_k)\right) \right],\]
where \(\mathbf{u}_k\) is the array of \(u\) values at all \(j\) (in the \(x\) direction) at a given \(k\) (\(y\) direction).
Example, 2D heat equation
Solve the 2D heat equation with zero Dirichlet boundaries and a unity initial condition in the domain center. We have \(u=u(t,x,y)\), \(\Gamma=1\), \(L_x=L_y=2\) and the domains \(x\) and \(y\) range from -1 to 1.
\[\frac{\partial u}{\partial t} = \Gamma\frac{\partial^2u}{\partial x^2} +\Gamma\frac{\partial^2u}{\partial x^2}.\]
Solutions using both psuedospectral and finite difference with a second-order central difference approximation are implemented. In Python, for \(N_x=N_y=21\), the psuedospectral method is considerably slower.
doFD = False
#--------- setup
Lx = 2.0
Ly = 2.0
Nx = 21
Ny = 21
G = 1.0
trun = 0.2
ntimes = 11
if doFD:
x = np.linspace(-1.0, 1.0, Nx)
y = np.linspace(-1.0, 1.0, Ny)
else:
x = cheb.chebpts2(Nx)[::-1]
y = cheb.chebpts2(Ny)[::-1]
X,Y = np.meshgrid(x,y)
X = X.T; Y = Y.T
u0 = np.zeros((Nx,Ny))
for j in range(Nx):
for k in range(Ny):
if x[j]>=-0.5 and x[j]<=0.5 and y[k]>=-0.5 and y[k]<=0.5 : u0[j,k] = 1.0
#--------- define PDE for psuedospectral
def rates_cheby(uu, t):
u = np.reshape(uu, (Nx,Ny)) # convert 1D array to 2D for eval
d2udx2 = np.empty((Nx, Ny))
d2udy2 = np.empty((Nx, Ny))
for k in range(Ny):
c = dct(u[:,k], type=1); c[0]*=0.5; c[-1]*=0.5
app = np.append(cheb.chebder(c, 2),[0,0]); app[0]*=2; app[-1]*=2
d2udx2[:,k] = idct(app, type=1)
for j in range(Nx):
c = dct(u[j,:], type=1); c[0]*=0.5; c[-1]*=0.5
app = np.append(cheb.chebder(c, 2),[0,0]); app[0]*=2; app[-1]*=2
d2udy2[j,:] = idct(app, type=1)
dudt = G*(d2udx2 + d2udy2)
dudt[:,0] = dudt[:,-1] = 0.0
dudt[0,:] = dudt[-1,:] = 0.0
return np.reshape(dudt, shape=Nx*Ny)
#--------- define PDE for finite difference
def rates_FD(uu, t):
u = np.reshape(uu, (Nx,Ny)) # convert 1D array to 2D for eval
i = np.arange(1,Nx-1); ip=i+1; im=i-1
j = np.arange(1,Ny-1); jp=j+1; jm=j-1
dx = x[1]-x[0]
dy = y[1]-y[0]
dudt = np.zeros((Nx,Ny))
dudt[IJ(i,j)] = G*( (u[IJ(im,j)] - 2.0*u[IJ(i,j)] + u[IJ(ip,j)])/dx**2 +
(u[IJ(i,jm)] - 2.0*u[IJ(i,j)] + u[IJ(i,jp)])/dy**2 )
dudt[:,0] = dudt[:,-1] = 0.0
dudt[0,:] = dudt[-1,:] = 0.0
return np.reshape(dudt, shape=Nx*Ny)
#--------- solve the ODE system (method of lines)
uu0 = np.reshape(u0, shape=Nx*Ny) # convert 2D array to 1D for solve
times = np.linspace(0,trun, ntimes)
rates = rates_FD if doFD else rates_cheby
uu = odeint(rates, uu0, times) # solve the problem
#--------- (plot solution at three times, not shown)u(0.2,0,0) = 0.328780332997131
Notes
Note that for explicit solvers as have been discussed in the examples, the stable time step size \(\tau\) for advective terms (first derivatives) is \(\tau\sim\Delta x\). For diffusive terms (second derivatives) \(\tau\sim\Delta x^2\). Fourier pseudospectral methods use a uniform grid so that \(\Delta x\sim 1/N\) and \(\tau\sim 1/N\) for advection and \(\tau\sim 1/N^2\) for diffusion. But Chebyshev interpolation uses a nonuniform grid with \(\Delta x = 1-\cos(\pi/N)\) for \(N+1\) points. The first two terms of a Taylor expansion of \(\cos(x)\) are \(1-x^2/2\), so that \(\Delta x\sim 1/N^2\), and \(\tau\sim 1/N^2\) for advection and \(\tau\sim 1/N^4\) for diffusion. This is a significant constraint. Boyd (“Chebyshev and Fourier Spectral Methods”) discusses using implicit or semi-implicit approaches.
In pseudospectral methods, the Fourier or Chebyshev basis functions are used to evaluate derivatives accurately. From that perspective, it is not so important whether the domain has periodic or non-periodic boundary conditions. For example, a uniform grid could be used on a problem with Dirichlet boundaries and derivatives computed using a Fourier basis if the domain were mirrored so as to make continuous periodic profiles appropriate for the Fourier basis. The accurate derivatives could then be used to advance a time-dependent problem.
Consider a vertical jet in a free stream with nonperiodic far-stream side boundaries. A Chebyshev pseudospectral method would use a grid with the highest clustering of points near the free-stream boundaries where nothing is happening. This could be mitigated by taking the domain as periodic since it will be periodically continuous (given the free stream on either side), and then derivatives computed using a uniform grid and a Fourier basis. In pipe flows, most of the action happens near the walls and a finer clustering near the wall may be desired.