clc
%Click the run bottom and refer to the comaand window
%These are the inputs that can be modified by the user
%A = nxn invertible square matrix
A=[8,12,1,7;1,2,3.001,6;57,9,1,4;1.047,108,98,5];
%n = number of equations
n=4;
%Bits_used_for_mantissa = number of bits used for mantissa in
%floating point representation
Bits_used_for_mantissa=23;

disp(sprintf('Adequacy in Solution of Simultaneous Linear Equations'))
disp(sprintf('2007 Fabian Farelo, Autar Kaw'))
disp(sprintf('University of South Florida')) 
disp(sprintf('United States of America'))
disp(sprintf('kaw@eng.usf.edu\n'))


disp(sprintf('NOTE: This worksheet demonstrates the use of Matlab to illustrate\nthree different techniques in finding the condition number of a\ncoefficient matrix [A]nxn.'))

%---------------------------------------------------------------
disp(sprintf('\n************************** Introduction ***********************'))

disp(sprintf('\nThe condition number allows one to quantify the accuracy in solution of\n[A] [X] = [C], where [A]nxn is an invertible square matrix, [X]nx1\nis the solution vector, and [C]nx1 is the right hand side array.\nMultiply the condition number by machine epsilon, and compare the result\nto 0.5 x 10-m to find out at least how many m significant digits\nare correct in solution.'))
disp(sprintf('\n        0.5 x 10^(-m) < Cond(A)*(machine epsilon)              (1.1)'))
disp(sprintf('\nThe following simulation uses three different techniques to determine\nthe condition number of coefficient matrix [A]nxn.'))
%---------------------------------------------------------------
disp(sprintf('\n**************************** Input Data ************************'))
disp(sprintf('\nBelow are the input parameters to begin the simulation. This is the only\nsection that requires user input. Once the values are entered, Matlab\nwill calculate the condition number using an exact method as well\nas two approximate methods.'))
disp(sprintf('\nInput Parameters:\n'))
disp(sprintf('A = nxn invertible square matrix'))
disp(sprintf('n = number of equations'))
disp(sprintf('Bits_used_for_mantissa = number of bits used for mantissa in floating\npoint representation\n\n'))
disp(sprintf('These are the defaults parameters, which can be changed at\nthe top of the worksheet'))
A
n
Bits_used_for_mantissa
%---------------------------------------------------------------
disp(sprintf('\n************ Calculating the Condition Number *****************'))
disp(sprintf('\nIn this section, three distinct methods are used to calculate the condition\nnumber of coefficient matrix [A]. Each has its own advantages while\nutilizing theorems that relate the norm of a matrix to the conditioning\nof the matrix.\n'))
disp(sprintf('\n======================================================================='))
%---------------------------------------------------------------
disp(sprintf('\n\nMethod 1: Finding the exact value of the Condition Number of a matrix'))
disp(sprintf('The following method finds the exact condition number of a square matrix.'))
disp(sprintf('The exact formula is given by\n'))

disp(sprintf('         Cond(A) = ||A||*||A-1||        (3.1.A)                    \n   '))        
 
disp(sprintf('Once the condition number is calculated, it can then be used to'))
disp(sprintf('solve for m, the number of significant digits that one can trust in\nthe solution.'))
disp(sprintf('Please note that, although this is the most direct method it'))
disp(sprintf('may not be practical in its computational time for higher order'))
disp(sprintf('matrices because this method requires calculation of the inverse of'))
disp(sprintf('coefficient matrix [A].  The problem in finding the inverse lies in'))
disp(sprintf('solving n sets of n equations which can be computationally intensive'))
disp(sprintf('for large coefficient matrices.'))

%Calculating  [A-1].
Ainverse=A^(-1);
%Calculating  ||A||.
NA=norm(A,inf);
%Calculating  ||A-1||.
NAInv=norm(Ainverse,inf);
format short g
disp(sprintf('\n|A|=%d',NA))
disp(sprintf('\n|A^(-1)|=%d',NAInv))
disp(sprintf('\nCond_A = ||A||*||A-1||\n'))
%Calculating the condition number of [A] knowing that Cond(A)=||A||*||A-1||.
Cond_A=NA*NAInv

