# Why not use the convolution theorem for explicit timestepping?

$$$$\frac{\partial C}{\partial t} + u\frac{\partial C}{\partial x} + v\frac{\partial C}{\partial y} = 0$$$$

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

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

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?

• What is u and v? If those are just constants, the equation can be easily solved numerically exactly by separation of the variables. Moreover, this is the same for any operator that can be diagonalized numerically, let it be with the FFT or not. Where's the catch? Most relevant PDEs don't fall under this class, and if they do, people already have applied this approach. Apr 7, 2022 at 21:02
• @davidhigh Yes, both u and v are constants. Are you saying you can give me an analytical solution by tackling it analytically? I'm curious about what such a solution would look like. Apr 7, 2022 at 23:44

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.

• Taking the Nth power of something has roughly log(N) cost. Apr 7, 2022 at 8:58
• In addition, the method that you describe is appallingly dissipative. If you did take a million time steps, the outcome would be worthless. So not a very valuable learning experience. Nevertheless, you do well to wonder about these things. Apr 7, 2022 at 17:18

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.

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.

• Oh, that's interesting. Could you get into more details regarding what is $M_x$ and $M_y$? I think you are saying they can be diagonal. I also don't quite understand $\partial_t C = (M_x + M_y) C$, because $\partial_t C$ and $C$ are scalars, whereas $M_x + M_y$ is a matrix. How can a (non-trivial) matrix be equal to a scalar? Apr 8, 2022 at 14:24
• ... and does not the method you described work for boundary conditions other than the periodic boundary? Apr 8, 2022 at 14:38
• C is a scalar function, which upon discretization in the chosen basis is expressed as a vector of expansion coefficients, while the linear operator becomes a matrix in this basis. These concepts are independent of the boubdary conditions. For non-periodic boundary conditions, however, you get different set of eigenfunctions (not Fourier modes). You can obtain these by diagonalizing the differential operator matrix in a basis which satisfies the boundary conditions you want. Apr 8, 2022 at 15:20