% Example of Matrix Iteratiions for Jacobi and Gauss-Seidel Methods

A=[4 3 0; 3 4 -1;0 -1 4];
b=[24 ;30; -24];
n=10;
z=zeros(size(b));
Exact = [3; 4; -5];

D = diag(diag(A));
L=tril(D-A);
U=triu(D-A);


% ------------------------------------------------
% Jacobi Iteration
PJ=D\(L+U);
c = D\b;
x(:,1)=z;
v(:,1)=x(:,1)-Exact;
res(:,1)=b-A*x(:,1);
for k = 2:n+1
   x(:,k) = PJ*x(:,k-1) + c;
   v(:,k)=x(:,k)-Exact;
   res(:,k)=b-A*x(:,k);
end

% Plot errors
figure(1)
hold off
plot(max(res),'b+-')
hold on
plot(max(v),'ro:')
hold off
title('Jacobi Iteration')
xlabel('Iteration')
ylabel('Error')
legend('Residual','Error')

% ------------------------------------------------
% Gauss-Seidel Iteration
PG = (D-L)\U;
c = (D-L)\b;
z=zeros(size(b));
x=zeros(size(b));
x(:,1)=z;
v(:,1)=x(:,1)-Exact;
res(:,1)=b-A*x(:,1);
for k = 2:n+1
   x(:,k) = PG*x(:,k-1) + c;
   v(:,k)=x(:,k)-Exact;
   res(:,k)=b-A*x(:,k);
end

% Plot errors
figure(2)
hold off
plot(max(res),'b+-')
hold on
plot(max(v),'ro:')
hold off
title('Gauss Seidel Iteration')
xlabel('Iteration')
ylabel('Error')
legend('Residual','Error')


