%% Chapter 7 Exercise Solutions
%
%% Exercise 7.1
% \begin{array}{rl}
%     4x_1 +x_2 -x_3 +x_4 &=14,\\
%     x_1 +4x_2 -x_3 -x_4 &=10,\\
%     -x_1 -x_2 +5x_3 +x_4 &=-15,\\
%     x_1 -x_2 +x_3 +3x_4 &=3.
% \end{array}
% The Jacobi method for this system is
% \begin{array}{rl}
%   x_1^{(k+1)} &= \frac{1}{4}(14 - x_2^{(k)} + x_3^{(k)} - x_4^{(k)}), \\
%   x_2^{(k+1)} &= \frac{1}{4}(10 - x_1^{(k)} + x_3^{(k)} + x_4^{(k)}), \\
%   x_3^{(k+1)} &= \frac{1}{5}(-15 + x_1^{(k)} + x_2^{(k)} - x_4^{(k)}), \\
%   x_4^{(k+1)} &= \frac{1}{3}(3 - x_1^{(k)} + x_2^{(k)} - x_3^{(k)}).
% \end{array}
% First iteration:
% \begin{array}{rl}
%   x_1^{(1)} &= \frac{1}{4}(14 - 0 + 0 - 0) = 3.5, \\
%   x_2^{(1)} &= \frac{1}{4}(10 - 0 + 0 + 0) = 2.5, \\
%   x_3^{(1)} &= \frac{1}{5}(-15 + 0 + 0 - 0) = -3, \\
%   x_4^{(1)} &= \frac{1}{3}(3 - 0 + 0 - 0) = 1.
% \end{array}
% Second iteration:
% \begin{array}{rl}
%   x_1^{(2)} &= \frac{1}{4}(14 - 2.5 - 3 - 1) = 1.875, \\
%   x_2^{(2)} &= \frac{1}{4}(10 - 3.5 - 3 + 1) = 1.125, \\
%   x_3^{(2)} &= \frac{1}{5}(-15 + 3.5 + 2.5 - 1) = -2, \\
%   x_4^{(2)} &= \frac{1}{3}(3 - 3.5 + 2.5 + 3) = 1.6667.
% \end{array}
%
%% Exercise 7.2
% The Gauss-Seidel method for this system is
% \begin{array}{rl} x_1^{(k+1)} &= \frac{1}{4}(14 - x_2^{(k)} + x_3^{(k)} - x_4^{(k)}), \\ x_2^{(k+1)} &= \frac{1}{4}(10 - x_1^{(k+1)} + x_3^{(k)} + x_4^{(k)}), \\ x_3^{(k+1)} &= \frac{1}{5}(-15 + x_1^{(k+1)} + x_2^{(k+1)} - x_4^{(k)}), \\ x_4^{(k+1)} &= \frac{1}{3}(3 - x_1^{(k+1)} + x_2^{(k+1)} - x_3^{(k+1)}).\end{array}
% First iteration:
% \begin{array}{rl}x_1^{(1)} &= \frac{1}{4}(14 - 0 + 0 - 0) = 3.5, \\ x_2^{(1)} &= \frac{1}{4}(10 - 3.5 + 0 + 0) = 1.625, \\ x_3^{(1)} &= \frac{1}{5}(-15 + 3.5 + 1.625 - 0) = -1.975, \\ x_4^{(1)} &= \frac{1}{3}(3 - 3.5 + 1.625 + 1.975) = 1.0333. \end{array}
% Second iteration:
% \begin{array}{rl} x_1^{(2)} &= \frac{1}{4}(14 - 1.625 - 1.975 - 1.0333) = 2.3417, \\ x_2^{(2)} &= \frac{1}{4}(10 - 2.3417 - 1.975 + 1.0333) = 1.6792, \\ x_3^{(2)} &= \frac{1}{5}(-15 + 2.3417 + 1.6792 - 1.0333) = -2.4025, \\ x_4^{(2)} &= \frac{1}{3}(3 - 2.3417 + 1.6792 + 2.4025) = 1.58. \end{array}
%
%% Exercise 7.3
% Checking that the coefficient matrix has all real eigenvalues.

A = [ 4, 1, -1, 1 ;
      1, 4, -1, -1 ;
     -1, -1, 5, 1 ;
      1, -1, 1, 3 ];

eig(A)

% All eigenvalues of $A$ are real so we can use
% \omega = 1 + \left( \frac{\rho(T_J)}{1 + \sqrt{1 - \rho(T_J)^2}} \right)^2.
% The iteration matrix for the Jacobi method is
% T_J = D^{-1} (L + U)

L = tril(A, -1);
U = triu(A, 1);
D = A - L - U;

TJ = -inv(D) * (L + U)

% Calculating \rho(T_J)

rhoTJ = max(abs(eig(TJ)))

