% BRIEF RECAP of October 8th, 1998 NEEP271 lecture
% written by Fabian Waleffe 98/10/08
%
% We talked about solving linear systems of equations
% e.g. 2 equations in 2 unknowns or
%      M equations in M unknowns
% 
% The basic idea to obtain the solution is GAUSSIAN ELIMINATION
% where one systematically eliminates 1 unknown from the equations
% down to the point where there is only one unknown left.
%
% Linear systems are generally written in the Matrix vector form:
%
%                     A x = b
% 
% where A is a matrix (i.e. a table of numbers), 
%       x and b are column vectors
% 
% for instance for the system
%   
%    5 x1 + 3 x2 = 2;   3 x1 + 2 x2 = 1 

A=[ 5  3 ; 3  2]

b=[2 ; 1]

% The solution by Gaussian Elimination is built-in Matlab's "\" command:

x= A\b    

% check:

A*x

% Here is a larger example:

A=rand(10)
b=[1:10]'
x=A\b

% check

A*x

% a better way to check the error is to compute a RELATIVE ERROR:

norm(A*x-b)/norm(b)

% "norm" measures the size of a vector, by default :
%  norm(v) = sum over i of v(i)^2  (i.e. sum of squares of each component)

% Matlab's "\" is pretty powerful, in particular it also works for
% OVERDETERMINED systems (more equations than unknowns), in which case
% it finds the BEST SOLUTION IN LEAST-SQUARES SENSE. Example:


A=[ 5  3 ; 3  2; 2 1]
b=[2 ; 1 ; 0]

% There are too many equations: 3 for only 2 unknowns
% so this system does not have an exact solution,
% yet, Matlab gives:

x=A\b

% This is the BEST SOLUTION IN A LEAST-SQUARES SENSE

A*x

norm(A*x-b)/norm(b)

disp('Try computing the relative error for other x''s ')

xbest=A\b;

x=xbest .* (1+0.1*randn(size(x)))  ,  norm(A*x-b)/norm(b)

disp('You''ll never get anything smaller than')

norm(A*xbest-b)/norm(b)