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

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');

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);
pause(0.01)
end

• Why should the solution split in half? That's a 1d effect. Your visualization shows the 2d surface on which the wave travels edge on. Rotate the figure a bit and you will probably see how you have a wave that is traveling radially outward, as one should expect. Commented May 22, 2022 at 17:02
• @WolfgangBangerth i know it travels radially outward. But the dynamics of the shape does not behave similarly like 1D? Commented May 22, 2022 at 23:45
• @WolfgangBangerth even in 2D, as far as I remembered it should be spreading with the height equals half of the initial. Commented May 23, 2022 at 9:21
• No, the height must decrease with time and distance in 2d because the circumference of the wave increases with time and distance. Commented May 23, 2022 at 15:50

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!

• But the code is correct. The mistake is only in my writing in this site. Commented May 22, 2022 at 10:20