I am trying to solve the 2D Wave Equation numerically and plot the results.

Using the finite difference method, I have solved the 2D wave equation to get the following equation. When I try to plot the results, nothing happens.

The source function is working but plotting p(x,y,t) does not work. Here is the code:

%Numerically solve 2D Wave Equation
clear
clc
close all

%define constants
A = 2;
Omega = 2*pi/2.5;
StandardDev = 1;
speed = 1;

Nt = 100; %Number of temporal steps
Nx = 100; %Number of steps in x direction
Ny = 100; %Number of steps in y direction

tMin = 0;
tMax = 10;

xMin = 0;
xMax = 10;

xNot = xMax/2;

yMin = 0;
yMax = 10;

yNot = yMax/2;

t = linspace(tMin, tMax, Nt);
x = linspace(xMin, xMax, Nx);
y = linspace(yMin, yMax, Ny);

%Compute secondary parameters
DeltaT = t(2) - t(1);
DeltaX = x(2) - x(1);
DeltaY = y(2) - y(1);

%%Numerically solve PDE
%Initialize containers to hold all P values
f = zeros(Nx, Ny, Nt);
p = zeros(Nx, Ny, Nt);

%Set Initial Conditions
p(:,:,1) = 0;
p(:,:,2) = 0;

%Set Boundary Conditions
p(1,:,:) = 0;
p(Nx,:,:) = 0;
p(:,1,:) = 0;
p(:,Ny,:) = 0;

%Iterate through t values
for l=1:Nt
    %Iterate through x values
    for i=1:Nx
        %Iterate through y values
        for j=3:Ny
            f(i,j,l) = A*sin(Omega*t(l))*exp(-1/(2*StandardDev^2)*((x(i) - xNot)^2 + (y(j) - yNot)^2));
        end
    end
end

for l=3:Nt
    %Iterate through x values
    for i=2:Nx
        %Iterate through y values
        for j=2:Ny
            p(i,j,l) = 2*p(i,j,l-1) - p(i,j,l-2) + DeltaT^2*(f(i,j,l) + speed^2*((p(i+1,j,l-1) - 2*p(i,j,l-1) + p(i-1,j,l-1))/deltaX^2 + (p(i,j+1,l-1) - 2*p(i,j,l-1) + p(i,j-1,l-1))/deltaY^2));
        end
    end
end

%Plot the results
figure
for t=1:100
    hold on
    clf
    surf(x, y, p(:,:,t))
    zlim([-2 2])
    colormap winter(1)
    drawnow
end
hold off

Answer 1

[Nikhil Sonavane]

There are two spell errors in the following lines of code-

 p(i,j,l) = 2*p(i,j,l-1) - p(i,j,l-2) + DeltaT^2*(f(i,j,l) + speed^2*((p(i+1,j,l-1) - 2*p(i,j,l-1) + p(i-1,j,l-1))/deltaX^2 + (p(i,j+1,l-1) - 2*p(i,j,l-1) + p(i,j-1,l-1))/deltaY^2));

deltaX should be DeltaX and similarly for deltaY as well. Even after correcting this there are some indexing errors in this line of code. In my suggestion you should correct this logic before moving forward.