%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  第６回レポートのサンプルコード    2014/12/22    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n=4;
A=[2, -1,  0,  0; 
  -1,  2, -1,  0;
   0, -1,  2, -1;
   0,  0, -1,  2];
b = ones(n,1)  %  Vector b can also be inputed by:  b=[1,1,1,1]';
sol_x=A\b  % Exact solution.

D=diag(diag(A)); 
E=[0,  0,  0,  0; 
  -1,  0,  0,  0;
   0, -1,  0,  0;
   0,  0, -1,  0];
F=A-D-E;

%% EとFの他の作り方
% E=-diag(ones(n-1,1),-1);
% F=-diag(ones(n-1,1),1);

%% Jacobi method
M=D; N=A-M;

%% SOR Method (if w=1, it is the Gauss-Seidel method)
% w=1.71;
% M=D/w+E; N=(w-1)/w*D+F;

x=zeros(n,1); %Initial value of x.
errlist=zeros(100,1);
k=1;
while k<=100 
   last_x = x;
   x = -M\(N*last_x-b); 
   %%In case of Jacobi method, the above 1 line can be replaced by the following 3 lines:
   %for i=1:n
   %  x(i) = (b(i) - N(i,:)*last_x)/A(i,i);
   %end    
   errlist(k) = norm(x-sol_x)/norm(sol_x);
   k = k+1;
   if norm(x-last_x) < 1E-10
     break       
   end
end
k = k-1;

printf('Iteration stops after %d steps.\n', k);
printf('The error \|x^{k}-x^{k-1}\| is %15.10f.\n', norm(x-last_x));

hold off; %To clear the current figure
plot(1:k,log(errlist(1:k)), '-');
grid on
legend('Jacobi method');
