ECE356S – Linear Systems and Control –

EXPERIMENT 1

Introduction to MATLAB and Simulink for Systems Analysis

1 Purpose

The purpose of this experiment 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.

2 Preparation

Before you go to your lab session, read through carefully the description of this laboratory provided in

the following pages. Useful terms and commands are highlighted in bold font, while laboratory work is

written in bold italics. Identify the parts you have to do in the lab.

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,

A = [1 2 3; 4 5 6; 7 8 9]

produces the matrix

A =

1 2 3

4 5 6

7 8 9

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.

I. Inverse of a Matrix and Solution of Linear Equations:

MATLAB provides a command to compute the inverse of an n × 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 x = A\b. The MATLAB operation “\” is called mldivide or

matrix left divide. This is illustrated in the following calculation.

(a) Save the following code to a file and name it inv matrix.m. This is a script file which you can

1

run in MATLAB.

n = 500;

Q = orth(randn(n,n));

d = logspace(0,-10,n);

A = Q*diag(d)*Q’;

x = randn(n,1);

b = A*x;

tic, y = inv(A)*b; toc

err = norm(y-x)

res = norm(A*y-b)

pause

tic, z = A\b; toc

err = norm(z-x)

res = norm(A*z-b)

Consult the MATLAB documentation to understand what each line means. Then

run “inv matrix” in MATLAB and note the outputs.

(b) Similarly, there is an operation called mrdivide or matrix right divide, A/B, which has the

effect of AB−1 when B is invertible. Let A = randn(4, 4). Compute the inverse of A using

inv, mldivide, and mrdivide. Show your results to your TA. Explain how you would

compare the accuracy of the three methods. For square matrices with low dimension, all

three calculations should give comparable accuracy.

II. Eigenvalues and Eigenvectors

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.

(a) Let A be given by

A =

7 2 −3

4 6 −4

5 2 −1

We will use this A from item (a) to (e) in this part. Use the MATLAB command [V, D] =

eig(A) to determine the eigenvalues and eigenvectors of A.

(b) The command eig(A) gives eigenvectors which are normalized to have norm 1. If you compute the eigenvectors by hand, you are more likely to come up with eigenvectors which are not

normalized. For the matrix in II(a), you can use the command eig(A,’nobalance’)

to produce more “readily recognizable” eigenvectors. Recall that eigenvectors are determined up to a scalar multiple. From the results of eig(A,’nobalance’), determine 3

eigenvectors of A whose entries are all integers.

(c) 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 (a).

(d) The eigenvalue-eigenvector equations can be written as one matrix equation

AV = V D,

2

with D a diagonal matrix consisting of the eigenvalues of A. From the results of IV(a),

determine norm(AV-VD). Show more significant digits in MATLAB using the command format long. From this determine the exact values of the eigenvalues and

eigenvectors (up to a scalar multiple). Now verify that AV − V D = 0.

(e) 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 II(a) can be diagonalized.

(f) 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

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. Try using eig on this A. Does the resulting [V,D] satisfy

AV = V D? Is V invertible?

(g) When A cannot be diagonalized, one can transform it to a form called the Jordan (canonical)

form. The MATLAB command is jordan. For the matrix

A =

0 4 3

0 20 16

0 −25 −20

find its Jordan form.

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 (1)

y = Cx + Du, (2)

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

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

MATLAB provides tools to solve such linear systems. We illustrate some of these tools using a second

order system.

3

Consider a second order system described by the differential equation

y¨ + 2 ˙y + 4y = 4u

The transfer function from the input u to the output y for this system is given by

G(s) = 4

s

2 + 2s + 4

(3)

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 (4)

y =

1 0

x (5)

MATLAB provides a data object which conveniently describes the state space system (1)-(2). It is created

by the command sys=ss(A,B,C,D).

(a) Create the sys object corresponding to the second order system described by (4)-(5).

(b) 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.

(c) To determine the response due only to initial conditions, it would be convenient to create a second

sys object using sys init=ss(A,0,C,0) (0 is a matrix of appropriate dimension). Use the command

[Y,T,X]=initial(sys init,x0) to determine the response to initial conditions when x0 =

0

1

. Again plot the state responses.

(d) In general, the system may have both nonzero initial conditions as well as inputs. For example, the

commands

t=0:0.01:20;

u=sin(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.

(e) If you are given a state space description, the command ss2tf generates the numerator and denominator polynomials of the transfer function associated with the state space system. Use the ss2tf

command, followed by the tf command to produce the transfer function for the systems

described by (4)-(5). Verify that it is the same as (3).

4