Solve Poisson Equation (Dirichlet boundary condition) via Jacobi and Gauss-Seidel-Iteration
65 views (last 30 days)
Show older comments
Hello everyone,
I would like to solve the Poisson Equation with Dirichlet boundary condition in Matlab with the Jacobi- and the Gauss-Seidel Iteration. After I completed running the iterations for some easy matrices, I would like to solve the Poisson Equation with f(i,j)=-4 (as the unknown b in Ax=b) and boundary conditions phi(x,y)=x^2+y^2. The length of the subintervals in (0,1)x(0,1) should be 1/32 (1/N). That means we have 31 subintervals. I tried the following Matlab code to generate a 2D figure of the solutions (Jacobi) but it does not work. I know there are a lot of mistakes in here but I would be very grateful for any help!
% Define domain
a = 0; b = 1;
c = 0; d = 1;
tol=10^(-6);
% Define grid sizes
N = 32; % number of sub-intervals %number of points=33
hx = (b-a)/N; % length of sub-intervals in x-axis
hy = (d-c)/N; % length of sub-intervals in y-axis
%define u
u=zeros((N-1)^2,1); %(N-1)^2 unknown u(i,j)
s2 = 4*ones(N-1,1);
s= -1*ones(N-2,1);
B = diag(s2,0) + diag(s,1) + diag(s,-1);
% Sparse matrix B
B = sparse(B);
% Build sparse identity matrix
I = speye(N-1);
A = kron(B,I) + kron(I,B);
%f(i,j)=-4; %ersetzt b
%phi(a,b) ersetzt x
h2=1/(N*N); %h^2=h2
f(i,j)=-4;
phi(i,j)=i^2+j^2;
k=0;
b=-4*ones((N-1)^2,1);
err=Inf;
while err>tol
uold=u;
for i=1:(N-1)^2
s=0;
for j=1:(N-1)^2
if j~=i
s=s+A(i,j)*uold(j);
end
end
u(i)=(1/A(i,i))*(b(i)-s);
end
k=k+1;
err=norm(uold-u);
end
% Generate 2D arrays of grids
[X,Y] = meshgrid(a:hx:b,c:hy:d);
0 Comments
Answers (1)
Kavya Vuriti
on 7 Aug 2019
You can try creating two variables i and j which have linearly spaced elements between 0 to 1 with the width of subinterval as 1/32 and generate 2D arrays of i,j values using meshgrid function.
i=linspace(a,b,N); %linearly spaced elements between a=0 and b=1 with N=32
j=linspace(c,d,N);
[I,J] = meshgrid(i,j);
Also f(i,j) can be given as function handle with boundary conditions as shown below:
f=@(i,j) -4;
phi=I^2+J^2;
I can see from the code that the error doesn’t converge. Try checking the while loop Jacobi implementation so that the error converges. To generate a 2D figure, you may plot u obtained in each iteration with respect to iteration number.
0 Comments
See Also
Categories
Find more on Geometry and Mesh in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!