Numerical solution of $2D$ wave equation using Fourier transform and finite differences

Question

This is the $2$-dimensional wave equation

$$ u_{tt} = u_{xx} + u_{yy} $$

with initial condition $u(x,y,0)=f(x,y)$ and $u_{t}(x,y,0) = 0$.

The inverse Fourier transform used is

$$ u(x,y,t) = \iint \hat{u}(\omega_{x}, \omega_{y}, t)e^{(2\pi) \omega_{x}i x} e^{(2\pi) \omega_{y} i y}d\omega_{x} d\omega_{y} $$

this is also the version that is used by Matlab's FFT. Applying this to wave equation we get

$$ \iint \hat{u}_{tt}e^{(2\pi) \omega_{x}i x} e^{(2\pi) \omega_{y}i y}d\omega_{x} d\omega_{y} = (2 \pi)^{2} (\omega_{x}^{2} + \omega_{y}^{2}) \iint \hat{u} e^{(2\pi) \omega_{x} i x} e^{(2\pi) \omega_{y} i y}d\omega_{x} d\omega_{y} $$ so the Wave equation in frequency space is:

$$ \hat{u}_{tt} = -(2\pi)^{2} (\omega_{x}^{2} + \omega_{y}^{2}) \hat{u} $$

Although this gives exact solution very easy, I tried to solve it numerically using Matlab.

The numerical method I used is finite difference:

$$ \hat{u}_{tt} \approx \frac{\hat{u}_{t}(t) - \hat{u}_{t}(t-\Delta t)}{ \Delta t} $$ $$ \approx \frac{\hat{u}(t) - 2 \hat{u}(t- \Delta t) + \hat{u}(t - 2\Delta t)}{ (\Delta t)^{2}} $$

And for the initial condition: $\hat{u}(\Delta t) = \hat{u}(0)$.

I did it successfully in 1D case, $\hat{u}(\omega,t)_{tt} = (2\pi \omega)^{2}\hat{u}(\omega,t)$, using similar formulation. But I get something wrong in the solution in 2D case, the solution does not blow up but the behavior does not match the exact solution (using Gaussian curve as initial condition, then the wave should split with it's height equals half of the initial height, instead it decays drastically). See animation below:

---

Matlab code

clear;
tic;

dx = 0.005; dy = 0.005;
xmax = 1; xmin = -1;
ymax = 1; ymin = -1;
periodx = xmax-xmin; periody = ymax-ymin;
x = [(xmin):dx:(xmax)];
y = [(ymin):dy:(ymax)];
nx = length(x);
ny = length(y);
multiplierx = (1)/periodx;
multipliery = (1)/periody;
[X, Y] = meshgrid(x,y);

if (mod(nx,2) == 0)
    w_right = [0:1:((nx/2)-1)];
    w_left = [(-nx/2):1:-1];
    w = [w_right, w_left];
    wx = multiplierx*w;
else  
    w_right = [0:1:(((nx-1)/2))];
    w_left = [(-(nx-1)/2):1:-1];
    w = [w_right, w_left];
    wx = multiplierx*w;
end

if (mod(ny,2) == 0)
w_right = [0:1:((ny/2)-1)];
w_left = [(-ny/2):1:-1];
w = [w_right, w_left];
wy = multipliery*w;
else  
w_right = [0:1:(((ny-1)/2))];
w_left = [(-(ny-1)/2):1:-1];
w = [w_right, w_left];
wy = multipliery*w;
end

[Wx, Wy] = meshgrid(wx, wy);

f = @(x,y) 0.5*exp(-(x.^2 + y.^2)/0.1);
Z = f(X,Y);

dt = 0.005;
t = [0:dt:2]; nt = length(t);

u(:,:, 1) = f(X, Y);
u(:,:, 2) = f(X, Y);

uhat(:,:,1) = fft2(u(:,:, 1));
uhat(:,:,2) = fft2(u(:,:, 2));


surfl(u(:,:,1));
zlim([-1,1]);
view([0, 0]);
colormap(winter);
xlabel('x'); ylabel('y');
shading interp;

for i = [3:1:nt]

  uhat(:,:,i) = (uhat(:,:,i-1) + (uhat(:, :, i-1) - uhat(:, :, i-2)))./(1 + (((dt*2*pi)^2)*(Wx.^2 + Wy.^2)));
  u(:,:,i) = real(ifft2(uhat(:,:,i)));
  uhat(:,:,i) = fft2(u(:,:,i));
  disp(i);  
end 

toc;

for i = [1:1:nt]  
  surfl(u(:,:,i));
  zlim([-1,1]);
  view([0, 0]);
  colormap(winter);
  shading interp;
  pause(0.01)
end