%% Simulation 5 Companion

% clc % clear command window
% clear all % clear all variables

%% Matrix Computations in MATLAB
% Goal: Execute matrix computations using MATLAB.
% Let us execute a few matrix computations using the matrix.
% A = \left[\begin{array}{ccc}1&0&0\\2&1&-2\\3&2&1\end{array}\right]
% First let us define A as an array.
% Reminder: uncomment the start of the line and hit run section.

% A = [1 0 0; 2 1 -2; 3 2 1] % define 3 x 3 matrix A 

%% Compute Powers
% Next let us calculate A^2.
% Note: we do not need to use pointwise .^2 becuase for the first time we are actually doing matrix exponentiation ^2.

% A^2 % A squared

%% Compute Determinant
% How about the determinant of A?

% det(A) % determinant of A

%% Compute Inverse
% And the inverse of A?

% inv(A) % inverse of A

%% Construct Identity Matrix
% Let us construct the 3 \times 3 identity matrix.

% I = eye(3) % 3 x 3 identity matrix

% If we multiply A by I we should obtain A. Let us see if this is true.
%% Compute Matrix Multiplication
% Note: again do not need to use pointwise .* becuase for the first time we are actually doing matrix multiplication *.

% A % print A
% A*I % compute I*A

% Yay they are equal. :)
% We can also check using the isequal command.

% isequal(A,A*I) % checks if A equals A*I: outputs 1 if true, 0 if false

% Confirmation has been made.
%% Solve Linear System Ax = b
% Let us try to solve the equation:
% A\mathbf{x} = \mathbf{b} where \mathbf{b} = \left[\begin{array}{c}2\\-3\\12\end{array}\right]
% First we define \mathbf{b} as a column vector.

% b = [-2; -3; 12] % define column vector b

% Now let us calculate the solution to A\mathbf{x} = \mathbf{b} by "dividing–on–the–left" by A using the left-division operator \.

% x = A\b % define x to be the solution to Ax = b

% To be extra safe let us verify that A\mathbf{x} does in fact equal \mathbf{b}.

% b % print b
% A*x % compute A*x
% isequal(A*x,b) % checks if A*x equals b: outputs 1 if true, 0 if false

% Confirmation has been made.
%% Find Eigenvalues and Eigenvectors
% We can obtain the eigenvalues of A quite easily using eig.

% eig(A) % outputs a column vector of eigenvalues of A

% Alternatively I can obtain eigenvectors and eigenvalues together by indicating that the the ouput is an array [evec, eval].

% [M, D] = eig(A) % M is a matrix whose columns are the eigenvectors, and D is a diagonal matrix of corresponding eigenvalues

% For example the 2nd eigenvector and corresponding eigenvalue are located as:

% M(:,2) % second column of matrix M whose columns are eigenvectors
% D(2,2) % entry in position (2,2) of matrix D which is the diagonal matrix of eigenvalues

% To be extra safe let us verify A\mathbf{v} = \lambda\mathbf{v} for this choice of eigenvector \mathbf{v} and eigenvalues \lambda.

% A*M(:,2) % compute A(2nd evec)
% D(2,2)*M(:,2) % (2nd eval)(2nd evec)
% isequal(A*M(:,2),D(2,2)*M(:,2)) % checks if A*M(:,2) equals D(2,2)*M(:,2): outputs 1 if true, 0 if false

% Yay! ... Wait. They look equal. Why is isequal outputting 0 for false? This is due to errors due to rounding. Rounding is the bane of our existence.
% In this case it would be more appropriate to use isapprox which assesses whether the arrays are approximately equal.

