Spectral Methods

Examples of spectral methods

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

\[\left(\cfrac{\partial }{\partial x^{2}}+\cfrac{\partial }{\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:

\[\begin{split}\begin{array}{c} f=\sum a_{j,k}e^{\mathbf{i}(jx+ky)}\\ g=\sum b_{j,k}e^{\mathbf{i}(jx+ky)} \end{array}\end{split}\]

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

\[\begin{split}\begin{array}{l} \cfrac{\partial f}{\partial x}=\sum j\mathbf{i}a_{j,k}e^{\mathbf{i}(jx+ky)}\\ \cfrac{\partial^{2} f}{\partial x^{2}}=-\sum j^{2}a_{j,k}e^{\mathbf{i}(jx+ky)}\\ \cfrac{\partial f}{\partial y}=\sum k\mathbf{i}a_{j,k}e^{\mathbf{i}(jx+ky)}\\ \cfrac{\partial ^{2}f}{\partial y^{2}}=-\sum k^{2}a_{j,k}e^{\mathbf{i}(jx+ky)}\\ \end{array}\end{split}\]

\[\begin{split}\begin{array}{l} -\sum (j^{2}+k^{2})a_{j,k}e^{\mathbf{i}(jx+ky)}=\sum b_{j,k}e^{\mathbf{i}(jx+ky)}\\ -(j^{2}+k^{2})a_{j,k}=b_{j,k} \end{array}\end{split}\]

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

\[\begin{array}{l} a_{j,k}=-\cfrac{b_{j,k}}{j^{2}+k^{2}} \end{array}\]

Example 1 of spectral methods

Consider the following PDE with constant \(a\)

\[\cfrac{\partial u}{\partial t}=a\cfrac{\partial u}{\partial x}\]

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^{\mathbf{i}2\pi n\cfrac{x}{L}}=\text{cos}(2\pi n\cfrac{x}{L})+\mathbf{i}\text{sin}(2\pi n\cfrac{x}{L})\), where \(\mathbf{i}\) is the imaginary number and \(L\) is the domian 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(x,t)=\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t)e^{\mathbf{i}{2\pi}\cfrac{nx}{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=\cfrac{L}{N}\), and \(x_{j}=j\Delta x = j\cfrac{L}{N}\). Let \(u_{j}=u(x_{j})\),

Then the inverse discrete Fourier transform(IDFT) of \(\hat{u}_{n}\) evaluated at grid points \(x_{j}\) and denoted \(F^{-1}\), is given by

\[\begin{split}\begin{array}{l} \displaystyle u(x_{j},t)=\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t)e^{\mathbf{i}{2\pi}\cfrac{nx_{j}}{L}}\\ \displaystyle u(x_{j},t)=u_{j}(t)=\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t)e^{\mathbf{i}{2\pi}\cfrac{nj\Delta x}{L}}\\ \displaystyle u(x_{j},t)=u_{j}(t)=\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t)e^{\mathbf{i}{2\pi}\cfrac{nj}{L}\cfrac{L}{N}}\\ \displaystyle u(x_{j},t)=u_{j}(t)=\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t)e^{\mathbf{i}{2\pi}\cfrac{nj}{N}}\\ \end{array}\end{split}\]

\[\begin{split}u(x_{j},t)=u_{j}(t)=F^{-1}(\hat{u}_{n}(t))=\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t)e^{\mathbf{i}{2\pi}\cfrac{nj}{N}},j=0,\cdots,n-1.\\\end{split}\]

or

\[\begin{split}u_{j}=F^{-1}(\hat{u}_{n})=\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}e^{\mathbf{i}{2\pi}\cfrac{nj}{N}},j=0,\cdots,n-1.\\\end{split}\]

The corresponding discrete Fourier transform(DFT) is given by

\[\hat{u}_{n}=F(u_{j})=\cfrac{1}{N}\sum_{j=0}^{N-1}u_{j}e^{-\mathbf{i}{2\pi}\cfrac{nj}{N}},n=-\cfrac{N}{2}+1,\cdots,\cfrac{N}{2}\]

