%richardson iteration for approximating the solution to Ax = b
%compares actual error to error bound

A = [1,0,1/2;0,1,1;1,0,1];
b = [1;2;3];
x_0 = [1;1;1];        %initial guess of solution
I = eye(3);           %identity matrix
N = 10;               %number of iterations

%first we do richardson iteration to solve Ax = b

x = x_0;
for i=1:N
    x = (I-A)*x+b;
end

disp('estimate of the solution is...')
x

disp('exact solution is')
x_exact = linsolve(A,b)

%now we check error bound

[P,D] = eig(I-A);          %this gives diagonalization I-A = P*D*P^(-1)
r = max(abs(diag(D)));     %maximum eigenvalue of I-A in absolute value
kappa_P = max(sqrt(eig(transpose(P)*P)))/min(sqrt(eig(transpose(P)*P)));

disp('theoretical error bound is...')
r^N*kappa_P*norm(x_0-x_exact)

disp('actual error is...')
norm(x-x_exact)            %Euclidean distance between vectors x, x_exact





