Spectral Methods for Hyperbolic Problems

Introduction

What are Spectral and Pseudospectral Methods?

Instead of discretising the domain, we express functions in terms of a truncated expansion over a particular basis. The choice of basis leads to advantages. For example, the Fourier basis, $\phi_n(x) = \exp(in\xi)$ are eigenvectors of the differential operator. Other commonly chosen bases are Chebyshev polynmials and Legendre polynomials.

Domain discretisation approximates a function up to first, or at most polynomial. In comparison, the spectral approximation of analytic functions as well as their derivatives converge exponentially.

Since multiplication in the spectral space is more expensive than the transform itself, pseudospectral methods perform multiplication in the original space.

What are Hyperbolic Problems?

Euler equations, shallow water equations, wave equations are some examples of hyperbolic systems of partial differential equations. Unlike parabolic/elliptic PDEs, which have inherent regularity, hyperbolic problems exhibit a finite speed of propogation and wave-like solutions. These solutions tend to develop discontinuities even when the initial condition is smooth. In the context of hyperbolic problems, discontinuities are often called shocks.

This report focuses on two simple, one-dimensional, hyperbolic problems to test the methods on:

Both are defined on the domain $[0,2\pi]$ and assumed to have periodic boundary conditions.

The solution to the linear advection equation is a simple translation of the initial condition in the positive $x$ direction. Thus the solution at $t=1$ matches perfectly with the initial condition.

Burgers’ equation, on the other hand, has a complicated solution. One can think of particles at $x$ travelling with the speed of $u(x)$ toward the positive $x$ direction.

Implementation

The code implementing the techniques and generating the plots is available at github.com/aadi-bh/spectral4hyp. It uses the Fourier pseudospectral method with RK4 time stepping and dealiasing to solve both these PDES.

Exact solutions are computed using a Lax-Friedrich scheme on a fine mesh.

A Note About Plots

Unless specified, every plot shows the initial condition, computed/processed solution, and exact solution on the left, with the pointwise error on the right.

Aliasing Errors

When sampling a function at only finitely many points, our coefficients are bound to suffer from aliasing errors.

Aliasing

These aliasing errors and their nonlinear mixing contributes to error.

Gibbs Phenomenon

Suppose $u_N$ is the truncated Fourier expansion of $u$. Then \[ ||u - u_N ||_{L^2[0, 2\pi]} = 2\pi\sum_{{|n|> N}} |\hat u_n|^2 . \] Thus it’s clear that the rate of convergence in $L^2$ is dependent on the decay of the Fourier coefficients.

If $u$ is only piecewise smooth, then \[ u_N\left(x_0 + \frac{2z}{2N+1}\right) \sim \frac12[u(x_0^+) + u(x_0^-)]) + \frac{1}{\pi}[u(x_0^+) - u(x_0^-)]\textrm{Si}(z) \] tells us that there is only linear pointwise convergence at points away from the discontinuity.
$\textrm{Si}(z)\to \frac\pi 2$ as $z\to\infty\implies$ a complete loss of convergence as $\bm x$ approaches the discontinuity.

Smooth case

Continuous case

Discontinuous case

Implications of the Gibbs Phenomenon

The lack of regularity/dissipation in hyperbolic problems makes tbe formation of discontinuities inevitable. Such discontinuities are not captured well by spectral/pseudospectral methods. This is a source of concern, and hence these methods have traditionally not been applied to such problems.

“With the situation being even more complex and the details partially unresolved for the discrete expansions, the Gibbs phenomenon and the loss of fast global convergence is often perceived as an argument against the use of spectral expansions for the representation of piecewise smooth functions and, ultimately, against the use of spectral methods for the solution of problems with discontinuous solutions.”

Filters

Oscillations are an indication of high-order modes. The Gibbs oscillations can be seen to be the cause of the slow decay of the Fourier coefficients, leading to slow convergence. Filters essentially force better convergence by reducing the amplitude of high-order modes.

