david/stabilized_spectral_methods
The Exact Solution Operator
We start with the heat equation \[u_t = \kappa u_xx\] Introducing the Fourier transform of the initial condition: \[\hat{u}_0(\xi) = \frac{1}{\sqrt{2\pi}} \int_{x=0}^{1} u(x,0) e^{-i\xi x} dx\] The exact solution is \[u(x,t) = \frac{1}{\sqrt{2\pi}} \int_{\xi=-\infty}^\infty \hat{u}_0(\xi) e^{i\xi x-\xi^2t}\] We will consider the domain \(0\le x\le 1\), with grid points \(x_j = jh, \ j=0,1,\dots,m-1\) where \(h=1/m\). Defining \(\hat{U}(\xi)\) as the discrete Fourier transform of the initial condition: \[\hat{U}(\xi) = \frac{h}{\sqrt{2\pi}} \sum_{j=0}^{m-1} U_j^0 e^{-i\xi x_j}\]
The exact solution is \[U^n_j = \sum \hat{U}(\xi) e^{-\xi^2 t}\]
In the following, it is convenient to represent the discrete Fourier transform and its inverse as matrix multiplication; we use \(R\) to denote this matrix. Thus \(\hat{U}=RU\) and \(U=R^{-1}\hat{U}\). We can write the exact solution as \[U^1 = R^{-1} [-\kappa \xi^2 \Delta t] R U^0.\] Here \([f(\xi)]\) denotes the diagonal matrix whose \(k\)th entry is \(f(\xi_k)\). In other words, the solution operator amounts to multiplying each Fourier mode by the factor \(-\kappa \xi^2 \Delta t\).
In fact, we can construct the exact solution operator in similar manner for any linear evolution PDE...
Hence a straightforward spectral semi-discretization is \[U'(t) = R^{-1} [-\kappa \xi^2] R U\] In fact this is exact. Supposing that we use the explicit Euler method in time, we obtain the full discretization
\[\begin{align*} U^{n+1} & = U^n + \Delta tR^{-1} [-\kappa \xi^2] R U^n \\ & = R^{-1} [1-\Delta t\kappa \xi^2 ] R U^n\end{align*}\]
This leads to a severe timestep restriction. We now consider methods of the same type, but where a different set of eigenvalues are used.
Let us take the opposite route, discretizing first in time. We will again use the explicit Euler method, but we will require that our full discretization be exact: \(U^n_j = e^{i\xi j h - \kappa \xi^2 t}.\) Inserting this gives
\[\begin{align*} \frac{U^{n+1}_j - U^n_j}{\Delta t} & = f(\xi) U^n_j \\ \frac{e^{-\kappa \xi^2\Delta t}-1}{\Delta t} U^n_j & = f(\xi) U^n_j \\ \frac{e^{-\kappa \xi^2\Delta t}-1}{\Delta t} & = f(\xi).\end{align*}\]
This leads us to the method \[U^{n+1} = U^n + \Delta tR^{-1}[f(\xi)]R U^n.\] From one point of view, this method is quite remarkable – it solves the (stiff) diffusion equation using the explicit Euler method, but is unconditionally stable! In fact, we have just reconstructed the exact solution operator. Notice that our choice of \(f(\xi)\) yields eigenvalues that lie (as they must) on the boundary of the explicit Euler stability region.
Again, we can construct the exact solution operator for any linear evolution equation using (say) Euler’s method and the dispersion relation...
Next, we consider the advection-diffusion equation \[u_t = -v u_x + \kappa u_{xx}.\] Using the approach above, we can construct the exact discrete solution: \[U'(t) = R^{-1} [(-iv\xi - \kappa \xi^2)\Delta t] R U.\] Hence the obvious semi-discretization is (written in an operator-split form) \[U'(t) = R [-iv\xi] R^{-1} U + R [-\kappa \xi^2] R^{-1} U,\] leading to (with Euler’s method)
\[\begin{align*} U^{n+1} & = U^n + \Delta tR [-iv\xi] R^{-1} U^n + \Delta tR[-\kappa \xi^2] R^{-1} U^n \\ & = R [1-\Delta t(iv\xi -\kappa \xi^2)] R^{-1} U^n\end{align*}\]
On the other hand, if we combine our discretizations of the advection equation and the heat equation above in an operator splitting approach, we have \[U^{n+1} = U^n + R \left[\frac{e^{-iv\xi}-1}{\Delta t}\right] R^{-1} U^n + R \left[\frac{e^{-\kappa \xi^2}-1}{\Delta t}\right] R^{-1} U^n.\]
(compare the eigenvalues)
Burger’s equation
We now develop an explicit discretization for the dissipative Burger’s equation that is convergent for a fixed CFL number.
