Spectral method (FFT scheme)

With the methods that we have used so far, we have been interested mainly in calculating the new value of the function in the same space as the solution exists in. For example, the diffusion equation has been solved in each example above in position space. However, there are times when it can be useful to perform the evolution of a function in some other space. Spectral methods use these ideas by Fourier transforming the problem and performing the calculation in Fourier space, it being easier in some cases for the solution to be calculated there than in position space. The discussion here follows that in A First Course in Computational Physics by P. DeVries.

The process that spectral methods use is to take an evolving function in some space, transform it into Fourier space, calculate the evolution of the function in Fourier space, and then transform back to the original space by taking the inverse Fourier transform. This process is then iterated for each time step of the algorithm.

Let's describe this process in more detail by considering the diffusion equation^3.10:

$$\frac{\partial }{\partial t} A(x,t) = \kappa \frac{\partial ^2}{\partial x^2} A(x,t).$$ (3.42)

The Fourier transform of $A(x,t)$ is

$$\alpha(k,t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty A(x,t) e^{-ikx} dx.$$ (3.43)

This has transformed the function $A$ from $x$-space to $k$-space. Note that I will use the terms $x$-space, position space and inverse space interchangeably, as well as the terms $k$-space, momentum space and Fourier space. We can now rewrite $A(x,t)$ as

$$A(x,t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty \alpha(k,t) e^{+ikx} dk.$$ (3.44)

These transformations between $x$ and $k$ are occurring for the same instance in time. Equation (3.44) allows us to rewrite Equation (3.42) as

$$\frac{\partial }{\partial t} \frac{1}{\sqrt{2\pi}} \int_{-\infty}... ...rtial x^2} \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty \alpha(k,t) e^{+ikx} dk.$$ (3.45)

On the left-hand side of this equation, the partial derivative operates only on $\alpha(k,t)$ since it is the only thing dependent upon $t$. On the right-hand side, the partial derivative acts only on the exponential, since this is the only thing dependent upon $x$. We take the derivatives through the integral signs and act them on the things affected.

$$\frac{\partial }{\partial t} \int_{-\infty}^\infty \alpha(k,t) e^{+ikx} dk = \kappa \frac{\partial ^2}{\partial x^2} \int_{-\infty}^\infty \alpha(k,t) e^{+ikx} dk,$$ (3.46)

$$\Rightarrow \int_{-\infty}^\infty \frac{\partial }{\partial t} \alpha(k,t) e^{+ikx} dk = \kappa \int_{-\infty}^\infty \alpha(k,t) \frac{\partial ^2}{\partial x^2} e^{+ikx} dk,$$ (3.47)

$$\Rightarrow \int_{-\infty}^\infty \frac{\partial }{\partial t} \alpha(k,t) e^{+ikx} dk = \kappa \int_{-\infty}^\infty \alpha(k,t) \bigl[ -k^2 \bigr] e^{+ikx} dk.$$ (3.48)

We now multiply both sides of this equation by $e^{-ik'x}$. This is is just a mathematical trick, since we can ``see'' a possibility of constructing a Dirac delta function out of an integral of exponentials. Anyway, to carry on with the explanation, we then integrate both sides over $x$, giving

$$\int_{-\infty}^\infty \frac{\partial }{\partial t} \alpha(k,t) e^{+ikx} dk = \kappa \int_{-\infty}^\infty \alpha(k,t) \bigl[ -k^2 \bigr] e^{+ikx} dk,$$ (3.49)

$$\Rightarrow \int_{-\infty}^\infty e^{-ik'x} \int_{-\infty}^\infty \frac{\partial }{\partial t} \alpha(k,t) e^{+ikx} dk dx = \kappa \int_{-\infty}^\infty e^{-ik'x} \int_{-\infty}^\infty \alpha(k,t) \bigl[ -k^2 \bigr] e^{+ikx} dk dx.$$ (3.50)

We now exchange the order of integration to give

$$\int_{-\infty}^\infty \int_{-\infty}^\infty \frac{\partial }{\partial t} \alpha(k,t) e^{-ik'x} e^{+ikx} dx dk = \kappa \int_{-\infty}^\infty \int_{-\infty}^\infty \alpha(k,t) \bigl[ -k^2 \bigr] e^{-ik'x} e^{+ikx} dx dk,$$ (3.51)

$$\Rightarrow \int_{-\infty}^\infty \int_{-\infty}^\infty \frac{\partial }{\partial t} \alpha(k,t) e^{+i(k-k')x} dx dk = \kappa \int_{-\infty}^\infty \int_{-\infty}^\infty \alpha(k,t) \bigl[ -k^2 \bigr] e^{+i(k-k')x} dx dk.$$ (3.52)

We note that

$$\int_{-\infty}^\infty e^{ik(k-k')x} dx = 2 \pi\delta(k-k'),$$ (3.53)

and performing the integral over $k$, this function will pick out all values of $k$ that are equal to $k'$ and we get the differential equation

$$\frac{\partial }{\partial t} \alpha(k,t) = -k^2 \kappa \alpha(k,t),$$ (3.54)

where we have done the trivial relabeling of $k'$ to $k$. Note that we now only have one first order partial derivative, instead of a first order and a second order partial derivative--this has simplified the calculation rather a lot. Equation (3.54) has the exact analytical solution

$$\alpha(k,t) = e^{-k^2 \kappa t} \alpha(k,0),$$ (3.55)

where $\alpha(k,0)$ is calculated from the initial distribution of $A = A(x,0)$ via the Fourier transform

$$\alpha(k,0) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty A(x,0) e^{-ikx} dx.$$ (3.56)

This has completely solved the problem, and we don't really need to use a computer to do the time evolution for us. however, it gives us a technique useful in more general and more difficult problems. Once we have $\alpha(k,t)$ we just take the inverse Fourier transform to obtain the desired solution:

$$A(x,t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty \alpha(k,t) e^{+ikx} dk,$$ (3.57)

$$= \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty e^{-k^2 \kappa t} \alpha(k,0) e^{+ikx} dk.$$ (3.58)

Overall, we have three steps:

1. Transform the original differential equation and initial conditions into Fourier space.
2. Solve the (hopefully simpler) equation in transform space.
3. Perform the inverse Fourier transform to get the solution in the original variables.

The only problem now is to know what $k$ values to use once the solution is in Fourier space. This is quite a problem with this method as such a choice depends upon the implementation of the Fast Fourier Transform (FFT) in the language or with the system you are using. In the rest of this discussion I shall assume that the implementation is Matlab's fft. Imagine that we have $n_x$ grid points in space, and we assume that they are spaced equally at intervals of $\Delta x$. The range of $x$ points is therefore $X = (nx-1) \Delta x$. Since we only have a finite number of $x$ points^3.11 each one must map to a point in $k$-space, and so the number of $k$ points is equal to $n_x$. In Matlab's implementation of the FFT, the points in $k$ are separated by an amount $\Delta k$ which is related to the parameters in $x$ by

$$\Delta k = \frac{1}{X} = \frac{1}{(n_x-1) \Delta x}.$$ (3.59)

In the problems solved in this course, the $x$ values lie roughly equally about $x=0$. We therefore set up our $k$ values in a similar manner, where we have a set $k=0$ (zero momentum) position and successively larger values of $k$ on either side, each separated by $\Delta k$ until a sufficient number of points is obtained. More explicitly this is

$$k = \cdots, -\Delta k, 0, \Delta k, \cdots$$ (3.60)

where the number of points in $k$ is equal to the number of points in $x$. So, if $n_x$ were equal to 5, with $\Delta x = 1$, then we would have the $x$ points:

$$x = -2, -1, 0, 1, 2,$$ (3.61)

the range of $x$ points is $4( = 2 - -2)$, so $\Delta k = \frac{1}{4}$ and hence we have the $k$ points:

$$k = -\frac{1}{2}, -\frac{1}{4}, 0, \frac{1}{4}, \frac{1}{2}.$$ (3.62)

This technique seems pretty tricky, but it seems to be fairly specific to certain problems. What if we want to solve a differential equation that is more complicated than the diffusion equation, as is often the case in real life? Well, a possible answer to this is to use what is known as the pseudo-spectral method (also known as the split-step FFT scheme), which we discuss in the next section.

Footnotes

... equation^3.10

I told you you'd see a lot of the diffusion equation!

... points^3.11

Note that we are using a Fast Fourier Transform (FFT) here, which is derived from the Discrete Fourier Transform (DFT). These transforms are, however, not Fourier transforms in the strict sense since they are discretised finite data-set approximations of a continuous set of data with infinite extent. Although they are related remember that they aren't the same thing. I will throw the terms around here a bit, but just keep this point in mind.