Introduction

to Matlab Prof. Dr. Y. Samim ÜNLÜSOY

Mech. Eng. Dept.,
Prof. Dr. Y. Samim Middle East Technical
ME 513 Vehicle Dynamics 1
Ünlüsoy University, Ankara
 MATLAB [Link]
MATLAB

Double-click
Load or save files here

stored
Watch your Write your
values here commands here

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 2

 INTRODUCTION TO MATLAB
„ Matrices
„ In MATLAB, a matrix is a rectangular array of numbers.
‹ Scalars are 1-by-1 matrices.
‹ Vectors are matrices with only one row or column.
„ MATLAB has other ways of storing both numeric and
nonnumeric data; but in the beginning, it is usually best to
think of everything as a matrix.
„ Where other programming languages work with numbers
one at a time, MATLAB allows you to work with entire
matrices quickly and easily.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 3

 HOW TO ENTER MATRICES
„ You can enter matrices into MATLAB in
several different ways:
‹ Enter an explicit list of elements.
‹ Load matrices from external data files.
‹ Generate matrices using built-in
functions.
‹ Create matrices with your own
functions in M-files.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 4

 Entering a Matrix
„ You have only to follow a few basic
conventions:
‹ Separate the elements of a row
with blanks or commas.
‹ Use a semicolon, ; ,to indicate the
end of each row.
‹ Surround the entire list of
elements with square brackets, [ ].

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 5

 Entering A Matrix
„ Start by entering the matrix as a list of its
elements. To enter, simply type in the
Command Window
‹ A = [16 3 2 13; 5 10 11 8; 9 6 7 12;
4 15 14 1]
„ Once you have entered the matrix, it is
automatically remembered in the MATLAB
workspace. You can refer to it simply as A.
Now that you have A in the workspace.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 6

 MatLab Command Window
>> A = [16 3 2 13; 5 10 11 8; 9 6 7 12;
4 15 14 1]

A=

16 3 2 13
5 10 11 8
9 6 7 12
4 15 14 1

>>
Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 7
 MatLab Command Window
Double click here to see
what is stored in A.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 8

 Transpose of a Matrix
„ There are two possibilities to transpose a
matrix.
‹ Use an apostrophe or single quote, '.
‹ Type transpose(A), where A is the
matrix transpose of which is sought.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 9

A= >> A'
MatLab ans =
Command
16 3 2 13
5 10 11 8
9 6
4 15
7 12
14 1
Window 16 5
3 10
9 4
6 15
2 11 7 14
13 8 12 1

>> transpose(A)

ans =

16 5 9 4
3 10 6 15
2 11 7 14
13 8 12 1

>>

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 10

 Inverse of a Matrix

„ inv(X) is the inverse of the square

matrix X.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 11

C= >> C

16 3 2 13 C=
5 10 11 8
9 6 15 12 16 3 2 13
4 15 14 1 5 10 11 8
9 6 15 12
4 15 14 1

MatLab >> B=inv(C)

Command
Window
C=

0.0760 -0.1985 0.0417 0.1005

0.0196 0.1544 -0.1250 0.0098
-0.0417 -0.1250 0.1250 0.0417
-0.0147 0.2279 -0.0417 -0.1324

>>

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 12

 Subscripts
„ The element in row i and column j of A is
denoted by A(i,j).
„ For example, A(4,2) is the number in the
fourth row and second column.
„ So it is possible to compute the sum of the
elements in the fourth column of A by
typing
‹ A(1,4) + A(2,4) + A(3,4) + A(4,4)

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 13

 „ The expression
Colon Operator
1:10
is a row vector containing the integers from 1 to 10
1 2 3 4 5 6 7 8 9 10

„ To obtain nonunit spacing, specify an increment.

100:-7:50
is
100 93 86 79 72 65 58 51
and
0:pi/4:pi
is
0 0.7854 1.5708 2.3562 3.1416

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 14

>> 1:10

ans = MatLab
1 2 3 4 5 6 7 8 9 10 Command
>> 100:-7:50
Window
ans =

100 93 86 79 72 65 58 51

>> 0:pi/4:pi

ans =

0 0.7854 1.5708 2.3562 3.1416

>>

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 15

 Colon Operator
„ Subscript expressions involving colons refer to
portions of a matrix.
A(1:k,j)
is the first k elements of the jth column of A. So
sum(A(1:4,4))
computes the sum of the fourth column.
The colon by itself refers to all the elements in a row
or column of a matrix and “end” refers to the last row
or column. So
sum(A(:,end))
computes the sum of the elements in the last column
of A.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 16

