Represent the solution $u(x,t)$ as $$u(x,t) = \sum_n \hat{u}_n(t)\phi_n(x),$$ where $\phi_n$ are some basis functions and $\hat{u}_n$ are corresponding weighting factors.
In a spectral method we take $\phi_n$ as $e^{2\pi inx/L}=\cos(2\pi nx/L)+i\sin(2\pi nx/L)$, where $i$ is the imaginary number and $L$ is the domain size.
Here, $\phi_n$ are periodic functions and we assume periodic boundary conditions.
This gives $u$ as a discrete inverse Fourier transform, and we take a (truncated) summation from $n=-N/2+1$ to $n=N/2$.
This results in
$$u = \sum_{n=-N/2+1}^{N/2} \hat{u}_n e^{2\pi inx/L}.$$In particular, consider a grid of N points with uniform spacing $\Delta x$. If the points are indexed by $j$, starting at $j=0$, we have $\Delta x=L/N$, and $x_j=j\Delta x=jL/N$. Let $u_j\equiv u(x_j)$.
The exponent of $e$ in the IDFT (and DFT) is an angle in the complex plane, and the multiplicative factor $n$ in the exponent effectively steps the angle around the circle in the complex plane. The figure below illustrates this for an even and an odd number of points.
The vertical dashed lines connect points with conjugate symmetry (the imaginary parts are the same in magnitude but have opposite signs). If $u$ is real, then $\hat{u}$ will have conjugate symmetry so that the IDFT of $\hat{u}$ will give cancellation of the imaginary parts, yielding the real-valued $u$.
For numerical solution, it is important that the ordering of the indicies $n$ match the ordering provided by the fast Fourier transform (and its corresponding inverse).
This is illustrated below in Python for five points. This shows that the points are ordered starting at $n=0$ and traversing counter-clockwise: $$n=0,\,1,\,2,\,-2,\,-1.$$ For general $N$, we can write (Python)
n = np.arange(N)
n[int(N/2)+1:]-= N
from scipy.fft import fft, ifft
import numpy as np
u = np.random.rand(5)
uhat = fft(u)
for i in uhat:
print(f"{i.real:.4f} + {i.imag:.4f}j")
3.2301 + -0.0000j -1.0486 + -0.2433j -0.2637 + -0.2236j -0.2637 + 0.2236j -1.0486 + 0.2433j
The inverse discrete Fouier transform (IDFT) of $\hat{u}_n$ evaluated at grid points $x_j$ and denoted $F^{-1}$, is given by $$u_j = F^{-1}(\hat{u}_n):\phantom{xx} u_j = \sum_{n=-N/2+1}^{N/2} \hat{u}_n e^{2\pi inj/N},\;\;\;j=0,\ldots,N-1.$$ The corresponding discrete Fourier transform (DFT) is given by $$\hat{u}_n=F(u_j):\phantom{xx} \hat{u}_n = \frac{1}{N}\sum_{j=0}^{N-1} u_j e^{-2\pi inj/N},\;\;\;n=-\frac{N}{2}+1,\ldots,\frac{N}{2}.$$
These are evaluated using packaged fast Fourier transform (FFT) routines (and the corresponding fast inverse).
Now we will insert the green equation into the PDE. Note that $\hat{u}$ is not a function of position, and $d e^{2\pi inx/L}/d x=(2\pi in/L)e^{2\pi inx/L}.$ Furthermore, the summation commutes with differentiation. The result, omitting the bounds of the summation is $$\sum_n e^{2\pi inx/L}\frac{d\hat{u}_n}{d t} = a\sum_n (2\pi in/L)\hat{u}_ne^{2\pi inx/L}.$$ This equation can be rewritten as $$\sum_n e^{2\pi inx/L}\left(\frac{d\hat{u}_n}{d t} - a\cdot(2\pi in/L)\hat{u}_n\right) = 0.$$ The functions $e^{2\pi inx/L}$ are orthogonal over the domain $L$, so each coefficient of these functions (which are independent of $x$) is zero, hence
$$\frac{d\hat{u}_n}{d t} = a\cdot(2\pi in/L)\hat{u}_n.$$ Given some initial condition for $u_j$, we can use the DFT to compute the initial $\hat{u}_n$. This ODE can then be solved for $\hat{u}_n(t)$. Then $u_j(t)$ can be computed from $\hat{u}_n(t)$ using the IDFT. Note, for this particular problem, the ODE has a simple analytic solution: $$\hat{u}_n(t) = \hat{u}_n(0)e^{(a\cdot 2\pi in/L)t}.$$Consider the same PDE, but with $a=a(x)$ instead of constant $a$. The previous approach won't work: in the development, we assumed $\hat{u}$ was independent of $x$, but with $a(x)$ the previous solution for $\hat{u}$ is a function of $x$, contradicting the assumption. Put another way, in the blue equation above, we have transformed from the physical space to the spectral space, but $a(x)$ remains in the physical space.
Consider again the PDE. While it remains linear with $a=a(x)$, the right-hand-side involves a product of functions of $x$, which is difficult to handle with with Fourier transforms. Similar issues arise from nonlinear terms like $u\partial u/\partial x$.
Here, we solve the problem in the physical space, but evaluate derivatives in the spectral space. On the right-hand-side of the PDE we write $u$ in terms of it's discrete Fourier transform,
\begin{align} \frac{\partial u}{\partial t} &= a\frac{\partial u}{\partial x}, \\ &= a\frac{\partial}{\partial x}\left(\sum_{n=N/2+1}^{N/2}\hat{u}_ne^{2\pi inx/L}\right), \\ &= a\sum_{n=N/2+1}^{N/2}\hat{u}_n\frac{2\pi in}{L}e^{2\pi inx/L}. \end{align}At this point, the pseudo-spectral method (also called the collocation method) evaluates this equation at grid points $x_j$: \begin{align} \frac{du_j}{d t} &= a\sum_{n=N/2+1}^{N/2}\hat{u}_n\frac{2\pi in}{L}e^{2\pi inj/N}. \end{align} We swap the order of the terms $\hat{u}_n$ and $2\pi in/L$, and write $\hat{u}_n$ in terms of $u_j$ using $\hat{u}_n=F(u_j)$: \begin{align} \frac{du_j}{dt} &= a\sum_{n=N/2+1}^{N/2}\left(\frac{2\pi in}{L}F(u_j)\right)e^{2\pi inj/N}. \end{align} But, the right-hand-side of this equation is just the IDFT of the term in parenthesis:
\begin{align} \frac{du_j}{dt} &= aF^{-1}\left(\frac{2\pi in}{L}F(u_j)\right). \end{align} Note that $F(u_j)$ is an array indexed by $n$, which multiplies $2\pi in/L$ termwise. Similarly, $F^{-1}(\cdot)$ is an array indexed by $j$. The array $n$ varies from $-\frac{N}{2}+1$ to $\frac{N}{2}$.Solve the PDE with the approach of Case 2, but with constant a=1.
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import odeint
L = 10.0
nx = 8
a = 1.0
ntflow = 1 # # flow-through times; result looks perfect for ntflow integer for any nx
# but requires more points to look good for non-integer ntflow
tend = ntflow*L/a
#---------- solution domain, initial condition
dx = L/nx # not L/(nx-1)
x = np.linspace(0.0, L-dx, nx)
u0 = 1.0/(np.cosh(2.0*(x-5.0)))
#---------- exact solution with more points for plotting
xx = np.linspace(0.0, L, 1000)
uu = 1.0/(np.cosh(2.0*(xx-5.0)))
#---------- solve the problem and plot
def rhsf(u, t):
N = len(u)
n = np.arange(N); n[int(N/2)+1:]-= N
i = 1j # 1j is the imaginary number in python; copy as i for clarity
return -a*ifft(2*np.pi*i*n/L*fft(u)).real
t = np.linspace(0,tend,11)
y = odeint(rhsf, u0, t)
fig,ax=plt.subplots()
ax.plot(xx,uu, '-', color='red', lw=2)
ax.plot(x,y[-1,:], 'o', color='gray', lw=1)
plt.rc('font', size=14)
plt.gca().set_xlim([0,10])
plt.xlabel('x')
plt.ylabel(r'u(x)');
plt.legend(['exact', 'spectral'], frameon=False);
Solve the viscous Burgers' equation: $$\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.
The exact solution (from wikipedia) is given by $$u(x,t) = -2\nu\frac{\partial}{\partial x}\ln\left\{(4\pi\nu t)^{-1/2}\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.
See also this link
import sympy as sp
ν, t, L, x, xp = sp.symbols('ν, t, L, x, xp')
#ex = sp.integrate(sp.cos(2*sp.pi*x/L)+2, (x,0,xp))/2/ν
ex = L/(4*sp.pi*ν)*sp.sin(2*sp.pi*xp/L) + xp/ν
ex = -2*ν*sp.diff(sp.ln(1/sp.sqrt(4*sp.pi*ν*t)*
sp.Integral(sp.exp(-(x-xp)**2/4/ν/t - ex),
(xp,-sp.oo, sp.oo))), x)
ex = ex.subs([(ν,0.1),(t,10.0),(L,10.0)]) # make sure these are the same as below
f = sp.lambdify(x, ex)
xe = np.linspace(0,10,100)
ue = np.empty(len(xe))
for i in range(len(xe)):
ue[i] = f(xe[i])
def spectral(nx):
L = 10.0
ν = 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 rhsf(u, t):
N = len(u)
n = np.arange(N); n[int(N/2)+1:]-= N
return -u*ifft(2*np.pi*1j*n/L*fft(u)).real - ν*ifft((2*np.pi*n/L)**2*fft(u)).real
t = np.linspace(0,tend,11)
u = odeint(rhsf, u0, t)
return x, u, u0
x,u, u0 = spectral(32)
fig,ax=plt.subplots()
ax.plot(x,u0, ':', color='blue', lw=1)
ax.plot(xe,ue, '-', color='red', lw=2)
ax.plot(x,u[-1,:], 'o', color='gray', lw=1)
plt.rc('font', size=14)
plt.gca().set_xlim([0,10])
plt.xlabel('x')
plt.ylabel(r'u');
plt.legend(['initial', 'exact', 'spectral'], frameon=False);
def FD(nx, upwind=False):
L = 10.0
ν = 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 rhsf(u, t):
if upwind:
return -u*(u - u[im])/dx + ν/dx/dx*(u[im] - 2*u + u[ip]) # upwind on the advective term: 1st order
else:
return -u*(u[ip] - u[im])/2/dx + ν/dx/dx*(u[im] - 2*u + u[ip]) # central difference: 2nd order
t = np.linspace(0,tend,11)
u = odeint(rhsf, u0, t)
return x, u, u0
nxs = np.array([2,4,6,8,10,14,20,30,45,60,80,128,256,512])
errSP = np.empty(len(nxs))
errFD = np.empty(len(nxs))
errFDUW = np.empty(len(nxs))
for i,nx in enumerate(nxs):
x,usp,u0 = spectral(nx)
x,ufd,u0 = FD(nx, upwind=False)
x,ufduw,u0 = FD(nx, upwind=True)
ue = np.empty(len(x))
for j in range(len(ue)):
ue[j] = f(x[j])
errSP[i] = np.linalg.norm((usp[-1,:]-ue)/ue)/nx
errFD[i] = np.linalg.norm((ufd[-1,:]-ue)/ue)/nx
errFDUW[i] = np.linalg.norm((ufduw[-1,:]-ue)/ue)/nx
#---------- Spectral, and show exponential convergence: fit err=a*exp(b*nx)
plt.loglog(nxs,errSP, 'bo', label="spectral")
b_loga = np.polyfit(nxs[2:11], np.log(errSP[2:11]), 1)
a,b = np.exp(b_loga[1]), b_loga[0]
xf = np.logspace(0,2,100)
plt.loglog(xf, a*np.exp(b*xf), ':', color='blue', label="");
#---------- FD, and show power law convergence: fit err=a*nx^b
plt.loglog(nxs,errFD, 'gs', label="FD central")
b_loga = np.polyfit(np.log(nxs[-4:]), np.log(errFD[-4:]), 1)
a,b = np.exp(b_loga[1]), b_loga[0]
xf = np.logspace(0,3,100)
plt.loglog(xf, a*xf**b, ':', color='green', label="");
plt.loglog(nxs,errFDUW, 'r^', label="FD upwind")
b_loga = np.polyfit(np.log(nxs[-4:]), np.log(errFDUW[-4:]), 1)
a,b = np.exp(b_loga[1]), b_loga[0]
xf = np.logspace(0,3,100)
plt.loglog(xf, a*xf**b, ':', color='red', label="");
#----------
plt.xlabel('# points')
plt.ylabel('Average Relative Error')
plt.xlim([1,1000])
plt.ylim([1E-10,1])
plt.legend(frameon=False);