Numerically computing the advection equation

Question

I am trying to write a program to compute the advection equation.

$$u_t +u_x = 0$$

I use the spectral method for the spatial derivative $u_x$ and the leapfrog method for the time derivative $u_t$.

That will give the an equation for the next time step:

$$u_{j}^{k+1} = u_{j}^{k-1} - \frac{\Delta t}{n}[f^{-1}(iN \pi f(U_{j}^{k}))]$$

where

$N = -n,...,n-1$, $f = $ fourier transform, $f^{-1} = $ inverse fourier transform, $n$ = number of discretized points.

Because I need two initial values (vectors) I use ode45 to compute the solution $u_{j}^{k+1}$ at time $\Delta t$ and having already the solution $u_{j}^{k}$ at $t = 0$ we can use those two vectors to start the algorithm.

I get wrong results, but I can't see what I am doing wrong. This is a plot

Maybe someone can see what I am doing wrong.

CODE:

w = 0.2;
A = 1;
n = 2^8;
deltat = 2*A/n;
h = 2*A/n;   %step size for N
x = -A:h:A-h; %Discretization of N
tm = 0:deltat:2;
u = @(x) exp(-(x/w).^2);
N = -n/2:n/2-1;   

sol = zeros(length(tm),length(x)); % rows are space / columns are time;

u0 = u(x); % initial condition at t = 0
[t, ksol] = ode45(@f,[0 deltat], u0);
u1 = ksol(end,:);
sol(2,:) = u1;
sol(1,:) = u0;

for k = 3:length(tm)+1
    f_hat = fft(sol(k-1,:),n);    % fft transform
    umiddle = deltat/n * ifft(1i * pi * N * f_hat); %% inverse fft
    sol(k,:) = sol(k-2,:) - umiddle;

end

mesh(abs(real(sol(2:end,:))))

Function:

  function sol = f(~,u)
    n = 2^8;
    N = -n/2:n/2-1;   
    D = diag(1i*N*pi); % diagonal matrix (* each u0 element)
    sol = -ifft(D*fft(u,n));
    end

Answer 1

I see several issues:

• The DFT computed with fft puts the zero mode at the beginning of the array, and if you want to compute the derivative, it is necessary to apply fftshift/ifftshift to the array N to make sure the derivative is correct. It is easy to see for yourself what the correct expression is by working it out with pen and paper, and see also the documentation for fftshift.
• You divide the derivative term by $n$, but I don't see why this is correct: the highest frequency mode in the solution on a mesh of size $N$ will be $u(x) = \cos(\pi x N/(2A))$, which has derivative of order $N$. Yet the way you compute the derivative, the highest frequency mode will have a derivative $\leq 1$, which doesn't seem right.
• I think you also lost a factor of two in the derivative term, because the approximation to $u_t$ should be $(u^{k+1}-u^{k-1})/(2\delta t)$.
• Your $\delta t$ is too large. If we analyze the stability of your scheme in the von Neumann sense, substituting the eigenfunctions $u(t,x) = e^{\mathrm{i}(\lambda t-\mu x)}$, the dispersion relation is $\mu\,\delta t = \sin(\lambda\,\delta t)$. The highest frequency mode on a mesh on $[-A,A]$ of $N$ points is $(-1)^j = \cos(\pi N\frac{(x/A)+1}{2})$, so $\delta t$ should be at most $2A/(\pi N)$ for stability.

Here are the changes, and this gives me the expected result of a wave propagating at speed 1.

w = 0.2;
A = 1;
n = 2^8;
deltat = 2*A/(pi*1.001*n); %%% Try changing 1.001 to 0.999
h = 2*A/n;   %step size for N
x = -A:h:A-h; %Discretization of N
tm = 0:deltat:2;
u = @(x) exp(-(x/w).^2);
N = -n/2:n/2-1;

sol = zeros(length(tm),length(x)); % rows are space / columns are time;

u0 = u(x); % initial condition at t = 0
%[t, ksol] = ode45(@f,[0 deltat], u0);
%u1 = ksol(end,:);
sol(1,:) = u0;
u1 = u0 + deltat * ifft(-1i*pi*ifftshift(N).*fft(sol(1,:)));
sol(2,:) = u1;

for k = 3:length(tm)+1
    f_hat = fft(sol(k-1,:));    % fft transform
    umiddle = ifft(-1i * pi * ifftshift(N) .* f_hat); %% inverse fft
    sol(k,:) = sol(k-2,:) + 2*deltat * umiddle;
    u1  = sol(k-1,:);
end

[meshX,meshT] = meshgrid(x, tm);
figure;
pcolor(meshT, meshX, real(sol(2:end,:)));
shading flat;