%Calculating machine epsilon:
Machine_epsilon=(2^(-Bits_used_for_mantissa));

%Calculating the number of significant digits that one can trust in solution:
r=Machine_epsilon*Cond_A;
sig_digits=floor(eval((solve('0.5*10^(-m)=r','m'))));
Trust_digits=max(0,sig_digits)

%---------------------------------------------------------------------
disp(sprintf('\n======================================================================='))
disp(sprintf('Method 2: Finding an approximate value of the condition number'))
disp(sprintf('The following numerical method finds the condition number of [A] using\nthe inequality'))

disp(sprintf('\nCond(A) >  [||DX|| / ||X+DX||]/[  ||DC|| / ||C|| ]        (3.2.A)\n'))
disp(sprintf('However, the value (X + DX) is equivalent to X\". The inequality then\nbecomes:'))

disp(sprintf('\nCond(A) >  [||DX|| / ||X\"||] / [||DC|| / ||C||]           (3.2.B)\n'))


disp(sprintf('where ||DX|| / ||X\"|| is the relative change in the norm\nof the solution vector and ||DC|| / ||C|| is the relative change\nin the norm of the right hand side vector. The ratio between\nthese two values quantifies the conditioning of a system of equations,\ndemonstrating the accuracy in solution.'))
disp(sprintf('That is, for any small change made in the right hand side array,'))
disp(sprintf('the resulting change in the solution vector will govern how accurate the\nsystem'))
disp(sprintf('\b is and therefore how many significant digits one can trust in the\nsolution of a'))
disp(sprintf('\b system of simultaneous linear equations.'))

disp(sprintf('\nIn this method, the condition number is calculated using')) 
disp(sprintf('equation (3.2.B) by first conducting the following steps:'))
disp(sprintf('\n1)  A right hand side vector [C] is chosen such that the solution\nvector equals 1. (i.e. [X] = [1,1,...,1]).'))
disp(sprintf('\n2) A new, unbiased right hand side vector [C\"] is then generated\nby Matlab. This is done by adding a random positive or negative value to\neach element of the old right hand side vector.'))
disp(sprintf('\n3) The new right hand side vector [C\"] can then be used to\ncalculate a new solution vector [X\"].'))
disp(sprintf('\n\nBy creating a small relative change in the right hand side\narray (i.e., ||DC||/||C|| < 1), the magnitude of the condition number\nwill be largely influenced by the relative error in the solution vector,\ndemonstrating how accurate the solution actually is.'))
disp(sprintf('\n\nNOTE: Each time the worksheet is executed, Matlab will generate\na new condition number. The user should run the worksheet several\ntimes to see if the estimate of the condition number approaches the exact\nvalue given in Method 1. The greatest of the generated values\nwill be the most accurate approximation and will never exceed the\ntrue condition number due to the above inequality, Eq.(3.2.B.)'))

disp(sprintf('Parameter names:'))
disp(sprintf('RHS = old right hand side array [RHS]'))
disp(sprintf('RHS1 = new right hand side array [RHS"]'))
disp(sprintf('X = old solution vector [X]'))
disp(sprintf('X1 =  new solution vector [X"]'))


