%% Assignment #1, Problem 5

%% Solving linear system (i)
% Gaussian elimination fails, 
% it does not give  matrix U with nonzero diagonal elements.
% Hence there is either no solution or infinitely many solutions.
%
% Note that the backslash command in machine arithmetic gives a warning "matrix is singular to working precision". But it does not return any solution (despite the fact that there are solutions). 
%
% Using symbolic computation in Matlab we obtain a particular solution xpart, and a basis V for the null space (one vector). Hence the solution is given by a line xpart+t*V with arbitrary t.

A = [2 -1 -1; -1 2 -1; -1 -1 2]
[L,U,p] = lu(A,'vector'); U            % matrix U from Gaussian elimination
b = [1;-3;2]
% Using Matlab with machine arithmetic
x = A\b                                % Matlab gives NaN (not a number)
% Using symbolic computation:
xpart = sym(A)\sym(b)                  % symbolic \ gives one solution and
                                       % warning that more solutions exist
V = null(sym(A))                       % basis for null space

%% Plotting planes for linear system (i)
% Here the three planes have an intersection which is a line,
% corresponding to the solutions xpart+t*V we found before.

T1 = diag( b(1)./A(1,:) )              % three points satisfying eq.1
T2 = diag( b(2)./A(2,:) )              % three points satisfying eq.2
T3 = diag( b(3)./A(3,:) )              % three points satisfying eq.3

fillpoints(stretch(T1,9),'y'); hold on % plot plane for eq.1
fillpoints(stretch(T2,3),'g')          % plot plane for eq.2
fillpoints(stretch(T3,5),'c');         % plot plane for eq.3

plotpoints(xpart-3*V,xpart+4*V,'r','LineWidth',3);  % plot line of solutions
hold off
nice3d; alpha(0.8); view(24,10); title('red line is solution set')

%% Solving linear system (ii)
% Gaussian elimination gives a matrix U with nonzero diagonal elements,
% hence there is a unique solution.

A = [2 -1 -1; -1 4 -1; -1 -1 2]
[L,U,p] = lu(A,'vector'); U            % matrix U from Gaussian elimination
b = [1;-3;2]
x = A\b                                % unique solution

%% Plotting planes for linear system (ii)
% The three planes intersect in a single point x
% which is the solution of the linear system we found before.

T1 = diag( b(1)./A(1,:) )              % three points satisfying eq.1
T2 = diag( b(2)./A(2,:) )              % three points satisfying eq.2
T3 = diag( b(3)./A(3,:) )              % three points satisfying eq.3

fillpoints(stretch(T1,9),'y'); hold on % plot plane for eq.1
fillpoints(stretch(T2,3),'g')          % plot plane for eq.2
fillpoints(stretch(T3,5),'c');         % plot plane for eq.3
plotpoints(x,'o'); label(x,'x')        % plot solution point
hold off
nice3d; alpha(0.75); view(-135,20); title('Three planes intersect at solution point x')