A= >> sum(A(1:4,4))
MatLab
Command
16 3 2 13
5 10 11 8 ans =
9 6
4 15
7 12
14 1
Window
34

>> sum(A(:,end))

ans =

34

>>

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 17

 Matrix Multiplication and Division
„ C = A*B
is the linear algebraic product of the matrices
A and B. For nonscalar A and B, the number
of columns of A must equal the number of
rows of B. A scalar can multiply a matrix of
any size.
„ .* Array multiplication.
A.*B is the element-by-element product of the
arrays A and B. A and B must have the same
size, unless one of them is a scalar.
Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 18
A= >> A*B
MatLab Command
Window
16 3 ans =
5 10

B= 154 195
175 70 16*7+3*14=154
7 12 16*12+3*1=195
14 1 >> A.*B 5*7+10*14=175
5*12+10*1=70
ans =

112 36 16*7=112
70 10 3*12=36
5*14=70
>> 10*1=10

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 19

 Matrix Multiplication and Division

„ / matrix right division.

B/A is roughly the same as B*inv(A). More
precisely, B/A = (A'\B')'.
A./B is the matrix with elements A(i,j)/B(i,j). A
and B must have the same size, unless one of
them is a scalar.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 20

 Matrix Multiplication and Division

„ \ : matrix left division.

If A is a square matrix, A\B is roughly the
same as inv(A)*B, except it is computed in a
different way. If A is an n-by-n matrix and B is
a column vector with n components, or a
matrix with several such columns, then X =
A\B is the solution to the equation AX = B
computed by Gaussian elimination.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 21

A= >> B/A
MatLab
Command
16 3 ans =
5 10

B= 0.0690 1.1793 Window

0.9310 -0.1793
7 12
14 1 >> B*inv(A)

ans =

0.0690 1.1793
0.9310 -0.1793

>>

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 22

A= >> A./B
MatLab
Command
16 3 ans =
5 10

B= 2.2857 0.2500 Window

0.3571 10.0000
7 12
14 1 >> A\B >> inv(A)*B

ans = ans =
0.1931 0.8069
1.3034 -0.3034 0.1931 0.8069
1.3034 -0.3034
>> >>

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 23

 Eigenvalues and Eigenvectors
„ E = eig(X) is a vector containing the
eigenvalues of a square matrix X.

„ [V,D] = eig(X) produces a diagonal

matrix D of eigenvalues and a full matrix
V whose columns are the corresponding
eigenvectors so that
X*V = V*D.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 24

A= >> A

16 3 2 13 MatLab A=
5 10 11 8 Command
9 6
4 15
15 12
14 1 Window 16 3 2 13
5 10 11 8
9 6 15 12
4 15 14 1

>> G=eig(A)

G=

36.3739
-7.1081
8.0000
4.7341

>>

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 25

 Eigenvalues and Eigenvectors
„ E = eig(A,B) is a vector containing the
generalized eigenvalues of square matrices A
and B.

„ [V,D] = eig(A,B) produces a diagonal matrix D

of generalized eigenvalues and a full matrix V
whose columns are the corresponding
eigenvectors so that

A*V = B*V*D.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 26

 Variables
„ MATLAB does not require any type
declarations or dimension statements. A
new variable name is automatically created
and appropriate amount of storage is
allocated. If the variable already exists,
MATLAB changes its contents and, if
necessary, allocates new storage.
„ Variable names consist of a letter, followed
by any number of letters, digits, or
underscores.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 27

 Variables
„ MATLAB uses only the first 31 characters
of a variable name.
„ MATLAB is case sensitive; it distinguishes
between uppercase and lowercase letters.
A and a are not the same variable.
„ To view the matrix assigned to any variable,
simply enter the variable name and press
enter.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 28

 Functions
„ MATLAB provides a large number of
standard elementary mathematical
functions, such as
‹ abs,
‹ sqrt,

‹ exp, and

‹ sin, cos, tan etc..

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 29

 Advanced Mathematical Functions
„ Taking the square root or logarithm of a
negative number is not an error; the
appropriate complex result is produced
automatically.
„ MATLAB also provides many more
advanced mathematical functions,
including Bessel and gamma functions.
Most of these functions accept complex
arguments.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 30

 Help on Available Functions
„ For a list of the elementary mathematical
functions, type
help elfun
„ For a list of more advanced mathematical
and matrix functions, type
help specfun
help elmat

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 31

 Built-in and m-file Functions
„ Some of the functions, like sqrt and sin,
are built-in. They are part of the MATLAB
core so they are very efficient, but the
computational details are not readily
accessible.
„ Other functions, like gamma and sinh, are
implemented in m-files. You can see the
code and even modify it if you want.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 32

 PLOTS >>help plot
