Parallel Spectral Numerical Methods/Finding Derivatives using Fourier Spectral Methods

Finding Derivatives using Fourier Spectral Methods edit

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,  . Secondly, Matlab does not take full advantage of real input data. The DFT of real data   satisfies the symmetry property  , 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 edit

Let   be the discrete approximation of a function   which is sampled at the   discrete points   where  . Now take the Fourier Transform of   so that   The Fourier transform of   can then be computed from   as shown below:






( 1)

where   if   is odd. More details can be found in Trefethen[1].

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.





( A)

A Matlab program that uses Fourier transformations to compute the first two derivatives of   over  .
% 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
legend('df/dx','d^2f/dx^2','Fourier df/dx','Fourier d^2f/dx^2');

Exercises edit

  1. Let   be the Fourier series representation of a function  . Explain why   provided the series converges.

  2. [2] Consider the linear KdV equation   with periodic boundary conditions for   and the initial data




  3. Using separation of variables, show that the “solution” is   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[3], this solution has a fractal structure at times that are an irrational multiple of   and a quantized structure at times that are rational multiples of  . The Matlab program in listing B 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.





( B)

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
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']);

Notes edit

  1. Trefethen (2000)
  2. This question was prompted by an REU and UROP project due to Sudarshan Balakrishan which is available at

  3. Olver (2010)

References edit

Olver, P.J. (2010). "Dispersive Quantization". American Mathematical Monthly. 117: 599–610. {{cite journal}}: Cite has empty unknown parameter: |coauthors= (help)

Trefethen, L.N. (2000). Spectral Methods in Matlab. SIAM.