Can't understand a simple wave equation matlab code

Question

I'm trying to figure out how to draw a wave equation progress in a 2D graph with Matlab.

I found this piece of code which effectively draw a 2D wave placing a droplet in the middle of the graph (I almost fully commented it to simplify things) and then letting it expanding till the border, then bouncing back (how can this code do that? I don't see any code taking the graph dimensions into account) and forth

The code should be pretty simple to someone who knows how a wave equation work, but to me.. it's just dark.

What is the C1 constant? And how does the "un" line succeed into evolving the wave with the convolution? I can't find theory about that nor anything that explains this to me

c=1; % speed of wave;
dx=1; % space step;
dt=0.05; % time step;


szx=200;
szy=200; % size of the drawing area

tm=3000; % time
k=0.002; % decay factor
dsz=3; % droplet size
da=0.07; % droplet amplitude


x=0:dx:(szx-dx); % Generate a x vector from 0 to 200 with step 1
y=0:dx:(szy-dx); % Generate a y vector from 0 to 200 with step 1
t=0:dt:tm; % time, generate a time vector from 0 to 3000 with step 0.05

[X,Y] = meshgrid(x,y); % Create base planes for the graph

Lx=length(x);
Ly=length(y);

u=zeros(Ly,Lx); % initial value (a null matrix for u)
uo=u; % previous = curent => velocties =0

% Initialize the drawing and the axis
close all;
hf=figure;
ha=axes;
hi=imagesc(x,y,u);
set(ha,'clim',[-1 1]); % set colors


D=[0 1 0; 1 -4 1; 0 1 0]; % 2d laplace operator, this will be used in the convolution (probably a laplace 2d operator substitute?)

 % Kdt = decay factor * dt (timestep), that's how much the wave decays per time step
 kdt=k*dt;
 % !!!! this c1 is constant, but I can't figure out what it is
 c1=dt^2*c^2/dx^2;

 % droplet as gaussian, the initial droplet is nothing else but a gaussian in a 2D matrix
 xd=-2*dsz:dx:2*dsz;
 yd=-2*dsz:dx:2*dsz;
 [Xd,Yd] = meshgrid(xd,yd);
 Zd=-da*exp(-(Xd/dsz).^2-(Yd/dsz).^2);


 % Used to count droplets
 One_single_droplet = 0;


 for tt=t
     % !!!!!! Calculate wave equation evolution -> this is exactly where I can't figure out what does this mean
    un=(2-kdt)*u+(kdt-1)*uo+c1*conv2(u,D,'same');

    uo=u; % current become old
    u=un; % new become current

       % Draw the wave updating the graph matrix values (u)
       set(hi,'Cdata',u);
       drawnow;

    % droplets, just one
    if One_single_droplet == 0
        x0d= 100;
        y0d= 100; % droplet center

        % Place the droplet centered on x0d and y0d: adds to the u matrix (the wave function) the gaussian 2d droplet
        u(y0d-2*dsz:y0d+2*dsz,x0d-2*dsz:x0d+2*dsz)=...
            u(y0d-2*dsz:y0d+2*dsz,x0d-2*dsz:x0d+2*dsz)+Zd;

        One_single_droplet = 1;
    end

 end

Answer 1

It is appears to be a standard finite difference discretization via method of lines of the wave equation

$u_{tt} + c^2 \nabla^2u = 0, \quad (x,y)\in[0,L_x]\times[0,L_y], \quad t\geq0$

Standard finite difference approximations can be used to formulate a discrete problem by introducing the formally second-order accurate approximations

$u_{tt}\approx \frac{u^{n+1}_{ij}-2 u_{ij}^n + u^{n-1}_{ij}}{\Delta t^2}$

and

$\nabla^2u = u_{xx} + u_{yy} \approx \frac{u_{i+1,j}^n - 2u_{ij}^n + u_{i-1,j}}{\Delta x^2} + \frac{u_{i,j+1}^n - 2u_{ij}^n + u_{i,j-1}}{\Delta y^2}$

where $u_{ij}^n = u(x_i,y_j,t_n)$ expressed the small-amplitude surface elevation at node $(x_i,y_j)$ at time $t_n$. If a uniform grid $x_i=hi$, $i=0,1,...,N_x$ and $y_j=jh$, $j=0,1,...,N_y$ with $h=\tfrac{L_x}{N_x}=\tfrac{L_y}{N_y}$ is introduced and with time defined as $t_n = n\Delta t$, $n=0,1,...$, then these approximations can be substituted for the continuous derivatives and then rewritten for the discrete equation to be expressed in the form

$u_{ij}^{n+1} = 2u_{ij}^n - u_{ij}^{n-1} + c^2\frac{\Delta t^2}{h^2}(u_{i-1,j}^n + u_{i+1,j}^n - 4 u_{ij}^n + u_{i,j-1}^n + u_{i,j+1}^n)$

which provides the basis for locally (explicitly) updating the solution for each grid point in the spatial domain. Remark, that the update requires storage of solutions at two previous steps in time.

This is a discrete scheme for the wave equation without the damping term of the form $−ku_t$ included on the right hand side. The discretization of $u_t$ can be done explicitly using a standard backward difference operations.

Answer 2

We could start with the equations and explain the code, but I'll go the other way because it's a useful methodology for figuring out what code is doing. First let's consider

conv2(u,D,'same');

which applies the difference stencil D to u. If you are not familiar with conv2(), the operation (expressed using $0$-based indexing for D) is

$$ u\star D = \sum_{k=-1}^1 \ \sum_{l=-1}^1 u_{i+k,j+l} D_{k,l} $$

For constant stencils, this representation is much more efficient than storing a sparse matrix or writing your own code to apply the difference stencil. You can experiment easily in 1D to understand conv and conv2.

> conv([0,0,1,0,0], [1,2,3], 'same')
ans =

   0   1   2   3   0

This difference stencil is multiplied by c1 = c^2*dt^2/dx^2 which happens to be the CFL number squared. For the moment, we can observe that conv2(u,D,'same')/dx^2 is the centered difference approximation to the Laplacian $\nabla^2 u$. Let's call this sub-expression Lap_u.

Now for the rest of the expression

un = (2-kdt)*u + (kdt-1)*uo + c^2*dt^2*Lap_u;

Let's abuse syntax by moving the part not depending on kdt = k*dt to the left side

(un - 2*u + uo)/dt^2 = -k*(u-uo)/dt + c^2 * Lap_u;

The left side is a discrete centered difference approximation to $u_{tt}$ and (u-uo)/dt is a first-order approximation to $u_t$, so the overall expression is approximating

$$ u_{tt} = -k u_t + c^2 \nabla^2 u . $$

This is the Leapfrog method for the wave propagation term and forward Euler for diffusion. This becomes more clear by writing the equation as a first order system

$$ \begin{align} u_t &= v \\ v_t &= -k v + c^2 \nabla^2 u \end{align} $$

where the discrete scheme at time index $n$ is

$$ \begin{align} \frac{u^n - u^{n-1}}{\Delta t} &= v^{n-1/2} \\ \frac{v^{n+1/2} - v^{n - 1/2}}{\Delta t} &= -k v^{n - 1/2} + c^2 \nabla^2 u^{n} . \end{align} $$

The Leapfrog method is symplectic (energy conserving in this case), but the even and odd time levels $n$ are decoupled, so some decay (or other filtering) is necessary to prevent the solution from oscillating in time. You should be able to observe this in the code by taking $k \to 0$.