Matlab Program to Use FDM Scheme Assignment Solution.

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