Spectral Methods

Case 1

Consider the following PDE with constant $a$: $$\frac{\partial u}{\partial t}=a\frac{\partial u}{\partial x}.$$

Represent the solution $u(x,t)$ as $$u(x,t)=\sum _n{\hat{u}}_n(t){\varphi }_n(x),$$ where ${\varphi }_n$ are some basis functions and ${\hat{u}}_n$ are corresponding weighting factors.

In a spectral method we take ${\varphi }_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.

• Note that $r{e}^{i\theta }=r\cos (\theta )+ir\sin (\theta )$ so that polar coordinates $r,,\theta$ in the complex plane map to Cartesian components $x=r\cos (\theta )$ and $y=r\sin (\theta )$.

Here, ${\varphi }_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 $\mathrm{\Delta }x$. If the points are indexed by $j$, starting at $j=0$, we have $\mathrm{\Delta }x=L/N$, and $x_j=j\mathrm{\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$.

Conjugate symmetry is illustrated below in Python for six points.

from scipy.fft import fft, ifft, fftshift
import numpy as np

u = np.random.rand(6)
uhat = ifft(u)
for i in uhat:
    print(f"{i.real:.4f} + {i.imag:.4f}i")

0.3354 + -0.0000i
0.0999 + -0.1277i
-0.0076 + -0.0542i
0.1084 + -0.0000i
-0.0076 + 0.0542i
0.0999 + 0.1277i

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).

The points are ordered starting at $n=0$ and traversing counter-clockwise: $$n=0,,1,,2,,3,,-2,,-1.$$ For general $N$, we can write (Python)

n = np.arange(N)
n[int(N/2)+1:]-= N

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,\dots ,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,\dots ,\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}/dx=(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}{dt}=a\sum _n(2\pi in/L){\hat{u}}_n{e}^{2\pi inx/L}.$$ This equation can be rewritten as $$\sum _n{e}^{2\pi inx/L}(\frac{d{\hat{u}}_n}{dt}-a\cdot (2\pi in/L){\hat{u}}_n)=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}{dt}=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}.$$

Case 2

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$.

Psuedo-spectral method

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 its discrete Fourier transform,

At this point, the pseudo-spectral method (also called the collocation method) evaluates this equation at grid points $x_j$: 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)$: But, the right-hand-side of this equation is just the IDFT of the term in parenthesis: 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}$.

Example 1

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);

Example 2

Solve the viscous Burgers’ equation: $$\frac{\partial u}{\partial t}=-u\frac{\partial u}{\partial x}+\nu \frac{{\partial }^{2}u}{\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 (from wikipedia ) is given by

$$u(x,t)=-2\nu \frac{\partial }{\partial x}\ln [(4\pi \nu t{)}^{-1/2}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\exp [-\frac{(x-x^{\prime }{)}^2}{4\nu t}-\frac{1}{2\nu }{\int }_0^{x^{\prime }}{u}_0(x^{\prime \prime })d{x}^{\prime \prime }]d{x}^{\prime }].$$

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])

Spectral solution

$$\frac{d{u}_j}{dt}=-u{F}^{-1}(\frac{2\pi in}{L}F(u_j))+\nu F^{-1}({\left(\frac{2\pi in}{L}\right)}^{2}F(u_j)),$$

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             # n[int(N/2):]-=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

Plot results

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);

Finite difference solver

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

Compare spectral and finite difference on multiple grids

• 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).

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);