„ PLOT(X,Y) plots vector Y versus vector X.
„ If X or Y is a matrix, then the vector is plotted
versus the rows or columns of the matrix,
whichever line up.
„ If X is a scalar and Y is a vector, length(Y)
disconnected points are plotted.
PLOT(Y) plots the columns of Y versus their index.
If Y is complex, PLOT(Y) is equivalent to
PLOT(real(Y),imag(Y)).
In all other uses of PLOT, the imaginary part is
ignored.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 33

 PLOTS >>help plot
„ Various line types, plot symbols and colors may be obtained with
PLOT(X,Y,S) where S is a character string made from one element
from any or all the following 3 columns:

b blue . point - solid

g green o circle : dotted
r red x x-mark -. dashdot
c cyan + plus -- dashed
m magenta * star
y yellow s square
k black d diamond

„ For example, PLOT(X,Y,'c+:') plots a cyan dotted line with a plus at

each data point; PLOT(X,Y,'bd') plots blue diamond at each data
point but does not draw any line.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 34

 PLOTS >> help plot
„ PLOT(X1,Y1,S1,X2,Y2,S2,X3,Y3,S3,...)
combines the plots defined by the (X,Y,S)
triples, where the X's and Y's are vectors or
matrices and the S's are strings.

„ For example, PLOT(X,Y,'y-',X,Y,'go') plots

the data twice, with a solid yellow line
interpolating green circles at the data
points.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 35

 PLOTS >> help plot
„ The PLOT command, if no color is specified,
makes automatic use of the colors specified by the
axes ColorOrder property. The default ColorOrder
is listed in the table above for color systems where
the default is blue for one line, and for multiple
lines, to cycle through the first six colors in the
table. For monochrome systems, PLOT cycles
over the axes LineStyleOrder property.

„ See also SEMILOGX, SEMILOGY, LOGLOG,

PLOTYY, GRID, CLF, CLC, TITLE, XLABEL,
YLABEL, AXIS, AXES, HOLD, COLORDEF,
LEGEND, SUBPLOT, STEM.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 36

 Solution of Linear Ordinary Differential
Equations with Constant Coefficients
„ >> help tf

„ tf : Creation of transfer functions or conversion to transfer

function.

‹ sys = tf(num, den)

„ creates a continuous-time transfer function sys with numerator
num and denominator den. The output sys is a tf object.
‹ s = tf('s')

„ specifies the transfer function H(s) = s (Laplace variable).

„ You can then specify transfer functions directly as rational

expressions in s, e.g.,
‹ s = tf('s'); H = (s+1)/(s^2+3*s+1)

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 37

 Time Response for LTI Systems : lsim

Matlab command lsim (sys, u, t) is used to

simulate time response of LTI models to
arbitrary inputs described by u.
The time vector t consists of regularly
spaced time samples and u is a matrix with
as many columns as inputs and whose i-th
row specifies the input value at time t(i).

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 38

 Time Response for LTI Systems : lsim
lsim(sys,u,t,x0) specifies the initial state vector x0
at time t(1) (for state-space models only). x0 is set
to zero when omitted.
lsim(sys1, sys2,..., u, t, x0) simulates the response
of multiple LTI models sys1, sys2,... on a single
plot.
You can also specify a color, line style, and
marker for each system, as in
lsim(sys1,'r',sys2,'y--',sys3,'gx',u,t).

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 39

 Time Response for LTI Systems : lsim
For example,
s = tf('s'); H = (s+1)/(s^2+3*s+1)
t = 0:0.01:10;
u = sin(t);
lsim(H,u,t) ;
simulates the response of a single-input model
represented by the transfer function H to the input
u(t)=sin(t)
for 10 seconds by 0.01 second intervals.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 40

 Time Response for LTI Systems : lsim
Linear Simulation Results
1

0.8

0.6

0.4

0.2
Amplitude

-0.2

-0.4

-0.6

-0.8

-1
0 1 2 3 4 5 6 7 8 9 10
Time (sec)

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 41

 m-files
„ There are two types of m-files :
‹ Script m-files, and
‹ Function m-Files
„ M-files can be scripts that simply execute a
series of MATLAB statements, or they can be
functions that also accept arguments and
produce output.
„ You create M-files using a text editor, then use
them as you would any other MATLAB
function or command.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 42

 m-files
„ Script M-Files
‹ Do not accept input arguments
or return output arguments.
‹ Operate on data in the
workspace.
‹ Useful for automating a series
of steps you need to perform
many times.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 43

 Script M-Files
