Spectral method

Here we presume an understanding of basic multivariate calculus and Fourier series. If ${\displaystyle g(x,y)}$ is a known, complex-valued function of two real variables, and g is periodic in x and y (that is, ${\displaystyle g(x,y)=g(x+2\pi ,y)=g(x,y+2\pi )}$) then we are interested in finding a function f(x,y) so that

${\displaystyle \left({\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial y^{2}}}\right)f(x,y)=g(x,y)\quad {\text{for all }}x,y}$

where the expression on the left denotes the second partial derivatives of f in x and y, respectively. This is the Poisson equation, and can be physically interpreted as some sort of heat conduction problem, or a problem in potential theory, among other possibilities.

If we write f and g in Fourier series:

${\displaystyle f=:\sum a_{j,k}e^{i(jx+ky)}}$

${\displaystyle g=:\sum b_{j,k}e^{i(jx+ky)}}$

and substitute into the differential equation, we obtain this equation:

${\displaystyle \sum -a_{j,k}(j^{2}+k^{2})e^{i(jx+ky)}=\sum b_{j,k}e^{i(jx+ky)}}$

We have exchanged partial differentiation with an infinite sum, which is legitimate if we assume for instance that f has a continuous second derivative. By the uniqueness theorem for Fourier expansions, we must then equate the Fourier coefficients term by term, giving

(*) ${\displaystyle a_{j,k}=-{\frac {b_{j,k}}{j^{2}+k^{2}}}}$

which is an explicit formula for the Fourier coefficients a_j,k.

With periodic boundary conditions, the Poisson equation possesses a solution only if b_0,0 = 0. Therefore, we can freely choose a_0,0 which will be equal to the mean of the resolution. This corresponds to choosing the integration constant.

To turn this into an algorithm, only finitely many frequencies are solved for. This introduces an error which can be shown to be proportional to ${\displaystyle h^{n}}$, where ${\displaystyle h:=1/n}$ and ${\displaystyle n}$ is the highest frequency treated.

1. Compute the Fourier transform (b_j,k) of g.
2. Compute the Fourier transform (a_j,k) of f via the formula (*).
3. Compute f by taking an inverse Fourier transform of (a_j,k).

Since we're only interested in a finite window of frequencies (of size n, say) this can be done using a fast Fourier transform algorithm. Therefore, globally the algorithm runs in time O(n log n).

We wish to solve the forced, transient, nonlinear Burgers' equation using a spectral approach.

Given ${\displaystyle u(x,0)}$ on the periodic domain ${\displaystyle x\in \left[0,2\pi \right)}$, find ${\displaystyle u\in {\mathcal {U}}}$ such that

${\displaystyle \partial _{t}u+u\partial _{x}u=\rho \partial _{xx}u+f\quad \forall x\in \left[0,2\pi \right),\forall t>0}$

where ρ is the viscosity coefficient. In weak conservative form this becomes

${\displaystyle \left\langle \partial _{t}u,v\right\rangle =\left\langle \partial _{x}\left(-{\frac {1}{2}}u^{2}+\rho \partial _{x}u\right),v\right\rangle +\left\langle f,v\right\rangle \quad \forall v\in {\mathcal {V}},\forall t>0}$

where ${\displaystyle \langle f,g\rangle$ :=\int _{0}^{2\pi }f(x){\overline {g(x)}}\,dx} following inner product notation. Integrating by parts and using periodicity grants

${\displaystyle \langle \partial _{t}u,v\rangle =\left\langle {\frac {1}{2}}u^{2}-\rho \partial _{x}u,\partial _{x}v\right\rangle +\left\langle f,v\right\rangle \quad \forall v\in {\mathcal {V}},\forall t>0.}$

To apply the Fourier-Galerkin method, choose both

${\displaystyle {\mathcal {U}}^{N}:=\left\{u:u(x,t)=\sum _{k=-N/2}^{N/2-1}{\hat {u}}_{k}(t)e^{ikx}\right\}}$

and

${\displaystyle {\mathcal {V}}^{N}:=\operatorname {span} \left\{e^{ikx}:k\in -N/2,\dots ,N/2-1\right\}}$

