4.13 How to numerically solve Poisson PDE on 2D using Jacobi iteration method?

Problem: Solve \(\bigtriangledown ^2 u= f(x,y)\) on 2D using Jacobi method. Assume \(f(x,y)=-\mathrm {e}^{-(x - 0.25)^2 - (y - 0.6)^2}\). Use mesh grid norm and relative residual to stop the iterations.

Mathematica

n = 21; (*number of grid points*) 
uNew = Table[0, {n}, {n}]; (*solution grid*) 
residual = uOld = uNew; 
 
h = 1/(n - 1); (*spacing*) 
c = 0.1;  (*controls the tolerance*) 
epsilon = c*h^2; (*error tolerance*) 
 
grid = N@Table[{i*h,j*h},{i,0,n-1},{j,0,n-1}]; 
 
(*the force function*) 
f[x_,y_] := -Exp[-(x - 0.25)^2 - (y - 0.6)^2]; 
 
(*evaluate force on the grid*) 
force = Map[f[#[[1]],#[[2]]]&,grid,{2}]; 
normF = Sqrt[h]*Norm[Flatten@force,2];(*force norm*) 
converged = False; 
 
k = 0; (*iteration counter*) 
While[Not[converged], 
 k++; 
 For[i = 2, i < n, i++, 
  For[j = 2, j < n, j++, 
   uNew[[i, j]] = (1/4)*(uOld[[i-1,j]]+uOld[[i+1,j]] 
    + uOld[[i,j-1]]+uOld[[i,j+1]]-h^2*force[[i,j]]); 
 
   residual[[i,j]] = force[[i, j]]-1/h^2*(uOld[[i-1,j]]+ 
    uOld[[i+1,j]]+uOld[[i,j-1]]+uOld[[i,j+1]]-4*uOld[[i,j]]) 
   ] 
 ]; 
 uOld = uNew; (*updated at end*) 
 normR = Sqrt[h] Norm[Flatten@residual, 2]; 
 If[normR/normF < epsilon, converged = True] 
]; 
 
(*plot solution*) 
ListPlot3D[uOld, PlotRange -> All, ImagePadding -> 30, 
  ImageSize -> 400, Mesh -> {n, n}, 
  AxesLabel -> {"x", "y", "u(x,y)"}, 
  PlotRangePadding -> 0, BaseStyle -> 12, 
  PlotLabel -> Row[{"Converged in ", k, " iterations", 
            " ,epsilon=", epsilon, " ,h=", h}]]
 

Matlab

close all; clear all; 
n = 21; % grid points, including internal points 
c = 0.1; % constant of tolerance 
myforce = @(X,Y) -exp( -(X-0.25).^2 - (Y-0.6).^2 ); 
 
h     = 1/(n-1);        % grid spacing 
tol   = c * h^2;        % tolerance 
res   = zeros(n,n);     % storage for residual 
u     = res;            % storage for solution 
uNew  = u; 
k     = 0;                     % loop counter 
[X,Y] = meshgrid(0:h:1,0:h:1); % coordinates 
f     = myforce(X,Y); 
normf = norm( f(:),2); % find norm 
 
%Now start the iteration solver, stop when 
%relative residual < tolerance 
figure; 
i     = 2:n-1; 
j     = 2:n-1;      %the iterations vectorized 
done  = false; 
 
while ~done 
    k = k+1; 
 
    uNew(i,j) = (1/4)*( u(i-1,j) + u(i+1,j) + ... 
                  u(i,j-1) + u(i,j+1) - h^2 * f(i,j) ); 
    res(i,j) = f(i,j) - (1/h^2)*( u(i-1,j) + ... 
           u(i+1,j) + u(i,j-1) + u(i,j+1) - 4*u(i,j) ); 
 
    %uncomment to see it update as it runs 
 
    %mesh(X,Y,u);  %hold on; 
    %title(sprintf('solution at iteration %d',k)); 
    %zlim([0,.07]); 
    %drawnow; 
 
    if norm(res(:),2)/normf < tol 
        done = true; 
    end; 
 
    u = uNew; 
end 
 
mesh(X,Y,u); 
title(sprintf('solution at iteration %d',k))