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