61 views (last 30 days)

Show older comments

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.

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

Start Hunting!