Instructions
Requirements and Specifications
Source Code
clc, clear all, close all
%% Parameters
n = 199;
p = 199;
T = 20;
a = 50;
b = 50;
beta = 4;
Nx = n+1;
Ny = p+1;
m = 2000;
dt = T/m;
dx = a/(n+1);
dy = b/(p+1);
%% Create matrix U that will store the solution
U = zeros(n+1, p+1, m);
% Boundary conditions
U(1,:,:) = 0;
U(n+1,:,:) = 0;
U(:,1,:) = 0;
U(:,p+1,:) = 0;
%% Initial conditions
% Left side of the pyramid
for i = 1:Nx/2
i = round(i);
j = round(1 + ((Ny-2)/(Nx-2))*(i-1));
val = 2/Nx *i;
U(i,j,:) = val;
end
% Right side of the pyramid
for i = Nx/2:Nx
i = round(i);
j = round(Ny/2 + ((2-Ny)/Nx)*(i - Nx/2));
val = 2-2*i/Nx;
U(i,j,:) = val;
end
%%
%% Solve
for k = 2:m-1
for i = 2:n
for j = 2:p
dUx = (U(i+1,j,k) - 2*U(i,j,k) + U(i-1,j,k))/dx^2;
dUy = (U(i,j+1,k) - 2*U(i,j,k) + U(i,j-1,k))/dy^2;
U(i,j,k+1) = 2*U(i,j,k) - U(i,j,k-1) + (beta*dt^2)*(dUx + dUy);
end
end
end
%% Now create animated plot
figure
x = linspace(0, a, n+1);
y = linspace(0, b, p+1);
for k = 1:m
if mod(k,15) == 0
ti = (k-1)*dt;
M = squeeze(U(:,:,k)');
surf(x,y,M), view(2);
shading interp
title(sprintf("t = %.2f s", ti));
xlabel('X')
ylabel('Y')
zlim([-1,1])
colorbar
drawnow
end
end