Consider the PDE

\[\begin{split}\begin{array}{l} \displaystyle u(x,t)=\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t)e^{\mathbf{i}{2\pi}\cfrac{nx}{L}}\\ \displaystyle \cfrac{\partial u(x,t)}{\partial t}=\sum_{n=-N/2+1}^{N/2}\cfrac{\partial\hat{u}_{n}(t)}{\partial t} e^{\mathbf{i}{2\pi}\cfrac{nx}{L}}\\ \displaystyle \cfrac{\partial u(x,t)}{\partial x}={\mathbf{i}{2\pi}\cfrac{n}{L}}\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t) e^{\mathbf{i}{2\pi}\cfrac{nx}{L}}\\ \end{array}\end{split}\]

then

\[\displaystyle \sum_{n=-N/2+1}^{N/2}\cfrac{\partial\hat{u}_{n}(t)}{\partial t} e^{\mathbf{i}{2\pi}\cfrac{nx}{L}} =a{\mathbf{i}{2\pi}\cfrac{n}{L}}\sum_{n=-N/2+1}^{N/2}\hat{u}_{n}(t) e^{\mathbf{i}{2\pi}\cfrac{nx}{L}}\]

\[\cfrac{\partial\hat{u}_{n}(t)}{\partial t}=a{\mathbf{i}{2\pi}\cfrac{n}{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 ten be solved for \(\hat{u}_{n}(t)\). Then \({u}_{j}(t)\) can be computed form \(\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^{\left(a\mathbf{i}{2\pi}\cfrac{n}{L}\right)t}\]

Example 2 of spectral methods

Consider the following PDE:

\[\cfrac{\partial ^{2}u}{\partial x^{2}} +\cfrac{\partial ^{2}u}{\partial y^{2}}=f(x,y)\]

Represent the solution \(u(x,y)\) as

\[u(x,y)=\sum_{m,n}\hat{u}_{m,n}\phi_{m}(x)\phi_{n}(y)\]

where \(\phi_{m}\) and \(\phi_{n}\) are some basis functions and \(\hat{u}_{m,n}\) are corresponding weighting factors.

In a spectral method we take

\[\phi_{m}(x)=e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}=\text{cos}(2\pi m\cfrac{x}{L_{x}})+\mathbf{i}\text{sin}(2\pi m\cfrac{x}{L_{x}})\]

\[\phi_{n}(y)=e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}=\text{cos}(2\pi n\cfrac{y}{L_{y}})+\mathbf{i}\text{sin}(2\pi n\cfrac{y}{L_{y}})\]

where \(\mathbf{i}\) is the imaginary number. \(L_{x}\) and \(L_{y}\) are \(x\) and \(y\) direction Length. then

\[u(x,y)=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\]

\[u(x,y)=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi(m\cfrac{x}{L_{x}}+n\cfrac{y}{L_{y}}) }\]

Let

\[\begin{split}\begin{array}{c} \Delta x = \cfrac{L_{x}}{M}， x_{i}=i\Delta x= i\cfrac{L_{x}}{M} \\ \Delta y = \cfrac{L_{y}}{N}， y_{j}=j\Delta y= j\cfrac{L_{y}}{N} \\ \end{array}\end{split}\]

Then the inverse discrete Fourier transform(IDFT):

\[\begin{split}\begin{array}{l} \displaystyle u(x,y)=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi(m\cfrac{x}{L_{x}}+n\cfrac{y}{L_{y}}) }\\ \displaystyle u(x_{i},y_{j})=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi(m\cfrac{x_{i}}{L_{x}}+n\cfrac{y_{j}}{L_{y}}) }\\ \displaystyle u(x_{i},y_{j})=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi(m\cfrac{i\Delta x}{L_{x}}+n\cfrac{j\Delta y}{L_{y}}) }\\ \displaystyle u(x_{i},y_{j})=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi(m\cfrac{iL_{x}/M}{L_{x}}+n\cfrac{jL_{y}/N}{L_{y}}) }\\ \displaystyle u(x_{i},y_{j})=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi(\cfrac{im}{M}+\cfrac{jn}{N}) }\\ \end{array}\end{split}\]

