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