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:

enter image description here

Matlab code


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;
    w_right = [0:1:(((nx-1)/2))];
    w_left = [(-(nx-1)/2):1:-1];
    w = [w_right, w_left];
    wx = multiplierx*w;

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;
w_right = [0:1:(((ny-1)/2))];
w_left = [(-(ny-1)/2):1:-1];
w = [w_right, w_left];
wy = multipliery*w;

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

view([0, 0]);
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));


for i = [1:1:nt]  
  view([0, 0]);
  shading interp;
  • 2
    $\begingroup$ 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. $\endgroup$ Commented May 22, 2022 at 17:02
  • $\begingroup$ @WolfgangBangerth i know it travels radially outward. But the dynamics of the shape does not behave similarly like 1D? $\endgroup$
    – Redsbefall
    Commented May 22, 2022 at 23:45
  • $\begingroup$ @WolfgangBangerth even in 2D, as far as I remembered it should be spreading with the height equals half of the initial. $\endgroup$
    – Redsbefall
    Commented May 23, 2022 at 9:21
  • 1
    $\begingroup$ No, the height must decrease with time and distance in 2d because the circumference of the wave increases with time and distance. $\endgroup$ Commented May 23, 2022 at 15:50

1 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!

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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