% 
% Jacobi, G-S 和 SOR 的演示
%

% 以三对角线性方程组
clear all;
close all;

n = 128*8;
A = 2*eye(n) - diag(ones(n-1,1),1) - diag(ones(n-1,1),-1);
b = eye(n,1);

D = diag(diag(A));  % A 的对角线部分
L = D - tril(A);    % -L 为 A 的严格下三角部分
U = D - triu(A);    % -U 为 A 的严格上三角部分

x0 = zeros(n,1);  % 迭代初始值
itermax = 500;    % 迭代步数

semilogy(0,1,'b+', 0,1,'ro', 0,1,'kd');  
axis([-2,110,1e-3,2])
legend('Jacobi', 'GS', 'SOR');
title('k=0');    
hold on;
%pause(1)


x1 = x0; x2 = x0; x3 = x0;
for k = 1 : itermax    
    title(['k=',int2str(k)]);
    
    % Jacobi
    x1 = D \ ( (L+U)*x1 + b );
    relres_J = norm(b - A*x1) / norm(b);  
    semilogy(k, relres_J,'b+');
    legend('Jacobi', 'GS', 'SOR'); 
    
    % G-S
    x2 = (D-L) \ ( U*x2 + b );
    relres_G = norm(b - A*x2) / norm(b);  
    semilogy(k, relres_G,'ro');    
    legend('Jacobi', 'GS', 'SOR'); 
    
    % SOR
    omega = 1.5;
    x3 = (D-omega*L) \ ( ((1-omega)*D + omega*U)*x3 + omega*b );
    relres = norm(b - A*x3) / norm(b);  
    semilogy(k, relres,'kd');    
    legend('Jacobi', 'GS', 'SOR');
    % pause(0.1)
   
end