1. $\sigma(\eta) = 0$ for $|\eta| \ge 1$.
2. $\sigma(0) = 1$.
3. $\sigma^{(m)}(0) = \sigma^{(m)}(1) = 0 \quad\forall m\in [1, \ldots, q-1]$.

To apply the filter, we simply process the Fourier coefficients as follows: \[ \mathcal P_N u_N(x) = u_N^\sigma(x) = \sum_{n=0}^N \sigma\left(\frac nN\right)\ \hat u_n\ \phi_n(x), \] where $\phi_n$ are the Fourier basis functions, $\hat u_n$ are the Fourier coefficients, and $\sigma$ is the chosen filter.

Some common filters

Filtration Results

Suppose that $u(x)$ is piece-wise $C^{2p}[0, 2\pi]$ with a single discontinuity at $x~=~\xi.$
If $\mathcal P_N u_N(x)$ is the filtered approximation of the true solution $u(x)$ then we have \[ |u(x) - \mathcal P_Nu_N(x)| \le C\frac{1}{N^{p-1}d(x, \xi)^{p-1}}K(u) + C\frac{\sqrt{N}}{N^{2p}}||u^{(p)}||_{L^2} \] where $d(x,\xi)$ measures the distance of $x$ from $\xi$. Thus, we can recover $(p-1)$-order convergence if

1. $x$ is away from the discontinuity $\xi$.
2. The piecewise continuous function $u$ is regular enough.
3. Order (regularity) of the filter $\sigma$ is sufficiently large.

If $u$ is piecewise analytic then we can increase $p$ with $N$.

We first plot solutions of the linear advection equation at $t=1$, with the initial condition having a discontinuity at $x=\pi$.

No filter: $\sigma \equiv 1$

Step Cut-off: $\sigma(\eta) = \mathbb 1_{[-0.5, 0.5]}$.

Cesáro filter: $\sigma(\eta) = 1 - \eta$

Raised Cosine Filter: $\sigma(\eta) = \frac12\left(1+\cos(\pi \eta)\right)$

Lanczos filter: $\sigma(\eta) = \frac{\sin(\pi \eta)}{\pi \eta}$

Exponential filter: $\sigma(\eta) = \exp(-\alpha \eta^{2p})$

We now apply a filter to Burgers’ equation and plot the solution along with the exact at $t=1.3$.

Burgers' equation, no filter

Burgers' equation, exponential filter ($2p=4$)

Post Processing

The reader should take note that the techniques discussed are only post processing the solution that is computed with purely spectral or pseudospectral methods. This makes sense because altering how one solves for the solution is not as straightforward as merely altering the interpretation of the solution.
That said, the action of filters is similar to that of artificial viscosity. Hence they can stabilise numerical methods like viscosity can, but without introducing stiffness.

With the right choice of a filter for the problem at hand, we are able to remove the effects of the Gibbs phenomenon at points away from the discontinuity. The accuracy of filtered approximations decreases as we approach the point of discontinuity, as is apparent in the plots above. Unsurprisingly, the bound on the pointwise error blows up as $x\to\xi$.

“We assume that we know the first $2N+1$ expansion coefficients of a piecewise analytic function and that the function is known to have a discontinuity at $x=\xi$. The aim is to recover the value of $u$ at any point in the interval $[0,2\pi]$.”

Gibbs Complimentary Basis

Filtering essentially treats Gibbs oscillations as noise, reducing them in a nice manner. It turns out that rather than destroying the information present within the Fourier coefficients, if it is redistributed, one actually recovers an exponentially convergent series to a piecewise analytic function. The convergence holds even at points arbitrarily close to the discontinuity.

The trick is finding a new basis with certain properties. Then the truncated global expansion is re-expanded in this new basis. The catch is that this method relies on knowing the location of the discontinuity.

