4
$\begingroup$

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

plot comparing different methods

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – user7257 Jul 11 '14 at 13:50
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I hadn't really thought about stability varying for each independent k equation, but that makes sense, interesting! I might be missing the point, but couldn't either option be used with an implicit timestepping scheme; if so is there then a preference for either option purely in terms of the spacial discretisation? $\endgroup$ – Hemmer Jul 10 '14 at 19:19
  • $\begingroup$ 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$. $\endgroup$ – Geoff Oxberry Jul 10 '14 at 20:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.