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

Answer 1

Correction in the expression

It appears that complex iota $i$ has not been included in the exponents in the expression for the inverse Fourier transform. The correct expression is:

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

This would result in the following wave equation in frequency space, which is a second order ODE:

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

Numerical methodology

Instead of directly discretizing the second order derivative in time, it would be numerically wise to convert the second order ODE into a system of first order ODEs beforehand as:

Take $\hat{v} = \hat{u}_t$. This leads to the the system of first order ODEs as follows

\begin{align} \hat{u}_t &= \hat{v} \\ \hat{v}_t &= -(2\pi)^{2} (\omega_{x}^{2} + \omega_{y}^{2}) \hat{u} \end{align}

The above system can be solved for $\hat{u}$ and $\hat{v}$ subject to the initial conditions \begin{align} \hat{u}(\omega_x,\omega_y,0) &= \iint f(x,y) e^{-i(2\pi) \omega_{x}x} e^{-i(2\pi) \omega_{y}y}dx dy \\ \hat{v}(\omega_x,\omega_y,0) &= 0 \end{align}

In matrix form,the system could be written as \begin{align} \begin{bmatrix} \hat{u}_t \\ \hat{v}_t \end{bmatrix} &= \underbrace{\begin{bmatrix} 0 & 1 \\ -(2\pi)^{2} (\omega_{x}^{2} + \omega_{y}^{2}) & 0 \end{bmatrix}}_{\mathbf{A}} \underbrace{\begin{bmatrix} \hat{u} \\ \hat{v} \end{bmatrix}}_{\mathbf{\hat{U}}} \\ \implies \mathbf{\hat{U}}_t & = \mathbf{A} \mathbf{\hat{U}} \\ & \text{where } \quad \mathbf{\hat{U}}(\omega_x, \omega_y, 0) = \begin{bmatrix} \hat{u}(\omega_x,\omega_y,0) \\ \hat{v}(\omega_x,\omega_y,0) \end{bmatrix} \end{align}

You could use any standard finite difference scheme to discretize the first order derivatives in $\hat{u}$ and $\hat{v}$, such as forward Euler scheme, Runge-Kutta scheme, as per the desired accuracy.

For example, if you use forward Euler scheme, the first order derivative would be approximated as \begin{align} \mathbf{\hat{U}}_t \approx \frac{\mathbf{\hat{U}}(t+\Delta t) - \mathbf{\hat{U}}(t)}{\Delta t} \end{align}

Hope this helps!