or

\[\begin{split}u_{i,j}=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi(\cfrac{im}{M}+\cfrac{jn}{N}) }\\\end{split}\]

and

\[i=0,1,\cdots,M-1; \quad j=0,1,\cdots ,N-1\]

The corresponding discrete Fourier transform(DFT) is given by

\[\begin{split}\hat{u}_{m,n}=\cfrac{1}{MN} \sum_{i=0}^{M-1}\sum_{j=0}^{N-1}u_{i,j}e^{-\mathbf{i}2\pi(\cfrac{im}{M}+\cfrac{jn}{N}) }\\\end{split}\]

where

\[\begin{split}\begin{array}{l} m=-\cfrac{M}{2}+1,\cdots,\cfrac{M}{2}\\ n=-\cfrac{N}{2}+1,\cdots,\cfrac{N}{2}\\ \end{array}\end{split}\]

Consider the PDE:

\[\begin{split}\begin{array}{l} \displaystyle u(x,y)=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\\ \displaystyle \cfrac{\partial u(x,y)}{\partial x} ={\mathbf{i}2\pi m\cfrac{1}{L_{x}}}\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\\ \displaystyle \cfrac{\partial ^{2}u(x,y)}{\partial x^{2}} =\left\{\mathbf{i}2\pi m\cfrac{1}{L_{x}}\right\}^{2}\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\\ \displaystyle \cfrac{\partial u(x,y)}{\partial y} ={\mathbf{i}2\pi n\cfrac{1}{L_{y}}}\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\\ \displaystyle \cfrac{\partial ^{2}u(x,y)}{\partial y^{2}} =\left\{\mathbf{i}2\pi n\cfrac{1}{L_{y}}\right\}^{2}\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{u}_{m,n}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\\ \end{array}\end{split}\]

and

\[\begin{split}f(x,y)=\sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}\hat{f}_{m,n}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\\\end{split}\]

then

\[\begin{split}\begin{array}{l} \cfrac{\partial ^{2}u}{\partial x^{2}} +\cfrac{\partial ^{2}u}{\partial y^{2}}=f(x,y)\\ \displaystyle \Rightarrow \sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\hat{u}_{m,n}\left\{\left(\mathbf{i}2\pi m\cfrac{1}{L_{x}}\right )^2+\left(\mathbf{i}2\pi n\cfrac{1}{L_{y}}\right )^2\right\}\\ \displaystyle = \sum_{m=-M/2+1}^{M/2}\sum_{n=-N/2+1}^{N/2}e^{\mathbf{i}2\pi m\cfrac{x}{L_{x}}}e^{\mathbf{i}2\pi n\cfrac{y}{L_{y}}}\hat{f}_{m,n} \end{array}\end{split}\]

\[\hat{u}_{m,n}\left\{\left(\mathbf{i}2\pi m\cfrac{1}{L_{x}}\right )^2+\left(\mathbf{i}2\pi n\cfrac{1}{L_{y}}\right )^2\right\}=\hat{f}_{m,n}\]

\[\begin{split}\begin{array}{l} -\hat{u}_{m,n}\left\{\left(\cfrac{2\pi m}{L_{x}}\right )^2+\left(\cfrac{2\pi n}{L_{y}}\right )^2\right\}=\hat{f}_{m,n}\\ \hat{u}_{m,n} =-\cfrac{\hat{f}_{m,n}}{\left(\cfrac{2\pi m}{L_{x}}\right )^2+\left(\cfrac{2\pi n}{L_{y}}\right )^2}\\ \end{array}\end{split}\]