% So \rho(T_J) = 0.6054, therefore
% \omega = 1 + \left( \frac{0.6054}{1 + \sqrt{1-0.6054^2}} \right)^2 = 1.1136.
% The SOR method for this system is
% \begin{array}{rl}
%   x_1^{(k+1)} &= (1 - \omega)x_1^{(k)} + \frac{\omega}{4}(14 - x_2^{(k)} + x_3^{(k)} - x_4^{(k)}), \\
%   x_2^{(k+1)} &= (1 - \omega)x_2^{(k)} + \frac{\omega}{4}(10 - x_1^{(k+1)} + x_3^{(k)} + x_4^{(k)}), \\
%   x_3^{(k+1)} &= (1 - \omega)x_3^{(k)} + \frac{\omega}{5}(-15 + x_1^{(k+1)} + x_2^{(k+1)} - x_4^{(k)}), \\
%   x_4^{(k+1)} &= (1 - \omega)x_4^{(k)} + \frac{\omega}{3}(3 - x_1^{(k+1)} + x_2^{(k+1)} - x_3^{(k+1)}).
% \end{array}
%
% First iteration:
% \begin{array}{rl}
%   x_1^{(1)} &= (1 - 1.1136)(0) + \frac{1.1136}{4}(14 - 0 + 0 - 0) = 3.8977, \\
%   x_2^{(1)} &= (1 - 1.1136)(0) + \frac{1.1136}{4}(10 - 3.8977 + 0 + 0) = 1.6989, \\
%   x_3^{(1)} &= (1 - 1.1136)(0) + \frac{1.1136}{5}(-15 + 3.8977 + 1.6989 - 0) = -2.0944, \\
%   x_4^{(1)} &= (1 - 1.1136)(0) + \frac{1.1136}{3}(3 - 3.8977 + 1.6989 + 2.0944) = 1.0749.
% \end{array}
% Second iteration:
% \begin{array}{rl}
%   x_1^{(2)} &= (1 - 1.1136)(3.8977) + \frac{1.1136}{4}(14 - 1.6989 - 2.0944 + 1.0749) = 2.0994, \\
%   x_2^{(2)} &= (1 - 1.1136)(1.6989) + \frac{1.1136}{4}(10 - 2.0994 - 2.0994 - 1.0749) = 1.7227, \\
%   x_3^{(2)} &= (1 - 1.1136)(-2.0944) + \frac{1.1136}{5}(-15 + 2.0994 + 1.7227 - 1.0749) = -2.4910, \\
%   x_4^{(2)} &= (1 - 1.1136)(1.0749) + \frac{1.1136}{3}(3 - 2.0994 + 1.7227 + 2.4910) = 1.7763.
% \end{array}
%
%% Exercise 7.5
% The iteration matrices for the Gauss-Seidel and SOR methods are
% \begin{array}{rl}
%   T_{GS} &= -(L+D)^{-1} U, \\
%   T_{SOR} &= (D+\omega L)^{-1} ((1-\omega )D-\omega U).
% \end{array}

omega = 1.1136;
TGS = -inv(L + D) * U
TSOR = inv(D + omega * L) * ((1 - omega) * D - omega * U)
rhoTGS = max(abs(eig(TGS)))
rhoTSOR = max(abs(eig(TSOR)))

% The spectral radii of the iteration matrices for the three methods are
% \begin{array}{rl}
%   \rho(T_J) &= 0.6054, \\
%   \rho(T_{GS}) &= 0.3553, \\
%   \rho(T_{SOR}) &= 0.1968.
% \end{array}
% So the SOR method should converge the fastest.
%
%% Exercise 7.5

function [x, k] = jacobi(A, b, tol)

n = length(b);
x = zeros(n, 1);
maxiter = 100;
for k = 0 : 100
    xold = x;
    for i = 1 : n
        sum_ = 0;
        for j = 1 : n
            if j ~= i
                sum_ = sum_ + A(i,j) * xold(j);
            end
        end
        x(i) = (b(i) - sum_) / A(i,i);
    end
    r = b - A * x;
    if max(abs(r)) < tol
        break
    end
end

end
function [x, k] = gauss_seidel(A, b, tol)

n = length(b);
x = zeros(n, 1);
maxiter = 100;
for k = 0 : 100
    for i = 1 : n
        sum_ = 0;
        for j = 1 : n
            if j ~= i
                sum_ = sum_ + A(i,j) * x(j);
            end
        end
        x(i) = (b(i) - sum_) / A(i,i);
    end
    r = b - A * x;
    if max(abs(r)) < tol
        break
    end
end

end

function [x, k] = sor(A, b, omega, tol)

n = length(b);
x = zeros(n, 1);
for k = 0 : 100
    for i = 1 : n
        sum_ = 0;
        for j = 1 : n
            if j ~= i
                sum_ = sum_ + A(i,j) * x(j);
            end
        end
        x(i) = (1 - omega) * x(i) + omega * (b(i) - sum_) / A(i,i);
    end
    r = b - A * x;
    if max(abs(r)) < tol
        break
    end
end

end

% Define system
A = [ 4, 1, -1, 1 ;
      1, 4, -1, -1 ;
     -1, -1, 5, 1 ;
      1, -1, 1, 3 ];
b = [ 14 ; 10 ; -15 ; 3];

% Solve system
tol = 1e-6;
[x, jacobi_iterations] = jacobi(A, b, tol)
[x, GS_iterations] = gauss_seidel(A, b, tol)
[x, SOR_iterations] = sor(A, b, omega, tol)

%
%% Exercise 7.6
% The Jacobi iteration matrix for this system is
% T_J = -D^{-1}(L + U)  = \pmatrix{ 2 & 0 \cr 0 & 2 }^{-1} \pmatrix{ 0 & 1 \cr \alpha & 0 } = \pmatrix{ \frac{1}{2} & 0 \cr 0 & \frac{1}{2} } \pmatrix{ 0 & 1 \cr \alpha & 0} = \pmatrix{ 0 & \frac{1}{2} \cr \frac{1}{2}\alpha & 0}.
% Calculating the eigenvalues of T_J
% 0 =
% \left| \matrix{ -\lambda & \frac{1}{2} \cr \frac{1}{2}\alpha & -\lambda } \right|
%   = \lambda^2 - \frac{1}{4} \alpha
%
% therefore \lambda = \pm\frac{1}{2}\sqrt{\alpha}. For convergence \rho(T_J) <1 so \frac{1}{2}\sqrt{\alpha} < 1 therefore \alpha < \sqrt{2}.