Computing Fourier representation of space dependent advection operator via FFT

Consider the following equation on the circle:

$$\dfrac{\partial p(x,t)}{\partial t} = a(x)\dfrac{\partial p(x,t)}{\partial x} \equiv L(p) \enspace ,$$

where $L$ is the operator acting on $p(x,t)$.

Now, I would like to create a matrix of this operator $L$ using FFT/IFFT in Matlab. Note that since our domain is the circle, $x\in[0,2\pi]$

Taking Fourier transform, we get $\hat L(p)= \left[a(x)\dfrac{\partial p(x,t)}{\partial x}\right]^{\wedge}= \left[ a(x) (D_x \hat{p})^{V}\right]^{\wedge}$, where $\wedge$ and $V$ are FFT and IFFT operations, and $D_x$ is the $x$-derivative matrix in Fourier (formed by multiplying the $k$th fourier coefficient by $ik$).

My question is: From the expression I was able to solve the PDE using an ODE solver such as ODE45. But, I am interested in extracting the operator for $L$ to be used for other purposes. Hence, I am trying to extract matrix $\hat{L}$ in the following equation

$$\dfrac{\partial \hat{p}(x,t)}{\partial t}=\hat{L}\hat{p}$$

Any ideas ?

• If you have something which computes the matrix-vector product, you can compute the $j$th matrix column by column by feeding in the vector $e_j$ in, where $e_j$ is zero except for the $j$th entry. – Jesse Chan Aug 18 '15 at 18:45

The FFT algorithm computes the discrete Fourier transform (DFT), which transforms $x_{1:n}$ to $\hat x_{1:n}$ by $$\hat x_k = \sum_{j=1}^{n} x_j \, e^{-2\pi i (j-1) (k-1)/n}.$$ The DFT is a linear transform, as is clear from its formula, meaning that it computes multiplication by a matrix $F$: $$\hat x = Fx, \qquad F = \big( e^{-2\pi i (j-1) (k-1)/n} \big)_{k,j}.$$ IFFT then implements multiplication by $F^{-1}$.
So when you are computing $(a(x)(D_x \hat p)^\vee)^\wedge$, you are computing $$\hat L\hat p = F A F^{-1} D_x \hat p,$$ where $D_x$ multiplies by $i k$ (a diagonal matrix), and $A$ multiplies by $a(x)$ (a diagonal matrix). The only thing to be done here is to construct $F$ explicitly, and to form the matrix product explicitly. There may sometimes be a difference in normalization conventions and such, so $F$ has to be consistent with what the FFT actually computes.
Notice also that $FAF^{-1}$ here is like an application of the change-of-basis formula for the matrix $A$ from space basis to frequency basis.