---
title: 8-0. Fourier approximations
tags: Spectral method
---

**Textbook**
1. [Spectral method in Matlab](https://people.maths.ox.ac.uk/trefethen/spectral.html) by *Lloyd N. Trefethen*
2. [Finite difference and Spectral methods](https://people.maths.ox.ac.uk/trefethen/pdetext.html) by *Lloyd N. Trefethen*
3. A practical guide to pseudospectral method by *Bengt Fornberg*

**Suggested reading**
1. [Notes on FFT-based differentiation](https://math.mit.edu/~stevenj/fft-deriv.pdf)

---

## Example of FFT in `Matlab`

### Example 1: FFT and trigonometric functions

Consider $f(x) = \sin(x)$ and its corresponding Fourier transform:

N = 10;
x = linspace(0, 2*pi, N+1);
x = x(1:N);

f = sin(x);
fft(f)

The result shows

0 -5i 0 0 0 0 0 0 0 5i

**Remark:** One can see easily that the Fourier modes are scaled by the number of points, so if one double the points, the result of the mode should be doubled.

#### Exercise 1:
* Try the following grid and see what you found in the Fourier coefficients of $\sin(x)$

N = 20;
h = 2*pi/N;
x = (1:N)*h;

---

### Example 2: Differentiation with FFT

**Remark** The treatments for the 1st derivaative and the 2nd derivative are different at the highest Fourier mode $(k=M/2)$.

% grid points
  N = 24; h = 2*pi/N; x = h*(0:N-1)';

% v: the function
% vp1: dv/dx
% vp2: d^2v/dx^2
  v = exp(sin(x));
  vp1 = cos(x).*v;
  vp2 = cos(x).*vp1 -sin(x).*v;
  v_hat = fft(v);
  
% weights: Fourier modes for odd deriv.
  kk1 = [0:N/2-1 0 -N/2+1:-1]';
  
% w1: spectral approximation of dv/dx
  w1_hat =  1i*kk1.*v_hat;
  w1 = real(ifft(w1_hat)); 

  error = norm(w1-vp1,inf);
  disp(['max error for 1st deriv. = ' num2str(error)])
  
% weights: Fourier modes for even deriv.
  kk2 = [0:N/2 -N/2+1:-1]';
  
% w2: spectral approximation of d^2v/dx^2
  w2_hat =  -(kk2.^2).* v_hat;
  w2 = real(ifft(w2_hat)); 

  error = norm(w2-vp2,inf);
  disp(['max error for 2nd deriv. = ' num2str(error)])

The result shows

max error for 1st deriv. = 9.5757e-13
max error for 2nd deriv. = 2.0521e-12

#### Exercise 2:
* Try to use `kk1` as Fourier modes to compute the 2nd derivative and see how the result changes
* Use Matlab command `semilogy(abs(v_hat))` to determine how many grid points are required
* Change number fo grid points `N` to see how the result changes

---

### Example 3: Spatial shift

Given
$$
u(x) = \sum_k \hat{u}_k e^{ikx},
$$
it is easy to see that a spatial shift to the function can be written as
$$
u(x-c) = \sum_k \hat{u}_k e^{ik(x-c)} = \sum_k e^{-ikc}\hat{u}_k e^{ikx},
$$
where $c$ is a constant.

% grid points
  N = 24; h = 2*pi/N; x = h*(0:N-1)';

% v: the function
  v = exp(sin(x));
  v_hat = fft(v);
  
% weights: Fourier modes for odd deriv.
  kk = [0:N/2 -N/2+1:-1]';

% c: spatial shift
  c = 1;
  vc_hat = exp(-1i*kk*c).*v_hat;
  vc = real(ifft(vc_hat));

% Compare the original one and the shifted one
  figure
  plot(x, v, x, vc)

---

### Example 4: Adding points

Given $N$ sample points of a function $u(x)$, we use Fourier transform to obtain the Fourier coefficients and interpolate the function as
$$
u(x) \approx u_N(x) = \sum^{N/2}_{-N/2+1} \hat{u}_k e^{ikx}.
$$
Suppose we want to evaluate $2N$ sample points from $u_N(x)$, we can rewrite the equation as
$$
u_{2N}(x) = \sum^{N}_{-N+1} \tilde{u}_k e^{ikx},
$$
where
$$
\tilde{u}_k = \hat{u}_k, \quad k=-N/2+1,\cdots,N/2-1,
$$
and
$$
\tilde{u}_k = \frac{1}{2}\hat{u}_k, \quad k=-N/2,N/2,
$$
and $\tilde{u}_k = 0$, for all the rest of the modes.

**Remark** The implementation in `Matlab` requires a bit careful as the multiplication scalar in FFT and iFFT are defined differently.

% grid points
  N = 24; h = 2*pi/N; x = h*(0:N-1)';

% v: the function
  v = exp(sin(x));
  v_hat = fft(v);

% new points
  M = 2*N;
  h2 = 2*pi/M; x2 = h2*(0:M-1)';
  v_tilde = zeros(M,1);
  v_tilde(1:N/2) = v_hat(1:N/2);
  v_tilde(N/2+1) = 0.5*v_hat(N/2+1);
  v_tilde(M-N/2+1) = 0.5*v_hat(N/2+1);
  v_tilde(M-N/2+2:M) = v_hat(N/2+2:N);
% 
  v_tilde = v_tilde*M/N;
  
  plot(x, v, '*-', x2, real(ifft(v_tilde)), 'ro-')

---

### Example 5: Aliasing error visualization

We use fine grid to evaluate the "exact" Fourier modes and see how well the approximation on the modes is for coarse grid.

% f: the function
  f = @(x) exp(sin(x)) + sin(15*x);

% coarse grid points
  N = 20; h = 2*pi/N; x = h*(0:N-1)';

% vc: sampling of the function at coarse grids
  vc = f(x);
  vc_hat = fft(vc);

% fine grid points
  M = 60; h2 = 2*pi/M; y = h2*(0:M-1)';

% vf: sampling of the function at fine grids
  vf = f(y);
  vf_hat = fft(vf);
  
% Comparing the two
  clf
  semilogy(abs(vf_hat(1:M/2+1)/M), 'o')
  hold on
  semilogy(abs(vc_hat(1:N/2+1)/N), '*')
  semilogy(abs(vc_hat(1:N/2+1)/N-vf_hat(1:N/2+1)/M), '+')
  legend('exact', 'approximate', 'difference')

---

### Example 6: Aliasing error

Consider taking derivative of the square of a function
$$
\frac{d}{dx} u^2(x) = 2u(x)u'(x)
$$
Let's see the difference in the following two approaches

% grid points
  N = 24; h = 2*pi/N; x = h*(0:N-1)';

% v: the function
  v = exp(sin(x));
  v_hat = fft(v);
  
% ve: the exact solution of (v^2)'
  ve = exp(2*sin(x)).*(2*cos(x));

% v2: v(x)^2
  v2 = v.^2;
  v2_hat = fft(v2);

% weights: Fourier modes for odd deriv.
  kk1 = [0:N/2-1 0 -N/2+1:-1]';

% w2: spectral approximation of (v^2)_x
  w2_hat =  1i*kk1.*v2_hat;
  w2 = real(ifft(w2_hat));

  error = norm(w2-ve,inf);
  disp(['max error for (v^2)_x = ' num2str(error)])


% w1: spectral approximation of v_x
  w1_hat =  1i*kk1.*v_hat;
  w1 = real(ifft(w1_hat));

  error = norm(2*v.*w1-ve,inf);
  disp(['max error for 2*v*v_x = ' num2str(error)])

The result shows

max error for (v^2)_x = 8.2602e-09
max error for 2*v*v_x = 2.883e-12

#### Exercise 6:
* Consider $u(x)=\sin(3x)$, using $10$ sample points to approximate this function. Do we have enough sample points to approximate the function accurately?
* Now consider $v(x)=u(x)^2$, using $10$ sample points to approximate this function. Do we have enough sample points to approximate the function accurately? What's the magnitude of the error?

---

### Example 7: De-Aliasing by zero-padding

% grid points
  N = 24; h = 2*pi/N; x = h*(0:N-1)';

% v: the function
  v = exp(sin(x));
  v_hat = fft(v);

% 3/2-rule
  M = ceil((3/2)*N);
  h2 = 2*pi/M; x2 = h2*(0:M-1)';
  v_tilde = zeros(M,1);
  v_tilde(1:N/2) = v_hat(1:N/2);
  v_tilde(N/2+1) = 0.5*v_hat(N/2+1);
  v_tilde(M-N/2+1) = 0.5*v_hat(N/2+1);
  v_tilde(M-N/2+2:M) = v_hat(N/2+2:N);
% 
  v_tilde = v_tilde*M/N;
  vm = real(ifft(v_tilde));
  vm2_tilde = fft(vm.^2);
  
% v2_hat: Fourier modes of v^2
  v2_hat = zeros(N,1);
  v2_hat(1:N/2) = vm2_tilde(1:N/2);
  v2_hat(N/2+2:N) = vm2_tilde(M-N/2+2:M);
  v2_hat = v2_hat*N/M;
  
% ve: the exact solution of (v^2)'
  ve = exp(2*sin(x)).*(2*cos(x));

% weights: Fourier modes for odd deriv.
  kk1 = [0:N/2-1 0 -N/2+1:-1]';

% w2: spectral approximation of (v^2)_x
  w2_hat =  1i*kk1.*v2_hat;
  w2 = real(ifft(w2_hat));

  error = norm(w2-ve,inf);
  disp(['max error for (v^2)_x = ' num2str(error)])

% w1: spectral approximation of v_x
  w1_hat =  1i*kk1.*v_hat;
  w1 = real(ifft(w1_hat));

  error = norm(2*v.*w1-ve,inf);
  disp(['max error for 2*v*v_x = ' num2str(error)])

The result shows

max error for (v^2)_x = 4.4851e-09
max error for 2*v*v_x = 2.883e-12

It looks like the process improves a bit but not much, and it also looks like expand the whole computation as $2u(x)u'(x)$ is much metter.

**But**, let's take another look at the error by reading the Fourier coefficients (continue from the above Matlab code):

% grid points
  M = 60; h2 = 2*pi/M; y = h2*(0:M-1)';

% v2p: (v^2)_x
  v2p = exp(2*sin(y)).*cos(y)*2;
  v2p_hat = fft(v2p);

  semilogy(abs(v2p_hat(1:M/2+1)/M), 'o'); hold on
%
  semilogy(abs(w2_hat(1:N/2+1)/N), '*')
  semilogy(abs(w2_hat(1:N/2+1)/N-v2p_hat(1:N/2+1)/M), '*')
%
  u = 2*v.*w1;
  u_hat = fft(u);
  semilogy(abs(u_hat(1:N/2+1)/N), '+')
  semilogy(abs(u_hat(1:N/2+1)/N-v2p_hat(1:N/2+1)/M), '+')
%
  legend('exact', 'de-alias(DA)', 'exact-DA', 'full expand(FP)', 'exact-FP')

One easily find that the Fourier coefficients for the full-expanded one is much larger comparing to the de-aliased one.