Can't understand a simple wave equation matlab code

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

-
This question is more suited to scicomp.stackexchange.com or maybe even SO. – Inquest Feb 11 '12 at 6:50

1 Answer

This is a discretized version of the wave equation that is used by iterating. So, let us consider

$$\frac{1}{c^2}\frac{\partial^2u}{\partial t^2}-\frac{\partial^2u}{\partial x^2}-\frac{\partial^2u}{\partial y^2}=0.$$

Now, you can discretize the derivatives in the following way

$$\frac{\partial u}{\partial t}=\frac{u(i,k,l)-u(i-1,k,l)}{\Delta t}$$

$$\frac{\partial u}{\partial x}=\frac{u(i,k,l)-u(i,k-1,l)}{\Delta x}$$

$$\frac{\partial u}{\partial y}=\frac{u(i,k,l)-u(i,k,l-1)}{\Delta x}$$

Then just repeat them to get second derivatives and you will get

$$\frac{\partial^2 u}{\partial t^2}=\frac{u(i+1,k,l)-2u(i,k,l)+u(i-1,k,l)}{(\Delta t)^2}$$

$$\frac{\partial^2 u}{\partial x^2}=\frac{u(i,k+1,l)-2u(i,k,l)+u(i,k-1,l)}{(\Delta x)^2}$$

$$\frac{\partial^2 u}{\partial y^2}=\frac{u(i,k,l+1)-2u(i,k,l)+u(i,k,l-1)}{(\Delta x)^2}.$$

Now, you can write an iteration scheme to solve your equation by substitution of these discrete derivatives as

$$\frac{u(i+1,k,l)-2u(i,k,l)+u(i-1,k,l)}{c^2(\Delta t)^2}=\frac{u(i,k+1,l)-2u(i,k,l)+u(i,k-1,l)}{(\Delta x)^2}+\frac{u(i,k+1,l)-2u(i,k,l)+u(i,k-1,l)}{(\Delta x)^2}$$

and so

$$u(i+1,k,l)=2u(i,k,l)-u(i-1,k,l)+c_1\left[u(i,k+1,l)-4u(i,k,l)+u(i,k-1,l)+u(i,k+1,l)+u(i,k-1,l)\right].$$

As you see from this scheme, the discretized Laplace operator has the coefficients of the matrix D in the matlab code and $c_1$ is exactly the one used in the code. Laplace operator is computed through conv2 that multiplies two matrices as it should. But note that a decay term has been added (those multiplied by k) and you must use for $u(0,k,l)$ the initial Gaussian wave-packet to start iteration.

This matlab code gives a fine example of application of an iteration scheme to solve a pde.

-
You're right! The 2-kdt and kdt+1 terms seems to balance the compromise between gain/decay, the rest is finally clear. Thank you very much – Paul Feb 11 '12 at 22:12