Pseudo-Spectral Methods

I was previously affiliated with the University of Toronto, where I wrote some pages for the math department's wiki on how to get pseudo-spectral algorithms up and running for nonlinear wave equations such as Nonlinear Schrödinger (NLS) and Korteweg - de Vries (KdV).  Much of this was based on Trefethen's Spectral Methods in MATLAB.

While that's still an excellent book that I recommend to people who are starting out, it has a very serious flaw in its discussion of nonlinear problems.  Trefethen's p27.m example, which simulates KdV solitons does not include a de-aliasing step.  Early in my career, this led me to think that de-aliasing was unnecessary.  That was, until, I started looking at problems with cubic nonlinearities; the p27.m example is only quadratic.

Consider focusing, cubic NLS on the torus, \(\mathbb{T}\),
\[
i u_t + u_{xx} + |u|^2 u = 0, \quad u: \mathbb{T} \times \mathbb{R} \to \mathbb{C}.
\]
When posed on the whole real line, instead of the torus, it admits soliton solutions,
\[
\Phi(x,t) = e^{i \lambda^2 t} \sqrt{2} \lambda {\rm sech}(\lambda x), \quad \lambda \neq 0.
\]
Since they are exponentially localized, if we take \( \lambda \) sufficiently large, we can approximately simulate it on the finite domain size.

Modifying p27.m for NLS, we have:

% Set discretization parameters
N = 512;
dt = .1/N;
x = (2*pi/N)*(-N/2:N/2-1)';
% Set initial condition
lam = 10;
u = sqrt(2) * lam * sech(lam * x);
v = fft(u);

% Construct integration operators
k = [0:N/2-1 0 -N/2+1:-1]';
E = exp(-.5i * dt * k.^2);
E2 = E.^2;
g = 1i * dt;

% Set simulation parameters
tmax = 5; nplt = floor((tmax/50)/dt); nmax = round(tmax/dt);
udata = u; tdata = 0;

% Integrate using the Integrating Factor method
for n = 1:nmax
  t = n*dt;
  a = g.*fft(abs( ifft(     v   ) ).^2 .* ifft(v));
  b = g.*fft(abs( ifft(E.*(v+a/2)) ).^2 .* ifft(E.*(v+a/2)) );
  c = g.*fft(abs( ifft(E.*v + b/2) ).^2 .* ifft(E.*v + b/2) );
  d = g.*fft(abs( ifft(E2.*v+E.*c) ).^2 .* ifft(E2.*v+E.*c));
  v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
  if mod(n,nplt) == 0
    u = ifft(v);
    udata = [udata u]; tdata = [tdata t];
  end
end

% Plot and save
figure;
for n = 1:length(tdata)
  plot(x, abs(udata(:,n)),...
    x, real(udata(:,n)), x, imag(udata(:,n)),...
    x,abs(udata(:,1)), 'linewidth',2);
  xlim([-pi, pi]);
  ylim([-16,16]);

  set(gca,'fontsize',16);
  xlabel('x');
  legend('|u|', 'Re. u', 'Im. u', '|u_0|');
  title(sprintf(' t = %g', tdata(n)));
  drawnow;
  saveas(gcf,sprintf('frame%i.png', n));
end

This generates the following output:

Notice the appearance of very high frequency oscillations.  These are aliasing errors. De-aliasing them, by zeroing out the upper half of the frequencies in Fourier, produces the following, desired, output:

This is the result of the following code.  Notice that I doubled the number of grid points to compensate for the de-aliasing.

% Set discretization parameters
N = 1024;
dt = .1/N;
x = (2*pi/N)*(-N/2:N/2-1)';
% Set initial condition
lam = 10;
u = sqrt(2) * lam * sech(lam * x);
v = fft(u);
v(N/4+2:end - N/4+1) = 0;

% Construct integration operators
k = [0:N/2-1 0 -N/2+1:-1]';
E = exp(-.5i * dt * k.^2);
E2 = E.^2;
g = 1i * dt;

% Set simulation parameters
tmax = 5; nplt = floor((tmax/50)/dt); nmax = round(tmax/dt);
udata = u; tdata = 0;

% Integrate using the Integrating Factor method
for n = 1:nmax
  t = n*dt;
  a = g.*fft(abs( ifft(     v   ) ).^2 .* ifft(v));
  a(N/4+2:end - N/4+1) = 0;
  b = g.*fft(abs( ifft(E.*(v+a/2)) ).^2 .* ifft(E.*(v+a/2)) );
  b(N/4+2:end - N/4+1) = 0;
  c = g.*fft(abs( ifft(E.*v + b/2) ).^2 .* ifft(E.*v + b/2) );
  c(N/4+2:end - N/4+1) = 0;
  d = g.*fft(abs( ifft(E2.*v+E.*c) ).^2 .* ifft(E2.*v+E.*c));
  d(N/4+2:end - N/4+1) = 0;
  v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
  if mod(n,nplt) == 0
    u = ifft(v);
    udata = [udata u]; tdata = [tdata t];
  end
end

% Plot and save
figure;
for n = 1:length(tdata)
  plot(x, abs(udata(:,n)),...
    x, real(udata(:,n)), x, imag(udata(:,n))...
    , 'linewidth',2);
  xlim([-pi, pi]);
  ylim([-16,16]);

  set(gca,'fontsize',16);
  xlabel('x');
  legend('|u|', 'Re. u', 'Im. u');
  title(sprintf(' t = %g', tdata(n)));
  drawnow;
  saveas(gcf,sprintf('frame%i.png', n));
end

While I still recommend Trefethen's book to people, I would supplement it with another text that addresses the aliasing problem. Two books I like are Boyd's, Chebyshev and Fourier Spectral Methods, and Canuto et al's, Spectral Methods: Fundamentals in Single Domains.