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 ?

  • 1
    $\begingroup$ 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. $\endgroup$
    – 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.


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.