Solve Poisson Equation (Dirichlet boundary condition) via Jacobi and Gauss-Seidel-Iteration
37 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Kavya Vuriti
2019-8-7
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Error Detection and Correction 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!