Parallel Spectral Numerical Methods/Finding Derivatives using Fourier Spectral Methods

From Wikibooks, open books for an open world
< Parallel Spectral Numerical Methods
Jump to: navigation, search

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, 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[edit]

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 \qquad \text{where} \quad k \in {\frac{-n}{2}+1,...\frac{n}{2}}.
The Fourier transform of \frac{\partial^2 u_j}{\partial x^2} can then be computed from \hat{u}_k as shown below:

\text{FFT}(\frac{\partial^\nu u_j}{\partial x^\nu}) \equiv (ik)^\nu \hat{u}_k

 

 

 

 

( 1)

where \hat{u}_{n/2}=0 if \nu 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 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[edit]

  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. [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\quad\text{ if }\quad 0< x\leq\pi

    and

     u(x,0)= 1\quad\text{ if }\quad \pi<x\leq2\pi

    1. Using separation of variables, show that the “solution” is u(t,x)=\frac{1}{2}-\frac{2}{\pi}\sum_{j=0}^{\infty}\frac{\sin((2j+1)x-(2j+1)^3t)}{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.

    2. As explained by Olver[3], 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 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.

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

    Notes[edit]

    1. Trefethen (2000)
    2. This question was prompted by an REU and UROP project due to Sudarshan Balakrishan which is available at http://www.lsa.umich.edu/math/undergrad/researchandcareeropportunities/infinresearch/pastreuprojects.

    3. Olver (2010)

    References[edit]

    Olver, P.J. (2010). "Dispersive Quantization". American Mathematical Monthly 117: 599-610. 

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