2-Dim wave equation: Finite Differences

Consider the 2-Dim wave equation
it can be expressed using a system of finite differences using the relations between position, velocity and acceleration: and
The two equations can be written as (semi-implicit Euler Method):
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
© Numerical Factory 2020