%section 8.2 computer problem #2

A = [7,1,-1,2;1,8,0,-2;-1,0,4,-1;2,-2,-1,6];
b = [3;-5;4;-3];
n = length(b);

x_0 = [0;0;0;0];         %initial guess
x_exact = [1;-1;1;-1];   %exact solution

%first we do Jacobi iteration

Q = zeros(n,n);
for i=1:n
    Q(i,i) = A(i,i);  %Q = diagonal matrix with same diagonal entries as A
end

x = x_0;
ct = 0;               %iteration count
while max(abs(x-x_exact)) > 0.0001
    ct = ct+1;
    y = (Q-A)*x + b;
    for j=1:n
        x(j) = y(j)/Q(j,j);   %solve Qx = y
    end
end

disp('Jacobi solution is...')
x
disp('number of Jacobi iterations needed is...')
ct

%now we do Gauss-Seidel iteration

I = eye(n);

x = x_0;
ct = 0;               %iteration count
while max(abs(x-x_exact)) > 0.0001
    ct = ct+1;
    for i=1:n
        x(i) = b(i);
        for j=1:i-1
            x(i) = x(i)-A(i,j)*x(j);
        end
        for j=i+1:n
            x(i) = x(i)-A(i,j)*x(j);
        end
        x(i) = x(i)/A(i,i);
    end
end

disp('Gauss-Seidel solution is...')
x
disp('number of Gauss-Seidel iterations needed is...')
ct

%finally we do SOR iteration

w = 1.1;              %acceleration factor

x = x_0;
ct = 0;               %iteration count
while max(abs(x-x_exact)) > 0.0001
    ct = ct+1;
    for i=1:n
        x(i) = x(i)*A(i,i)*(1-w)/w;
        x(i) = x(i)+b(i);
        for j=1:i-1
            x(i) = x(i)-A(i,j)*x(j);
        end
        for j=i+1:n
            x(i) = x(i)-A(i,j)*x(j);
        end
        x(i) = w*x(i)/A(i,i);
    end
end

disp('SOR solution is...')
x
disp('number of SOR iterations needed is...')
ct
