% IAP2007 Introduction to MATLAB: Calculus, Linear Algebra, Differential Equations
% Instructor: Violeta Ivanova, violeta@mit.edu

% Example 1: Linear Systems, Eigenvalues & Eigenvectors: RCL Circuit
% This exercise is based on Prof. Stephen Hall's Signals and Systems lectures 
% for 16.01-16.04 Unified Engineering.

% NOTE: consult this session's handout for illustrations (slides 21-13).

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

% A. STATE EQUATION

% A1. For state equation x' = Ax, and (column) state vector x = (v; i), 
% create the matrix A for the following system:
% dv/dt = -1/2*v - 2*i
% di/dt =  1/2*v - 3*i

A = [-1/2   -2;     1/2     -3]

% A2. Define a (column) vector x0 for initial conditions: v(0) = 2 and i(0) = 1

x0 = [2; 1] 

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

% B. EIGENVALUES & EIGENVECTORS

% B1.  Compute eigenvalues and eigenvectors of the matrix A
% Hint: use function EIG with option 'nobalance' for non-symmetric A

[eigvecs, eigvals] = eig(A, 'nobalance')

% B2. Compute a column vector lambda including the two eigenvalues
% Hint: use function DIAG.

lambda = diag(eigvals)

% B3. Define the two column eigvectors eigvec1 and eigvec2:

eigvec1 = eigvecs(:, 1)
eigvec2 = eigvecs(:, 2)

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

% C. SOLUTION: v(t) and i(t)
      
% C1. Compute the coefficients' vector A from the initial conditions x0
% Hint: Use operator \ to solve equation x0 = eigenvectors * A for A.

a = eigvecs \ x0

% C2. Compute v(t) and i(t) for 0 < t < 5s
% Hint: the solution has an exponential form - see handout for formula.

t = [0 : 0.1 : 5];
v = a(1) * eigvec1(1) * exp(lambda(1)*t) + a(2) * eigvec2(1) * exp(lambda(2)*t);
i = a(1) * eigvec1(2) * exp(lambda(1)*t) + a(2) * eigvec2(2) * exp(lambda(2)*t);

% C3. Plot v(t) and i(t) vs. time t on the same plot. 

pv = plot(t, v);
hold on
pi = plot(t, i, 'r-');
title('RCL CIRCUIT');
xlabel('Time');
legend('Voltage', 'Current');
hold off