Chapter 7
Fourier spectral methods
Context: We now develop the ideas for solving partial differential equations outlined in the previous chapters. In this chapter, we specifically describe solution strategies using the Fourier basis. This leads to a solution method that belongs to the class of spectral methods.
7.1 Differential operators
To solve differential equations, we now use exactly the same methods that we developed in the previous chapter: Minimizing the residual using the Galerkin method. Our residual now has the general form \begin {equation} R(x,y,z,\ldots ; a_0, a_1, \ldots , a_N) = \mathcal {L} u_N(x,y,z,\ldots ) - f(x,y,z,\ldots ), \end {equation} where the unknown function \(u_N\) is represented here as a series expansion into a certain basis \(\varphi _n(x,y,z)\). The Galerkin method requires \begin {equation} (\varphi _n, R) = 0 \end {equation} for each \(n\).
We now discuss the Fourier basis for periodic functions on \(x\in [0,L]\) in one dimension, \begin {equation} \varphi _n(x) = \exp (i q_n x) \label {eq:fourierbasis2} \end {equation} with \(q_n = 2\pi n/L\). The operator \(\mathcal {L}\) can contain any differential operations that act on the basis functions, for example \begin {align} \frac {d}{dx} \varphi _n(x) &= iq_n \varphi _n(x) \\ \frac {d^2}{dx^2} \varphi _n(x) &= -q_n^2 \varphi _n(x). \end {align}
The derivatives of the (Fourier) basis functions result in the same basis function and an algebraic factor. It is also said that the basis functions diagonalize the differential operator. (This will be different for the finite elements discussed in the next chapter).
This property is particularly useful because, at least for linear differential equations, the residual becomes a trivial series expansion again and we can easily determine the coefficients using the orthogonality of the basis.
7.2 Poisson equation in one dimension
We use the (one-dimensional) Poisson equation, \begin {equation} \nabla ^2 \Phi \equiv \frac {\dif ^2 \Phi }{\dif x^2} = - \frac {\rho }{\varepsilon }, \label {eq:poisson-1d} \end {equation} as a demonstrator. Here \(\rho \) is a charge density and \(\Phi \) is the electrostatic potential. The residual is therefore \begin {equation} R(x)=\frac {\dif ^2 \Phi }{\dif x^2} + \frac {\rho }{\varepsilon }, \label {eq:poisson-1d-res} \end {equation} and the solution of Eq. \eqref{eq:poisson-1d} is given by \(R(x)=0\).
Formally, we now write the potential as the series expansion \begin {equation} \Phi (x) \approx \Phi _N(x) = \sum _{n=-(N-1)/2}^{(N-1)/2} a_n \varphi _n(x), \end {equation} whereby we will not explicitly specify the summation limits in the following. We also expand the right-hand side of Eq. \eqref{eq:poisson-1d} into a series with the same basis functions, \begin {equation} \rho _N(x) = \sum _{n=-(N-1)/2}^{(N-1)/2} b_n \varphi _n(x). \label {eq:seriesphi} \end {equation} Substituting this into Eq. \eqref{eq:poisson-1d-res} we obtain \begin {equation} R_N(x) = - \sum _n a_n q_n^2 \varphi _n(x) + \frac {1}{\varepsilon } \sum _n b_n \varphi _n(x). \end {equation} We now multiply this from the left by the basis functions, \((\varphi _k, R_N)\) (Galerkin method) and, due to the orthogonality of the basis functions, we obtain the equations \begin {equation} (\varphi _k, R_N) = - L q_k^2 a_k + L b_k/\varepsilon . \end {equation} The factor \(L\) appears because the basis functions are not normalized. The condition \((\varphi _k, R_N)=0\) leads to \(a_k = b_k/(q_k^2 \varepsilon )\). The approximate solution of the Poisson equation is thus given by \begin {equation} \Phi _N(x) = \sum _n \frac {b_n}{q_n^2 \varepsilon } \varphi _n(x). \label {eq:discrpoissonfouriersol} \end {equation} This is the Fourier series of the solution.
7.3 Transition to the Fourier transform
The Fourier basis Eq. \eqref{eq:fourierbasis2} is periodic on a finite domain of length \(L\). If we let the length \(L\) go to infinity, we get a formulation for non-periodic functions. This leads directly to the Fourier transform.
We write the series expansion as \begin {equation} \Phi _N(x) = \sum _{n=-N}^N a_n \varphi _n(x) = \sum _{n=-N}^N a_n \exp \left ( i q_n x \right ) = \sum _{n=-N}^N \frac {\Delta q}{2\pi }\,\tilde {\Phi }(q_n) \exp \left ( i q_n x \right ) \label {eq:fouriertrafo1} \end {equation} with \(\Delta q = q_{n+1}-q_n = 2\pi /L\) and rescaled coefficients \(\tilde {\Phi }(q_n)=L a_n\). Here, only the factor \(1=L \Delta q/2\pi \) was inserted on the right-hand side of Eq. \eqref{eq:fouriertrafo1}. This now helps to form the limits \(L\to \infty \) and \(N\to \infty \). In this case, \(\Delta q \to dq\) and the sum becomes the integral. This yields \begin {equation} \Phi (x) = \int _{-\infty }^\infty \frac {\dif q}{2\pi }\,\tilde {\Phi }(q) \exp \left ( i q x \right ), \label {eq:fouriertrafo2} \end {equation} the inverse Fourier transform.
The (forward) transform is obtained via a similar argument. We now know that \begin {equation} \tilde {\Phi }(q_n) = L a_n = L \frac {(\varphi _n, \Phi _N)}{(\varphi _n, \varphi _n)} = (\varphi _n, \Phi _N) = \int _0^L \dif x \, \Phi _N(x) \exp \left ( -i q_n x \right ). \end {equation} In the limiting case \(L\to \infty \) and \(N\to \infty \) this becomes \begin {equation} \tilde {\Phi }(q) = \int _{-\infty }^\infty \dif x \, \Phi (x) \exp \left ( -i q x \right ), \label {eq:fouriertrafo3} \end {equation} of the Fourier transform. The Fourier transform is useful to obtain analytical solutions for partial differential equations on infinite domains.
Note: A tilde \(\tilde {f}(q)\) denotes the Fourier transform of a function \(f(x)\). The Fourier transform is a function of the wave vector \(q\). In contrast, from the Fourier series we obtain countable coefficients \(a_n\). The reason for this is the periodicity of the function.
7.4 Poisson equation in multiple dimensions
Similar to how we constructed an approximate solution for a differential equation using a series expansion, we can use the approach Eq. \eqref{eq:fouriertrafo2} to obtain analytical solutions. In this section, this is demonstrated using the Poisson equation in three dimensions.
In three dimensions, the Poisson equation is \begin {equation} \nabla ^2 \Phi \equiv \frac {\partial ^2 \Phi }{\partial x^2} + \frac {\partial ^2 \Phi }{\partial y^2} + \frac {\partial ^2 \Phi }{\partial z^2} = - \frac {\rho }{\varepsilon }. \label {eq:poisson-3d} \end {equation} In contrast to Eq. \eqref{eq:poisson-1d}, the partial derivative \(\partial \) now appears here because \(\Phi (x,y,z)\) depends on three variables (the Cartesian coordinates).
The generalization of the Fourier basis and thus also of the Fourier transform to three dimensions is trivial. A basis is obtained by multiplying basis functions in the Cartesian directions (\(x\), \(y\) and \(z\)). Usually, you now need three indices for the coefficients, which denote the basis in \(x\), \(y\) and \(z\) respectively. The result is a series expansion \begin {equation} \begin {split} \Phi _{NMO}(x,y,z) =& \sum _{n=-(N-1)/2}^{(N-1)/2} \sum _{m=-(M-1)/2}^{(M-1)/2} \sum _{o=-(O-1)/2}^{(O-1)/2} a_{nmo} \varphi _n(x) \varphi _m(y) \varphi _o(z) \\ \equiv & \sum _{n=-(N-1)/2}^{(N-1)/2} \sum _{m=-(M-1)/2}^{(M-1)/2} \sum _{o=-(O-1)/2}^{(O-1)/2} a_{nmo} \varphi _{nmo}(x,y,z) \end {split} \end {equation} with (possibly different) truncation orders \(N\), \(M\) and \(O\). The basis set is given here by the functions \(\varphi _{nmo}(x,y,z)=\varphi _n(x)\varphi _m(y)\varphi _o(z)\). Orthogonality of this basis set is trivially derived from the orthogonality of the one-dimensional basis functions \(\varphi _n(x)\). The generalization of the Fourier transform follows directly from this. The Fourier inverse transform is written as \begin {equation} \Phi (x,y,z) = \int _{-\infty }^\infty \frac {\dif ^3 q}{(2\pi )^3}\,\tilde {\Phi }(q_x, q_y, q_z) \exp \left ( i q_x x + i q_y y + i q_z z\right ), \label {eq:fouriertrafo3d} \end {equation} where the Fourier transform \(\tilde {\Phi }\) now naturally depends on three wave vectors \(q_x\), \(q_y\) and \(q_z\). The differential operator \(\dif ^3 q=\dif q_x \dif q_y \dif q_z\) is a shorthand notation for three-dimensional integration.
We can now insert Eq. \eqref{eq:fouriertrafo3d} into the PDGL Eq. \eqref{eq:poisson-3d} and obtain \begin {equation} R(\v {r}) = \int _{-\infty }^\infty \frac {\dif ^3 q}{(2\pi )^3}\, \left [ \left (-q_x^2 - q_y^2 - q_z^2\right ) \tilde {\Phi }(\v {q}) + \frac {\tilde {\rho }(\v {q})}{\varepsilon } \right ] \exp \left ( i \v {q}\cdot \v {r} \right ) = 0 \end {equation} with \(\v {r}=(x,y,z)\) and \(\v {q}=(q_x,q_y,q_z)\). This equation must be fulfilled for every \(x,y,z\) and therefore the argument of the integration must disappear, i.e. \begin {equation} -q^2 \tilde {\Phi }(\v {q}) + \frac {\tilde {\rho }(\v {q})}{\varepsilon } = 0. \label {eq:fourierpoisson2d} \end {equation}
Note: An alternative argument is obtained by writing the Fourier transform of \(R(x,y,z)\): \begin {equation} R(q_x', q_y', q_z') = \int \dif ^3 r \, R(\v {r}) \exp \left ( -i q_x x - i q_y y - i q_z z \right ). \end {equation} This contains terms of the form \begin {equation} \int _{-\infty }^\infty \dif x\, \exp \left (i (q_x - q_x') x\right ) = 2\pi \delta (q_x - q_x'), \end {equation} which are expressions of the orthogonality of the basis functions. Since the basis functions are now “parameterized” with a continuous \(q_x\) (instead of a discrete \(n\)), a Dirac \(\delta \) function is obtained instead of the Kronecker \(\delta \) in the orthogonality relation.
Equation \eqref{eq:fourierpoisson2d} can easily be solved analytically. This yields \begin {equation} \tilde {\Phi }(\v {q}) = \frac {\tilde {\rho }(\v {q})}{\varepsilon q^2} \label {eq:fourierpoissonsol} \end {equation} with \(q=|\v {q}|\). This is equivalent to solving Eq. \eqref{eq:discrpoissonfouriersol} for the Poisson equation on a periodic domain. The difficulty now lies in evaluating the back and forth transformation for a given \(\rho (x,y,z)\).
Example: As an example, we now consider the solution for a point charge \(Q\) at the origin, \begin {equation} \rho (x,y,z) = Q \delta (x) \delta (y) \delta (z). \end {equation} The Fourier transform of the charge density \(\rho \) is obtained from Eq. \eqref{eq:fouriertrafo3}, \begin {equation} \tilde {\rho }(q_x,q_y,q_z) = Q. \end {equation} I.e. the Fourier transform of the electrostatic potential is given by (see Eq. \eqref{eq:fourierpoissonsol}) \begin {equation} \tilde {\Phi }(\v {q}) = \frac {Q}{\varepsilon q^2}, \end {equation} and thus the representation in real space is \begin {equation} \begin {split} \Phi (\v {r}) =& \int _{-\infty }^\infty \frac {\dif ^3 q}{(2\pi )^3}\, \frac {Q}{\varepsilon q^2} \exp \left ( i \v {q}\cdot \v {r} \right ) \\ =& \frac {Q}{(2\pi )^3 \varepsilon } \int _0^\infty \dif q \int _0^{2\pi } \dif \phi \int _{-1}^1 \dif (\cos \theta ) \, \exp \left ( i q r \cos \theta \right ) \end {split} \end {equation} where \(\dif ^3 q = q^2 \dif q \dif \phi \dif (\cos \theta )\) with azimuth angle \(\phi \) and elevation angle \(\theta \), was used (see also Fig. 7.1). We require here (without limiting the generality) that \(\v {r}\) points in the direction of the zenith.
One obtains \begin {equation} \begin {split} \Phi (\v {r}) =& \frac {Q}{(2\pi )^2 \varepsilon } \int _0^\infty \dif q \int _{-1}^1 \dif (\cos \theta ) \, \exp \left ( i q r \cos \theta \right ) \\ =& \frac {Q}{(2\pi )^2 \varepsilon } \int _0^\infty \dif q\, \frac {\exp ( i q r ) - \exp (-iqr)}{iqr} \\ =& \frac {Q}{(2\pi )^2 \varepsilon } \int _{-\infty }^\infty \dif q\, \frac {\sin q r}{qr} \\ =& \frac {Q}{4\pi \varepsilon r}, \end {split} \end {equation} where \(\int \dif x\,\sin x/x=\pi \) was used. This is the known solution for the electrostatic potential of a point charge. It is also called the fundamental solution or Green’s function of the (three-dimensional) Poisson equation.
