Wave propagation simulation with Matlab

Multiple sources inside a cavity

Given a bounded cavity, four sources are emitting together and synchronously. Interference phenomenom, that is wave superpositions resulting in a global wave either of greater amplitude (contructive interference) or same or lower amplitude (destructive interference, waves anihilate each others).
See and observe the following animation:

Simulation and computationnal aspects rely on finite difference method, see for example (french) this page and also for example (also french) introduction to Matlab programming.
Computation, simulation, and animation is performed using Matlab, from the following script:

% New beginning... clean...
clear all;close all
tic;% Starting clock

% wave propagation speed
c=8;

% Computation grid
Lx=100;Ly=100;
T=20;

deltax=Lx/150;
deltay=Ly/150;
deltat=sqrt(deltax^2+deltay^2)/(2*c);

gax=(c*deltat)^2/(deltax^2);
gay=(c*deltat)^2/(deltay^2);

xx=[0:deltax:Lx];nx=length(xx);
yy=[0:deltay:Ly];ny=length(yy);
tt=[0:deltat:T];nt=length(tt);

Reflection_parois=1; % 0 for absorbing boundary

sx=round(nx/4);sy=round(ny/4);
sx2=nx-sx;sy2=sy;
sx3=sx;sy3=ny-sy;
sx4=sx2;sy4=sy3;
sx=[sx sx2 sx3 sx4];sy=[sy sy2 sy3 sy4];


sfrq=1; % Source frequency


Mat=zeros(nx,ny,nt);
% Mat Indexes, without boundaries
i=[2:nx-1];j=[2:ny-1];

Mat(i,j,2)=Mat(i,j,1)+...
    gax/2*(Mat(i-1,j,1)+Mat(i+1,j,1)-2*Mat(i,j,1))+...
    gay/2*(Mat(i,j-1,1)+Mat(i,j+1,1)-2*Mat(i,j,1));

for k=3:nt-1
    Mat(i,j,k+1)=2*Mat(i,j,k)-Mat(i,j,k-1)+...
                gax*(Mat(i-1,j,k)+Mat(i+1,j,k)-2*Mat(i,j,k))+...
                gay*(Mat(i,j-1,k)+Mat(i,j+1,k)-2*Mat(i,j,k));
    
    Mat(sx,sy,k+1)=sin(2*pi*sfrq*k*deltat);
    
    if (Reflection_parois==0)
        Mat(1,:,k+1)=Mat(2,:,k);
        Mat(nx,:,k+1)=Mat(nx-1,:,k);
        Mat(:,1,k+1)=Mat(:,2,k);
        Mat(:,ny,k+1)=Mat(:,ny-1,k);
    end
    
end

disp(' ')
disp(['Fin du calcul en ' num2str(toc) ' s'])


ff=figure(1);clf;whitebg('w');
Film=[];
for k=1:1:nt
   pcolor(squeeze(Mat(:,:,k)));
   axis off;shading interp;caxis([-0.5 0.5]);
   Film=[Film getframe(ff)];
   %pause(0.1)
end

movie2avi(Film,'Film.avi')