Spectral Methods in Numerical Analysis
Method
For numerical solution of partial differential equations, instead of directly solving the PDE itself, one assumes a candidate solution as a Fourier series and circumvents the problem by obtaining a system of ordinary differential equations for the Fourier coefficients. This technique is called the spectral method.
Description
Regardless of which PDE we consider, under the assumption that the solution admits a Fourier expansion, we start from the idea that for a sufficiently large $N$ the solution is approximated by the Fourier series $u$: $$ u ( t , x ) = \mathcal{F} (u) \approx \sum_{k = -N}^{N} c_{k} (t) \exp \left( i k x \right) $$ Looking at what is known and unknown here, $\exp \left( i k x \right)$ corresponds to the $k$-th Fourier coefficient and can be computed, while $x$ is the value at grid points and is always known. In contrast, the $c_{k}$ at index $k$, which depends on time $c_{k} (t)$, is unknown. If these satisfy the governing equation expressed as an ODE system as follows, then they can be solved relatively straightforwardly using a numerical solver. $$ \begin{align*} {\frac{ d c_{-N} }{ dt }} =& f_{-N} \left( c_{-N} , \cdots , c_{N} \right) \\ \vdots & \\ {\frac{ d c_{N} }{ dt }} =& f_{N} \left( c_{-N} , \cdots , c_{N} \right) \\ \implies \dot{\mathbf{c}} =& \mathbf{f} \left( \mathbf{c} \right) \end{align*} $$
Now everything unknown can be summarized as the operator $\mathbf{f}$ that represents the system and the initial value $\mathbf{c} (0)$. As an example, consider solving the Kuramoto–Sivashinsky equation (Kuramoto–Sivashinsky equation) with a spectral method.
The equation starts as follows. Personally, I think this equation is an excellent example for studying spectral methods because it contains more than one linear term as well as a nonlinear term. $$ {\frac{ \partial u }{ \partial t }} + {\frac{ \partial }{ \partial x }} \left( {\frac{ 1 }{ 2 }} u^{2} \right) + {\frac{ \partial^{2} u }{ \partial x^{2} }} + {\frac{ \partial^{4} u }{ \partial x^{4} }} = 0 $$
Linear term 1
First, suppose the nonlinear term $u^{2}$ is absent. Assuming $u = \sum c \exp$ and substituting directly gives $$ \begin{align*} & {\frac{ d }{ d t }} \sum_{k = -N}^{N} c_{k} (t) \exp \left( i k x \right) \\ =& - {\frac{ \partial^{2} }{ \partial x^{2} }} \sum_{k = -N}^{N} c_{k} (t) \exp \left( i k x \right) \\ & - {\frac{ \partial^{4} }{ \partial x^{4} }} \sum_{k = -N}^{N} c_{k} (t) \exp \left( i k x \right) \end{align*} $$
Let us carry out all the partial differentiations. $$ \begin{align*} & \sum_{k = -N}^{N} \dot{c}_{k} (t) \exp \left( i k x \right) \\ =& - \sum_{k = -N}^{N} \left( ik \right)^{2} c_{k} (t) \exp \left( i k x \right) \\ & - \sum_{k = -N}^{N} \left( ik \right)^{4} c_{k} (t) \exp \left( i k x \right) \end{align*} $$
Orthogonality of the Fourier basis: The two sets $\left\{ e^{inx} \right\}_{n=-\infty}^\infty$ and $\left\{ \cos nx\ \right\}_{n=0}^\infty \cup \left\{ \sin nx \right\}_{n=1}^\infty$ form a $L^{2}(-\pi,\ \pi)$ orthonormal basis. Also, $\left\{ \cos nx \right\}_{n=0}^{\infty}$ and $\left\{ \sin nx \right\}_{n=1}^{\infty}$ are an orthonormal basis of $L^{2}(0,\ \pi)$.
Although the expression is written as a sum, if the KS equation had no nonlinear term, by the orthogonality of the Fourier basis one can solve separately for the $k$-th mode as follows. $$ \dot{c}_{k} (t) = k^{2} c_{k} (t) - k^{4} c_{k} (t) $$
Mathematically, no matter how many dimensions this has, it is not a difficult differential equation; however, from a numerical-analysis perspective the term $k^{4}$ can be quite burdensome. For example, if it is on the order of $N = 32$, it becomes a stiff problem in the sense of $k^{4} = 1,048,576 \approx 10^{6}$. Of course, this depends on the domain size, but it is natural that the fourth-power factor coming from the spatial derivatives in the Fourier expansion leaves a significant burden.
Intuitively, no matter what clever trick exists, the difficult PDE will not simply reduce to an ODE without cost. The load that would originally have been applied when approximating derivatives (as in finite difference methods) is effectively transferred into large per-step costs (e.g., of Rosenbrock methods) in time integration.
Nonlinear term
Properties of the Fourier transform: (c) Suppose $f$ is continuous and piecewise smooth (piecewise smooth). Then $$ \mathcal{F} \left[ f^{\prime} \right] (\xi) = i \xi \mathcal{F} f (\xi) $$
For the $k$-th Fourier mode, the nonlinear contribution $\left( u^{2} / 2 \right)_{x}$ must be reflected on the right-hand side as follows. $$ \mathcal{F} \left[ {\frac{ \partial }{ \partial x }} \left( {\frac{ 1 }{ 2 }} u^{2} \right) \right] (k) = i k \mathcal{F} \left[ {\frac{ 1 }{ 2 }} u^{2} \right] (k) $$ Therefore, the differential equation including the nonlinear term for the $k$-th mode is written as follows. Although $u$ appears to remain on the right-hand side, in practice the computation is performed in the frequency domain via the Fourier transform $\mathcal{F}$. $$ \dot{c}_{k} (t) = - ik \mathcal{F} \left[ {\frac{ 1 }{ 2 }} u^{2} \right] (k) + \left[ k^{2} - k^{4} \right] c_{k} (t) $$
Code
The video above shows an implementation in Julia that solves the Kuramoto–Sivashinsky equation for $L = 22$ with $64$ grid points and time step $h = 0.01$ over a duration of $t \in [0, 200]$ using a spectral method.
using DifferentialEquations, FFTW, SparseArrays, LinearAlgebra, Random, Plots, Arpack
function ks_spectral_rhs!(du_hat, u_hat, p, t)
k, Lk = p.k, p.Lk
u = real.(ifft(u_hat))
u22 = (u .^ 2) / 2
nonlinear_term = -im .* (k .* fft(u22))
linear_term = Lk .* u_hat
du_hat .= linear_term .+ nonlinear_term
return nothing
end
Looking at the right-hand side of the ODE system, one can easily verify that it matches the formulae discussed above. Since we must update the information for $u$ in the frequency domain, a Fourier transform is applied at every time step.
L = 22.0; Q = 64; dt = 0.01; T = 200
tspan = 0:dt:T
x = (0:Q-1) .* (L/Q)
k = vcat(collect(0:Q÷2-1), 0, collect(-Q÷2+1:-1)) .* (2pi / L)
Lk = k .^ 2 .- k .^ 4
Random.seed!(0)
prob = ODEProblem(ks_spectral_rhs!, fft(randn(Q)), extrema(tspan), (; k, Lk))
@time sol = solve(prob, ROS2(autodiff=false), saveat=collect(tspan), reltol=1e-4, abstol=1e-6);
u = stack(real.(ifft.(sol.u)))
In fact, a 2-stage Rosenbrock method was used to integrate the stiff ODEs. The inverse Fourier transform is used only once after all time-stepping is finished.
_p1 = heatmap(u, xlabel = "t", ylabel = "grid", xticks = [])
@time anime = @animate for tk in 1:100:size(u, 2)
p1 = vline(_p1, [tk], label = :none, color = :white)
p2 = plot(u[:, tk], label = :none, ylims = [-3, 3], color = :black, xlabel = "grid", ylabel = "u")
plot(p1, p2, layout = (2, 1), size = [600, 600])
end
mp4(anime, "ks_pde.mp4")
This is just visualization code to create the animation.