Solve Poisson Equation (Dirichlet boundary condition) via Jacobi and Gauss-Seidel-Iteration

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);

Answer 1

[Kavya Vuriti]

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.