Solve Poisson's Equation in Matlab

Solution Poisson's Equation In a Rectangular Domain Using 5-Point Stencil

Poisson's equation in a square domain $0 \leq x,y \leq 1$ with source term given by $q(x)=100*\sin(\pi x)*\sin(\pi y)$. The boundaries of the domain are maintained at $u=0$. Determine the solution using finite difference method.

Solution:

Discretizing domain uniformly with step size h and k along x and y direction. Using five point stencil, equation at each node is given by $$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}+\frac{u_{i,j+1}-2u{i,j}+u{i,j-1}}{k^{2}}=-q(i,j)$$ $$u_{i,j}=\frac{1}{\frac{2}{h^{2}}+\frac{2}{k^{2}}}(\frac{u_{i+1,j}+u_{i-1,j}}{h^{2}}+\frac{u_{i,j+1}+u{i,j-1}}{k^{2}}+q(i,j))$$ where $q(i,j)=100*\sin(\pi x_{i})*\sin(\pi y_{j})$ Here a MATLAB code is given below. Choosen $m=n=81$, these are the grid spacing.

clear;
close;
clc;

n = 81; m = 81;
x = 0:1/(n-1):1;
y = 0:1/(m-1):1;
h = 1/(n-1); k = 1/(m-1);
[Y,X] = meshgrid(y,x);
q = 100*sin(pi*X).*sin(pi*Y);
u = zeros(n);
u0 = u;
residual = 1;
% Gauss-Seidel Method
while residual>1e-6
    for i=2:n-1
        for j=2:n-1
            u(i,j) = ((u(i+1,j) + u(i-1,j))/h^2 + ...
                (u(i,j+1) + u(i,j-1))/k^2 + q(i,j))/(2/h^2 + 2/k^2);
        end
    end
    residual = max(max(abs(u0-u)));
    u0=u;
end
subplot(2,2,1)
contour(X,Y,u0,"ShowText","on");
xlabel("X")
ylabel("Y")
title("Contour Plot")
subplot(2,2,2)
surf(X,Y,u0);
xlabel("X")
ylabel("Y")
zlabel("u(X,Y)")
title("Surface Plot")
subplot(2,2,3)
% map = [0 0 0.3
%     0 0 0.4
%     0 0 0.5
%     0 0 0.6
%     0 0 0.8
%     0 0 1.0];
% colormap(map)
%set(gcf,'Position',[100 100 500 500])
imagesc(u0);
colorbar
drawnow
xlabel("X")
ylabel("Y")
title("Filled Contour Plot")

Published with MATLAB® R2020b

• Matlab Official website: Visit Matlab.com
• Matlab Documentation Center: Visit Matlab Documentation Center