echo on
%%%%% laplace_jacobi.m %%%%%

disp(sprintf('**** Exampel of Jacobi iteration for Laplace eq. ****\n'));

%Two dimensional 3x3 elliptic boundary problem (Laplace equation)
%with boundary values 0 (left, bottom right) and 1 up.
%Natural ordering of nodes: 1 left bottom, ..., 3 right bottom,
%4 left second row (from bottom) etc. 9 right up.
%Laplace matrix L is composed of sub-matrices T and -I where
format short;
b=[0 0 0 0 0 0 1 1 1]'; %boundary conditions (upper edge = 1)
T=diag([4 4 4],0)+ diag([-1 -1],1)+ diag([-1 -1],-1)
I=eye(3)
Z=zeros(3,3)
LA=[T -I Z; -I T -I; Z -I T;]
%All diagonal elements are the same. Matrix LA can be saved by diagonals,
%or not at all because their structure is known.
%The system dimension (number of inner points, elements of vector x) is 9.
%Use an array 5x5 for the representation of the system domain with 
%initial and boundary conditions. 
%Ordering of nodes: (2,2) left up, ..., (2,4) right up,
% (3,2) left second row (from bottom) etc. (3,3) right bottom.
%Using Jacobi iteration and boundary conditions
U=[ones(1,5);zeros(4,5)]	% copy for the solution
%Initial approximation + boundary conditions
%D=[ones(1,5);0, [4 4 4], 0; 0, [4 4 4], 0;0, [4 4 4], 0;zeros(1,5);]
D=[ones(1,5);0, [0 0 0], 0; 0, [0 0 0], 0;0, [0 0 0], 0;zeros(1,5);]
%flops(0);
%first approximation for inner points:
for i=2:4
  for j=2:4
    U(i,j)=(D(i-1,j)+D(i,j-1)+D(i+1,j)+D(i,j+1))/4
  end
end
D=U;
%other steps
echo off;
for k=1:18    %number of additional steps
   for i=2:4
     for j=2:4
       U(i,j)=(D(i-1,j)+D(i,j-1)+D(i+1,j)+D(i,j+1))/4;
     end
   end
  D=U;
  k=k+1;
end
echo on;
U	%approximate solution after 19 steps (accurate in three digits)
%flops;
x=LA\b	%exact solution