2
$\begingroup$

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

enter image description here

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
$\endgroup$
8
$\begingroup$

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;
$\endgroup$
  • $\begingroup$ Thank you for such an explanatory answer. I am not understanding what you mean in your second bullet point. If you have time could you elaborate or reference me to a site or pdf. Thanks $\endgroup$ – DoubleOseven Apr 5 '16 at 17:16
  • $\begingroup$ @DoubleOseven Perhaps something like this book? people.maths.ox.ac.uk/trefethen/spectral.html $\endgroup$ – Kirill Apr 5 '16 at 17:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.