2-Dim wave equation: Finite Differences

it can be expressed using a system of finite differences using the relations between position, velocity and acceleration: and

We consider the case of a 30x30 grid with different coefficients: . The terms of the first equation right hand side can be considered as forces because they modify the velocity.

clearvars;

%-------------------------------------------

% Initialize variables

%-------------------------------------------

dt = 0.01; %time Step

m=30; %files

n=m; %cols

kcx=150; %xkc=dt^2*c^2/dx^2;

kcy=240; %ykc=dt^2*c^2/dy^2;

pos = zeros(m,n,3);

vel = zeros(m,n,3);

forces = zeros(m,n,3);

x=linspace(-2,2,n); %dx=4/m

y=linspace(1,2,m); %dy=1/n

[X,Y] = meshgrid(x,y);

pos(:,:,1)=X;

pos(:,:,2)=Y;

%--------------------

% Time evolution

%--------------------

numIter= 600;

for step=1:numIter

%-----------------------------------------------------

% compute acceleration values (forces) that modify velocity

%-----------------------------------------------------

for i=1:m

for j=1:n

iright = min(j+1,n);

ileft = max(j-1,1);

ibottom = max(i-1,1);

iTop = min(i+1,m);

for k=1:3

forces(i,j,k) = kcx*(pos(i,iright,k)-2*pos(i,j,k)+pos(i,ileft,k))+...

kcy*(pos(ibottom,j,k)-2*pos(i,j,k)+pos(iTop,j,k));

end

end

end

%-----------------------------------------

% Generate a water drop in the center

%-----------------------------------------

if step < m

forces(12,12,3) = 5;

end

%------------------------------------

%Solver: Semi-Implicit Euler Method

%------------------------------------

vel = vel + forces*dt;

pos = pos + vel*dt;

%------------------------------------

% Boundary conditions (fixed borders)

%------------------------------------

pos(1,:,1)=X(1,:);

pos(1,:,2)=Y(1,:);

pos(1,:,3)=0;

pos(end,:,1)=X(end,:);

pos(end,:,2)=Y(end,:);

pos(end,:,3)=0;

pos(:,1,1)=X(:,1);

pos(:,1,2)=Y(:,1);

pos(:,1,3)=0;

pos(:,end,1)=X(:,end);

pos(:,end,2)=Y(:,end);

pos(:,end,3)=0;

%------------------------------------

% Draw the surface

%------------------------------------

surf(pos(:,:,1),pos(:,:,2),pos(:,:,3));

view(-64,36);

axis([-2 2 0.5 2 -0.2 0.2]);

colormap('winter');

drawnow;

%pause(0.005);

end