%Defining RHS as a vector.
RHS=ones(1,n);
%Defining RHS1 as a vector.
RHS1=ones(1,n);
%Defining X as a vector.
X=ones(1,n);
%Defining X1 as a vector.
X1=ones(1,n);
disp(sprintf('Step 1: Assigning values to [RHS] so that [X] equals [1,..,1].')) 
RHS=sum(A');
RHS
disp(sprintf('Solving for [X].'))
%X=eval(solve('A*X=RHS','X'))
h=horzcat(A,RHS');
re=rref(h);
X=re(1:n,n+1)
%---------------------------------------------------
disp(sprintf('\nStep 2: Generating the new, unbiased right hand side\nvector [RHS"] by adding a random + value between -0.001*RHS[i]\nand +0.001*RHS[i] to create a small DC value:'))
for i=1:n
% Randomizing the sign of the value added to each element of new right hand side vector.
sign_val=floor(1.5+rand(1));
RHS1(i)=RHS(i)+(-1)^sign_val*0.001*RHS(i)*rand(1);
end
RHS1
disp(sprintf('Step 3: Calculating the new solution vector [X"].'))
h=horzcat(A,RHS1');
re=rref(h);
X1=re(1:n,n+1)
%--------------
disp(sprintf('\nNotice the small difference between the values of the old'))
disp(sprintf('and new right hand side vectors. Now look at the new solution'))
disp(sprintf('vector [X1]. Do the values deviate much from the old solution'))
disp(sprintf('vector [X]? If not, then any small changes that are made in the'))
disp(sprintf('right hand side array do not affect the accuracy of the solution '))
disp(sprintf('vector, and the solution can therefore be trusted. Otherwise, '))
disp(sprintf('the system of equations is ill-conditioned. Below, the condition '))
disp(sprintf('number of [A] is calculated by determining the values required by '))
disp(sprintf('Equation (3.2.B).\n'))



%Calculating the relative change in the norm of the solution vector.
Num_Cond=norm(X1-X,inf)/norm(X1,inf);
%Calculating the relative change in the norm of the right hand side array.
Den_Cond=norm(RHS1-RHS,inf)/norm(RHS,inf);
format short g
disp(sprintf('\n|X1-X|/|X1|=%d',Num_Cond))
disp(sprintf('\n|RHS1-RHS|/|RHS1|=%d',Den_Cond))
disp(sprintf('\nCond_AA =[|X1-X|/|X1|]*[|RHS1-RHS|/|RHS1|]'))
%Calculating the condition number.
Cond_AA=Num_Cond/Den_Cond

   

disp(sprintf('\n======================================================================='))
disp(sprintf('\nMethod 3: A second technique in approximating the condition number'))

disp(sprintf('\nThe following numerical technique is simpler than Method 2 as')) 
disp(sprintf('it requires fewer calculations. This method utilizes the theorem'))
disp(sprintf('\n     Cond(A) >    ||A|| ||X||/||C||            (3.3.A)'))                                                                                                                                                                            
disp(sprintf('\nThe proof is as follows:'))
disp(sprintf('Let '))
disp(sprintf('[X] = [A]^(-1)* [C]'))                                                                                                                                                                                                                                           

disp(sprintf('\nApplying norm properties, this equation becomes')) 
disp(sprintf('||X|| < ||A-1|| ||C||                               '))                                                                                                                                                                                                           
disp(sprintf('\nMultiplying ||A|| to both sides gives'))
disp(sprintf('\n||A|| ||X|| < Cond(A) ||C||             '))                                                                                                                                                                                                          
disp(sprintf('\nby the definition of condition number.'))
disp(sprintf('\nAnd division by ||C|| results in the final inequality.'))
disp(sprintf('In this technique, however, the right hand side values of')) 
disp(sprintf('[C] matrix are chosen to equal [+1,  +1,  +1] with signs generated')) 
disp(sprintf('randomly. This will result in ||C|| = 1, minimizing the number of'))
disp(sprintf('calculations required to solve for Cond(A). Therefore, the calculation')) 
disp(sprintf('is reduced to:'))
disp(sprintf('\n      Cond(A) > ||A|| ||X||                        (3.3.B) '))                                                                                  
disp(sprintf('\nAgain, the user should re-execute the worksheet to see if the greatest'))
disp(sprintf('approximate condition number approaches the exact condition number '))
disp(sprintf('defined in Method 1.'))




%Randomly choosing a (+1) or (-1) value for each element of the right hand side array. 
for i=1:n
   sign_val=floor(1.5+rand(1));
   RHS(i)=(-1)^sign_val;
end
%Solving for the solution vector [X].
h=horzcat(A,RHS');
re=rref(h);
X=re(1:n,n+1)
%Calculating the condition number.
disp(sprintf('Cond_AAA = ||A|| ||X||'))   
Cond_AAA=norm(A,inf)*norm(X,inf)
