Two approaches to solving diffusion equation in Fourier space

Question

I want to numerically solve the diffusion equation $\partial_t u = D \partial_x^2 u$ in Fourier space, and can think of multiple ways to do it.

Setup

Option 1

Differentiating $u$ twice in Fourier space brings down two factors of $ik$ from the exponential $e^{i k x}$.

$$\partial_t \tilde{u}_k = D (ik)^2 \tilde{u}_k$$

which can be discretised, for example, like:

$$\frac{\tilde{u}^{n+1} - \tilde{u}^n}{\Delta t} = -Dk^2\left(\frac{\tilde{u}^{n+1} + \tilde{u}^{n}}{2} \right)$$

Option 2

Taking the second order finite difference equation in real space:

$$\partial_x^2 u \approx \frac{u_{i+1} - 2 u_i + u_{i-1}}{\Delta x^2}$$

where $u_{i+1}$ and $u_i$ are $\Delta x $ apart. Fourier transforming this, remembering that phase shift of $\Delta x$ is given by $e^{\Delta x k}$ we get:

$$\partial_x^2 \tilde{u}_k \approx \frac{e^{+\Delta x k} -2 + e^{-\Delta x k} }{\Delta x^2} \tilde{u}_k = 2 \frac{\cos(\Delta x k) - 1}{\Delta x^2} \tilde{u}_k\qquad(1)$$

Applying again to the diffusion equation:

$$\frac{\tilde{u}^{n+1}_k - \tilde{u}^n_k}{\Delta t} = D\frac{2 (\cos(\Delta x k)-1)}{\Delta x^2} \left(\frac{\tilde{u}^{n+1}_k + \tilde{u}^{n}_k}{2} \right)$$

I came across Option 2 in a few places, for example, in this paper.

Question

Which of these (if either) should I be using, and why? I'm guessing the discrete nature of the FFT is important, as the approaches are equivalent as $\Delta x \to 0$, seen by Taylor expanding Eq. (1):

$$2 \frac{\cos(\Delta x k) - 1}{\Delta x^2} = 2 \frac{1 - \frac{(\Delta x k)^2}{2} + \mathcal{O}(\Delta x^4) - 1}{\Delta x^2} = -k^2 + \mathcal{O}(\Delta x^2)$$

Any help is greatly appreciated!

Update

Thinking more about this, the properties of the DFT are actually important. The maximum $k$ is $k_{max} =2\pi (N/2)$ where $\Delta x = L / N$. So as $\Delta x \to 0$, $k_{max} \to \infty$. Therefore the two methods are different for large $k$ at any $\Delta x$:

This looks like the latter option is applying some smoothing at high $k$.

Answer 1

Typically, the first approach is used. The range of frequencies used will limit an explicit time integrator, and is well-known for finite-difference discretizations. The relationship to the finite difference case follows by noting that $D/(\Delta{x}^{2})$ is a frequency, so for low frequencies, you can take large time steps, but for high frequencies, stability limits the time step. Since the Fourier basis consists of eigenfunctions of the diffusion operator, another way to look at it is that the Fourier-transformed spectral discretization yields a collection of scalar test equations that can be integrated independently (it diagonalizes the operator), and then you can use standard arguments from numerical stability theory for ODEs to high frequencies will limit your time step due to stability. You should probably consider an implicit time integrator for the high-frequency modes, unless you are willing to use small time steps due to this stability limitation.