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.
cme_assignment8_b 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")
- Matlab Official website: Visit Matlab.com
- Matlab Documentation Center: Visit Matlab Documentation Center
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")
- Matlab Official website: Visit Matlab.com
- Matlab Documentation Center: Visit Matlab Documentation Center
Comments
Post a Comment
If you have any doubts, please let me know!