%uses condition number to estimate number of sig digits lost solving Ax = b

format long

A = [1,1/2,1/3,1/4;1/2,1/3,1/4,1/5;1/3,1/4,1/5,1/6;1/4,1/5,1/6,1/7];
b = [0;1/2;1/3;1/4];

eig_values = eig(A);                 %eigenvalues of A
sing_values = sqrt(eig(transpose(A)*A));   %eigenvalues of A^T A

%the condition number can be computed from 
%eigenvalues of A when A is symmetric (or normal)

%if A is not symmetric, it must be computed from singular values
%singular values are square roots of eigenvalues of A^T A

disp('the condition number computed from eigenvalues of A^T A is...')
kappa = max(sing_values)/min(sing_values)

disp('the condition number computed from eigenvalues of A is...')
kappa = max(abs(eig_values))/min(abs(eig_values))

disp('the expected number of sig digits lost in solution to Ax = b is...')
floor(log10(kappa))   %log_10(kappa) rounded down to nearest integer

disp('the solution to Ax = b is...')
linsolve(A,b)
