### 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[#[],#[]]&,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)) 