Consider a one degree of freedom mass-spring-damper
system representing a road vehicle travelling on a
rough road surface.
x
Obtain the motion of the vehicle
m in terms of the vehicle body mass
displacement y and its derivatives.
c k
z
Assume that the road surface profile z
is sinusoidal with a frequency of 5 rad/s
and an amplitude of 1 cm.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 44

 m
x Script M-Files
c k
The equation of motion is given as :
z
 + cx + kx = cz + kz
mx
The transfer function is then given by :

X(s) cs + k
G(s) = =
Z(s) ms2 + cs + k

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 45

 Script M-Files
% March 2004
% Response of a vehicle travelling on a
% sinusoidal road surface profile.
% -------------------------------------------------
% Define Parameters and Variables
% m : vehicle mass [kg]
% c : total damping coefficient of suspension dampers [m/s]
% k : total spring constant of suspension springs [kN/m]
% z : vehicle body displacement [m]
% -------------------------------------------------
% Enter Data Values
m=1260;
c=4000;
k=60000;
.....

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 46

..... Script M-Files
% ---------------------------------------------------------
% Generate the time array for the solution
t=0:0.01:10; % time array (from 0 to 10 seconds with 0.01
second increments)
% -----------------------------------------------
% Generate the road surface profile
A=0.01; % road surface profile amplitude [m]
omega=5; % road surface profile frequency [rad/s]
z=A*sin(omega*t); % road surface profile [m]
% ----------------------------------------
% Plot the road surface profile
figure(1)
plot(t,z)
title(‘Sinusoidal Road Surface Profile - Amplitude : 1 cm, Frequency : 5
rad/s’);
xlabel(‘Time[s]’);ylabel(‘y[m]’);
.....

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 47

.....
Script M-Files
% -------------------------------------------
% Generate the transfer functions
s=tf(‘s’); % Declare s as the Laplace variable
% from z to y
hz=(c*s+k)/(m*s^2+c*s+k);
% from z to ydot (velocity)
hzdot=s*hz;
% from z to ydoubledot (acceleration)
hzddot=s*hzdot;
% -----------------------------------------------
% Plot the displacement, velocity, and acceleration of
% the vehicle body
figure(2)
% Plot of displacement – note that hz is the TF from z to y
lsim(hz,z,t)
.....

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 48

 Script M-Files
.....
figure(3)
% Plot of velocity – note that hzdot is the TF from z to
ydot
lsim(hzdot,z,t)
.....

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 49

 Script M-Files
- Click on “File”, then “New”, and “M-file”.
- A window to enter m-script file will open.
- Write or copy Matlab statements given in the
example.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 50

 Script M-Files
- Click “here” to save
and immediately
run the m-script
file.

- Click “here” just to

save the m-script
file.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 51

 Sinusoidal Road Surface Profile - Amplitude : 1 cm, Frequency : 5 rad/s
0.01

0.008

0.006

Script
m-files
0.004

0.002
y [m]

-0.002
This is what
-0.004
you get after
-0.006 running the m-
-0.008
script file with
program
-0.01
0 1 2 3 4 5
Time [s]
6 7 8 9 10 defaults !

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 52

 Sinusoidal Road Surface Profile - Amplitude : 1 cm, Frequency : 5 rad/s
0.02

Script
0.015

0.01

m-files
ydot [m/s]

0.005

-0.005
You can
-0.01 modify the
plots using the
-0.015
edit facility on
-0.02 0 1 2 3 4 5 6 7 8 9 10 the Matlab plot
Time [s] (sec)

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 53

 Sinusoidal Road Surface Profile - Amplitude : 1 cm, Frequency : 5 rad/s

0.1

Script
0.08

0.06

m-files
ydoubledot [m/s2]

0.04

0.02

-0.02
You can
-0.04 modify it
-0.06 using the edit
-0.08
facility on the
Matlab plot
-0.1
0 1 2 3 4 5 6 7 8 9 10
Time [s] (sec)

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 54

 m-files
„ Function m-Files
‹ Can accept input arguments
and return output arguments.
‹ Internal variables are local to
the function by default.
‹ Useful for extending the
MATLAB language for your
application.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 55

 Function m-files
function [Fx,Fy]=tireforces(munom,Fx,Fz,alpha,slip,vwl)
.....
.....
.....

Call from the main program body (where

munoml,Fxfl,Fzfl,alphafl,slipfl,vwlfl are assigned values for which
Fxfl,Fyfl are to be calculated) :

[Fxfl,Fyfl]=tireforces(munoml,Fxfl,Fzfl,alphafl,slipfl,vwlfl);
and values of Fxfl,Fyfl calculated in the function will be returned.

Prof. Dr. Y. Samim Ünlüsoy ME 304 Control Systems 56