Description
1 Purpose
The purpose of this laboratory is to introduce you to MATLAB and Simulink as tools
for systems analysis. MATLAB is used primarily for numerical computations, while
Simulink provides a block diagram style representation of systems for analysis, simulation, and design. Please attach a single MATLAB script for each separate task. This
will not only simplify the work of the TAs in grading your work but will also help you
with building the habit of having scripts that automate your work. Furthermore save
these files so that they can be a useful reference for you in the future assignment and
courses. At the end of this file you will find a list with all the files you are required to
submit.
2 Matlab scripts
Read the MATLAB guide about how to create scripts here!
3 Linear Algebra
MATLAB is heavily based on linear algebra for doing its computations. In MATLAB,
matrices are entered row by row. The end of a row is indicated by a “;” or a carriage
return. For example:
1 A = [ 1 2 3 ; 4 5 6 ; 7 8 9 ]
produces the matrix:
A =
1 2 3
4 5 6
7 8 9
(1)
Individual rows and columns can be extracted from A. For example, A(:, 1) and
A(1, 🙂 give the first column and row of A, respectively.
• How would you select the third column?
1
• Without trying on MATLAB, what do you think would be the result of A(2,2)?
Check the result on MATLAB.
3.1 Inverse of a Matrix and Solution of Linear Equations:
You will create and upload a unique script for all the calculations you operate
in this section 3.1. Name this file ex 3 1.m and add comments where required.
MATLAB provides a command to compute the inverse of an n x n matrix A, namely
inv(A). However, if the objective is to solve a linear system of equations Ax = b, there is
a better alternative than using x = inv(A) * b, namely A\b. The MATLAB operation
“\” is called mldivide or matrix left divide. This is illustrated in the following calculation.
1. Save the following code to a file and name it inv matrix.m. This is a script file
which you can run in MATLAB. You can use the lines below as an example
of how your submitted files should look like.
1 n = 5 0 0; % dimen si on o f the ma trix
2 Q = o r t h (randn (n , n) ) ;
3 d = l o g s p a c e (0 , −10 ,n) ;
4 A = Q∗ di a g (d) ∗Q’ ;
5 x = randn (n , 1 ) ; % the s o l u t i o n o f the equ a ti on , known t o
compare methods
6 b = A∗x ;
7 ti c , y = in v (A) ∗b ; t o c % s o l v e e q u a ti o n with in v and time
8 e r r i n v = norm(y−x ) % compute e r r o r
9 r e s i n v = norm(A∗y−b)
10 pause
11 ti c , z = A\b ; t o c % s o l v e e q u a ti o n with ml di vi d e and time
12 e r r ml = norm( z−x ) % compute e r r o r
13 res m ; = norm(A∗z−b)
Consult the MATLAB documentation to understand what each line
means. Then run inv matrix in MATLAB and note the outputs.
2. Similarly, there is an operation called mrdivide or matrix right divide, A/B, which
has the effect of solving the equation XA = B with respect to X. This method
can be used only when A and B have the same number of columns. In this exercise
we will compute the solution of the equation AX = B using all the three comands
inv, mldivide, and mrdivide, in a similar way as above.
Define A = randn(4, 4) and also X = randn(4, 4). Find the solution to the
equation AX = B using the commands inv, mldivide, and mrdivide. In order
to use mrdivide successfully, notice that solving AX = B with respect to X is the
same as XT AT = BT were AT denotes the transpose of A.
Be sure that the code used for the comparison of the three results is included in
the submitted file ex 3 1.m. For square matrices with low dimension, all three
calculations should give comparable accuracy.
2
3.2 Eigenvalues and Eigenvectors
You will create and upload a unique script for all the calculations you operate
in this section 3.2. Name this file ex 3 2.m and add comments where required.
Of great importance in the analysis of linear systems are the notions of eigenvalues and
eigenvectors. A scalar λ is an eigenvalue and a non-zero vector x is an eigenvector of a
linear transformation A if Ax = λx. If A is a matrix, then eigenvalues and eigenvectors
are only defined if A is square.
1. Let A be given by
A =
7 2 −3
4 6 −4
5 2 −1
(2)
We will use this A for items 1 to 4 in this part. Use the MATLAB command
[V, D]=eig(A) to determine the eigenvalues and eigenvectors of A.
2. For manual computation of eigenvalues, usually you determine the characteristic
polynomial of A given by det(sI − A). Then you find the roots of the characteristic polynomial, i.e., find solutions of the characteristic equation det(sI − A) = 0.
Determine the characteristic polynomial of A using the command poly(A), and
determine the eigenvalues by applying the command roots to the resulting polynomial. Compare the answer with those from 1.
3. The eigenvalue-eigenvector equations can be written as one matrix equation
AV = V D, (3)
with D a diagonal matrix consisting of the eigenvalues of A. From the results
of point 1, determine norm(AV-VD). Show more significant digits in MATLAB
using the command format LONG before the expression you evaluate. From this
determine the exact values of the eigenvalues and eigenvectors (up to a scalar
multiple). Now verify that AV − V D = 0.
4. Recall that if the eigenvalues of A are distinct, the eigenvectors are linearly independent. In that case, the V matrix computed using eig is invertible and that
V
−1AV = D. This process is called diagonalizing A. Check that the matrix
in (2) can be diagonalized. Write one line in the comments justifying why it is
diagonalizable or not.
5. When A has repeated eigenvalues, it may not be possible to diagonalize A. Suppose
A has λ as a repeated eigenvalue, then det(λI − A) = 0 and the number of times
λ repeats as roots is called its algebraic multiplicity. Thus the following matrix
A =
1 1
0 1
(4)
3
has 1 as its only eigenvalue with algebraic multiplicity 2. Determine by hand
calculation the eigenvector corresponding to the eigenvalue 1, and verify there is
only one independent eigenvector. Report this eigenvector in the comments. Try
using eig on this A. Does the resulting [V,D] satisfy AV = V D? Is V invertible?
Report the answers in the comments.
6. When A cannot be diagonalized, one can transform it to a form called the Jordan
(canonical) form. The MATLAB command is jordan. Find the Jordan form for
the matrix:
A =
0 4 3
0 20 16
0 −25 −20
(5)
4 Ordinary Differential Equations and Transfer Functions
Control systems studied in this course are modelled by in the time domain by state
equation of the form:
x˙ = Ax + Bu
y = Cx + Du (6)
or by input/output models of the form
y
(n)
(t)+a1y
(n−1)(t)+a2y
(n−2)(t)+· · ·+any(t) = b0u
(m)
(t)+· · ·+bmu(t), with m ≤ n
(7)
In the state equation, the transfer matrix from u to y is given by:
G(s) = C(sI − A)
−1B + D
In this course, we focus on single-input single-output (SISO) systems, i.e., u and y are
scalar-valued. In this case, the transfer function G(s) is a scalar-valued proper rational
function. In the case of an input/output model, the transfer function G(s) is given by
G(s) = b0s
m + b1s
m−1 + · · · + bm
s
n + a1s
n−1 + · · · + an
(8)
MATLAB provides tools to analyze such linear systems. We illustrate some of these
tools using a second order system.
4.1 ODE and transfer functions
You will create and upload a unique script for all the calculations you operate
in this section 4.1. Name this file ex 4 1.m and add comments where required.
You will need to upload the required plots as well. Consider a second order system
described by the differential equation
y¨ + 2 ˙y + 4y = 4u
4
The transfer function from the input u to the output y for this system is given by
G(s) = 4
s
2 + 2s + 4
(9)
To model this as a state equation, let
x(t) =
y(t)
y˙(t)
Then the state x satisfies the differential equation
x˙ =
0 1
−4 −2
x +
0
4
u (10)
y =
1 0
x (11)
MATLAB provides a data object which conveniently describes the state space system
(6). It is created by the command sys=ss(A,B,C,D).
1. Create the sys object corresponding to the second order system described by (10)-
(11).
2. One can study the response of the second order system due to particular inputs
such as a step (to get the step response), or the response due to nonzero initial
conditions with no input, or generally with nonzero initial conditions and inputs.
Here, determine the step response using the command [Y,T,X]=step(sys), where
sys is the object you created using the ss command. Use plot(T,X) to plot the
trajectories of both states. Save the plot1 as .png file and name it plot 4 2.png.
3. Use the command [Y, T, X]= initial (sys, x0) to determine the response to
initial condition x0 =
0
1
and zero input. Again plot the state responses and
save it as a .png file and name it plot 4 3.png.
4. In general, the system may have both nonzero initial conditions as well as inputs.
For example, the commands:
1 t = 0 : 0 . 0 1 : 2 0 ;
2 u=s i n ( t ) ;
define an input vector u whose components are given by the values of sin(t) at
various times specified by the time vector t. Use [Y,t,X]=lsim(sys,u,t,x0) to
produce the response of the second order system due to the initial condition and
sinusoidal input defined above. Plot the state trajectories and save it as a .png
file and name it plot 4 4.png.
1 We highly discourage using the screenshot function of your OS. Use instead the MATLAB builtin commands for exporting plots as images. For example MATLAB documentation provides useful
information. A simple way to save your plot is (in the plot window): File -> Export as. In MATLAB
you can control virtually everything about your plots, try for example to add a meaningful legend,
it will be useful for the next labs. If you have problems contact your lab TAs.
5
