Parallel Spectral Numerical Methods/Finding Derivatives using Fourier Spectral Methods

Finding Derivatives using Fourier Spectral Methods

Spectral methods are a class of numerical techniques that often utilize the FFT. Spectral methods can be implemented easily in Matlab, but there are some conventions to note. First note that Matlab’s “fft” and “ifft” functions store wave numbers in a different order than has been used so far. The wave numbers in Matlab and in most other FFT packages are ordered, $0,1,...,\frac{n}{2},-\frac{n}{2}+1,-\frac{n}{2}+2,...,-1$. Secondly, Matlab does not take full advantage of real input data. The DFT of real data $u$ satisfies the symmetry property $\hat{u}(-k)=\hat{u}(k)$, so it is only necessary to compute half of the wave numbers. Matlab’s ``fft" command does not take full advantage of this property and wastes memory storing both the positive and negative wave numbers. Third, spectral accuracy (exponential decay of the magnitude of the Fourier coefficients) is better for smooth functions, so where possible be sure your initial conditions are smooth – When using a Fourier spectral method this requires that your initial conditions are periodic.

Taking a Derivative in Fourier Space

Let $u_j$ be the discrete approximation of a function $u(x)$ which is sampled at the $n$ discrete points $x_j\in h,2h,...,ih,..,2\pi -h,2\pi$ where $h=2\pi /n$. Now take the Fourier Transform of $u_j$ so that $\text{FFT}(u_j)\equiv {\hat{u}}_k\text{where}k\in \frac{-n}{2}+1,...\frac{n}{2}.$ The Fourier transform of $\frac{{\partial }^{\nu }{u}_j}{\partial x^{\nu }}$ can then be computed from ${\hat{u}}_k$ as shown below:

where ${\hat{u}}_{n/2}=0$ if $\nu$ is odd. More details can be found in Trefethen.

Thus, differentiation in real space becomes multiplication in Fourier space. We can then take the inverse fast Fourier Transform (IFFT) to yield a solution in real space. In the next section we will use this technique to implement forward Euler and backward Euler timestepping schemes to compute solutions for several PDEs.

Template:NumBlk A Matlab program that uses Fourier transformations to compute the first two derivatives of $f(x)={\sin }^2\left(wx\right)$ over $(1,1+2\pi ]$.

% Approximates the derivative of a periodic function f(x) using Fourier
% transforms.  Requires a linear discretization and a domain (a,b] such
% that f(x) = f(x+b-a)

% domain
a = 1;
b = 1 + pi/2;
N = 100;
dx = (b-a)/N;
x = a + dx*(0:N-1);

% function
w = 2;
f = sin(w*x).^2;

% exact derivatives
dfdx = 2*w*sin(w*x).*cos(w*x);
d2fdx2 = 4*w^2*cos(w*x).^2 - 2*w^2;

% fourier derivatives
Nx = size(x,2);
k = 2*pi/(b-a)*[0:Nx/2-1 0 -Nx/2+1:-1];
dFdx = ifft(1i*k.*fft(f));
d2Fdx2 = ifft(-k.^2.*fft(f));

% graph result
figure;
plot(x,dfdx,'r-',x,d2fdx2,'g-',x,dFdx,'k:',x,d2Fdx2,'b:','LineWidth',2);
legend('df/dx','d^2f/dx^2','Fourier df/dx','Fourier d^2f/dx^2');

Exercises

1. Let $u(x)=\sum _k{\hat{u}}_k\exp (ikx)$ be the Fourier series representation of a function $u(x)$. Explain why $\frac{{\mathrm{d}}^{\nu }u}{\mathrm{d}{x}^{\nu }}=\sum (ik{)}^{\nu }{\hat{u}}_k\exp (ikx),$ provided the series converges.

2. Consider the linear KdV equation $u_t+u_{xxx}=0$ with periodic boundary conditions for $x\in (0,2\pi ]$ and the initial data

   $u(x,0)=0\text{ if }0<x\le \pi$

   and

   $u(x,0)=1\text{ if }\pi <x\le 2\pi$

3. Using separation of variables, show that the “solution” is $u(t,x)=\frac{1}{2}-\frac{2}{\pi }\sum _{j=0}^{\mathrm{\infty }}\frac{\sin ((2j+1)x-(2j+1{)}^{3}t)}{2j+1}.$ Quotation marks are used because the expression for the solution that is given does not converge when differentiated either once in time or twice in space.

4. As explained by Olver, this solution has a fractal structure at times that are an irrational multiple of $\pi$ and a quantized structure at times that are rational multiples of $\pi$. The Matlab program in listing Template:EquationNote uses the Fast Fourier transform to find a solution to the linearized KdV equation. Explain how this program finds a solution to the linearized KdV equation.

5. Compare the numerical solution produced by the Matlab program with the analytical solution. Try to determine which is more accurate and see if you can find evidence or an explanation to support your suggestions.

Template:NumBlk Matlab Program Code download

% This program computes the solution to the linearly dispersive
% wave equation using the Fast Fourier Transform

N = 512; % Number of grid points.
h = 2*pi/N; % Size of each grid.
x = h*(1:N); % Variable x as an array.
t = .05*pi; % Time to plot solution at
dt = .001; % Appropriate time step.
u0 = zeros(1,N); % Array to hold initial data
u0(N/2+1:N)= ones(1,N/2); % Defining the initial data
k=(1i*[0:N/2-1 0 -N/2+1:-1]); % Fourier wavenumbers
k3=k.^3;
u=ifft(exp(k3*t).*fft(u0)); % Calculate the solution
plot(x,u,'r-'); % Plot the solution
xlabel x; ylabel u; % Label the axes of the graphs
title(['Time ',num2str(t/(2*pi)),' \pi']);