Why not use the convolution theorem for explicit timestepping?

Question

Consider the advection equation

\begin{equation} \frac{\partial C}{\partial t} + u\frac{\partial C}{\partial x} + v\frac{\partial C}{\partial y} = 0 \end{equation}

I want to do a forward time, center space finite difference approximation with the Lax method and a periodic boundary

\begin{equation} C_{i, j}^{n+1} = \frac{1}{4}(C_{i-1, j}^{n}+C_{i+1, j}^{n}+C_{i, j-1}^{n}+C_{i, j+1}^{n}) - \frac{\Delta t}{2\Delta x}\left[u(C_{i+1, j}^{n}-C_{i-1, j}^{n})+v(C_{i, j+1}^{n}-C_{i, j-1}^{n})\right] \end{equation}

This is essentially an image convolution operation with a 3-by-3 kernel, and it seems obvious that we can apply the convolution theorem to optimize the algorithm with FFT. With this trick, the computation time barely depends on the number of timesteps $N$ (because I just need to calculate the $N$-th power of the FFT-transformed kernel and multiply it with the FFT-transformed data matrix representing the initial condition), so calculating millions of steps forward takes roughly the same time as a single step. However, this is not a popular method in practice, which is odd given its tremendous speedup.

What am I missing? What is the catch?

Answer 1

This is a linear PDE, and so while this technique works here, it would not work for any nonlinear PDE. Often times when people are solving these equations it is to get experience with common solution techniques (i.e. in a class setting) and this would not be applicable to more general scenarios. Also, taking the $N^{th}$ power of the matrix should be roughly the same cost as doing the explicit time stepping as the explicit time stepping over N iterations is literally applying the time-evolution matrix $N$ times.

Answer 2

For the linear, constant coefficient advection equation on a torus, one can simply use the exact solution. So there are no "popular" numerical methods for this problem.

Answer 3

Reharding the question in the comment: After discretization you get a system of the form

$ \partial_t C = (M_x + M_y) C $

where $M_x$ and $M_y$ are commuting matrices obtained by discretizing the time-independent linear system operator. In your case, these are just the differential operators, but they can also be more general operators (e.g. $x^2 \partial_x^2$).

The solution to your initial value problem then is given by

\begin{align} C(t,x,y) &= \exp(t(M_x + M_y)) C(t=0,x,y)\\ &= \exp(t M_x) \exp(tM_y) C(t=0,x,y)\\ \end{align}

If you now are able to evaluate the matrix exponentials and apply them to the initial value function, you obtain the exact solution within the chosen discretization.

In general, matrix exponentiation takes $O(N^3)$ using e.g. scaling-and-squaring or an eigendecomposition, and thus exact time-propagation is possible only for small dimensions/basis sizes, or systems where a separation of the dimensions is possible (as im your case). Therefore, for larger and/or coupled systems, one has use Krylov subspace methods or the whole bunch of other (approximate) propagation methods, which usually only scale linearly in the system dimension.

In your case, however, things are more easy, as you can simply choose Fourier eigenfunctions for the discretization. The resulting matrices then are diagonal and can thus be readily exponated by exponentation of the diagonal entries. This is the numerical parallel to the fact mentioned by @DavidKetcheson.

Where's the catch? Nowhere. People just tend to use the tools they have to each and any equation, and thus often apply general integrators to cases that can be solved much more easily. But you'll also won't revolutionize things if you take more care.