% isapprox(A*M(:,2),D(2,2)*M(:,2),"tight") % outputs array of 1s and 0s: 1 for each entry of the arrays that are approximately equal
% the specification "tight" means that approximate equality occurs when values are within 10^(-12) of each other (that's close!) 

% We can see that this outputs a 3 x 1 array of 1s. This means every entry of the 3 x 1 vectors A*M(:,2) and D(2,2)*M(:,2) are approximately equal.
%% Plot Eigenvalues on the Complex Plane
% Now that we have our eigenvalues we can plot them in the complex plane.
% First we need to be able to obtain the real and imaginary parts of our eigenvalues.

% eigenvalues = eig(A) % column vector of eigenvalues of A
% real(eigenvalues) % column vector of real parts of eigenvalues
% imag(eigenvalues) % column vector of imaginary parts of eigenvalues

% Next we plot the imaginary parts against the real parts.

% figure % new figure window
% plot(real(eigenvalues),imag(eigenvalues),'x','Linewidth',2) % plot imag(eigenvalues) against real(eigenvalues) and mark points using thick 'x's
% grid on
% xlabel("Re(z)")
% ylabel("Im(z)")
% title("eigenvalues of A")

%
%% Systems of 1st–Order ODEs
% Goal: Solve and plot the solutions for 0 \leq t \leq 12 to the system:
% \left\{\begin{array}{l}y_1'=-2y_1+y_2+1-\cos 4t\\y_2' = y_1-2y_2+y_3\\y_3'=y_2-2y_3\end{array}\right. with initial values y_1(0) = y_2(0) = y_3(0) = 0
% We have actually already solved such systems in Simulation04 wherein we converted 2nd–order ODEs into systems of 1st–order ODEs.
% Our system is in standard form:
% \mathbf{y}' = \mathbf{f}(t,\mathbf{y}) with \mathbf{y}(0) = \mathbf{y}_0
% where:
% \mathbf{f}\left(t,\mathbf{y}\right) = \left[\begin{array}{c}-2y_1+y_2+1-\cos 4t\\ y_1-2y_2+y_3\\y_2-2y_3\end{array}\right] and \mathbf{y_0} = \left[\begin{array}{c}0\\0\\0\end{array}\right]
% We will use ode45 once more. Here is a reminder. You should not uncomment anything in this box.

%==========================================================================
%
% ode45: Numerically solves the initial value problem dy/dt = f(t,y) with
% y(0) = y0.
%
%    [t,y] = ode45(f,tspan,y0)
%
%--------------------------------------------------------------------------
%
% ------
% INPUT:
% ------
%   f : (function_handle) dy/dt = f(t,y)
%   tspan = [t0 tf] numerically solves from t0 to tf
%   y0 : (number/vector) initial value/s
%
% ------
% OUTPUT:
% ------
%   t : (column vector) t-values in numerical solution 
%   y : (column vector) values y(t) in numerical solution
%
%==========================================================================

% Let us define:
% \mathbf{f}\left(t,\mathbf{y}\right) = \left[\begin{array}{c}-2y_1+y_2+1-\cos 4t\\ y_1-2y_2+y_3\\y_2-2y_3\end{array}\right]
% as a function handle for a column–vector of functions applied to the variable t and the column vector y = [y(1);y(2);y(3)].

% f = @(t,y) [-2.*y(1)+y(2)+1-cos(4.*t); y(1)-2.*y(2)+y(3); y(2)-2.*y(3)]

% Next we define 0 \leq t \leq 12.

% tspan = [0 12]

% And next we define the initial values column vector:
% \mathbf{y}(0) = \left[\begin{array}{c}0\\0\\0\end{array}\right]

% y0 = [0;0;0]

% We are ready to solve the differential equation.

% [t,y] = ode45(f,tspan,y0)

% Recall that y has multiple columns here. The ith column y(:,i) are the values y_i(t).
% Let us plot the solution y_1(t), y_2(t), y_3(t) onto a single set of axes.

% figure % new figure window
% plot(t,y(:,1)) % plot y_1(t) against t
% hold on 
% plot(t,y(:,2)) % plot y_2(t) against t
% hold on 
% plot(t,y(:,3)) % plot y_3(t) against t
% hold off
% xlabel("t") % label the horizontal axis as t
% legend("y_1(t)","y_2(t)","y_3(t)") % label the graphs in the order they were plotted
% title("Graph of y1, y2, and y3 versus t.") % indicate that this is y1, y2, and y3 versus t