% 16.01 / 16.02 Unified Engineering, Fall 2006
% Introduction to MATLAB: Matrix Computations
% Instructor: Violeta Ivanova, MIT Academic Computing, violeta@mit.edu

% EXERCISE 2 System of Linear Equations
% MATLAB vectors, matrices, matrix operations, curve fitting, eigenvalues, plotting. 

% The coordinates (x1,y1), (x2,y2), (x3,y3) and (x4,y4) of four points are
% given: (1, 54), (5, 82), (8, 40), and (10, 10).
% Find a cubic equation y(x) that passes through the four points i.e.:
% y = c1*x^3 + c2*x^2 + c3*x + c4

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% A. SOLVE A SYSTEM OF LINEAR EQUATIONS FOR C1, C2, C3, AND C4.

% A1. Define column vectors X and Y from the points' coordinates 
% X = (x1, x2, x3, x4) and Y = (y1, y2, y3, y4)

% One way to create a column vector is by transposing a row vector:
X = [1 5 8 10]; X=X'

% The simplest way is to enter the elements separated by semi-colon.
Y = [54; 82; 40; 10]

% A2. Define matrix A from X
% If C = (c1, c2, c3, c4) is the vector we want to solve for, 
% we need to define a matrix from the vector X, so that Y = AC.

% Method 1:
% Note the formatting of the matrix below - new line denotes end of a row and
% extra spaces are ignored i.e. you can format your code to be easy to read.
% This method is too long and is given mainly to illustrate the format. 

A = [X(1)^3     X(1)^2     X(1)    1
     X(2)^3     X(2)^2     X(2)    1
     X(3)^3     X(3)^2     X(3)    1
     X(4)^3     X(4)^2     X(4)    1]  

% Method 2:
% You can use this method even if you do not know the size of X.
A = [X.^3   X.^2   X.^1    X.^0]

% Method 3:
% Hint: this method uses the built-in function ONES.
% If you do not know the size of X, you can use: ones(length(X),1);
A = [X.^3   X.^2   X   ones(4,1)]

% A3. Solve Y = AC for C

% Method 1 (not recommended) Using the inverse matrix of A, i.e. c = A^(-1)*Y
Ainv = A^(-1)
C = Ainv * Y

% Method 2 (recommended) Using Gaussian Elimination 
% Hint: this method applies the \ operator
C = A \ Y

% Method 3 Using Row Reduction of [A|Y]
% Hint: This method applies the built-in function RREF. First, create a new 
% matrix AY in which matrix A forms the first four collumns and vector Y 
% forms the last collumn. Then use function RREF. 
AY = [A Y]
R = rref(AY)
C = R(:,5)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% B. COMPUTE A QUADRATIC AND A STRAIGHT LINE FIT USING POLYFIT

% B1. Quadratic fit
C2 = polyfit(X, Y, 2);

% B2 Straight line fit
C1 = polyfit(X, Y, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% C. PLOT THE CUBIC, QUADRATIC, STRAIGHT LINE AND POINTS

% C1. Compute a cubic curve, a quadratic curve, and a straight line 
% Hint: use built-in function POLYVAL
x = [0 : 0.1 : 10];
y = polyval(C, x);
y1 = polyval(C1, x);
y2 = polyval(C2, x);

% C2. Plot the cubic curve, quadratic curve, straight line, and points on
% on the same figure
p = plot (x, y, 'r-');
hold on;
p2 = plot (x, y2, 'g-');
hold on;
p1 = plot (x, y1, 'k-');
hold on;
P = plot (X, Y, 'bo');
hold off;
xlabel('X', 'FontSize', 12);
ylabel('Y', 'FontSize', 12);
title('CURVE FITTING TO POINTS', ...
       'FontSize', 12, ...
       'FontWeight', 'bold');
legend ('Cubic Fit', 'Quadratic Fit', 'Straight Line Fit', 'Data Points');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% D. EIGENVVALUES AND EIGENVECTORS

% An eigenvalue and eingenvector of a square matrix A are a scalar lambda
% and a nonzero vector V such that A*lambda = A*V.
 
% D1. Compute the eigenvalues and eigenvectors of A
% Hint: use built-in function EIG.

% With two output parameters, e.g. S and Lambda, EIG computes the eigenvectors 
% in the columns of S and stores the eigenvalues along the diagonal of Lambda.

[S, Lambda] = eig(A)

% D2. Check solution
Sinv = S^(-1)
B = S*Lambda*Sinv

% B should be the same as A
A