Contents

7. Pseudo-spectral method and aliasing errors

As we have seen, the evaluation of the nonlinear term in equation (6.2) pertaining to the Galerkin discretisation of the Burgers equation is computationally of order \(\mathcal{O}(K^2)\). It is therefore desirable to find more efficient ways to compute this sum. The most imporant way to do this is via Fourier transforms leading to an order \(\mathcal{O}(K \log K)\) for the same operation which is significantly less than \(\mathcal{O}(K^2)\) for large \(K\).

To illustrate such an algorithm, consider the product \(w_j=w(x_j)\), \(j=1,\ldots,N\) of two grid functions \(u_j=u(x_j)\) and \(v_j=v(x_j)\) in physical space

(7.1)\[ \begin{equation} w_j = u_j \cdot v_j \ . \end{equation} \]

We proceed by transforming to spectral space with the help of

\[ \begin{equation} u_j = \sum_{|k|\leq K} \hat{u}_k e^{\mathrm{i} k x_j} \ , \ \ \ \hat{u}_k = \frac{1}{N}\sum_{j=1}^N u_j e^{-\mathrm{i}kx_j} \end{equation} \]

and the discrete form of the orthogonality relation

\[ \begin{equation} \frac{1}{N}\sum_{j=1}^N e^{\mathrm{i}k x_j} e^{\mathrm{i}m x_j} = \delta_{k,-m+n\cdot N} \ . \end{equation} \]

In the last equation, \(n\in \mathbb{Z}\) leads to an arbitrary multiple of \(N\) which is closely related to aliasing errors (see further down).
Finally one obtains

(7.2)\[\begin{split} \begin{equation} \hat{w}_l = \frac{1}{N} \sum_{|k|\leq K} \sum_{|m|\leq K} \hat{u}_k \hat{v}_m \sum_{j=1}^N e^{\mathrm{i} k x_j} e^{\mathrm{i}(m-l)x_j} = \sum_{\substack{k=-K\\ l=k+m+n\cdot N\\ |m|\leq K}}^K \hat{u}_k \hat{v}_{m} \ , \end{equation} \end{split}\]

i.e. the well-known result that a multiplication in physical space corresponds to a convolution in spectral space. To put it into other words, the fairly expensive evaluation of a convolution in spectral space, equation (7.2), is equivalent to the direct evaluation of a pointwise multiplication in physical space, equation (7.1).

Comparing our latest result, equation (7.2), to equation (6.2) we see that an efficient evaluation of a nonlinear product of two variables given in spectral space, \(\hat{u}_k\) and \(\hat{v}_m=i m \hat{u}_m\) may be done via the following steps

  1. Transform \(\hat{u}_k\) and \(\hat{v}_m\) to physical space using FFT: \(u_j=\mathcal{F}^{-1}(\hat{u}_k)\), \(v_j=\mathcal{F}^{-1}(\hat{v}_m)\)

  2. Multiplication in physical space: \(w_j=u_j\cdot v_j\)

  3. Transform back to spectral space: \(\hat{w}_l=\mathcal{F}(w_j)\). This final results then reads:

(7.3)\[\begin{split} \begin{equation} \hat{w}_l = \sum _{\substack{k=-K\\ l=k+m + n\cdot N\\ |m|\leq K}}^K \mathrm{i} m \hat{u}_k \hat{u}_m = \sum _{\substack{k=-K\\ l=k+m + n\cdot N\\ |m|\leq K}}^K \hat{u}_k \hat{v}_m \ . \end{equation} \end{split}\]

Such an evaluation of the spectral convolution in physical space is usually termed pseudo-spectral evaluation of the nonlinear terms, since the calculations are not done exclusively in spectral space.

7.1. Aliasing errors

According to the sampling theorem, the highest wave number that can be represented as a grid function \(f_j\) with \(j=1,\ldots N=2K\) is \(K\). Higher wave numbers \(k_h>K\) are mapped to this representable region \(|k|\leq K\) via

(7.4)\[ \begin{equation} k=k_h-nN \ , \ \ n\in\mathbb{Z} \ , \ \ |k|\leq K \ . \end{equation} \]

An example of such a case is given in Figure Fig. 7.1 illustrating the misrepresentation of wave numbers when mapped on a grid function with not enough grid points.

../../_images/alias.png

Fig. 7.1 Illustration of aliasing errors, i.e. the mapping of higher wave numbers to lower ones. Here, \(N=8\) (number of points in physical space) and \(K=N/2=4\) (highest representable wavenumber), and the wave numbers of the two waves are \(k_{0}=9\) (black) and \(k_1=k_0-N=1\) (red).

7.2. Pseudo-spectral method

Such aliasing errors also appear when evaluating a convolution sum via Fourier transforms. To understand this, compare the nonlinear term in equations (6.2) and (7.3): The difference is the appearance of the term \(\ldots +n\cdot N\) in the expression computed via the FFT. For the computed result to be correct, measures have to be taken to avoid the inclusion of these additional parts (aliasing errors) in the sum.

Another way to look at aliasing errors is by realising that the non-linear evaluation in \(\hat{w}_l\) in equation (7.2) involves factors of \(\hat{v}_m\) with \(|m|=|l-k|\leq 2K\), i.e. wave numbers up to twice as high as representable on the grid. These wave numbers are thus mapped back to the represented wave-number space according to the above expression (7.4). Thus the sum will contain errors due to these spurious contributions, exactly the aliasing errors.

The additional conditions on the convolutional sum when using a pseudo-spectral evaluation of the nonlinear terms are however not directly possible to implement using the above algorithm for computing \(\hat{w}_l\) given in equation (7.2). There are however some possibilities to remove (or at least reduce) aliasing errors even in the pseudo-spectral case. The main options are, commonly termed dealiasing:

  • remove via phase shifting: In this method, the nonlinear terms are calculated twice with phase shifts, which can then be used to remove the spurious contributions.

  • apply a filter in spectral space before and after the transform to physical space which removes the highest \(1/3\) of the wavenumber content. This means that the actual resolution in physical space is only \(2N/3\). Therefore, this method is called the \(2/3\)-rule.

  • Another popular variant is the so-called 3/2-rule: The original grid in physical space is refined by a factor \(3/2\) to \(M=3/2N\) in every direction, and the nonlinear multiplications are then performed on this finer grid. Afterwards, the resulting product is transformed back to spectral space, cutting away the wave numbers with \(|k|>K\). Using the \(3/2\)-rule, the aliasing errors are completely removed for a quadratic nonlinearity as in the case of the Burgers equation, and thus the pseudo-spectral method becomes equivalent to a true Galerkin method in spectral space. As discussed in the literature, it turns out that explicitly removing aliasing errors yields better results than just incrasing the resolution (i.e. the filtering step is necessary).

For Chebyshev and Legendre methods based on Gaussian quadrature similar issues arise, which can be treated using overintegration. In this technique, the mesh is extended again by a factor of \(3/2\) prior to quadrature, to make the Gauss formula exact for the required order.

Note that aliasing errors are more severe for marginally resolved cases, as then the energy content in the highest resolved modes is typically larger. Alternative methods, such as finite-difference methods, are equally affected by aliasing errors, however their suppression is less obvious. Different formulations of the nonlinear term such as the skew-symmetric form, albeit analytically identical, have been shown to reduce aliasing errors, and are thus commonlly employed in particular for large-eddy simulations (LES).

6. Nonlinear problems 8. NB Numerical accuracy