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

4 views (last 30 days)
Richard Batelaan on 4 Jul 2020
Edited: Nikhil Sonavane on 14 Jul 2020
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

Nikhil Sonavane on 14 Jul 2020
Edited: Nikhil Sonavane on 14 Jul 2020
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.

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!