# Two approaches to solving diffusion equation in Fourier space

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$.

• You should get a copy of "Numerical Solutions of Partial Differential Equations: Finite Difference Methods" by GD Smith. It will help you understand stability analysis related to these schemes. Jul 11, 2014 at 13:50

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.
• I suppose either option could be used. Like I said, typically, the first approach is used, because it diagonalizes the diffusion operator in 1-D. The second approach does not necessarily diagonalize the resulting discrete approximation of the diffusion operator, although I'd think the two would agree in the limit as $\Delta{x} \rightarrow 0$. Jul 10, 2014 at 20:37