1. Orthonormality: \[ \left\langle\Phi_k^\lambda, \Phi_l^\lambda\right\rangle_\lambda = \delta_{kl}\quad\forall k,l. \]
2. Spectral convergence: The expansion of a function $g(\xi)$ which is analytic in $[-1, 1]$, in the basis $\Phi_k$ converges exponentially fast with $\lambda$: \[ \max_{\xi\in[-1, 1]}\left|g(\xi) - \sum_{k=0}^\lambda\langle g, \Phi_k^\lambda\rangle_\lambda \Phi_k^\lambda(\xi)\right| \le \exp(-q_1\lambda),\quad q_1 > 0. \]
3. The Gibbs condition: There exists a number $\beta < 1$ such that if $\lambda~=~\beta N$ then $\exists\alpha < 1$ suchthat for $k>N, l\le\lambda$: \[ \left|\left\langle \Phi_l^\lambda(\xi), \Psi_k(x_\xi) \right\rangle_\lambda\right|\ \max_{\xi\in[-1, 1]} |\Phi_l^\lambda(\xi)| \le \left(\frac{\alpha N}{k}\right)^\lambda. \] This means that high modes (large $k$) of the Fourier basis $\Psi_k$ on the lower modes (small $l$) in the basis $\Phi_l^\lambda$ is exponentially small in the interval $-1\le \xi \le 1$ for $\lambda$ proportional to $N$.

Resolving the Gibbs Phenomenon

Let \[\xi(x) = -1 + 2\left(\frac{x-a}{b-a}\right)\] be a map from $x\in[a,b]$ to $\xi\in[-1,1]$.

1. $|(f, \Psi_k)| \le C\quad\forall k$ and
2. $\lim_{N\to\infty} |f(x) - f_N(x)| = 0.$

Let ${\Phi_k^\lambda(\xi)}$ be a Gibbs complimentary basis to the family ${\Psi_k(x)}$, with $\lambda = \beta N$. Furthermore, assume \[ \langle f - f_N, \Phi_l^\lambda\rangle_\lambda = \sum_{k=N+1}^\infty (f, \Phi_l^\lambda)\ \langle\Phi_l^\lambda, \Psi_k\rangle_\lambda. \] Then \[ \max_{x\in[a,b]}\left| f(x) - \sum_{l=0}^\lambda\langle f_N, \Phi_l^\lambda\rangle\ \Phi_l^\lambda(\xi(x)) \right| \le \exp(-qN), \quad q > 0. \]

The knowledge of the interval of analyticity of the true solution is enough to construct an exponentially converging approximation to it, from just the knowledge of the original Fourier coefficients.

Gegenbauer Polynomials

Implicit in the earlier statement is the assumption that there not only exists a Gibbs complimentary basis for the Fourier basis, but also an explicit expression for it is available. Thankfully, Gegenbauer polynomials exist and they do the job.

Gegenbauer polynomials for $\lambda=3$

Results

We present the results of the re-expansions with the linear advection equation and the inviscid Burgers’ equation.

Linear advection equation with $N=16, \lambda=3$.

Linear advection equation with $N=64, \lambda=6$.

Burgers' equation with $N=32, \lambda=4$.

Burgers' equation with $N=64, \lambda=6$.

Burgers' equation with $N=64, \lambda=8$.

The question of whether high-order accuracy can be recovered in the case of non-linear equations is still open. The numerical results seem to confirm this fact.

Conclusions

Spectral methods have numerous advantages over traditional methods for solving partial differential equations. They are highly accurate and have superior phase properties. The above results show us how spectral methods can be combined with post processing to retain those properties, even when the traditional estimates and results do not apply, such as in the case of hyperbolic problems.

Unanswered Questions

1. In practice, how are shock locations computed?
2. Conservation and RH condition.

References

1. D. Gottlieb, J.S. Hesthaven, Spectral methods for hyperbolic problems, Journal of Computational and Applied Mathematics, Volume 128, Issues 1–2, 2001, Pages 83-131, ISSN 0377-0427,
   https://doi.org/10.1016/S0377-0427(00)00510-0.

2. Hesthaven, J., Gottlieb, S., & Gottlieb, D. (2007). Spectral Methods for Time-Dependent Problems (Cambridge Monographs on Applied and Computational Mathematics). Cambridge: Cambridge University Press.