where ${\displaystyle {\hat {u}}_{k}(t):={\frac {1}{2\pi }}\langle u(x,t),e^{ikx}\rangle }$. This reduces the problem to finding ${\displaystyle u\in {\mathcal {U}}^{N}}$ such that

${\displaystyle \langle \partial _{t}u,e^{ikx}\rangle =\left\langle {\frac {1}{2}}u^{2}-\rho \partial _{x}u,\partial _{x}e^{ikx}\right\rangle +\left\langle f,e^{ikx}\right\rangle \quad \forall k\in \left\{-N/2,\dots ,N/2-1\right\},\forall t>0.}$

Using the orthogonality relation ${\displaystyle \langle e^{ilx},e^{ikx}\rangle =2\pi \delta _{lk}}$ where ${\displaystyle \delta _{lk}}$ is the Kronecker delta, we simplify the above three terms for each ${\displaystyle k}$ to see

${\displaystyle {\begin{aligned}\left\langle \partial _{t}u,e^{ikx}\right\rangle &=\left\langle \partial _{t}\sum _{l}{\hat {u}}_{l}e^{ilx},e^{ikx}\right\rangle =\left\langle \sum _{l}\partial _{t}{\hat {u}}_{l}e^{ilx},e^{ikx}\right\rangle =2\pi \partial _{t}{\hat {u}}_{k},\\\left\langle f,e^{ikx}\right\rangle &=\left\langle \sum _{l}{\hat {f}}_{l}e^{ilx},e^{ikx}\right\rangle =2\pi {\hat {f}}_{k},{\text{ and}}\\\left\langle {\frac {1}{2}}u^{2}-\rho \partial _{x}u,\partial _{x}e^{ikx}\right\rangle &=\left\langle {\frac {1}{2}}\left(\sum _{p}{\hat {u}}_{p}e^{ipx}\right)\left(\sum _{q}{\hat {u}}_{q}e^{iqx}\right)-\rho \partial _{x}\sum _{l}{\hat {u}}_{l}e^{ilx},\partial _{x}e^{ikx}\right\rangle \\&=\left\langle {\frac {1}{2}}\sum _{p}\sum _{q}{\hat {u}}_{p}{\hat {u}}_{q}e^{i\left(p+q\right)x},ike^{ikx}\right\rangle -\left\langle \rho i\sum _{l}l{\hat {u}}_{l}e^{ilx},ike^{ikx}\right\rangle \\&=-{\frac {ik}{2}}\left\langle \sum _{p}\sum _{q}{\hat {u}}_{p}{\hat {u}}_{q}e^{i\left(p+q\right)x},e^{ikx}\right\rangle -\rho k\left\langle \sum _{l}l{\hat {u}}_{l}e^{ilx},e^{ikx}\right\rangle \\&=-i\pi k\sum _{p+q=k}{\hat {u}}_{p}{\hat {u}}_{q}-2\pi \rho {}k^{2}{\hat {u}}_{k}.\end{aligned}}}$

Assemble the three terms for each ${\displaystyle k}$ to obtain

${\displaystyle 2\pi \partial _{t}{\hat {u}}_{k}=-i\pi k\sum _{p+q=k}{\hat {u}}_{p}{\hat {u}}_{q}-2\pi \rho {}k^{2}{\hat {u}}_{k}+2\pi {\hat {f}}_{k}\quad k\in \left\{-N/2,\dots ,N/2-1\right\},\forall t>0.}$

Dividing through by ${\displaystyle 2\pi }$, we finally arrive at

${\displaystyle \partial _{t}{\hat {u}}_{k}=-{\frac {ik}{2}}\sum _{p+q=k}{\hat {u}}_{p}{\hat {u}}_{q}-\rho {}k^{2}{\hat {u}}_{k}+{\hat {f}}_{k}\quad k\in \left\{-N/2,\dots ,N/2-1\right\},\forall t>0.}$

With Fourier transformed initial conditions ${\displaystyle {\hat {u}}_{k}(0)}$ and forcing ${\displaystyle {\hat {f}}_{k}(t)}$, this coupled system of ordinary differential equations may be integrated in time (using, e.g., a Runge Kutta technique) to find a solution. The nonlinear term is a convolution, and there are several transform-based techniques for evaluating it efficiently. See the references by Boyd and Canuto et al. for more details.