MCEN30020

Systems Modelling & Analysis

Simon Illingworth

2021 semester 1
Contents

1 Mechanical and electrical systems 5

1.1 The mass-spring damper system . . . . . . . . . . 6
1.2 Rotating systems . . . . . . . . . . . . . . . . . . 9
1.3 Electrical circuits . . . . . . . . . . . . . . . . . . 10
1.4 A new device . . . . . . . . . . . . . . . . . . . . 14
1.5 Summary of elements . . . . . . . . . . . . . . . . 16
1.6 Links and further reading . . . . . . . . . . . . . 17
1.A Energy of the mass-spring system . . . . . . . . . 19

2 Linear state-space models 22

2.1 The state x . . . . . . . . . . . . . . . . . . . . . 22
2.2 The matrix exponential . . . . . . . . . . . . . . . 24
2.3 Eigendecomposition of eAt . . . . . . . . . . . . . 26
2.4 Stability and frequency of oscillation . . . . . . . 28
2.5 Complex eigenvalues come in pairs . . . . . . . . 29
2.6 Including inputs . . . . . . . . . . . . . . . . . . . 29
2.7 Including outputs . . . . . . . . . . . . . . . . . . 31
2.8 Links and further reading . . . . . . . . . . . . . 32

3 Linearization 33
3.1 Finding equilibria . . . . . . . . . . . . . . . . . . 34
3.2 Linearization . . . . . . . . . . . . . . . . . . . . 35
3.3 From one state to many states . . . . . . . . . . . 37
3.4 Linearization with inputs . . . . . . . . . . . . . . 38

1
Contents 2

3.5 Cart-pendulum example . . . . . . . . . . . . . . 39

3.6 Links and further reading . . . . . . . . . . . . . 43
3.A Cart-pendulum equations of motion . . . . . . . . 43

4 The Laplace transform 46

4.1 Life without the Laplace transform . . . . . . . . 46
4.2 Definition of the Laplace transform . . . . . . . . 47
4.3 Laplace transform of common functions . . . . . . 47
4.4 Properties of the Laplace transform . . . . . . . . 50
4.5 The inverse Laplace transform . . . . . . . . . . . 54
4.6 Partial-fraction expansions . . . . . . . . . . . . . 54
4.7 Understanding complex roots . . . . . . . . . . . 58
4.8 Matlab commands . . . . . . . . . . . . . . . . . 60
4.9 Links and further reading . . . . . . . . . . . . . 60

5 Second-order systems 61
5.1 Standard form . . . . . . . . . . . . . . . . . . . . 61
5.2 Laplace transforms . . . . . . . . . . . . . . . . . 62
5.3 Poles . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.4 Free response . . . . . . . . . . . . . . . . . . . . 64
5.5 Transfer function . . . . . . . . . . . . . . . . . . 68
5.6 Impulse response . . . . . . . . . . . . . . . . . . 69
5.7 Step response . . . . . . . . . . . . . . . . . . . . 70
5.8 Characterizing the underdamped step response . . 73
5.9 Links and further reading . . . . . . . . . . . . . 77

6 First-order systems 81
6.1 Cup of tea as a first-order system . . . . . . . . . 81
6.2 Laplace transforms . . . . . . . . . . . . . . . . . 82
6.3 Poles . . . . . . . . . . . . . . . . . . . . . . . . . 83
6.4 Free response . . . . . . . . . . . . . . . . . . . . 83
6.5 Transfer function . . . . . . . . . . . . . . . . . . 84
6.6 Impulse response . . . . . . . . . . . . . . . . . . 84
6.7 Step response . . . . . . . . . . . . . . . . . . . . 85
Contents 3

6.8 Matlab code . . . . . . . . . . . . . . . . . . . . . 85

7 Transfer functions and block diagrams 87

7.1 Summary: first- and second-order systems . . . . 87
7.2 Transfer functions . . . . . . . . . . . . . . . . . . 88
7.3 Initial-value & final-value theorems . . . . . . . . 91
7.4 Non-minimum-phase transfer functions . . . . . . 94
7.5 Block diagrams . . . . . . . . . . . . . . . . . . . 97
7.6 Links and further reading . . . . . . . . . . . . . 98

8 Fluid and thermal systems 99

8.1 Fluid systems . . . . . . . . . . . . . . . . . . . . 99
8.2 Thermal systems . . . . . . . . . . . . . . . . . . 107
8.3 Links and further reading . . . . . . . . . . . . . 110

9 Bode plots 111

9.1 Introduction . . . . . . . . . . . . . . . . . . . . . 111
9.2 The frequency response . . . . . . . . . . . . . . . 114
9.3 The Bode plot . . . . . . . . . . . . . . . . . . . . 117
9.4 Bode plots by hand . . . . . . . . . . . . . . . . . 118
9.5 Second-order system example . . . . . . . . . . . 121
9.6 Important properties of Bode plots . . . . . . . . 123
9.7 General procedure for plotting Bode plots . . . . 125
9.8 Bode plots of common transfer functions . . . . . 126
9.A Derivation of the frequency response . . . . . . . 132
9.B The decibel . . . . . . . . . . . . . . . . . . . . . 133

10 Nyquist plots 135

10.1 Nyquist plots . . . . . . . . . . . . . . . . . . . . 135
10.2 Nyquist plots of common transfer functions . . . . 136

11 Fourier analysis 142

11.1 Fourier series . . . . . . . . . . . . . . . . . . . . 142
11.2 Complex-exponential form . . . . . . . . . . . . . 144
Contents 4

11.3 Relationship between them . . . . . . . . . . . . . 145

11.4 Example: Square wave . . . . . . . . . . . . . . . 145
11.5 Application: Pulse Width Modulation (PWM) . . 149
11.6 Application: PWM of a DC motor . . . . . . . . 153
11.7 Links and further reading . . . . . . . . . . . . . 160
11.A The discrete Fourier transform (DFT) . . . . . . . 161
11.B The fast Fourier transform (FFT) . . . . . . . . . 162
11.C Sampling frequency . . . . . . . . . . . . . . . . . 163
11.D Examples . . . . . . . . . . . . . . . . . . . . . . 163
11.E Spectral leakage . . . . . . . . . . . . . . . . . . . 168
11.F Aliasing . . . . . . . . . . . . . . . . . . . . . . . 171
11.G Links and further reading . . . . . . . . . . . . . 173

12 Discrete-time systems 174

12.1 Motivation: Line-follower robot . . . . . . . . . . 174
12.2 Difference equations . . . . . . . . . . . . . . . . 175
12.3 The z transform . . . . . . . . . . . . . . . . . . . 176
12.4 Mapping the s-plane to the z-plane . . . . . . . . 180
12.5 Converting from s to z . . . . . . . . . . . . . . . 184
12.6 Frequency warping . . . . . . . . . . . . . . . . . 185
12.A Deriving the bilinear transform . . . . . . . . . . 187
12.B Relationship between ωa and ωd . . . . . . . . . . 189

13 State-space and Laplace connections 191

13.1 Introduction . . . . . . . . . . . . . . . . . . . . . 191
13.2 Autonomous system ẋ = Ax . . . . . . . . . . . . 191
13.3 Non-autonomous system ẋ = Ax + Bu . . . . . . 193
13.A Calculating B −1 for a matrix Bn×n . . . . . . . . 195

14 MDoF vibration 198

14.1 The two-mass case . . . . . . . . . . . . . . . . . 198
14.2 The three-mass case . . . . . . . . . . . . . . . . 199
14.3 The n-mass case . . . . . . . . . . . . . . . . . . . 199
14.4 Natural frequencies and mode shapes . . . . . . . 200
1 Mechanical and electrical systems

In this first chapter we will model mechanical and electrical sys-

tems. We will start by modelling a mechanical mass-spring-damper
system and then see that each of its elements (mass, spring,
damper) have direct analogues in an electrical circuit (capacitor,
inductor, resistor). A key feature that the two systems share is
that both can oscillate. We will also see that the analogy between
mechanical and electrical systems is not quite complete—the me-
chanical mass doesn’t act quite like a capacitor but rather like
a grounded capacitor. This will lead us to discuss a relatively
new mechanical device that was first proposed only about twenty
years ago: the inerter. The inerter is now used in Formula 1
car suspensions where it more commonly goes by the name of a
J-damper.

Figure 1.1: Mass-spring-damper system.

5
1.1. The mass-spring damper system 6

1.1 The mass-spring damper system

A standard mass-spring-damper system is shown in figure 1.1.
There are three ingredients for a general mass-spring-damper sys-
tem. All three resist motion in some way. The first ingredi-
ent is mass: the mass resists acceleration because of its inertia.
The second ingredient is stiffness, which comes in the form of
a spring. The spring resists displacement because it doesn’t like
to be stretched or compressed away from its original length. The
third ingredient is damping: the damper resists velocity by pro-
viding a retarding force that is always in opposition to the mass’s
velocity. (Other models of damping exist, but this is by far the
most common, and the one we will use throughout.)
There is also a fourth ingredient. The system might be affected by
external forcing of some kind. By external we mean that the force
does its own thing rather than being produced by (for example) a
spring, which produces a force in response to what is happening
to it.

1.1.1 Without damping

Mass and stiffness are the bare minimum required for mechanical
vibration (i.e. harmonic motion). The governing equation in this
simplest case is
mq̈(t) + kq(t) = 0. (1.1)
You might remember that there is a natural frequency associated
with (1.1):
r
k
ωn = (1.2)
m
One simple way to find this is to substitute into (1.1) a trial
solution of q(t) = sin(ωt) and then solve for ω. The mass
p therefore
naturally oscillates with a natural frequency of ωn = k/m.
1.1. The mass-spring damper system 7

1.1.2 Why does it oscillate?

Why does (1.1) lead to oscillation? There are two ways to think
about this. The first is that the maths tells us so—the solution to
(1.1) is a sine wave, or a cosine wave, or any linear combination
of the two:
q(t) = A sin ωn t + B cos ωn t. (1.3)
The second way is to think in terms of energy. The system has
two elements: a mass and a spring. Both of them store energy.
The mass stores kinetic energy because of its velocity:
1
Emass (t) = mq̇ 2 (t).
2
The spring stores potential energy when it is stretched or com-
pressed:
1
Espring (t) = kq 2 (t).
2
The total energy of the system remains constant:
Etotal (t) = Emass (t) + Espring (t) = constant.
(To understand why, see § 1.A.) However, how much of that energy
is held by the mass and how much is held by the spring varies in
time. The mass and the spring constantly exchange energy with
each other. Thus we have two energy storage elements (a mass
and a spring) and a way for them to exchange energy (the mass
can do work on the spring; the spring can do work on the mass).
This constant exchange of energy between the two elements leads
to oscillation.

1.1.3 Including damping

Now let us also include the effect of damping. The governing
equation is then
mq̈(t) + cq̇(t) + kq(t) = 0, (1.4)
1.1. The mass-spring damper system 8

where c is the damping coefficient with units of (kg s−1 ).

This model says that the damper applies a force which opposes
velocity. (This is similar to the way the spring opposes displace-
ment.) This is easier to see if we rearrange (1.4) so that it looks
like F = ma:
mq̈(t) = −cq̇(t) − kq(t). (1.5)
The damper does not store any energy. In fact the only thing
the damper can do with energy is dissipate it and the rate at
which it does so is given by cq̇ 2 (to see where this comes from, see
§ 1.A). So now we have three elements, two of which (the mass
and spring) store energy and one of which (the damper) dissipates
energy.

1.1.4 Including forcing

The final ingredient is to allow for an external force to be applied
to the mass. This means an additional term on the right-hand
side of (1.5)
mq̈(t) = −cq̇(t) − kq(t) + f (t)
which we usually write as
mq̈(t) + cq̇(t) + kq(t) = f (t) (1.6)

Thinking in terms of energy again, the force f (t) acts as an ad-

ditional source of energy. Just think of pushing someone on a
swing. If you push them at the right moment then you will in-
crease their overall energy and cause them to swing higher. But if
you push them at the wrong time then, even though you are do-
ing work when you push, your pushing will actually reduce their
overall energy and cause them to swing lower. Thus the force f (t)
can serve both to increase and decrease the total energy of the
system (the work it does on the system can be both positive and
negative).
1.2. Rotating systems 9

1.2 Rotating systems

Turbines, gear trains, motors and pumps all involve the rotation
of mass. For these systems we can also define three elements
in analogy with the translational mass-spring-damper system: i )
rotational mass; ii ) rotational stiffness; and iii ) rotational damp-
ing.
The mass m for translating systems is replaced by the moment
of inertia J for rotating systems. For a point mass rotating at a
radius r from the axis of rotation this is

J = mr2

and has units of (kg m2 ). The moment of inertia of n point masses

rotating about an axis is
n
X
J= rk2 mk
k=1

and if the mass is modelled as being distributed continuously then

this sum becomes an integral over the body:
Z
J = r2 dm.

The force f for translating systems is replaced by torque T for

rotating systems. In analogy with the mass, the relationship be-
tween the torque T and the angular velocity Ω of the rotating
mass is
dΩ
T =J .
dt
(This is like F = m dv
dt
.) The relationship between torque and
angular velocity for the torsional spring and for the torsional
damper also follow the same patterns as the translational case.
1.3. Electrical circuits 10

The translational spring stiffness k is replaced by the torsional

spring stiffness K and
Z
T = Kθ = K Ω dt.
R
(This is like F = k vdt.) The translational damping c is replaced
by the torsional damping B:
T = BΩ.
(Which is like F = cv.) We therefore have the following cor-
respondences between translating and rotating mechanical sys-
tems:
velocity (v) ←→ angular velocity (Ω)
force (f ) ←→ torque (T )
mass (m) ←→ rotational mass (J)
damper (c) ←→ rotational damper (B)
spring (k) ←→ rotational spring (K)

1.3 Electrical circuits

Electrical circuits can also be modelled as being composed of three
ideal elements: capacitors, inductors and resistors. For each el-
ement we will now write the current i in terms of the voltage v;
but we could just as well choose to write the voltage v in terms
of the current i instead.

1.3.1 Capacitor
A capacitor stores energy in an electric field. The simplest imple-
mentation of a capacitor consists of two pieces of conducting ma-
terial separated by a dielectric material1 . The charge on an ideal
1
a material in which an electric field can be established without a signifi-
cant flow of charge through it
1.3. Electrical circuits 11

capacitor is proportional to the voltage difference across it:

q = Cv,

where C is the capacitance with units of farads (f = amp · sec/v).

Therefore the current i (the rate of change of the charge) is
dv
i=C .
dt
An ideal capacitor provides zero resistance to the flow of charge.

1.3.2 Inductor
An inductor stores energy in a magnetic field. The standard ex-
ample of an inductor is a helical coil of conducting wire. When
current flows through a conducting structure, a magnetic field is
established in the space or material around the structure. If the
current varies as a function of time then the magnetic field will
also vary. Lenz’s law then says that this changing magnetic field
will generate voltage differences in the conducting structure which
oppose the changing current. This property—whereby a voltage
difference is set up that opposes the rate of change of current
flow—is called inductance and the governing equation is
1 t
Z
i= v dt + i0 ,
L 0
where L is the inductance with units of henries (h = v · sec/amp).
An ideal inductor provides zero resistance to the flow of charge.

1.3.3 Resistor
A resistor dissipates energy. The simplest form of a resistor is
a thin metal wire. Unlike capacitors and inductors it doesn’t
store energy because it involves no electric-field or magnetic-field
effects.
1.3. Electrical circuits 12

All materials show some resistance to the flow of electric charge.

Materials with small resistance are called conductors and mate-
rials with large resistance are called insulators. An ideal resistor
satisfies
1
i = v.
R
A resistor does receive electrical power equal to
1 2
P = v .
R
It just doesn’t keep (i.e. store) any of it—it is all dissipated. This
means the energy is converted into other forms (heat) which can-
not easily be recovered.

1.3.4 RLC series circuit

Let’s now see what happens when we connect these three elements
together in series in a circuit, which is shown in figure 1.2. To
analyse the circuit we will need to use two laws:
1. Kirchoff ’s current law, which says that the sum of all cur-
rents entering into a circuit node is zero:
X
ik = 0. (1.7)
k

2. Kirchoff ’s voltage law, which says that the sum of all volt-
age drops around a circuit loop is zero:
X
vdrop = 0. (1.8)
elements

Applying these laws to the circuit in figure 1.2 we get

Z
di 1
L + Ri(t) + i(t) dt = vs (t). (1.9)
dt C
1.3. Electrical circuits 13

Figure 1.2: Series RLC circuit.

This is like the equation for the mechanical mass-spring damper
(1.4). In fact if we write the mechanical system (1.4) in terms of
velocity v(t) = q̇(t) instead of displacement then we get
Z
dv
m + cv(t) + k v(t) dt = f (t) (1.10)
dt
which is of exactly the same form as equation (1.9).

1.3.5 Mechanical-electrical analogy

By comparing equations (1.9) and (1.10) we can make the follow-
ing analogy between the mass-spring-damper (figure 1.1) and the
series RLC circuit (figure 1.2):
velocity (v) ←→ current (i)
force (f ) ←→ voltage (i)
mass (m) ←→ inductor (L)
damper (c) ←→ resistor (R)
spring (k) ←→ capacitor (1/C)
This is called the velocity-current analogy. This is the analogy
that is most intuitive and comes most naturally by comparing the
mass-spring-damper system of figure 1.1 to the electrical circuit of
figure 1.2. Notice that mass m is like inductance L; that damping
c is like resistance R; but that the spring stiffness k is like the
inverse of capacitance, 1/C.
1.4. A new device 14

1.4 A new device

I couldn’t resist telling you about the inerter. But to do so we
need to say a bit more about the analogy between mechanical and
electrical systems. It turns out that there is a second analogy in
which the correspondences are (velocity ↔ voltage) and (force ↔
current). This is called the velocity-voltage analogy and is the
most natural for introducing the inerter. To get ourselves there,
it is a good idea to first think carefully about the mass.

1.4.1 Something wrong with the mass

Notice something that is true of the mechanical and electrical
elements we have considered. For all of them except the mass,
what matters is the difference in something. Take for example
the spring: it has two ends (two ‘terminals’) and what matters
is the difference in displacement q between them. If I move both
ends of the spring by the same amount then no restoring force is
generated because all I did was move the spring as a rigid body
without compressing or stretching it. Therefore what R matters for
the spring is the difference in displacement, q = v dt, across it.
The same is true of a damper (difference in velocity v across it),
the capacitor (difference in dv/dt across it), the resistor
R (difference
in voltage v across it) and the inductor (difference in v dt across
it).
Notice that in every case there is what we might call an across
variable. We call it ‘across’ because it varies across the element.
For the mechanical system it is the velocity (or its derivative or
its integral). For the electrical system it is the voltage (or its
derivative or its integral).
The mass is different though. What matters for the mass is its
acceleration (dv/dt) rather than some difference in acceleration.
It seems to have only one terminal rather than two. The reason
1.4. A new device 15

for this is actually quite subtle: one of its terminals is the frame of
reference we have chosen (the ‘mechanical ground’). So the mass
behaves differently in an important way.
We can also define a ‘through’ variable’, so-called because it doesn’t
change across an element but instead is constant through it. For
the mechanical system the through variable is force. For the elec-
trical system the through variable is current.
To summarize, in the velocity-voltage analogy we have the follow-
ing correspondence:

velocity (v) ←→ voltage (v)

force (f ) ←→ current (i)
mass (m) ←→ capacitor (C)
damper (c) ←→ resistor (1/R)
spring (k) ←→ inductor (1/L)

But we now know that the mass does not behave exactly like a
capacitor. Instead it behaves like a grounded capacitor. This is
because one of the mass’s terminals is connected to the mechanical
ground, which is the frame of reference.
Is it possible to come up with a device that behaves like a mass but
has two terminals instead of one? In other words, is it possible to
come up with the true mechanical analogue of the capacitor?

1.4.2 The inerter

The answer, of course, is yes. Such a device was developed by
Professor Malcolm Smith at the University of Cambridge and goes
by the name of the inerter. The inerter is now used in Formula
1 car suspensions, where it more commonly goes by the name of
a J-damper. An inerter can be realized by a rack sliding in a
cylinder which drives a flywheel, as shown in figure 1.3.
 1.5. Summary of elements
The Inerter Concept and Its Application 16
Malcolm C. Smith

One method of realisation

rack pinions

terminal 2 gear flywheel terminal 1

r1 = 1.3:
Figure radius An
of rack pinion γ of=theradius
implementation of gyration
inerter. Takenoffrom
flywheel
here.
r2 = radius of gear wheel m = mass of the flywheel
There are two main advantages toα1the
r 3 = radius of flywheel pinion = γ/r 3 and α2 = r2 /r1
inerter:
1. A broader Equation rangeofofmotion:
mechanical F systems can
= (mα12 α 2 be realized when
2 ) (v̇2 − v̇1 )
compared to using a standard mass. (You can probably imag-
(Assumes mass of gears, housing etc is negligible.)
ine that the electrical circuits one could generate would be
limited
SICE Annual if any
Conference, Fukui,capacitors always had to be grounded, and this18
Japan, 4 August 2003

is the situation we have for mass.)

2. The effective mass of an inerter (its ‘inertance’) can be sig-
nificantly greater than the mass of the inerter itself. This is
because the effective mass (the inertance) is determined by the
moments of inertia of the gear and flywheel.

1.5 Summary of elements

A summary of the elements we have considered in this chapter
are shown in table 1.1. Note three important things:
• The rows have been arranged according to the velocity-voltage
analogy which, remember, is the most appropriate for talking
about the inerter. But we could have chosen the velocity-
current analogy instead.
• For the mass and the rotational mass I have indicated in the
1.6. Links and further reading 17

symbols that one of their two terminals is ‘grounded’. We would

not normally do this in, say, a free-body diagram, and I do not
expect you to do it. I have done it here to highlight that these
two elements differ from all other elements in having only one
‘free’ terminal.
• In each equation when I write velocity (or voltage) v, you should
read this as a difference in velocity (or voltage) across the ele-
ment. So for the spring, for example, we really have
Z
F = k (v2 − v1 ) dt

where v1 , v2 are the velocities at each of the two terminals.

1.6 Links and further reading

• Astrom & Murray § 2.1 Modelling concepts
• SDoF resonance vibration test
• Presentation on the inerter by Prof. Malcolm Smith
• Article on the inerter (J-damper)
• Keio University inerter
• Analogous electrical and mechanical systems
• Manual transmission, how it works?
1.6. Links and further reading 18

Name Symbol Equation Energy/power

Mass F = m dv
dt
E = 12 mv 2
11 2
R
Spring F = k v dt E= 2k
F

Damper F = cv P = cv 2

Rot. mass T = J dΩ
dt
E = 12 JΩ2

1 1 2
R
Rot. spring T =K Ω dt E= 2K
T

Rot. damper T = BΩ P = BΩ2

Capacitor i = C dv
dt
E = 12 Cv 2

i = L1 v dt E = 12 Li2
R
Inductor
1 1 2
Resistor i= R
v P = R
v

(Inerter) F = b dv
dt
E = 12 bv 2

Table 1.1: Summary of mechanical and electrical elements.

The electrical elements are arranged according to the velocity-
voltage analogy.
1.A. Energy of the mass-spring system 19

1.A Energy of the mass-spring system

1.A.1 Without damping
Substituting the general solution (1.3) into the expression for the
kinetic energy of the mass:
1
Emass = mq̇ 2 (t)
2
1
= mωn2 (A cos ωn t − B sin ωn t)2
2
1
= k(A2 cos2 (ωn t) + B 2 sin2 (ωn t) − 2AB sin(ωn t) cos(ωn t)).
2
Doing the same for the potential energy of the spring:
1
Espring = kq 2 (t)
2
1
= k(A sin ωn t + B cos ωn t)2
2
1
= k(A2 sin2 (ωn t) + B 2 cos2 (ωn t) + 2AB sin(ωn t) cos(ωn t)).
2
Adding them together to get the total energy and using the iden-
tity sin2 x + cos2 x = 1, we get
1 1
Etotal = k(A2 + B 2 ) = mωn2 (A2 + B 2 ), (1.11)
2 2
which is constant in time.
Another way to see the same thing, which also allows us to see
the effect of damping, is to start with the expression for the total
energy:
1 1
Etotal = mq̇ 2 (t) + kq 2 (t)
2 2
1.A. Energy of the mass-spring system 20

and then look at how it varies in time by differentiating it:

dEtotal 1 d 1 d
= m (v 2 ) + k (q 2 )
dt 2 dt 2 dt
dv dq
= mv + kq
dt dt
= mq̇ q̈ + kq q̇
= q̇(mq̈ + kq). (1.12)

The term in brackets in (1.12) is the governing equation of the

system, which we know from (1.1) is zero! Therefore
dEtotal
=0
dt
when there is no damping. This is the same as saying that Etotal
is constant in time and therefore amounts to the same thing as
(1.11).

1.A.2 Including damping

Now, what about when there is damping? In this case the gov-
erning equation becomes

mq̈(t) + cq̇(t) + kq(t) = 0

but it is still true that the total energy of the system is Etotal =
Emass + Espring (the damper stores no energy) and so we can still
write
dEtotal
= q̇(mq̈ + kq).
dt
Now to make use of our modified governing equation we can write
this as
dEtotal
= q̇(mq̈ + cq̇ + kq − cq̇)
dt
= −cq̇ 2 .
1.A. Energy of the mass-spring system 21

The damping constant satisfies c ≥ 0 and so this tells us that

dEtotal
≤0
dt
and so the system with damping loses energy.
2 Linear state-space models

How can we solve the types of linear dynamical systems that we

considered in chapter 1? There are two quite different ways to do
so and it is important that we look at both. The first option is to
use Laplace transform techniques. This will come later, starting
in chapter 4. The second option is to form a linear state-space
model of the system.1 That is the subject of this chapter.
We will now build up our understanding of linear state-space mod-
els piece-by-piece by first looking at the simplest case of the evo-
lution equation ẋ(t) = Ax(t); then at how to include the effect
of inputs; and finally at how to single out particular outputs of
interest.

2.1 The state x

A linear state-space model is of the form
ẋ(t) = Ax(t), (2.1)
where x ∈ Rn is called the system state and A is a matrix of
dimension n × n. But what is the state, x, and what does it
represent? The best way to see this is to look directly at an
example.
1
Notice that we specified that the state-space model is linear. That is
because nonlinear state-space models are also possible—we will look at them
a little bit in chapter § 3.

22
2.1. The state x 23

2.1.1 Example: second-order system

The damped mass-spring-damper system (1.4) that we saw in
chapter 1 is a second-order ODE:

mq̈(t) + cq̇(t) + kq(t) = 0. (2.2)

(We have set the forcing to be zero for now.) How would we
represent this system in the state-space form (2.1)? The first
thing to notice is that equation (2.2) is second-order (because the
highest time derivative is two) while the state-space form (2.1)
involves n equations, each of which is first-order (since the highest
derivative in it is one). So we somehow need to convert a second-
order equation into two first-order equations. We can do this by
doubling the number of equations. If we set the state to be
   
x1 (t) q(t)
x(t) = =
x2 (t) q̇(t)

then we can write down two equations

d
q(t) = q̇(t)
dt
d k c
q̇(t) = − q(t) − q̇(t)
dt m m
and the linear state-space model is
    
d q(t) 0 1 q(t)
= k . (2.3)
dt q̇(t) −m − mc q̇(t)

Notice that we used up one equation to state the obvious: the first
row of (2.3) just says that q̇(t) = q̇(t). This is the price we pay for
writing one second-order equation as two first-order equations. In
other words, one equation of second order became two equations,
each of first order. All of the dynamics of the mass-spring-damper
is contained in the second row of (2.3) (which is just the original
2.2. The matrix exponential 24

equation (2.2) divided through by the mass m). Writing this in

terms of x1 (t) and x2 (t) we get
    
d x1 (t) 0 1 x1 (t)
= k , (2.4)
dt x2 (t) −m − mc x2 (t)
which is our state-space model.

2.2 The matrix exponential

To solve equations like (2.1) and (2.4) we will need the matrix
exponential, which we will define shortly. Before that, let’s first
review the scalar exponential, eat .

2.2.1 The scalar exponential

You should already know that the equation
ẋ(t) = ax(t) (2.5)
with initial condition x(0) = x0 has solution
x(t) = eat x0 . (2.6)
We can check that this works by substituting the solution (2.6)
back into the governing equation (2.5):
dx d
= (eat x0 ) = aeat x0 = ax(t). (2.7)
dt dt
As a reminder, two important properties of the scalar exponential
are that
1. Its power-series representation is
1 1
eat = 1 + at + (at)2 + (at)3 + . . .
2! 3!
2. Its derivative is:
d at
e = aeat .
dt
2.2. The matrix exponential 25

2.2.2 The matrix exponential

What if instead of a scalar a we have a square matrix A like in
(2.1)? It turns out that the scalar exponential generalizes so that
we can define something similar for matrix equations. We are now
interested in solving the matrix equation:

ẋ(t) = A x(t) (2.8)

with initial condition x(0) = x0 , where x0 is now a vector. The

solution is
x(t) = eAt x0 (2.9)
where eA is called the matrix exponential.
Here are three important properties of the matrix exponential.
The first two are generalizations of properties of the scalar expo-
nential while the third is new.
1. Its power-series representation is
1 1
eAt = I + At + (At)2 + (At)3 + . . . (2.10)
2! 3!

2. Its derivative is
d At
e = AeAt . (2.11)
dt
3. It commutes with A (i.e. the order of their multiplication does
not matter):
AeAt = eAt A.

You can check the solution (2.9) by substituting it into the two
sides of (2.8) and using (2.11) for the derivative of eAt .
2.3. Eigendecomposition of eAt 26

2.3 Eigendecomposition of eAt

An eigenvalue (λi ) and eigenvector (vi ) pair of a square matrix A
satisfy
Avi = λi vi . (2.12)
If A is of dimension n × n then there will be n such pairs and
equation (2.12) will be true n times. These n equations can then
be written compactly as one big matrix equation:

AV = V Λ (2.13)

where  
λ1 0
 
| |
Λ=
 .. ,

V = v1 · · · vn  .
.
0 λn | |
That is, the matrix Λ contains on its diagonal the n eigenvalues
λi ; and the matrix V contains in its columns the n eigenvectors vi .
(If you’re not sure about the order of terms in (2.13), try writing
it out for a 2 × 2 matrix.)
Equation (2.13) can be rearranged for A:

A = V ΛV −1 (2.14)

Substituting (2.14) into the matrix exponential we get

−1 t
eAt = eV ΛV . (2.15)
−1
Can we write eV ΛV t in any other, simpler way? Let’s see what
−1
happens if we substitute eV ΛV t into the power series represen-
2.3. Eigendecomposition of eAt 27

tation of the matrix exponential (2.10):

−1 t t2
2
eV ΛV = I + tV ΛV −1 + V ΛV −1 + . . .
2!
t2
= I + tV ΛV −1 + V ΛV −1 V ΛV −1 + . . .


2!
t2
= V V −1 + tV ΛV −1 + V ΛΛV −1 + . . .


 2! 
1
= V I + Λt + (Λt) + . . . V −1
2
2!
Λt −1
=Ve V ,
where eΛt is a diagonal matrix containing scalar exponential terms
on its diagonal:  
eλ1 t 0
eΛt = 
 ... .

λn t
0 e
So the solution x(t) = eAt x0 can be written as
x(t) = eAt x0
= V eΛt V −1 x0
  λ1 t  
e 0 − uT1 −

| |
= v1 · · · vn   .. ..
 x0
  
.  .
| | 0 e λn t − uTn −
= (v1 uT1 x0 )eλ1 t + · · · + (vn uTn x0 )eλn t
| {z } | {z }
c1 c
  n 
| |
λ1 t   λn t  
= c1 e v1 + · · · + cn e vn . (2.16)
| |
In words (2.16) says that the solution in time is given by a linear
combination of the eigenvectors vi . The amount of each eigenvec-
tor that appears in the response is given by the constants ci (which
2.4. Stability and frequency of oscillation 28

themselves depend on the initial condition x0 ); and the way in

which each eigenvector vi evolves in time is determined by its
corresponding eigenvalue λi , giving terms of the form eλi t .

2.4 Stability and frequency of oscillation

In general a particular eigenvalue λi of A will be complex:

λi = σi ± jωi .

Each eigenvalue therefore gives rise to terms of the form

ci eλi t = ci e(σi +jωi )t = ci eσi t (cos ωi t + j sin ωi t)

and so we see that the real part, σi , gives information on the

exponential growth or decay of a given eigenvalue; and the imagi-
nary part, ωi , gives information on the frequency of oscillation of a
given eigenvalue. From this we can say the following things:
• a linear system is stable if all eigenvalues of A have a negative
real part
• a linear system is unstable if any of the eigenvalues of A has
a positive real part
• a linear system is marginally stable if any of the eigenvalues
of A has a real part of zero
• an eigenvalue λ = σ + jω with imaginary part ω will give rise
to oscillations at frequency ω
• if λ = σ + jω (with ω =
6 0) is an eigenvalue then its complex
conjugate λ∗ = σ − jω is also an eigenvalue
2.5. Complex eigenvalues come in pairs 29

2.5 Complex eigenvalues come in pairs

We have seen that the solution in terms of eigenvalues and eigen-
vectors is
   
| |
λ1 t   λn t  
x(t) = c1 e v1 + · · · + cn e vn . (2.17)
| |

We also know that the eigenvalues are in general complex. (The

coefficients ci and the eigenvectors vi will also be complex in gen-
eral.) But we also know that the solution in time must be real.
This appears to be a contradiction. How do we resolve it? The
answer is that any complex eigenvalues come in pairs. Therefore
individual terms in (2.17) can be complex, but they always com-
bine to give an overall solution x(t) that is real. We will say a bit
more on this in chapter 4.

2.6 Including inputs

Most of this chapter has been concerned with the equation ẋ(t) =
Ax(t). But what if there are also inputs? How do we include
them in the state-space model? By “inputs” I mean anything
that enters the system from outside and influences it in some
way. We can include the influence of inputs by adding an extra
term to equation (2.1):

ẋ(t) = Ax(t) + Bu(t), (2.18)

where u ∈ Rp is a vector of inputs and B is called the forcing

matrix and has dimensions of n × p.
2.6. Including inputs 30

2.6.1 Example: mass-spring-damper with forcing

Looking again at the example of the mass-spring-damper from
§ 2.1.1, this time including the forcing, we have

mq̈(t) + cq̇(t) + kq(t) = f (t) (2.19)

and the force f (t) is an input to the system. We can repeat

what we did in § 2.1.1: write (2.19) as two first-order equations
by setting the two states to be x1 = q and x2 = q̇. Since the
force f (t) is acting as the input to the system, we should also set
u(t) = f (t) so that we can use (2.18). Doing so we get the matrix
equation
      
d q(t) 0 1 q(t) 0
= k c + 1 u(t) (2.20)
dt q̇(t) − m
− m
q̇(t) m

1 T
and so the B matrix for this system is B = [0 m
] .

2.6.2 General solution with inputs

How do we solve (2.18) or (2.20)? In other words, how do we find
the solution in time when there are inputs? The solution is
Z t
At
x(t) = e x(0) + eA(t−τ ) Bu(τ ) dτ. (2.21)
0

(We will not derive this.) The solution (2.21) is not terribly easy
to use because the integral term can be quite complicated to eval-
uate. For that reason we will not use (2.21) in this course and
will instead use Laplace transform techniques when there are in-
puts. But (2.21) is included here for completeness. (It also helps
to motivate the Laplace transform techniques of chapter 4, since
(2.21) shows us what we would have to do without them!)
2.7. Including outputs 31

2.7 Including outputs

We are often more interested in some states than others. For ex-
ample, the state x for the mass-spring-damper in § 2.1.1 contains
both the displacement (x1 = q) and the velocity (x2 = q̇). But
suppose we only care about the displacement, q? For example,
we might be designing a feedback controller, and the controller
might make its control decisions based on the displacement q but
not the velocity q̇. How can we ‘pick out’ the displacement q from
the state x?
More generally, how could we pick out certain states that are of
particular interest? For example, our system might have eight
states, of which two are particularly important to us. The way
we can do this is by defining an output, y(t), such that
y(t) = Cx(t) + Du(t). (2.22)
Here y ∈ Rq is a vector of outputs and C and D are matrices.
In words, equation (2.22) says that the output y(t) is some linear
combination of the states x and some linear combination of the
inputs u. To see how this works, let’s look at a specific exam-
ple.

2.7.1 Example: mass-spring-damper with y = q

Looking again at the mass-spring-damper from § 2.6.1, the state-
space equation ẋ = Ax(t) + Bu(t) is
      
d x1 0 1 x1 0
= k c + 1 u. (2.23)
dt x2 − m − m x2 m

We need to choose the C and D matrices in (2.22) such that we

‘pick out’ the displacement q. This is achieved by
   
  q   x1
y= 1 0 = 1 0 , (2.24)
q̇ x2
2.8. Links and further reading 32

so that C = [1 0] and D = 0. Note that the D matrix is zero.

It turns out that this is often the case. So when would we need
it? The simplest possible example is when the output y is simply
some constant multiplied by the input u:

y(t) = Ku(t). (2.25)

Then D = K and the A, B and C matrices are not needed!2

Another example for which the D matrix is non-zero is the mass-
spring-damper system

mq̈(t) + cq̇(t) + kq(t) = f (t)

when we choose the output to be the acceleration, y(t) = q̈(t).

(We might guess that we need to do something a little different
here, since the acceleration is not found in either of the states.)
In this case the equation for the output is
 
 x1
+ m1 u
 k c
 
y = −m −m
x2
k c
so that C = [− m − m
] and D = [ m1 ].

2.8 Links and further reading

• Astrom & Murray § 2.2 State-space models
• DC motor

2
One way to think about this is that the order of equation (2.25) is zero
because the highest time derivative is zero. Then the A matrix is of dimension
0 × 0, which is another way of saying we don’t need it.
3 Linearization

Few physical elements display truly linear character-

istics. For example the relation between force on a
spring and displacement of the spring is always non-
linear to some degree. The relation between current
through a resistor and voltage drop across it also devi-
ates from a straight-line relation. However, if in each
case the relation is reasonably linear, then it will be
found that the system behaviour will be very close to
that obtained by assuming an ideal, linear physical
element, and the analytical simplification is so enor-
mous that we make linear assumptions wherever we
can possibly do so in good conscience.
Robert Cannon, Dynamics of Physical Systems, 1967
The first two chapters have dealt with linear systems and this
whole subject is concerned with linear systems: it could just as
well be called ‘Linear systems modelling & analysis’. But what
if we have a dynamical system that is nonlinear? That is the
subject of this chapter. We consider nonlinear systems in this
chapter only so that we can learn how to linearize them, so even
this chapter is concerned with linear systems.
Imagine a pendulum hanging freely from a pivot. This is a nonlin-
ear system because there is a (sin θ) term in the governing equa-
tion. The pendulum has two equilibria. The first equilibrium is

33
3.1. Finding equilibria 34

when the pendulum hangs vertically downwards. This equilib-

rium is stable because if we disturb the pendulum away from the
vertical, it will go back there after oscillating for a while. The
second equilibrium is when the pendulum balances perfectly up-
right. This is an equilibrium because all forces are in balance.
However, it is unstable because the slightest perturbation will
push it further away from vertical. In this chapter we look at two
things:
1. How to find equilibria of nonlinear systems; and
2. How to linearize a nonlinear system about a given equilibrium.
Together these will allow us to look at important characteristics
of linearized systems (including their stability and frequency of
oscillation) about an equilibrium.

3.1 Finding equilibria

Consider a general nonlinear state-space system:

ẋ = f (x). (3.1)

(There is no input for the moment.) The function f is in general

nonlinear. The state x is a vector so that writing out (3.1) in full
we have    
ẋ1 f1 (x1 , x2 , . . . , xn )
 ẋ2   f2 (x1 , x2 , . . . , xn ) 
 ..  =  . (3.2)
   
..
.  . 
ẋn fn (x1 , x2 , . . . , xn )
In words this says that the rate of change of the first state x1 is a
nonlinear function of all of the states (including itself); that the
rate of change of the second state x2 is a nonlinear function of all
of the states (including itself); and so on.
3.2. Linearization 35

The vector  
x1
x =  ...  (3.3)
 
xn
is an equilibrium if  
0
 .. 
f (x) =  .  . (3.4)
0
In words: the vector x is an equilibrium if, when we ‘f ’ it, we get
back a vector of zeros. This makes sense. The governing equation
(3.1) tells us that the time derivative of x is given by f (x). So
if f (x) gives back zero, then the time derivative of x is zero, and
the state does not change in time. In other words, we are at an
equilibrium.
Note that an equilibrium x is not necessarily unique. A dynami-
cal system of the form (3.1) can have zero, one, or many equilib-
ria.

3.2 Linearization
Once we have an equilibrium x, the next thing we are interested
in is linearizing the system about this equilibrium. To keep things
simple, let’s consider a dynamical system with a single-state to
start with:
ẋ = f (x).
An example function f (x) is plotted in figure 3.1.
Now consider a small perturbation δ away from the equilibrium
x so that
x=x+δ
3.2. Linearization 36

Figure 3.1: Function f (x) and its slope at equilibrium x.

then we can say that

δ =x−x
δ̇ = ẋ.

A Taylor series expansion about x then gives

δ̇ = ẋ = f (x)
∂f 1 ∂ 2f
≈ f (x) + (x − x) + (x − x)2 + . . .
∂x x 2! ∂x2 x
or
∂f 1 ∂ 2f
δ̇ ≈ δ+ δ2 + . . . (3.5)
∂x x 2 ∂x2 x
Here comes the key point: when we linearize, we retain only
the leading-order (i.e. linear) term so that we are left with sim-
ply:
∂f
δ̇ = δ (3.6)
∂x x
or
δ̇ = aδ (3.7)
with the constant a given by
∂f
a= .
∂x x
3.3. From one state to many states 37

Graphically, a represents the slope of f (x) evaluated at x (see

figure 3.1).

3.3 From one state to many states

How could we generalize what we did for a single state to many
states? Well, suppose we have a system with two states:

ẋ1 (t) = f1 (x1 (t), x2 (t)) (3.8)

ẋ2 (t) = f2 (x1 (t), x2 (t)). (3.9)

Pay careful attention to what this says: it says that ẋ1 (t) depends
on both x1 (t) and x2 (t); and that ẋ2 (t) depends on both x1 (t) and
x2 (t). How would we do the same Taylor series expansion for this
system? This time we will have two equations; and we need to
account for derivatives not only with respect to x1 but also with
respect to x2 . The first equation is

∂f1 1 ∂ 2 f1
δ̇1 ≈ f1 (x1 , x2 ) + δ1 + δ12 + . . .
∂x1 x 2 ∂x21 x
∂f1 1 ∂ 2 f1
+ δ2 + δ22 + . . .
∂x2 x 2 ∂x22 x

Retaining only linear terms and using f1 (x1 , x2 ) = 0:

∂f1 ∂f1
δ̇1 ≈ δ1 + δ2 .
∂x1 x ∂x2 x

Applying the same logic to the second equation we get

∂f2 ∂f2
δ̇2 ≈ δ1 + δ2 .
∂x1 x ∂x2 x
3.4. Linearization with inputs 38

These two equations can be put into a single matrix equation:

∂f1 ∂f1
 

δ̇ =  ∂x
 1 ∂x2 
∂f2 ∂f2  δ
∂x1 ∂x2 x
= A δ.
So that was the two-dimensional case. You might be able to
imagine how this generalizes to n dimensions. The A matrix be-
comes
∂f1 ∂f1
 
 ∂x1 · · · ∂xn 
 . .. .. 
A=  .
. .  .
.  (3.10)
 ∂fn ∂fn 
···
∂x1 ∂xn x
So A is the matrix of derivatives of f , all evaluated at the equilib-
rium x. The dimension of A is n×n because there are n functions
f1 . . . fn , and each function needs to be differentiated with respect
to each of the n states x1 . . . xn . This matrix is important enough
that it has a name: we call it the Jacobian of f . It represents the
linearization of the nonlinear function f evaluated at the equilib-
rium x.

3.4 Linearization with inputs

So far we have considered nonlinear systems without inputs:
ẋ(t) = f (x).
What if there are also inputs u? The way we can include inputs
is by writing
ẋ(t) = f (x, u) (3.11)
so that ẋ(t) now depends not only nonlinearly on x, but also
nonlinearly on u. Remember that x in (3.11) is a vector and this
3.5. Cart-pendulum example 39

is also true of the input u—there is more than one of them in

general. Linearizing (3.11) we get
ẋ(t) = Ax(t) + Bu(t), (3.12)
where A is given by
∂f1 ∂f1
 
 ∂x1 · · · ∂xn 
 . ... .. 
A=
 .. . 
 (3.13)
 ∂fn ∂fn 
···
∂x1 ∂xn x,u
and B is given by
∂f1 ∂f1
 
 ∂u1 · · ·
 . ∂um 
B =  .. . .. .. 
 .
.  (3.14)

 ∂fn ∂fn 
···
∂u1 ∂um x,u
Notice that the input u has shown up in two ways. First, we now
have to evaluate the derivatives of f with respect to u to get the
B matrix. Second, the value of u at the equilibrium, u, is required
in order to evaluate both A and B.

3.5 Cart-pendulum example

The cart-pendulum system is a nice example to look at a few
different parts of what we have learnt so far. In particular we can
look at i) forming a nonlinear state-space model; ii) finding its
equilibria; and iii) linearizing about its ‘up’ or ‘down’ equilibria.
In its down position we have something like a crane; in its up
position we have something like a balancing device such as a one-
wheeled electric skateboard, a hoverboard, or a SpaceX Falcon
Heavy booster in its landing phase.
3.5. Cart-pendulum example 40

se

M
i f

fee
i d

m
Figure 3.2: Cart-pendulum system.
The equations of motion for the cart-pendulum system are derived
in § 3.A. They are

m(ẍ cos θ + lθ̈ + g sin θ) = d (3.15)

(M + m)ẍ − ml(θ̇2 sin θ + θ̈ cos θ) = f + d cos θ. (3.16)

The first equation represents the equation of motion for the pen-
dulum and the second equation represents the equation of motion
for the cart. Notice that they are nonlinear. We want to put these
two equations in the form

ẋ = f (x, u).

We can do this by first grouping together the highest-order deriva-

tives on one side of the equations of motion:
      
m cos θ ml ẍ mg sin θ d
+ =
M + m ml cos θ θ̈ −mlθ̇2 sin θ f + d cos θ
or
   −1  
ẍ m cos θ ml d − mg sin θ
= .
θ̈ M + m ml cos θ f + d cos θ + mlθ̇2 sin θ
3.5. Cart-pendulum example 41

Expanding and simplifying we get1

−mlθ̇2 sin θ − mg cos θ sin θ − f

ẍ =
m cos2 θ − (M + m)
(m + M )g sin θ + mlθ̇2 sin θ cos θ + f cos θ − ( M
m
+ sin2 θ)d
θ̈ =
l(m cos2 θ − (M + m))
Now define the state x to be

x1 = x
x2 = θ
x3 = ẋ
x4 = θ̇,

and the nonlinear state-space model is

ẋ1 = x3
ẋ2 = x4
−mlx24 sin x2 − mg cos x2 sin x2 − f
ẋ3 =
m cos2 x2 − (M + m)
(m + M )g sin x2 + mlx24 sin x2 cos x2 + f cos x2 − ( M
m
+ sin2 x2 )d
ẋ4 =
l[m cos2 x2 − (M + m)]
This is complicated, and you can probably imagine that it’s not
very nice to linearize! But it’s not as bad as it looks for two
reasons:
• The first two equations are easy to linearize because they are
already linear and each only depends on one other state.
1
A reminder that the inverse of a 2 × 2 matrix is
 −1    
−1 a b 1 d −b 1 d −b
A = = = .
c d det(A) −c a ad − bc −c a
3.5. Cart-pendulum example 42

• The third and fourth equations are harder, but the only tricky
case is differentiating them with respect to x2 .
So out of the sixteen entries in the A matrix, many of them are
zero, and only two of them are difficult to evaluate. These are
∂f3 /∂x2 and ∂f4 /∂x2 .

3.5.1 Linearization about pendulum down

For the pendulum in the down position, θ = 0. The linearized
equations are

mẍ + mlθ̈ + mgθ = d (3.17)

(M + m)ẍ + mlθ̈ = f + d. (3.18)

The state-space model corresponding to these equations is

        
1 0 0 0 x 0 0 1 0 x 0 0  
0 1 0 0  d  θ  0
    0 0 1  θ  0
    0 f
 = +
0 0 m ml dt ẋ 0 −mg 0 0 ẋ 0 1 d
0 0 M + m ml θ̇ 0 0 0 0 θ̇ 1 1

3.5.2 Linearization about pendulum up

For the pendulum in the up position, θ = π = 180◦ .

mẍ − mlθ̈ + mgθ = −d (3.19)

(M + m)ẍ − mlθ̈ = f − d. (3.20)

The state-space model corresponding to these equations is

        
1 0 0 0 x 0 0 1 0 x 0 0  
 d  θ  = 0   θ +0 0  f
0 1 0 0     0 0 1    

0 0 m −ml dt ẋ 0 −mg 0 0 ẋ 0 −1 d
0 0 M + m −ml θ̇ 0 0 0 0 θ̇ 1 −1
3.6. Links and further reading 43

3.6 Links and further reading

• Astrom & Murray § 5.4
• Control of double inverted pendulum

3.A Cart-pendulum equations of motion

Here we derive the equations of motion for the cart-pendulum
system of § 3.5. Perhaps the key step in deriving these equations
is to realize that any motion of the pendulum affects the cart,
and any motion of the cart affects the pendulum. In other words,
the cart and the pendulum are coupled together. How do we
reflect this in the free-body diagrams? We can do so by including
internal forces Rx and Ry . They represent the forces applied to
the pendulum by the cart, and the forces applied to the cart by
the pendulum. Newton tells us that these have to be equal and
opposite: this means that, however we choose to set up the force
Rx for the pendulum, we had better make sure that its sign is
reversed when we draw it for the cart.

Figure 3.3: Free-body diagrams for the pendulum and the cart.
3.A. Cart-pendulum equations of motion 44

The free-body diagrams for the cart and for the pendulum are
shown in figure 3.3. Note that the coordinates of the cart are
(x, y) and the coordinates of the pendulum are (xp , yp ) with

xp = x + l sin θ
yp = −l cos θ.

3.A.1 Pendulum
First let’s consider the pendulum. Summing forces in the x-
direction:
d 2 xp
Rx + d cos θ = m
dt2
d2
= m 2 (x + l sin θ)
dt
d
= m (ẋ + lθ̇ cos θ)
dt
or
Rx = m(ẍ − lθ̇2 sin θ + lθ̈ cos θ) − d cos θ (3.21)
Summing forces in the y-direction:
d2 yp
Ry − mg + d sin θ = m
dt2
d2
= m 2 (−l cos θ)
dt
d
= m (lθ̇ sin θ)
dt
or
Ry = m(lθ̈ sin θ + lθ̇2 cos θ) + mg − d sin θ (3.22)
Taking moments about the pendulum’s point-mass, noting that
the moment of inertia of a point-mass is zero (I = 0):

Rx l cos θ + Ry l sin θ = I θ̈ = 0. (3.23)

3.A. Cart-pendulum equations of motion 45

Substituting (3.21) and (3.22) into (3.23) and after some rear-
ranging we get

m(ẍ cos θ + lθ̈ + g sin θ) = d (3.24)

3.A.2 Cart
Now for the cart. We only need consider the force balance in x
because the cart is supported in y and therefore cannot accelerate
in that direction. The force balance in the x-direction is

M ẍ = f − Rx .

We can use (3.21) again to eliminate the internal force Rx , giv-

ing

(M + m)ẍ − mlθ̇2 sin θ + mlθ̈ cos θ = f + d cos θ (3.25)

4 The Laplace transform

The subject of chapters 5–7 will be transfer functions of lin-

ear systems. To prepare ourselves, in this chapter we consider
Laplace transforms, which are essential for forming transfer func-
tions.

4.1 Life without the Laplace transform

By the end of an engineering degree you could be forgiven for not
wanting to see another Laplace transform table for a while. But
why do they come up so often? One way to understand that is
to ask what we would do without them. Suppose I give you the
equation of a second-order system such as:

mq̈(t) + cq̇(t) + kq(t) = f (t) (4.1)

and I ask you to find the solution in time, q(t), when the initial
condition is zero and the forcing is

f (t) = sin ω1 t + sin ω2 t. (4.2)

How would you do it? It turns out that this is a difficult problem
to solve in the time domain. One way to do it would be to convert
the governing equation (4.1) to a linear state-space model and
then solve for q(t) using (2.21). Converting to a state-space model
is straightforward enough but solving the integral in (2.21) will

46
4.2. Definition of the Laplace transform 47

be quite tricky in general. Another problem with this approach is

that it is not so clear what the effect of changing the input (4.2)
might be on the solution. We would just have to recalculate it
and see what happens.
Laplace transform techniques give us a way of analysing systems
with inputs like (4.1) in a general way, and it also provides, for
free, useful insights into what is going on with the system to aid
our understanding.

4.2 Definition of the Laplace transform

The Laplace transform X(s) of a function x(t) is defined as
Z ∞
X(s) = L{x(t)} = e−st x(t) dt. (4.3)
0

Here s is the Laplace variable which is in general complex:

s = σ + jω.

Notice that the integral goes from t = 0 to t = ∞. We therefore

implicitly assume that the function satisfies x(t) = 0 for t < 0
(If it is non-zero for t < 0 then the integral will ‘miss it’.) We
also have to assume that limt→∞ e−st x(t) = 0 so that the integral
converges.

4.3 Laplace transform of common functions

Here we will derive the Laplace transform of a select few functions
before compiling a list of all the most important Laplace transform
pairs in table 4.1.
4.3. Laplace transform of common functions 48

4.3.1 The exponential function e−at

The Laplace transform of e−at is:
Z ∞
−at
L(e ) = e−st e−at dt
Z0 ∞
= e−st−at dt
Z0 ∞
= e−(s+a)t dt
0
1 ∞
=− e−(s+a)t 0
s+a
1
=− (e−∞ − e0 )
s+a
so
1
L{e−at } = (4.4)
s+a
Notice that the integral above only converges if Re(s) > −a.
(Otherwise e(−s+a)t will blow up for t → ∞ rather than going to
zero.) The therefore say that the region of convergence (ROC) of
e−at is Re(s) > −a. We will not dwell on this any further (that
is for mathematicians to worry about) but I wanted to at least
mention it. When we learn a new tool, we should also know when
that tool can and cannot be used. The region of convergence
(ROC) for this and all other functions is given in table 4.1.

4.3.2 Unit step function u(t)

The unit step is sometimes denoted by u(t) and sometimes by
1(t). It is defined as
(
0 t<0
u(t) = 1(t) =
1 t > 0.
4.3. Laplace transform of common functions 49

Its Laplace transform is

Z ∞ ∞
1
L{u(t)} = e−st dt = − e−st
0 s 0
1
= − (e−∞ − e0 )
s
so
1
L{u(t)} = L{1(t)} = (4.5)
s

4.3.3 Impulse function δ(t)

The impulse function δ(t) is defined as
δ(t) = lim g(t)
→0

with (
1/ 0 ≤ t ≤ 
g(t) =
0 0.
An important property of the impulse function is that, for any
function f (t), Z ∞
δ(t − τ )f (t)dt = f (τ ),
0
so the integral ‘picks out’ the value of f (t) at t = τ . The Laplace
transform of the impulse function is
Z ∞
L{δ(t)} = e−st δ(t)dt
0
= e−st t=0
= e0
so
L{δ(t)} = 1 (4.6)
The impulse function δ(t) also goes by the name of the Dirac delta
function after Paul Dirac.
4.4. Properties of the Laplace transform 50

4.3.4 Common Laplace transform pairs

Table 4.1 contains Laplace transforms of some common functions.
Some points to note:
• RoC stands for Region of Convergence (see § 4.3.1).
• u(t) is the unit step function (see § 4.3.3).

4.4 Properties of the Laplace transform

Having looked at some common Laplace transform pairs, we now
look at some important properties of the Laplace transform.

4.4.1 Linearity/superposition
For two functions f1 (t) and f2 (t),
Z ∞
L{αf1 (t) + βf2 (t)} = e−st (αf1 (t) + βf2 (t))dt
0
Z ∞ Z ∞
−st
=α e f1 (t)dt + β e−st f2 (t)dt
0 0
= αL{f1 (t)} + βL{f2 (t)}
= αF1 (s) + βF2 (s)

so
L{αf1 (t) + βf2 (t)} = αF1 (s) + βF2 (s) (4.7)
Note that this does not apply to the product of functions, i.e.
L{f1 (t)f2 (t)} =
6 F1 (s)F2 (s)!
4.4. Properties of the Laplace transform 51

x(t) X(s) RoC

x(t)e−st dt
R
x(t) (definition)

δ(t) 1 all s
1
u(t) Re(s) > 0
s
1
e−at u(t) Re(s) > −a
s+a
m!
tm e−at u(t) Re(s) > −a
(s + a)m+1
s
cos(ω0 t) u(t) Re(s) > 0
s + ω02
2

ω0
sin(ω0 t) u(t) Re(s) > 0
s + ω02
2

s+a
e−at cos(ω0 t) u(t) Re(s) > −a
(s + a)2 + ω02
ω0
e−at sin(ω0 t) u(t) Re(s) > −a
(s + a)2 + ω02
Table 4.1: Common Laplace transform pairs.
4.4. Properties of the Laplace transform 52

4.4.2 Time-scale
What if we stretch time by a factor of a? The Laplace transform
of f (at) is
Z ∞
L{f (at)} = e−st f (at) dt
0

Now define a new time, τ = at. Then dτ = a dt and

Z ∞
1 −s(τ /a)
L{f (at)} = e f (τ ) dτ
0 a
1 ∞ −s(τ /a)
Z
= e f (τ ) dτ
a 0
so
1 s
L{f (at)} = F (4.8)
a a

4.4.3 Derivative
The Laplace transform of the derivative of x(t) is
Z ∞  
−st d
L{x(t)} = e x(t) dt.
0 dt
Rb Rb
Now we use integration by parts, a uv 0 dt = uv|ba − a u0 vdt, with
u = e−st and v 0 = dtd
x(t) to give
Z ∞   Z ∞
−st d −st ∞
e x(t) dt = e x(t) 0 − −se−st x(t) dt
0 dt
Z ∞0
= −x(0) + s e−st x(t) dt
0

so
L{ẋ(t)} = sX(s) − x(0) (4.9)
4.4. Properties of the Laplace transform 53

4.4.4 Second derivative

For the second derivative we can use the result for the first deriva-
tive twice:
 
d
L{ẍ(t)} = L ẋ(t)
dt
= sL{ẋ} − ẋ(0)
= s(sX(s) − x(0)) − ẋ(0)
= s2 X(s) − sx(0) − ẋ(0)
so
d2
 
L x(t) = s2 X(s) − sx(0) − ẋ(0) (4.10)
dt2

4.4.5 Higher-order derivatives

For the nth derivative of x(t) we can apply the rule for the first
derivative n times. This gives
 n 
d
L n
x(t) = L{x(n) (t)}
dt
= sn X(s) − sn−1 x(0) − sn−2 ẋ(0) − . . . − x(n−1) (0)

4.4.6 Integration
If we integrate x(t) and then take the Laplace transform we get
Z t 
1
L x(t) dt = X(s) (4.11)
0 s

4.4.7 Time-shift
If we time-shift x(t) by an amount of time τ then we get
L{x(t − τ )u(t − τ )} = e−τ s X(s) (4.12)
where u(t) is the step function.
4.5. The inverse Laplace transform 54

4.4.8 Table of properties

A summary of the important properties of the Laplace transform
are summarized in table 4.2. Note that:
• x(n) denotes the nth time derivative of x.
• u(t) is the unit step function.

4.5 The inverse Laplace transform

The inverse Laplace transform is defined as
Z σ+j∞
1
x(t) = X(s) est ds, (4.13)
2πj σ−j∞
where σ is a real number such that the contour path of integration
is in the region of convergence of X(s). This looks like a nasty
integral but there is good news—we won’t ever need to compute
it! Instead we can use the table of Laplace transforms and their
inverses in table 4.1. However, in order to use the tables we
must generally break up whatever Laplace transform we have into
smaller, bite-size chunks. For this we use partial fractions.

4.6 Partial-fraction expansions

We can use partial fraction expansions to break up a function
to make it easier to take its inverse Laplace transform. Start-
ing with a rational function G(s), the partial-fraction expansion
decomposes it into a sum of lower-order rational functions:
dm sm + dm−1 sm−1 + · · · + d0
G(s) =
cn sn + cn−1 sn−1 + · · · + c0
n
X bi
= .
i=1
(s − p i )
4.6. Partial-fraction expansions 55

time Laplace property

αx1 (t) + βx2 (t) αX1 (s) + βX2 (s) superposition

x(t − τ )u(t − τ ) e−sτ X(s) time shift

1 s
x(at) X time scaling
a a
e−at x(t) X(s + a) frequency shift
dx(t)
sX(s) − x(0) first derivative
dt
d2 x(t)
s2 X(s) − sx(0) − ẋ(0) second derivative
dt2
Pn
x(n) (t) sn X(s) − k=1 sn−k x(k−1) (0) nth derivative
Rt 1
0
x(τ )dτ X(s) integration
s
Table 4.2: Important properties of the Laplace transform.
4.6. Partial-fraction expansions 56

Where pi are called the poles of G(s) and the coefficients bi are
sometimes called the residues. (We assume for now that there are
no repeated poles.)

4.6.1 Distinct poles

If all poles pi are distinct and real, pi 6= pj for all i 6= j then the
partial fraction is of the form
n
X bi
G(s) = .
i=1
(s − pi )

1
Example: G(s) = s2 +s

We first factor the denominator:

1
G(s) = .
s(s + 1)

The poles pi are distinct. Using the ’cover up’ method we have

1
b1 = = 1.
s(s
A + 1)

s=0

and
1
b2 = XX = −1.
s
(s
+X1)
X s=−1

Finally, write the partial fraction form of G(s):

n
X bi
G(s) =
i=1
(s − pi )
b1 b2
= +
s s+1
1 1
= − .
s s+1
4.6. Partial-fraction expansions 57

4.6.2 Repeated poles

If any of the poles pi are repeated then the k-repeated poles gen-
erate terms of the form
b1 b2 bk
+ + ··· + (4.14)
(s − p) (s − p)2 (s − p)k
1
For example, G(s) = (s+1) 2 has two repeated poles and will there-

fore contain two terms in its partial fraction expansion:

1 b1 b2
G(s) = 2
= + .
(s + 1) (s + 1) (s + 1)2
How do we find these constants? This is best seen with an exam-
ple.

1
Example: G(s) = (s+2)2 (s+1)

We have
1 b1 b2 c1
G(s)) = 2
= + 2
+ .
(s + 2) (s + 1) (s + 2) (s + 2) (s + 1)
Multiply both sides through by the denominator of G(s):

1 = b1 (s + 2)(s + 1) + b2 (s + 1) + c1 (s + 2)2 .

We can find c1 easily by setting s = −1:

1 = c1 (−1 + 2)2 =⇒ c1 = 1

Setting s = −2 gives

1 = b2 (−2 + 1) =⇒ b2 = −1

Then setting s to any other number will give us b1 . Let’s go for

s = 0:
1 = 2b1 + b2 + 4c1 = 2b1 + 3 =⇒ b1 = −1
4.7. Understanding complex roots 58

and so we have
1 1 1
G(s) = − 2
−
(s + 1) (s + 2) (s + 2)

and the solution in time is

L−1 {G(s)} = e−t − e−2t − te−2t

The final term has a ‘t’ in there because (table 4.1)

tm −at
 
−1 1
L = e
(s + a)m+1 m!

which, for m = 1, gives

 
−1 1
L = te−at .
(s + a)2

4.6.3 Second-order terms in the denominator

If there are second-order terms in the denominator:
1 A Bs + C
G(s) = = + .
s(s2 + s + 10 s s2 + s + 10

4.7 Understanding complex roots

Let’s consider a term in the partial fraction expansion of the
form
b
T (s) =
s−p
with p = σ + jω and b = |b| ejφ , i.e. p and b are complex. Then
in the time domain this term gives

L−1 {T (s)} = be(σ+jω)t ,

4.7. Understanding complex roots 59

which is complex. Why is this complex when we expect the re-

sponse in time to be real? Because there exists another term
which is its complex conjugate:
b∗
T ∗ (s) = .
s − p∗
Examining the pair together:
b∗
 
−1 ∗ −1 b
L {T (s) + T (s)} = L +
s − p s − p∗
= be(σ+jω)t + b∗ e(σ−jω)t
= |b|eσt (ejφ ejωt + e−jφ e−jωt )
= |b|eσt (ej(φ+ωt) + e−j(φ+ωt) ).

Now let φ̃ = φ + ωt. By Euler’s formula ej φ̃ = cos φ̃ + j sin φ̃

then
L−1 (T (s) + T ∗ (s)) = |b|eσt (cos(φ̃) + j sin(φ̃) + cos(−φ̃) + j sin(−φ̃))
= 2|b|eσt cos(φ̃)
= 2|b|eσt cos(φ + ωt)
which is real. Alternatively we can see the same thing by using
Laplace transforms directly:
b b∗
T (s) + T ∗ (s) = +
s − p s − p∗
···
= 2
s − sp − sp∗ + pp∗
···
=
(s + a)2 + ωo2
for some constants a and ωo . Then in time we have
L−1 {T (s) + T ∗ (s)} = e−at (c1 sin ωo t + c2 cos ωo t)
for some constants c1 and c2 (see table 4.1).
4.8. Matlab commands 60

4.8 Matlab commands

• s = tf(’s’); to set up the Laplace variable
• residue(); for partial fractions

4.9 Links and further reading

• xkcd: eiθ
• Astrom & Murray § 8.5 Laplace transforms
5 Second-order systems

In this chapter we will apply Laplace-transform techniques to

analyse second-order systems. We will look at how to put second-
order systems into what is called standard form; and then look at
second-order systems’ free response (i.e. their response to initial
conditions) and their forced response (in particular, their response
to an impulsive input or to a step input).

5.1 Standard form

We have already seen two examples of second-order systems in
chapter 1. These were the mass-spring damper system:

mq̈(t) + cq̇(t) + kq(t) = f (t). (5.1)

and the RLC circuit:

We will now look at how to put both of these systems into the
same template equation which we call standard form. The ad-
vantage of standard form is that we can analyse different systems
in the same way. For example, a mass-spring-damper system and
an RLC circuit are very different physical systems and yet they
display similar behaviours: both can oscillate, for example. Stan-
dard form brings out these similarities very nicely. The standard
form equation is

q̈(t) + 2ζωn q̇(t) + ωn2 q(t) = Ku(t) (5.2)

61
5.2. Laplace transforms 62

where ωn is the undamped natural frequency, ζ is called the damp-

ing ratio, and K is a gain. Here q could be the displacement of
a mass-spring-damper, or the voltage in an RLC circuit, to give
two examples.
How can we put the mass-spring damper system into the standard
form (5.2)? We start with the governing equation (5.1) and divide
both sides by the mass m:
c k 1
q̈(t) + q̇(t) + q(t) = f (t).
m m m
Notice that this already looks a lot like the standard form (5.2).
In fact comparing the two equations, we can see straight away
that we must have
c k 1
2ζωn = ; ωn2 = ; K= .
m m m
Solving for ωn and ζ, we find the three equations
r
k c 1
ωn = ; ζ= √ ; K= .
m 2 km m
Notice that the undamped natural frequency, ωn , has come out
to be the same as that found in chapter 1 (eq. (1.2)) where we as-
sumed no damping—so calling it the undamped natural frequency
is well justified. Notice also that the damping ratio is zero when
there is no damping, i.e. ζ = 0 when c = 0.

5.2 Laplace transforms

We now move into the Laplace domain. Taking Laplace trans-
forms of (5.2), including terms corresponding to the initial con-
ditions, we have

(s2 + 2ζωn s + ωn2 )Q(s) = KF (s) + q0 s + q̇0 + 2ζωn q0 ,

5.3. Poles 63

or, rearranging for the displacement Q(s):

K q0 s + q̇0 + 2ζωn q0
Q(s) = F (s) + 2 . (5.3)
s2 2
+ 2ζωn s + ωn s + 2ζωn s + ωn2
The first term in (5.3) is called the transfer function. It tells
us how the displacement of the mass, Q(s), is related to the ex-
ternal forcing, F (s), in the Laplace domain. What about the
second term in (5.3)? This represents the response due to any
initial conditions—in particular the initial displacement, q0 , and
the initial velocity, q̇0 . Equation (5.3) tells us everything about
the system’s response.

5.3 Poles
The two terms on the right-hand side of (5.3) share the same
denominator, which suggests that it has some special significance.
Its significance is that its roots (its poles) give us the natural
frequencies of the system. The poles pi are therefore given by the
roots of the denominator, s2 + 2ζωn s + ωn2 , which are
p
p1,2 = −ζωn ± ωn ζ 2 − 1. (5.4)

These two poles can be real or complex depending on the value

of the damping ratio ζ. For 0 ≤ ζ < 1, the poles are complex and
oscillation occurs. We call this underdamped. For ζ > 1 the poles
are real and there is no oscillation. We call this overdamped. For
ζ = 1, the poles are real and repeated and there is no oscillation.
We call this critically damped.
Program poles.m plots the positions of the poles in the complex
plane for a range of damping ratios ζ. The results are plotted
in figure 5.1. The poles are purely imaginary for zero damping
(ζ = 0); follow a circular path for 0 < ζ < 1; and are real for
ζ ≥ 1.
5.4. Free response 64

1
0.8
0.6
0.4
0.2
)
1,2

0
Im(λ

−0.2
−0.4
−0.6
−0.8
−1

−1 −0.8 −0.6 −0.4 −0.2 0 0.2

Re(λ1,2)

% poles.m plot the poles of a second-order system as the damping ratio is varied

omn = 1; % UNDAMPED natural frequency

zeta = [0:0.01:3]'; % range of damping ratios (vector)

% initialize p1 & p2, the vectors containing the poles

% (p2 is the complex conjugate of p1)
p1 = zeros(length(zeta),1);
p2 = zeros(length(zeta),1);

% for each value of zeta, find the two poles

for i=1:length(zeta)
p1(i) = -zeta(i)*omn + omn*sqrt(zeta(i).ˆ2-1);
p2(i) = -zeta(i)*omn - omn*sqrt(zeta(i).ˆ2-1);
end

% plot
figure, hold on, grid on
plot(p1,'b'), plot(p2,'r')
% plot first entry (corresponding to zeta=0) with a cross
plot(p1(1),'bx'), plot(p2(1),'rx')
axis([-1.2*omn 0.2*omn -1.2*omn 1.2*omn])
xlabel('Re(p {1,2})'), ylabel('Im(p {1,2})')
set(gca,'Box','On','XMinorTick','On','YMinorTick','On');

Figure 5.1: Output from poles.m.

5.4 Free response
Let’s now use (5.3) to see how the system responds without any
forcing: that is, with F (s) = 0. This is called the system’s free
response. We could also call it the response to an initial condi-
5.4. Free response 65

tion since any response is due purely to the initial values of the
displacement q0 and the velocity q̇0 . Without forcing, F (s) = 0
we have
q0 s + q̇0 + 2ζωn q0
Q(s) = . (5.5)
s2 + 2ζωn s + ωn2

Write this as
c1 s + c0
Q(s) = (5.6)
s2 + 2ζωn s + ωn2

with c1 = q0 and c0 = q̇0 + 2ζωn q0 . We write (5.5) in this more

general form because it means that we can get the impulse re-
sponse from (5.6) ‘for free’, which will be useful in § 5.6.
Now we want to take inverse Laplace transforms of (5.6). This
will give us the solution in time. But we need to be careful because
the way in which we find the inverse Laplace transform depends
in an important way on the damping ratio ζ. We now look at
the three cases of underdamped (ζ < 1), overdamped (ζ > 1) and
critically damped (ζ = 1).

5.4.1 Underdamped ζ < 1

Underdamped means that the damping ratio satisfies ζ < 1. This
gives rise to a complex-conjugate pair of poles:
p
p1,2 = −ζωn ± ωn ζ 2 − 1
p
= −ζωn ± jωn 1 − ζ 2

The second line follows because ζ 2 −1 is negative and so its square-

root is imaginary.
The key to solving the underdamped case is to rearrange (5.6) in
such a way that we can use the Laplace transforms in table 4.1
5.4. Free response 66

directly. This requires a bit of manipulation:

c1 s + c0
Q(s) =
s2 + 2ζωn s + ωn2
c1 s + c0
=
(s + ζωn )2 + ωn2 − ζ 2 ωn2
c1 s + c0
= p
(s + ζωn ) + (ωn 1 − ζ 2 )2
2

c1 s + c0
=
(s + σ)2 + ωd2
c0 − σc1
 
s+σ ωd
= c1 2
+
2
(s + σ) + ωd ωd (s + σ)2 + ωd2
p
where ωd = ωn 1 − ζ 2 is the damped natural frequency and σ = ζωn .
Now substituting in c1 = q0 and c0 = q̇0 + 2ζωn q0 and taking in-
verse Laplace transforms we get

q̇0 + ζωn q0 −ζωn t

q(t) = q0 e−ζωn t cos ωd t + e sin ωd t. (5.7)
ωd

5.4.2 Critically damped ζ = 1

Critically damped means that ζ = 1. This gives rise to a pair of re-
peated poles since the denominator becomes s2 + 2ωn s + ωn2 = (s + ωn )2 :

p1 = −ωn
p2 = −ωn ,

and the response to an initial condition (5.6) becomes

c1 s + c0
Q(s) =
(s + ωn )2
A1 A2
= + .
s + ωn (s + ωn )2
5.4. Free response 67

Solving for A1 and A2 we find

A1 = c1
A2 = c0 − ωn c1 .
Using c1 = q0 and c0 = q̇0 + 2ζωn q0 and taking inverse Laplace
transforms we get

q(t) = q0 e−ωn t + (ωn q0 + q̇0 )te−ωn t (5.8)

Notice the ‘t’ in the second term. This appears because of the
repeated poles (see the Laplace transforms in table 4.1).

5.4.3 Overdamped ζ > 1

Overdamped means that ζ > 1. The two poles are real:
p
p1 = −ζωn − ωn ζ 2 − 1
p
p2 = −ζωn + ωn ζ 2 − 1
and the partial fraction expansion takes the form
c1 s + c0
Q(s) =
(s − p1 )(s − p2 )
A1 A2
= + .
s − p1 s − p2
Solving for A1 and A2 gives
c1 p 1 + c0
A1 =
p1 − p 2
c1 p 2 + c0
A2 = .
p2 − p 1
We could go further by substituting in values for p1 , p2 , c0 and
c1 , but we won’t bother. Taking inverse Laplace transforms we
get
q(t) = A1 ep1 t + A2 ep2 t (5.9)
5.5. Transfer function 68

with A1 , A2 , p1 , p2 as given above. The important point, which

holds irrespective of the precise values of A1 and A2 , is that there
is no oscillation, only exponential growth or decay.
Program secondFree.m plots the free response for a particular
initial condition and for a range of damping ratios. The results
are shown in figure 5.2. For zero damping, oscillation continues
indefinitely. For 0 < ζ < 1, the response is still oscillatory, but
also decays. For ζ ≥ 1, the response no longer oscillates, instead
showing pure decay.

5.5 Transfer function

As mentioned in § 5.2, equation (5.3) tells us everything about the
system’s response. That is, it tells us about the system’s response
to forcing F (s) and to initial conditions (q0 q̇0 ). We can learn
more about the system’s response to forcing (its forced response)
by setting the initial condition to zero. Then (5.3) becomes

K
Q(s) = F (s), (5.10)
s2 + 2ζωn s + ωn2
or
Q(s) = G(s)F (s).
G(s) is called the transfer function. For the second-order system
considered in this chapter it is given by

K
G(s) =
s2 + 2ζωn s + ωn2

The transfer function describes how an input F (s) is converted

into an output Q(s) by the system G(s)—all in the Laplace do-
main. We can represent this pictorially in the following way:
5.6. Impulse response 69

t
Io
G(s) takes as its input F (s) and gives as its output Q(s).

5.6 Impulse response

The impulse response means that we set the input f (t) = δ(t),
where δ(t) is impulse function (see § 4.3.3). The Laplace transform
of δ(t) is simply L{δ(t)} = 1 and so the forcing in the Laplace
domain is simply
F (s) = 1.
We can substitute this into the expression for the forced response
of Q(s) (5.10) to give

K
Q(s) = G(s)F (s) = G(s) · 1 = G(s) = .
s2 + 2ζωn s + ωn2

Now we have Q(s) but we want q(t). For this we need to take
inverse Laplace transforms but, as for the free response, we need
to be careful to treat the underdamped, critically-damped and
overdamped cases separately. The good thing about having used
the general form (5.6) is that we can now reuse it. Comparing
(5.6) to the expression for Q(s) above, we see that if we set

c0 = K
c1 = 0

then we can reuse the results for the free response for the impulse
response. Since this involves simply repeating the steps we have
already performed (with different values of c0 and c1 we state only
the final results here.
5.7. Step response 70

5.6.1 Underdamped
The underdamped impulse response is

K −ζωn t
q(t) = e sin(ωd t)
ωd

5.6.2 Critically damped

The critically damped impulse response is

q(t) = Kte−ωn t

5.6.3 Overdamped
The overdamped impulse response is

K
q(t) = (ep1 t − ep2 t )
p1 − p2

Program secondImpulse.m plots the impulse response for a range

of damping ratios. The results are shown in figure 5.3. For zero
damping, oscillation continues indefinitely. For 0 < ζ < 1, the re-
sponse is still oscillatory, but also decays. For ζ ≥ 1, the response
no longer oscillates, instead showing pure decay.

5.7 Step response

For the step response we set the forcing equal to the step function,
f (t) = u(t) = 1(t). That is, f (t) = 1 for t > 0 and f (t) = 0
otherwise. Taking the Laplace transform of f (t) we get
1
F (s) = .
s
5.7. Step response 71

We substitute this into the expression for the forced response of

Q(s) (5.10) to give

Q(s) = G(s)F (s)

K 1
= 2 2
×
s + 2ζωn s + ωn s
K
= .
s(s + 2ζωn s + ωn2 )
2

We will now spend some time on the step response for the under-
damped case and carefully characterize it. Before that, we first
consider the two other cases (overdamped and critically damped)
and get them out of the way.

5.7.1 Critically damped

The poles are repeated and real, p1 = p2 = −ωn and
A B1 B2
Q(s) = + + .
s s + ωn (s + ωn )2
Then
q(t) = A + B1 e−ωn t + B2 te−ωn t
with
K K K
A= , B1 = − , B2 = − .
ωn2 ωn2 ωn

5.7.2 Overdamped
The poles are real and distinct and
A B1 B2
Q(s) = + + .
s s − p1 s − p2
Then
q(t) = A + B1 ep1 t + B2 ep2 t
5.7. Step response 72

with
K K K
A= , B1 = , B2 = .
p 1 p2 p1 (p1 − p2 ) p2 (p2 − p1 )

5.7.3 Underdamped step response

In the underdamped case we have a complex-conjugate pair of
poles. Learning our lesson from the solution to the underdamped
free response, we use the same trick and write Q(s) as
K
Q(s) = (5.11)
s((s + σ)2 + ωd2 )
p
with σ = ζωn and ωd = ωn 1 − ζ 2 . Notice though that there
is an extra ‘s’ in the denominator because of the step input,
F (s) = 1/s. We choose the following partial fraction decompo-
sition so that we can re-use the results for the free response:
A c1 s + c0
Q(s) = + . (5.12)
s (s + σ)2 + ωd2
The first term is new but taking its inverse Laplace transform is
simple: L−1 (A/s) = A. The second term is the same as (5.6) so
we already know how to solve it. We just need to work out the
values of the constants c0 and c1 . To do so we equate equations
(5.11) and (5.12):
K A c1 s + c0
= + , (5.13)
s((s + σ)2 + ωd2 ) s (s + σ)2 + ωd2
multiply both sides by the denominator and rearrange:
K = (A + c1 )s2 + (2Aσ + c0 )s + A(σ 2 + ωd2 ). (5.14)
Now equating powers of s we can solve for A, c0 and c1 :
K −K
A= , c0 = −2Kζωn , c1 = (5.15)
ωn2 ωn2
5.8. Characterizing the underdamped step response 73

Substituting these into the expression for the free response and
tidying up we find
!!
K ζ
q(t) = 1 − e−ζωn t cos ωd t + p sin ωd t (5.16)
ωn2 1 − ζ2

We will spend some time to characterize the underdamped step

response (5.16). It will sometimes be convenient to write (5.16)
in a slightly different form by using the following identity:

cos(θ) + β sin(θ) = M sin(θ + φ),

p
with M = β 2 + 1 and φ = tan−1 ( β1 ). Applying this to (5.16)
we can rewrite the underdamped step response (5.16) as
!
K e−ζωn t
q(t) = 2 1− p sin(ωd t + φ) (5.17)
ωn 1 − ζ2
√ 
−1 1−ζ 2
with φ = tan ζ
. The advantage of this representation
is that we only have the sine term instead of both sine and cosine
terms.
Program stepResp.m plots the step response for a range of damp-
ing ratios. The results are shown in figure 5.4.

5.8 Characterizing the underdamped step re-

sponse
Let’s now characterize some important features of the under-
damped step response (5.17) shown in figure 5.4. We will look
at peak time; overshoot; rise time; and settling time.
5.8. Characterizing the underdamped step response 74

5.8.1 Peak time Tp

The peak time Tp is the time required to reach the first peak of
the overshoot. We can find it by taking the derivative of (5.16)
and setting it to zero:
   
K −σt σ −σt
q̇(t) = 2 σe cos ωd t + sin ωd t − e (−ωd sin ωd t + σ cos ωd t)
ωn ωd
K −σt σ 2
 
= 2e sin ωd t + ωd sin ωd t
ωn ωd
K −σt σ 2
 
= 2e + ωd sin ωd t.
ωn ωd

This is satisfied when

sin ωd Tp = 0
so the first non-zero Tp is

π
Tp =
ωd

5.8.2 Overshoot Mp
The overshoot Mp is the amount by which the maximum response
qmax overshoots the final value qfinal = limt→∞ q(t) = K/ωn2 . We
express the overshoot as a ratio:
qmax − qfinal
Mp = ,
qfinal
so that an overshoot of 0.05, for example, means that the max-
imum response qmax is 5 % larger than the final response qfinal .
The overshoot and the peak time are related: Mp tells us the
size of the peak response, while Tp tells us the time at which it
happens. Therefore to get the overshoot Mp we set t = Tp in the
5.8. Characterizing the underdamped step response 75

underdamped step response (5.16) and find a maximum response

of
K
qmax = 2 (1 + e−σπ/ωd ).
ωn
The overshoot is then
qmax − qfinal
Mp =
qfinal
= 1 + e−σπ/ωd − 1
− σπ
=e ω d

ζω
√n π
−
=e ωn 1−ζ 2 .
Simplifying:
√−ζπ
1−ζ 2
Mp = e (5.18)
Notice what this says: that the overshoot depends only on the
damping ratio ζ! If someone gives us a step response, we could
therefore estimate ζ from knowledge of Mp alone by rearranging
(5.18). (Remember that to get Mp we need to know both qmax
and qfinal .)

5.8.3 Rise time Tr

The rise time is the time required for the response to rise from
10% to 90%; or from 5% to 95%; or from 0% to 100% of its final
value. For underdamped second-order systems, the 0% to 100%
rise time is normally used.1 We therefore obtain the rise time Tr
by setting q(Tr ) = qfinal :
  
K −σTr σ
q(Tr ) = 2 1 − e cos ωd Tr + sin ωd Tr
ωn ωd
K
= qfinal = 2 .
ωn
1
(For overdamped systems, the 10% to 90% rise time is normally used.
5.8. Characterizing the underdamped step response 76

This requires that

σ
cos(ωd Tr ) + sin(ωd Tr ) = 0
ωd
or p
ωd 1 − ζ2
tan(ωd Tr ) = − = − . (5.19)
σ ζ
If we use (5.19) as it is then we will find a negative rise time!
This is just an artefact of the tan−1 function, which repeats with
a period of π. To get around this we use the fact that tan(x) =
tan(x + π). Using this (and the fact that tan(−x) = − tan(x)) we
get a final answer of
p !!
1 1 − ζ 2
Tr = π − tan−1
ωd ζ

5.8.4 δ-settling time Tδ

The δ-settling time is the time required to reach and stay within
δ of the final value qfinal . To simplify things a little, let’s assume
that K = ωn2 so that qfinal = 1. The question is then: After how
long does q(t) always lie in the range

1 − δ ≤ q(t) ≤ 1 + δ.

To answer we use the form (5.17) of the step response (with re-
member K = ωn2 ):

e−ζωn t
q(t) = 1 − p sin(ωd t + φ).
1 − ζ2
The two envelopes of the overall curve are then
e−ζωn t
1± p .
1 − ζ2
5.9. Links and further reading 77

(The sine wave part of the solution wiggles within this envelope.)
We want to set this envelope equal to (1±δ) and so we have

e−ζωn t
±δ = ± p
1 − ζ2

with the time set to t = Tδ , since this is the definition of Tδ . Sub-

stituting this in and rearranging gives a δ-settling time of
p
− ln(δ 1 − ζ 2 )
Tδ =
ζωn

5.9 Links and further reading

• Impact testing for modal analysis
• Create custom impulse responses in REAPER
• da Vinci robot stitches a grape back together
• High-speed robot hand
• Oscillating buildings
• What the phugoid is longitudinal stability?
5.9. Links and further reading 78

x(t) 0

-1
0 0.5 1 1.5 2 2.5 3 3.5 4
t
1
x(t)

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t
1
x(t)

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t

% secondFree.m: free response of second-order system

omn = 2*pi*1; % UNDAMPED natural frequency
t = [0:0.001:1] * 4; % time
u = zeros(length(t),1); % input is zero
x0=1; v0=0; % initial condition (displacement & velocity)
zeta = [0:0.05:1.2]; % vector of damping ratios
for i=1:length(zeta)
if zeta(i) < 1; % underdamped
omd = omn * sqrt(1-zeta(i)ˆ2); % damped natural frequency
x = (x0*cos(omd*t) + (v0+zeta(i)*omn*x0)*sin(omd*t)/omd) .* exp(-zeta(i)*omn*t);
subplot(3,1,1), hold on
elseif zeta(i) == 1 % critically damped
x = x0*exp(-omn*t) + (omn*x0 + v0)*t.*exp(-omn*t);
subplot(3,1,2), hold on
else % overdamped
p1 = -zeta(i)*omn - omn*sqrt(zeta(i)ˆ2-1); A1 = ((p1+2*zeta(i)*omn)*x0 + v0) / (p1-p2);
p2 = -zeta(i)*omn + omn*sqrt(zeta(i)ˆ2-1); A2 = ((p2+2*zeta(i)*omn)*x0 + v0) / (p2-p1);
x = A1*exp(p1*t) + A2*exp(p2*t);
subplot(3,1,3), hold on
end
plot(t, x)
% check solution is correct using state-space and lsim
%A = [0 1; -omnˆ2 -2*zeta(i)*omn]; B = [0 1]'; C = [1 0];
%y = lsim(ss(A, B, C, 0), u, t, [x0 v0]');
%plot(t, y, 'r--')
xlabel('t'), ylabel('x(t)'), set(gca,'Box','On','XMinorTick','On','YMinorTick','On');
end

Figure 5.2: secondFree.m with q0 = 1 and q̇0 = 0. Top: ζ < 1.

Middle: ζ = 1. Bottom: ζ > 1.
5.9. Links and further reading 79

0.2

x(t) 0

-0.2
0 0.5 1 1.5 2 2.5 3 3.5 4
t

0.05
x(t)

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t
0.05
x(t)

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t

% secondImpulse.m: plot response to an impulseive input

omn = 2*pi*1; % UNDAMPED natural frequency
K = 1; % gain
t = [0:0.001:1] * 4; % time
zeta = [0:0.05:1.2]; % vector of damping ratios
for i=1:length(zeta)
if zeta(i) < 1 % underdamped
omd = omn * sqrt(1-zeta(i)ˆ2); % damped natural frequency
x = K/omd * exp(-zeta(i)*omn*t).*sin(omd*t);
subplot(3,1,1), hold on
elseif zeta(i) == 1 % critically damped
x = K*t.*exp(-omn*t);
subplot(3,1,2), hold on
else % overdamped
p1 = -zeta(i)*omn - omn*sqrt(zeta(i)ˆ2-1);
p2 = -zeta(i)*omn + omn*sqrt(zeta(i)ˆ2-1);
x = K/(p1-p2) * (exp(p1*t) - exp(p2*t));
subplot(3,1,3), hold on
end
plot(t, x)
% check solution is correct using state-space and impulse()
%A = [0 1; -omnˆ2 -2*zeta(i)*omn]; B = [0 1]'; C = [1 0];
%y = impulse(ss(A,B,C,0), t);
%plot(t, y, 'r--')
xlabel('t'), ylabel('x(t)')
set(gca,'Box','On','XMinorTick','On','YMinorTick','On');
end

Figure 5.3: secondImpulse.m. Top: ζ < 1. Middle: ζ = 1.

Bottom: ζ > 1.
5.9. Links and further reading 80

x(t) 1

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t
1
x(t)

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t
1
x(t)

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5 4
t

% secondStep.m: plot response to a step input

omn = 2*pi*1; % UNDAMPED natural frequency
K = 1; % gain
t = [0:0.001:1] * 4; % time
zeta = [0:0.05:1.2]; % vector of damping ratios
for i=1:length(zeta)
if zeta(i) < 1 % underdamped
omd = omn * sqrt(1-zeta(i)ˆ2); % damped natural frequency
x = K/omnˆ2 * ( 1 - exp(-zeta(i)*omn*t) .* ...
(cos(omd*t) + (zeta(i)/sqrt(1-zeta(i)ˆ2))*sin(omd*t)) );
subplot(3,1,1), hold on
elseif zeta(i) == 1 % critically damped
A = K/omnˆ2; B1 = -K/omnˆ2; B2 = -K/omn;
x = A + B1*exp(-omn*t) + B2*t.*exp(-omn*t);
subplot(3,1,2), hold on
else % overdamped
p1 = -zeta(i)*omn - omn*sqrt(zeta(i)ˆ2-1);
p2 = -zeta(i)*omn + omn*sqrt(zeta(i)ˆ2-1);
A = K/(p1*p2); B1 = K/(p1*(p1-p2)); B2 = K/(p2*(p2-p1));
x = A + B1*exp(p1*t) + B2*exp(p2*t);
subplot(3,1,3), hold on
end
plot(t, x)
% check solution is correct using state-space and impulse()
%A = [0 1; -omnˆ2 -2*zeta(i)*omn]; B = [0 1]'; C = [1 0];
%y = step(ss(A,B,C,0), t);
%plot(t, y, 'r--')
xlabel('t'), ylabel('x(t)')
set(gca,'Box','On','XMinorTick','On','YMinorTick','On');
end

Figure 5.4: stepResp.m. Top: ζ < 1. Middle: ζ = 1. Bottom:

ζ > 1.
6 First-order systems

Having considered second-order systems we will now go back-

wards, so to speak, and consider first-order system. What is a
first-order system? To fix ideas we start by considering the spe-
cific example of a cup of tea.

6.1 Cup of tea as a first-order system

Imagine a cup of tea sitting on the kitchen bench that has just
been made. (You can replace the tea with any drink of your
choice. The important thing is that it is warmer than its sur-
roundings.) What will happen to it? Anyone can tell you that
it will cool down. If you have studied a certain amount of math-
ematics then you can also say something about the way it will
cool: in such a way that the rate of change of its temperature is
proportional to its temperature at that time (or rather the tem-
perature difference between the tea and its surroundings). The
differential equation that describes in mathematics what we just
said in words is
Ṫ (t) = −αT (t).
For this to make sense we need α > 0 so that the hot tea cools (or
a cold drink warms). It must be stressed that the temperature
T (t) here represents the difference in temperature between the
tea and it surroundings. I should have perhaps called it ∆T (t),
but that would have introduced ∆ signs everywhere. Now let’s

81
6.2. Laplace transforms 82

set the constant α to α = 1/τ , where τ is called the time constant

(the reason for doing this will become clear shortly):
1
Ṫ (t) = − T (t). (6.1)
τ
The solution to this equation, which you should already know,
is
T (t) = T0 e−t/τ . (6.2)
This can be verified by substituting (6.2) back into (6.1).
The cup of tea is a first-order system. A first-order system wants
to return to equilibrium and does so with exponential decay, as
described by equation (6.2). (This is true provided it is stable,
which requires τ > 0. If it is unstable then it will display expo-
nential growth instead of exponential decay.) The ‘first’ refers to
the fact that the highest time derivative is one, which also means
its transfer function has a single pole.
There is an important simplification that we made in the analysis
of the cup of tea which will apply to all of the first-order systems
we study. We did not consider any variation in the temperature
in different parts of the cup. Instead we set the temperature to
be a single number that varies in time but not in space. This
approximation will be fine provided that we are interested in the
temperature of the cup of tea as a whole—its average temperature.
If we also wanted to describe how the temperature at the top of
the cup differs from the temperature at the bottom of the cup,
then (6.1) would be an inadequate description of the tea.

6.2 Laplace transforms

Let’s see what we get when we take Laplace transforms of a first-
order system. For this consider a general first-order system,
τ ẏ(t) + y(t) = Ku(t).
6.3. Poles 83

The key difference from the cup of tea that we just considered
is that there is now an input u(t). (For the cup of tea, this
could be heating from a stove-top, for example.) Taking Laplace
transforms, including the effect of any initial conditions,

Y (s)(τ s + 1) = τ y0 + KU (s)

or
K τ y0
Y (s) = U (s) + . (6.3)
τs + 1 τs + 1
The first term is the transfer function and the second term repre-
sents the response to the initial condition y0 .

6.3 Poles
Both terms in (6.3) share the same denominator, which suggests
it has some special significance. The roots of this denominator,
τ s + 1 = 0, give the poles of the first-order system:

p = −1/τ.

If τ > 0 then the system is stable. If τ < 0 then the system is

unstable. The response of the system is always non-oscillatory
because the imaginary part of the single pole is always zero.

6.4 Free response

Taking the inverse Laplace transform of the second term in (6.3)
gives the response in time due to the initial condition y0 . This
gives
y(t) = y0 e−t/τ , (6.4)
which recovers the solution for the unforced cup of tea (6.2). This
is what we should expect—but we got there using Laplace trans-
forms this time.
6.5. Transfer function 84

6.4.1 Time constant τ

What is the significance of the time constant τ ? Substituting
t = τ into (6.4) we find

y(τ ) = e−1 y0 = 0.368 y0 .

So the time constant is the time taken for the free response to
decay to about 37 % of its initial value. There is nothing special
about the number 0.368 other than it is equal to e−1 . The time
constant is simply a characteristic time that allows us to compare
one first-order system to another one.

6.5 Transfer function

Setting the initial condition to zero in (6.3) we see that the trans-
fer function of a first-order system is

Y (s) K
G(s) = = (6.5)
U (s) τs + 1

6.6 Impulse response

Remember that for the impulse response we set u(t) = δ(t), the
inverse Laplace transform of which is simply L−1 {δ(t)} = 1. Sub-
stituting this into the general solution (6.3) (with y0 = 0):

K K/τ
Y (s) = ·1=
τs + 1 s + 1/τ

and taking inverse Laplace transforms:

K −t/τ
y(t) = e (6.6)
τ
6.7. Step response 85

6.7 Step response

The inverse Laplace transform of the step function is L−1 {u(t)} = 1/s.
Substituting this into the general solution (6.3) (with y0 = 0) we
get
K
Y (s) = .
s(τ s + 1)

The partial fraction decomposition is

A B
Y (s) = +
s τs + 1
with A = K and B = −Kτ :
K Kτ
Y (s) = −
s τs + 1
K K
= − .
s s + τ1

Taking inverse Laplace transforms, the solution in time is

y(t) = K(1 − e−t/τ ) (6.7)

6.8 Matlab code

Program firstResp.m plots the free response, the impulse re-
sponse and the step response for first-order systems with time
constants of τ = 1, 2, 3, 4, and 5. The results are shown in figure
6.1. For the impulse response, we have set α = τ for every case
so that the initial response is 1; For the step response, we have
set α = 1 for every case so that the final response is 1.
6.8. Matlab code 86

1
x(t)

0.5

0
0 2 4 6 8 10
t
1
x(t)

0.5

0
0 2 4 6 8 10
t
1
x(t)

0.5

0
0 2 4 6 8 10
t

%% firstResp.m: plot response of first-order system

t = [0:0.001:1] * 10; % time
tau = [1:1:5]; % range of time constants
for i=1:length(tau)
% free response
y0 = 1; % initial condition
y = y0 * exp(-t/tau(i));
subplot(3,1,1), hold on
plot(t, y), xlabel('t'), ylabel('y(t)')
set(gca,'Box','On','XMinorTick','On','YMinorTick','On');
% impulse response
K = tau(i); % so that initial response is 1
y = K/tau(i) * exp(-t/tau(i));
subplot(3,1,2), hold on
plot(t, y), xlabel('t'), ylabel('y(t)')
set(gca,'Box','On','XMinorTick','On','YMinorTick','On');
% step response
K = 1; % so that final response is 1
y = K * (1 - exp(-t/tau(i)));
subplot(3,1,3), hold on
plot(t, y), xlabel('t'), ylabel('y(t)')
set(gca,'Box','On','XMinorTick','On','YMinorTick','On');
end

Figure 6.1: Output from firstResp.m. Top: free response.

Middle: impulse response. Bottom: step response.
7 Transfer functions and block diagrams

Having looked at how to form the transfer function of first-order

and second-order systems, you can perhaps guess how we might
apply the procedure to more general nth -order systems. That is
the subject of this chapter. We will also introduce the concept of
a block diagram.

7.1 Summary: first- and second-order systems

Before we move on to the general case, let’s summarize the trans-
fer functions of first- and second-order systems that we considered
in the previous two chapters. The transfer function of a first-order
system is
K K/τ
G(s) = = ,
τs + 1 s+p
which has a single pole at p = −1/τ . The transfer function of a
second-order system is
K
G(s) = ,
s2 + 2ζωn s + ωn2

which has a pair of poles. These can be real or complex depending

on the value of the damping ratio ζ (see § 5.3).

87
7.2. Transfer functions 88

7.2 Transfer functions

In the general case we can write the differential equation of a
linear system (assuming only one input u and one output y) in
the form

an y (n) (t) + · · · + a1 ẏ(t) + a0 y(t) = bm u(m) (t) + · · · + b1 u̇(t) + b0 u(t),

(7.1)
(n) th (m)
where y is the n time derivative of y, and similarly u is the
mth time derivative of u. Notice that this expression allows for
time derivatives not only of the output y(t) but also of the input
u(t). (The highest time derivative of y(t) is n and the highest
time derivative of u(t) is m.) Now take Laplace transforms of
both sides of (7.1), taking into account the initial conditions for
both y(t) and u(t):

(an sn + · · · + a1 s + a0 )Y (s) = (bm sm + · · · + b1 s + b0 )U (s)

+ f (y0 , ẏ0 , . . . , u0 , u̇0 , . . . )

and solve for Y (s) to get

bm s m + · · · + b1 s + b0 f (y0 , ẏ0 , . . . , u0 , u̇0 , . . . )

Y (s) = U (s) +
a s + · · · + a1 s + a0
n a sn + · · · + a1 s + a0
|n {z } | n {z }
forced response free response

or
Y (s) = G(s)U (s) + Q(s).
If we assume that the initial conditions are all zero—which is
what we do when considering the transfer function—then we get
simply
Y (s) = G(s)U (s)
G(s) is the transfer function. It characterizes the forced re-
sponse of a linear system.
7.2. Transfer functions 89

7.2.1 Three ways to write the transfer function

There are three ways in which we can write the transfer function
G(s). The first is to write G(s) as the ratio of two polynomi-
als:
bm s m + · · · + b1 s + b0
G(s) = (7.2)
an s n + · · · + a1 s + a0
The second way, which we have seen previously in chapter 4, is
to use a partial fraction expansion:1
n
X Ai
G(s) = (7.3)
i=1
s − pi

This is what we did in chapters 5 & 6 to solve for the free re-
sponse, impulse response and step response. The constants Ai
are sometimes called the residues of G(s) and pi are the poles
of G(s).
The third way, which we haven’t seen before, is to factor the
numerator and denominator of (7.2) using their roots:
Qm
(s − zi )
G(s) = K Qni=1 (7.4)
i=1 (s − pi )

The m roots of the numerator, zi , are called the zeros of G(s).

This leads us on to the idea of poles and zeros.

7.2.2 Poles and zeros

• The roots of the denominator, pi , are called the poles of G(s).
We have already seen poles in chapters 5 & 6. Remember that
the real part of the poles is related to stability and that the
1
This assumes no repeated poles. If there are repeated poles then we
would need additional terms in (7.3) to account for them. See § 4.6.
7.2. Transfer functions 90

imaginary part of the poles is related to the frequency of oscil-

lation.
• The roots of the numerator, zi , are called the zeros of G(s).
These zeros are a new concept—we haven’t come across them
before. Where did they come from? Tracing back through
the equations, we can see that they come from the derivatives
of the input. In the first- and second-order systems we have
considered in earlier chapters, there were never any derivatives
of the input, and so there weren’t any zeros. We will go through
some examples of systems with zeros in lectures.

7.2.3 Relative degree

The relative degree of a transfer function G(s) is simply the num-
ber of poles minus the number of zeros. For a transfer function
with n poles and m zeros (like in equations (7.2) & (7.4)) the
relative degree, r, is therefore

r = n − m.

For example, the second-order system of chapter 5 has a relative

degree of r = 2 − 0 = 2; and the first-order system of chapter 6
has a relative degree of r = 1 − 0 = 1.

7.2.4 Dominant-pole approximation

The dominant-pole approximation is a method for approximating
a high-order transfer function with a (simpler) transfer function of
lower order. We can use it provided that the real part of the poles
we keep are much closer to the origin than the real part of the
poles we throw away. In particular we require that, roughly,

|σkeep | < 5|σdiscard |

7.3. Initial-value & final-value theorems 91

where σkeep is the real part of the poles we keep and σdiscard is the
real part of the poles we throw away. We will look at an example
in lectures.

7.3 Initial-value & final-value theorems

The initial-value and final-value theorems can be useful to char-
acterize a system’s response just after an input has been applied
(initial-value theorem) or long after an input has been applied
(final-value theorem). They can save us time by allowing us to
say things about a system’s step response (for example) without
calculating the step response in detail.

7.3.1 Final-value theorem (FVT)

We start with the final value theorem (FVT). It says that

lim y(t) = lim sY (s) (7.5)

t→∞ s→0

In words this says that the behaviour in the time domain for t
very large is related to the behaviour in the Laplace domain for
s very small.
We need to be careful when applying (7.5) because there are only
certain conditions under which it can be used. For (7.5) to make
sense, the limit limt→∞ y(t) needs to exist. This means that y(t)
must settle down to a definite value for t → ∞. But this won’t
be true for all systems. For this to be true, all poles of sY (s)
must lie in the left half-plane. If this is true then limt→∞ y(t)
exists. If instead sY (s) has poles on the imaginary axis or in the
right half plane then limt→∞ y(t) does not exist because y(t) will
contain oscillating or exponentially growing terms. Let’s look at
some examples where the final value theorem can and cannot be
applied.
7.3. Initial-value & final-value theorems 92

Example
Consider the step response (U (s) = 1/s) of
 
1 1
G(s) = 2 .
s + 0.2s + 1 s + 10

The final-value theorem (7.5) can be applied because all of the

poles of sY (s) = sG(s) 1s = G(s) are in the left half-plane. We
get
 
1 1 1
lim y(t) = lim sY (s) = lim s × × 2
t→∞ s→0 s→0 s s + 0.2s + 1 s + 10
 
1 1
= lim 2
s→0 s + 0.2s + 1 s + 10
1
=
10
so we know that the final value of this system’s step response is
1/10. Note that this is much easier than finding the system’s step
response and letting t → ∞, which is much more work!

Examples where FVT cannot be applied

Here are some examples where the final-value theorem cannot be
applied.
1. The impulse response (U (s) = 1) of
1
G(s) = .
s2 + 4
Applying the FVT seems to give limt→∞ y(t) = 0, so everything
seems to be fine. But the poles of G(s) are at ±j2, so the FVT
cannot be applied. (Try plotting this in Matlab and you will
see that it does not settle down to any final value.)
7.3. Initial-value & final-value theorems 93

2.
1
G(s) = with u(t) = sin ωt
s+1
The Laplace transform of the input is L{sin ωt} = s21+1 . Then
Y (s) has ‘poles’ (due to U (s)) at ±1, so the FVT cannot be
applied. (Again, try plotting in Matlab.)

7.3.2 Initial value theorem (IVT)

The initial value theorem (IVT) says that

y(0+ ) = lim sY (s) (7.6)

s→∞

In words this says that the behaviour in the time domain for t
very small is related to the behaviour in the Laplace domain for
s very large. Unlike the FVT, we can apply the IVT (7.6) even
when there are imaginary-axis or right half-plane poles.

Example
The impulse response of
1
G(s) = .
s2 + 0.2s + 1
Applying the IVT (7.6) we get
s
y(0+ ) = lim sY (s) = lim = 0.
s→∞ s→∞ s2 + 0.2s + 1
So we know that the system’s impulse response starts at zero.
This is much less work than finding the system’s impulse response
and letting t → 0!
7.4. Non-minimum-phase transfer functions 94

7.4 Non-minimum-phase transfer functions

We will now look at two examples of what are known as ‘non-
minimum-phase’ transfer functions. A full explanation of why
they are called this is a little complicated, so we will make do
by saying that ‘non-minimum-phase’ means that there is more
phase lag than there could be. This in turn means that these
systems tend to have some kind of delay in their dynamics and
this is important for both modelling and especially for control.
We will now look at two examples of non-minimum phase transfer
functions: i) a time delay and ii) a system with a right half-plane
zero.

7.4.1 Time delays

Imagine driving a car where the information coming to you through
the windscreen is delayed by one second. Or imagine that there
is a one-second delay between pressing the car’s brake pedal and
the car slowing down. This would be a difficult car to drive! Or
imagine that the shower is too cold, so you nudge the dial up a
little hotter, and then wait a bit to see if it did the trick. This
can take a while to get right. These are all examples of systems
with a time delay.
Now imagine the simplest possible case in which the output y(t)
is simply the input u(t) at an earlier time:
y(t) = u(t − τ ) (7.7)
where τ is the time delay. What is the Laplace transform of (7.7)?
We can use the time-shift property of the Laplace transform (table
4.2), giving
Y (s) = e−sτ U (s). (7.8)
Therefore the transfer function of a time delay is
G(s) = e−sτ
7.4. Non-minimum-phase transfer functions 95

7.4.2 Right half-plane zeros

Like poles, zeros can occur in the left half-plane or the right half-
plane. Consider the following transfer function:
1−s
G(s) = , (7.9)
(s + 2)(s + 3)

which has two LHP poles and one RHP zero. The system’s step
response is plotted in figure 7.1. Notice that the step response
is initially in the opposite direction to its final value. We can
understand why this happens by using the initial- and final–value
theorems. Using the FVT, the system’s final response is
1 1
lim y(t) = lim sY (s) = lim sG(s) = . (7.10)
t→∞ s→0 s→0 s 6
We would also like to know the initial value of the system’s deriva-
tive, ẏ(t) which, using the IVT (and L{ẏ(t)} = sY (s)) is

lim ẏ(t) = lim s · sY (s) = lim sG(s) = −1. (7.11)

t→0 s→∞ s→∞

The initial gradient of y(t) is therefore in the opposite direction

to its final value. So in the same way that a time delay makes us
wait an amount of time τ before anything happens, a right half-
plane zero makes us wait a certain amount of time before things
are going in the direction we want. (Since whichever direction we
want the system to go, it will go in the direction opposite to that
to start with.)
An example of a system with a right half-plane zero is a bicycle.
We will go through a simple example in lectures, but here I want to
include a quote from Wilbur Wright of the Wright brothers:
“I have asked dozens of bicycle riders how they turn
to the left. I have never found a single person who
stated all the facts correctly when first asked. They
7.4. Non-minimum-phase transfer functions 96

0.2

0.1
y(t)

-0.1
0 0.5 1 1.5 2 2.5 3 3.5 4
t

%% stepRHPzero.m: plot step response of a system with a RHP zero

clear all, close all
% set up the transfer function
s = tf('s');
G = (1-s) / ((s+2)*(s+3));
% calculate its step response
t = [0:0.01:4]; % time
y = step(G, t); % output
% plot
figure, plot(t, y)
set(gca,'Box','On','XMinorTick','On','YMinorTick','On');
xlabel('t'), ylabel('y(t)')

Figure 7.1: Step response of the transfer function (7.9), which

has a right half-plane zero.
almost invariably said that to turn to the left, they
turned the handlebar to the left and as a result made
a turn to the left. But on further questioning them,
some would agree that they first turned the handlebar
a little to the right, and then as the machine inclined
to the left they turned the handlebar to the left, and
as a result made the circle inclining inwardly.”
Since the bike’s initial response is in the wrong direction, the best
thing we can do to compensate, initially, is to steer the bike in
the opposite direction to where we want to go!
7.5. Block diagrams 97

7.5 Block diagrams

Block diagrams can be used to represent systems that have both
inputs and outputs. They are a conceptually simple way to
1. break a complicated system into smaller, simpler systems; or
2. combine systems together to form an overall system.
We will first look at some basic blocks; and then see how blocks
combine when they are connected together in series, in parallel,
and in feedback.

7.5.1 Basic blocks

We start with the three simplest blocks representing multiplica-
tion, differentiation, and integration, all shown below.
1. multiplier or gain:

L {ku(t)} = kU (s)

2. differentiator:

L {u̇(t)} = sU (s)

3. integrator:
Z 
1
L u(t) dt = U (s)
s

7.5.2 Connecting blocks

Imagine a DC motor driving a mechanical load of some kind. If we
have a transfer function for the DC motor and a transfer function
for the mechanical load, then block diagrams provide a way to
connect them together. Three common connections that we might
7.6. Links and further reading 98

want to perform are series connection; parallel connection; and

feedback connection. We will derive each of the three cases in
lectures.

1. Series connection

Example: DC motor and load.

2. Parallel connection

Example: combined effect of two loudspeakers at a certain point

in a room.

3. Feedback connection

Example: feedback control loop for an (open-loop unstable) seg-

way.

7.6 Links and further reading

• Why bicycles do not fall: Arend Schwab at TEDxDelft
• Bicycle dynamics at TU Delft Bicycles are complicated!
8 Fluid and thermal systems

In this chapter we model fluid and thermal systems.

8.1 Fluid systems

We will now consider fluid systems for which we are interested in
the level of liquid (water, say) in some part of the system.

Some examples of applications of such systems are

• water level regulation in boilers
• irrigation lakes
• cooling towers
• waste water tanks.

99
8.1. Fluid systems 100

8.1.1 Single-tank system

First consider a system for which there is a single tank. This is
shown in figure 8.1.

Figure 8.1: Single-tank system.

We assume that the fluid is incompressible. The parameters of

the system are listed below:
• m(t): mass of fluid in the tank;
• qi (t), qo (t): volume flow rate in, out
• h(t): liquid level in the tank
• A: cross-sectional area of the tank
• pa : ambient pressure (of the surrounding air)
• p: pressure at the bottom of the tank
• ρ: density of the fluid
• R: resistance of the fluid valve
8.1. Fluid systems 101

To find the governing equation for this system we need to apply

three principles:
1. Conservation of mass. The rate of change of mass in the
tank must be equal to the mass flow rate in minus the mass
flow rate out:
ṁ(t) = ṁi (t) − ṁo (t). (8.1)
The rate of change of mass is ṁ(t) = ρq(t) so (8.1) becomes

ρq(t) = ρqi (t) − ρqo (t)

or, since the flow is incompressible (ρ is constant),

q(t) = qi (t) − qo (t).

We also know that the volume flow rate q(t) is related to the
height of water in the tank:
d dh
q(t) = (Ah(t)) = A .
dt dt

2. Pressure change across the tank. The pressure due to a

column of water of height h is
weight ρgAh
∆p = p − pa = = = ρgh.
area A
Notice that the important thing is the pressure difference,
∆p = p − pa .
3. Fluid flow rate across the exit valve. For this we define
a resistance of the valve as the ratio of the pressure difference
across the valve and the mass flow-rate through it:
∆p p − pa
ṁo = ρqo (t) = = .
R R
8.1. Fluid systems 102

You might wonder how we would find this resistance R, to

which there are two things to say. First, we didn’t worry
about how the resistance of a circuit or the damping coeffi-
cient for a mass-spring-damper were found. We just assumed
that they could be and that we had access to the value. (This
whole course is really about how to model systems of elements
combined together, rather than the detailed modelling of the
individual elements themselves.) Second, if you do want to
know more on how R can be found for fluid systems, see the
additional material for this chapter.
Now we combine these three equation to get a single equation in
terms of two things: the pressure difference ∆p(t) and the volume
flow-rate into the tank qi (t). We get

q(t) = qi (t) − qo (t)

dh ∆p
A = qi (t) −
dt ρR
A d∆p ∆p
= qi (t) − . (8.2)
ρg dt ρR

Taking Laplace transforms of (8.2) and rearranging we get

∆P (s) ρg/A
= g (8.3)
Qi (s) s + RA

This gives the transfer function between the the volume flow-
rate in Qi (s) (the input) and the pressure difference ∆P (s) (the
output). If we are interested instead in the height of the water,
H(s), then we can use ∆p(t) = ρgh(t) to give

H(s) 1/A
= g (8.4)
Qi (s) s + RA
8.1. Fluid systems 103

8.1.2 Two-tank system

What if there is more than one tank? A two-tank system is shown
in figure 8.2. Notice now that the volume flow-rate in, qi , comes
into tank 1, and that the volume flow-rate out, qo , comes out of
tank 2. We still assume that the fluid is incompressible. We want
to find the dynamic equations for h1 (t) and h2 (t) in terms of the
volume flow-rate into tank 1, qi (t). To do so we go through the
same process as for the single-tank system, but this time there
are two tanks to consider.

Figure 8.2: Two-tank system.

Tank 1
1. Conservation of mass.

ṁ = ṁi − ṁo
d
(ρA1 h1 ) = ρA1 ḣ1 = ρqi − ρq1
dt
8.1. Fluid systems 104

2. Pressure difference between the tanks.

p1 = pa + ρgh1
p2 = pa + ρgh2
=⇒ p21 = p2 − p1 = ρg(h2 − h1 )

3. Fluid flow rate across the exit valve.

p1 − p2
ρq1 =
R1
ρgh1 − ρgh2
=
R1

Now we combine all equations and eliminate internal rates:

ρgh1 − ρgh2
ρA1 ḣ1 = ρqi −
R1
gh1 − gh2
A1 ḣ1 = qi −
R1
or
g g 1
ḣ1 (t) + h1 (t) − h2 (t) = qi (t) (8.5)
A1 R1 A1 R1 A1

Tank 2
Now we do the same thing for tank 2.
1. Conservation of mass.
ṁ = ṁi − ṁo
ρA2 ḣ2 = ρq1 − ρqo

2. Pressure change across the tank.

p2 = pa + ρgh2
p3 = p a
=⇒ ∆p = p2 − p3 = ρgh2
8.1. Fluid systems 105

3. Fluid flow rate across the exit valve.

∆p ρgh2
ρqo = =
R2 R2

Now we combine all dynamic equations and eliminate internal

rates:
ρgh1 − ρgh2 gh2
ρA2 ḣ2 = −
R1 R2
gh1 gh2 gh2
A2 ḣ2 = − −
R1 R1 R2
or  
g g g
ḣ2 − h1 + + h2 = 0 (8.6)
A2 R1 A2 R1 A2 R2

Combining them together

Representing the two equations (8.5) and (8.6) for the two tanks
in a single matrix equation we get
      
g g
h1  − A1 R1
d    A1 R1
 h1  1/A1 
 =  +  qi .
   
dt    g 
g g

   
h2 A2 R1
− A2 R1
+ A2 R2
h2 0

Notice that this is in the form of a state-space model:

   

h1 
d   h1 
  = A   + Bqi , (8.7)
 
dt    
h2 h2

where the states are the two heights h1 , h2 and the input is qi .
We want to find the two transfer functions H 1 (s)
Qi (s)
and H 2 (s)
Qi (s)
. To do
8.1. Fluid systems 106

this we can take Laplace transforms of the matrix equation (8.7)

(assuming zero initial conditions):
   
H1 (s) H1 (s)
s = A  + BQi (s)
   
 
   
H2 (s) H2 (s)
 
H1 (s)
(sI − A)   = BQi (s)
 
 
H2 (s)
 
H1 (s)
 = (sI − A)−1 BQi (s).
 

 
H2 (s)

Inverting the 2 × 2 matrix we get

   −1  
H1 (s) s + A gR − A1gR1  1/A1 
1 1
=  Qi (s)
     
   
     
H2 (s) − A2gR1 s + A2gR1 + g
A2 R2
0

or    
H 1 (s)
1 A2 R1 R2 s + g(R1 + R2 )
 
 Qi (s) 
=
 H (s)  D(s) 
  

2
Qi (s)
R2 g

with
1
D(s) = .
A1 A2 R1 R2 s2 + g(A2 R2 + A1 (R1 + R2 ))s + g 2
8.2. Thermal systems 107

There are two transfer functions and they share the same de-
nominator D(s), which is second-order. The transfer function for
H1 (s) also has a single zero while the transfer function for H2 (s)
does not have any zeros. We just did something that we have not
done before: we took Laplace transforms of a state-space model.
We will look at this again in more detail in chapter 13 when we
look at the connections between state-space models and transfer
functions.

8.2 Thermal systems

Now we turn to thermal systems. We have already considered a
thermal system in chapter 6—a cup of tea. There we stated as
part of our analysis that the cup of tea had a characteristic time
constant, τ . (See for example eq. (6.2).) But we did not perform
any further modelling to explain where that time constant might
come from. Here we will consider thermal systems in a little more
detail, although the models are still quite simple. Some possible
application of thermal systems are
• heaters;
• refrigerators;
• engines;
• computer-chip design (see figure 8.3);
• communication satellites.
The important parameters for a thermal system are
• Ti , To : temperature inside, outside the system of interest;
• qi , qo : heat flow rate in, out of the system;
• ct : thermal capacity of the system;
8.2. Thermal systems 108

Figure 8.3: Simple thermal model of a computer chip.

• Rt : thermal resistance of the system.
The key assumptions that we make in the analysis are that i) the
system can be considered as a lumped model so that we model the
temperature variation in time but not in space; and ii) that no
work is performed on the system so that any temperature increase
is due to the addition of heat only. To model a thermal system
we apply two principles:
1. Conservation of energy. The rate of change of internal en-
ergy U of the system is given by the heat flux into the system
qi minus the heat flux out of the system qo :
U̇ (t) = qi (t) − qo (t).
We can also relate the internal energy U to the temperature
T:
U (t) = mu(t) = (ρV )cT (t) = ct T (t)
where u is the specific energy (J kg−1 ); c is the specific heat
capacity (J kg−1 K −1 ); and ct is the thermal capacity (J K −1 ).
2. Heat flux across the boundary. We define the thermal
resistance as the ratio of the temperature difference and the
heat flux out of the system:
∆T Ti − To
qo = = (8.8)
Rt Rt
8.2. Thermal systems 109

Example
Find the dynamic equations for ∆T (t) in terms of the heat input
qi (t) for the thermal system shown below.

• Ti , To : temperature inside, outside;

• qi , qo : heat flux in, out;
• ct : thermal capacity;
• Rt : thermal resistance;

1. Conservation of energy. The internal energy U (t) is related

to the heat fluxes in and out by

U̇ (t) = qi (t) − qo (t) (8.9)

and U (t) is related to Ti (t) by

U̇ (t) = ct Ṫi (t) = ct ∆T (t) (8.10)

(since To does not vary in time).

2. Energy flux across the boundary. We can use the defini-
tion of the thermal resistance (8.8) to relate the energy flux
out of the system to the temperature difference across it:

∆T (t) Ti (t) − To (t)

qo (t) = = . (8.11)
Rt Rt

Now we combine the three dynamical equations (8.9–8.11) and

eliminate internal rates to give

∆T (t)
ct ∆Ṫ (t) = qi (t) − .
Rt
8.3. Links and further reading 110

Taking Laplace transforms and rearranging we find that the trans-

fer function is

∆T (s) Rt
=
Qi (s) ct Rt s + 1

which has a single pole at p = − ct1Rt .

8.3 Links and further reading

• Rubicon Water—irrigation canal automation
• Rubicon Water—ABC Landline (2007)
9 Bode plots

This chapter is all about Bode plots: what they are and how to
plot them.

9.1 Introduction
Bode plots characterize the response of linear dynamical systems
to harmonic (i.e. sinusoidal) inputs. To give some simple exam-
ples, we might want to characterize the response of a mass-spring-
damper system to harmonic forcing of different frequencies; or the
response of an RLC circuit to a harmonic voltage of different fre-
quencies. To give a more complicated example, we might want
to characterize the response of a room (such as a concert hall) to
sound sources of different harmonic frequencies. For example if a
violin or a guitar plays here, how will it sound over there, and how
does this vary with frequency? Or we might want to characterize
the response of an aircraft to gusts of different frequencies. What
is the frequency for which a passenger at the back of the plane
is most likely to spill their coffee (a small gust) or be injured (a
large gust)? Can we make these less likely?
Why is it important to characterize a system’s response to har-
monic forcing? There are several reasons:
• The inputs of interest are often harmonic (or composed of a
small number of harmonic frequencies).

111
9.1. Introduction 112

• Non-harmonic inputs can be decomposed into their Fourier se-

ries (coming up in chapter 11)—and it often makes a lot of
sense to do so.
• Characterizing the response to harmonic inputs can provide
good insights into the system dynamics.
• It is essential for some closed-loop control techniques (more on
this in ELEN90055 Control Systems).
So how might we characterize the response to harmonic inputs?
One option is to set the input to
u(t) = sin ωt
and then use Laplace transform techniques like we have done in
earlier chapters.1 Let’s try this approach for a specific exam-
ple.

Example
We choose the example of a first-order system that is forced at
the frequency ω = π rad/s:
1
G(s) = and u(t) = sin πt.
s+1
First we take the Laplace transform of u(t) (see table 4.1):
π
U (s) = L{u(t)} = .
s2 + π2
Then the Laplace transform of the output is
π
Y (s) = G(s) U (s) =
(s + 1)(s2 + π 2 )
A Bπ Cs
= + 2 2
+ 2 .
s+1 s +π s + π2
1
u(t) = cos ωt would work fine too.
9.1. Introduction 113

Now take inverse Laplace transforms to get

y(t) = Ae−t + B sin πt + C cos πt

with the three constants given by (derivation not shown)

π 1 −π
A= , B= , C= .
π2 +1 π2 +1 π2 +1
How does this output y(t) relate to the sinusoidal input u(t) we
put in? For this we can use

B sin θ + C cos θ = M sin(θ + φ)

√
with M = B 2 + C 2 and φ = tan−1 B C


. We find that
s 2  2
1 −π
M= +
π2 + 1 π2 + 1
s r
1 + π2 1
= 2 2
= 2
(π + 1) π +1

and
−π/(π 2 + 1)
 
−1
φ = tan
1/(π 2 + 1)
= tan−1 (−π)
= −1.263 rad = −72.3 deg

and so finally, after quite some work, we have

r
π −t 1
y(t) = 2+1
e + 2+1
sin (πt − 1.263) (9.1)
π
| {z } π
| {z }
transient response steady-state response

The input u(t) (black) and the output y(t) (blue) are shown be-
low.
9.2. The frequency response 114

0.5
u(t), y(t)

-0.5

-1
0 2 4 6 8 10
t
That was quite a lot of work! But don’t worry, we will find a better
way. What did we learn? The key insight from (9.1) is that the
response to a sinusoidal input is made up of two components: a
transient response and a steady-state response.
The transient response decays with time due to the e−t term.
(Notice that the system pole at p = −1 appears in the transient
response via this e−t term.)
The steady-state response does not decay in time. It has the same
frequency as the input (ω = π) but it has a different amplitude
and is phase-shifted.
Now for a better way to get to the same answer: the frequency
response.

9.2 The frequency response

The purpose of the frequency response is to characterize the re-
sponse of a linear system to harmonic inputs in terms of the gain
and phase of the response.
9.2. The frequency response 115

9.2.1 Using a complex input signal

We will now look at how a linear system responds to a harmonic
input. To do so we will use an input in time of

u(t) = ejωt .

We do this because it is the easiest way to see what happens to

the gain and to the phase (which we will define soon). But we
expect time-domain signals to be real. So what does it mean to
have a complex input? Remember the identity

ejωt = cos ωt + j sin ωt.

So by using an input of u(t) = ejωt , we put in simultaneously

to the system sin ωt and cos ωt. What’s more, we know how the
system responds to each individually because the response to the
real part will be real; and the response to the imaginary part will
be imaginary. (This happens because all of the coefficients in the
ODE are real.) In other words, we can keep track of the two
parts because the imaginary part of the response will be ‘tagged’
with a ‘j’, like money stolen from the bank that is stained with
dye.

9.2.2 Response to harmonic forcing

With the idea of a complex input not entirely weird, we will state
the main result. If we apply an input u(t) = eiωt to a system with
transfer function G(s) then its output y(t) will be
n
X Bi
y(t) = epi t + G(jω) ejωt . (9.2)
i=1
pi − jω
| {z } | {z }
transient response steady-state response

(I’ll spare you the details of the derivation but it is provided in

§ 9.A if you are interested.) This looks a little bit complicated so
9.2. The frequency response 116

let’s break it down. The response is composed of two parts: a

transient response and a steady-state response.
The transient response appears the more complicated of the two.
But provided that the system is stable—i.e. provided that all of
its poles satisfy Re(pi ) < 0—then it will decay to zero and so it is
not our main interest. Once the transient response has decayed
then we are left with the bit in which we are more interested: the
steady-state response,

y(t) = G(jω) ejωt (9.3)

or, since u(t) = ejωt ,

y(t) = G(jω) u(t) (9.4)

G(jω) is called the frequency response. It is complex and

so, like any complex number, we can represent it in terms of a
magnitude (gain) and an angle (phase):

G(jω) = |G(jω)| ej∠G(jω) . (9.5)

Therefore

y(t) = G(jω) ejωt = |G(jω)| ej(ωt+∠G(jω)) . (9.6)

Remember that the input we applied was u(t) = ejωt and so (9.6)
tells us two important things:
• that relative to the input u(t), the output y(t) is scaled by a
factor |G(jω)|. This is the gain of G at frequency ω.
• that relative to the input u(t), the output is phase-shifted by
an amount ∠G(jω). This is the phase of G at frequency ω.
So now we know what the frequency response is, and we know
how to define its gain and phase. We are now ready to tackle
Bode plots.
9.3. The Bode plot 117

9.3 The Bode plot

The Bode plot is really two plots:
1. A gain plot in which we plot 20 log10 |G(jω)|
We call this unit of measurement the decibel. For an expla-
nation and some history behind its use, see § 9.B.
2. A phase plot in which we plot ∠G(jω)

On the x-axis for both plots we plot log10 ω

So this is what we plot. But what does the Bode plot tell us?
The key point is that it tells us how the system responds to a
harmonic input in steady-state (i.e. once any transients have died
away).
How do we actually plot the Bode plot? That is, how do we find
the gain and phase of some transfer function G(s)? Let’s look at
a specific example of a first-order system:
1
G(s) = . (9.7)
s+a
This is the transfer function, G(s), but (9.1) tells us that for the
frequency response we need G(jω). We can obtain the frequency
response from the transfer function by setting

s = jω.

Let’s try that for the first order system—we just need to write jω
wherever we see an s. This gives
1
G(jω) = .
jω + a
It is clear that G(jω) is complex. Our task now is to find the
gain, |G(jω)|, and the phase, ∠G(jω).
9.4. Bode plots by hand 118

The gain |G(jω)| is

|1|
|G(jω)| =
|jω + a|
1
=√ .
ω 2 + a2
The phase ∠G(jω) is

∠G(jω) = ∠(1) − ∠(jω + a)

= −∠(jω + a)
ω 
= − arctan .
a

9.4 Bode plots by hand

You can probably imagine that Matlab would have a command
to plot Bode plots. So why would we bother plotting them by
hand? There are a couple of reasons.
1. It is no good plotting something if you have no idea what it
means. If plotting is just a matter of typing stuff and pressing
enter, then engineers might be replaced by monkeys sooner or
later.
2. Often we need to design a transfer function to meet certain
requirements. This is especially true in control. In order to
know what to modify, we need to understand how changes in
a transfer function will translate to changes in the Bode plot.
Matlab cannot tell us this—we need to use our brains!
When we plot a Bode plot by hand we only need to get it roughly
right. There are a few common strategies we will use for all trans-
fer functions to help in plotting their Bode plot by hand:
1. Consider what happens at low frequencies.
9.4. Bode plots by hand 119

2. Consider what happens at high frequencies.

3. Consider what happens around the frequency of any poles and
zeros.
4. Consider the gradient of the gain in different frequency ranges.
Let’s look again at the first-order system.

9.4.1 At low frequencies

When we say low: low compared to what? For the first-order
system (9.7) the only other frequency to which we can compare
is the frequency of the pole at s = −a. Low for this transfer
function therefore means ω  |a|.2 Then the gain becomes
1
|G(jω)| =
|a|

and the phase becomes

∠G(jω) = − arctan(0) = 0.

(Provided that a > 0. If a < 0 then ∠G(jω) = −180◦ .) So

we have found what happens to the gain and the phase at low
frequencies.

9.4.2 At high frequencies

Similarly, when we say high: high compared to what? ‘High’ for
the first-order transfer function (9.7) means ω  a. Then the
gain becomes
1
|G(jω)| ≈
ω
2
We only plot the Bode plot for positive frequencies, hence it is the mag-
nitude of a that matters.
9.4. Bode plots by hand 120

and the phase becomes

∠G(jω) ≈ −90◦ .

(Again, provided that a > 0. If a < 0 then ∠G(jω) = −270◦ .)

Notice that the phase settles to a constant final value but that
the gain does not. Instead the gain becomes progressively smaller
as ω increases.

9.4.3 Around the frequencies of any poles or zeros

There is a single pole at s = −a. Setting ω = |a| the gain be-
comes3
1 1
|G(jω)| = √ =√
a2 + a2 2|a|
and the phase is
a
∠G(jω) = − arctan = −45◦ .
a

One more thing: the straight-line approximations that we will go

through in lectures imply a gain of |G(jω)|
√ = 1/|a| at ω = |a|. But
we just saw that the true gain is 1/ 2|a|. Therefore
√ the straight-
line approximation
√ is wrong by a factor of 1/ 2. In decibels this
is 20 log10 (1/ 2) ≈ -3 dB. (This will make more sense when we
draw it!)

9.4.4 Slope of the gain

At low frequencies the gain is simply 1/|a|, which is independent
of ω. In other words, the gain is constant at low frequencies
and so the slope of the gain is zero or, using the right units, 0
dB/decade.
3
Again, we only plot the Bode plot for positive frequencies, hence it is the
magnitude of a that matters.
9.5. Second-order system example 121

At high frequencies the gain is 1/ω. The slope of this on a log-log

scale is -1. But remember that we plot 20 log |G(jω)| and so the
slope is -20 dB/decade.
1
The Bode plot of the first-order system G(s) = with a = 1
s+a
is shown in figure 9.4.

9.5 Second-order system example

Let’s try another example: a second-order system
1
G(s) =
s2 + 2ζωn s + ωn2

with a damping ratio 0 < ζ < 1. Setting s = jω, the frequency

response is
1
G(jω) = 2 . (9.8)
(ωn − ω 2 ) + j2ζωn ω
Let’s again consider what happens at low frequencies; at high
frequencies; at ω = ωn ; and to the slope of the gain in different
frequency ranges.

9.5.1 Low frequencies

For ω  ωn we have
1
G(jω) ≈
ωn2
since the ωn2 term in (9.8) dominates everything else. The gain
and the phase are therefore approximately
1
|G(jω)| ≈ and ∠G(jω) ≈ 0◦ .
ωn2
9.5. Second-order system example 122

9.5.2 High frequencies

For ω  ωn we have
1
G(jω) ≈
−ω 2
since now the ω 2 term in (9.8) dominates everything else. The
gain and the phase are therefore approximately
1
|G(jω)| ≈ and ∠G(jω) ≈ −180◦ .
ω2
Notice that the phase settles to a constant final value but that
the gain does not. Instead the gain becomes progressively smaller
as ω increases.

9.5.3 At ω = ωn
For ω = ωn we have
1 −j
G(jω) = 2
= .
j2ζωn 2ζωn2
The gain and the phase are therefore
1
|G(jωn )| = and ∠G(jωn ) = −90◦ . (9.9)
2ζωn2

9.5.4 Slope of the gain

At low frequencies the gain is 1/ωn2 , which is independent of ω.
The gain is therefore constant at low frequencies and its slope is
zero or, using the right units, 0 dB/decade.
At high frequencies the gain is 1/ω 2 . The slope of this on a log-log
scale is -2. But remember that we plot 20 log |G(jω)| and so the
slope is -40 dB/decade.
1
The Bode plot of the second-order system G(s) =
s2 + 2ζωn s + ωn2
is shown in figure 9.5.
9.6. Important properties of Bode plots 123

9.5.5 Resonant peak

In figure 9.5 we see that there is a clear peak in the gain plot
when ζ is small. But at what value of ω does this occur? The
expression for the gain for some general value of ω is
1
|G(jω)| = p .
(ωn − ω ) + (2ζωn ω)2
2 2

Finding the peak in the gain plot means maximizing this quantity.
Since the numerator is just 1, this is the same as minimizing the
quantity
[(ωn2 − ω 2 ) + (2ζωn ω)2 ]1/2 . (9.10)
To find the minimum value of (9.10) we differentiate it with re-
spect to ω and set it to zero. This gives a frequency of
p
ωpeak = ωn 1 − 2ζ 2 . (9.11)
This is neither the undamped
p natural frequency ωn nor the damped
natural frequency ωd = ωn 1 − ζ 2 ! But for small ζ it will be close
to both. Substituting (9.11) into the expression for the gain we
find a maximum gain of
1
Gmax = |G(jωpeak )| = p .
2ζωn2 1 − ζ 2
For small ζ this is very close to the value of the gain at ω = ωn
that we found in (9.9).
One final remark: from (9.11) we see that there
√ is no resonant
peak if the damping ratio is greater than 1/ 2 = 0.707 (since
then ωpeak becomes negative, which is not meaningful).

9.6 Important properties of Bode plots

There are two important properties of Bode plots that will be
helpful when we come to plot them: i) the additive property for
9.6. Important properties of Bode plots 124

the gain; and ii) the additive property for the phase, both of which
we will now explain.

9.6.1 Additive property for gain

What if we multiply together two transfer functions and plot the
gain of their product? Let

G(s) = G1 (s)G2 (s).

Then

20 log10 |G(s)| = 20 log10 |G1 (s)G2 (s)|

= 20 log10 |G1 (s)| + 20 log10 |G2 (s)|.

So to plot the gain of G1 (s)G2 (s) (in dB), we can simply add
together the gain of G1 (s) (in dB) and the gain of G2 (s) (in dB).
This relies critically on the fact that we plot the gain using a
logarithmic scale.

9.6.2 Additive property for phase

Similarly, what if we multiply together two transfer functions and
plot the phase of their product?

∠ G(s) = ∠(G1 (s)G2 (s))

= ∠ G1 (s) + ∠ G2 (s).

So to plot the phase of G1 (s)G2 (s), we can simply add together

the phase of G1 (s) and the phase of G2 (s). In other words, the
phase follows the same rule as the gain. Note that this works
because if
G1 (s) = r1 ejθ1 , G2 (s) = r2 ejθ2
then
∠G(s) = ∠r1 r2 ej(θ1 +θ2 ) = θ1 + θ2 .
9.7. General procedure for plotting Bode plots 125

9.7 General procedure for plotting Bode plots

Here is a general procedure for plotting Bode plots:
1. Factor the transfer function G(s) into the form
K(s + b1 ) · · · (s + bm )
.
sq (s + a1 ) · · · (s + an )(s2 + 2ζ1 ω1 s + ω12 ) · · · (s2 + 2ζr ωr s + ωr2 )
(Note that ai = −pi and bi = −zi .)
2. Define the ω-axis so that it starts two decades before the lowest
important frequency; and ends two decades after the highest
important frequency, i.e.
ωmin = 0.01 min(pi , zi , ωi )
ωmax = 100 max(pi , zi , ωi )

3. Mark all break points on the ω-axis: |p1 |, . . . , |pn |, ω1 , . . . , ωr

and |z1 |, . . . , |zm |.

Gain plot
1. If there are no zeros or poles at the origin, start at 20 log10 |G(j0)| dB
with a horizontal asymptote.
2. If the term 1/sq is present, start at 20 log10 |G(j0)sq |−20q log10 ωmin dB
with a −20q dB/dec asymptote.
3. At |zi |, add to the slope +20 dB/dec. At |pi |, add to the slope
-20 dB/dec. At ωi , add to the slope
√ -40 dB/dec. Around ωi ,
mark on a resonant peak if ζi < 1/ 2 ≈ 0.7.

Phase plot
1. No poles/zeros at the origin: if K > 0 start the graph at 0◦
with a horizontal asymptote; if K < 0 start the graph at −180◦
with a horizontal asymptote.
9.8. Bode plots of common transfer functions 126

2. If the term 1/sq is present: if K > 0 start the graph at −90q ◦ ;

if K < 0 start at (−180 − 90q)◦ with a horizontal asymptote.
3. At |zi |, add to slope +45◦ /dec from 0.1|zi | to 10|zi |. At |pi |, add
to slope −45◦ /dec from 0.1|pi | to 10|pi |. (The phase around ωi
depends strongly on the corresponding damping ratio ζi and so
no general statements can be made, other than that the phase
will change by −180◦ overall.)

9.8 Bode plots of common transfer functions

Here we collect together the Bode plots of some common transfer
functions. They have all been plotted in Matlab using the bode()
command. Since we have plotted them in Matlab we of course
had to choose specific parameters—for example, G(s) = 1/(s + 1)
rather than G(s) = a/(s + a). You should be able to relate the
parameters chosen to important points on the plots.
9.8. Bode plots of common transfer functions 127

20

15
Magnitude (dB)

10

-5
91
Phase (deg)

90.5

90

89.5

89
10 0 10 1
Frequency (rad/s)

Figure 9.1: G(s) = s: differentiator.

0
Magnitude (dB)

-5

-10

-15

-20
-89
Phase (deg)

-89.5

-90

-90.5

-91
10 0 10 1
Frequency (rad/s)

Figure 9.2: G(s) = 1s : integrator.

9.8. Bode plots of common transfer functions 128

40
Magnitude (dB)

30

20

10

0
90
Phase (deg)

45

0
10 -2 10 -1 10 0 10 1 10 2
Frequency (rad/s)

Figure 9.3: G(s) = s + 1: single zero.

0
Magnitude (dB)

-10

-20

-30

-40
0
Phase (deg)

-45

-90
10 -2 10 -1 10 0 10 1 10 2
Frequency (rad/s)

1
Figure 9.4: G(s) = s+1
: single pole.
9.8. Bode plots of common transfer functions 129

40

20
Magnitude (dB)

-20

-40

-60

0
Phase (deg)

-45

-90

-135

-180
10 -2 10 -1 10 0 10 1 10 2
Frequency (rad/s)

Figure 9.5: G(s) = s2 +2ζω1n s+ω2 for ωn = 1 and ζ ranging

n
between ζ = 0.01 and ζ = 1.

0
Magnitude (dB)

-20

-40

-60
0
Phase (deg)

-45

-90

-135

-180
10 -1 10 0 10 1 10 2
Frequency (rad/s)

32
Figure 9.6: G(s) = (s+3)2
: pair of real poles.
9.8. Bode plots of common transfer functions 130

0
Magnitude (dB)

-5

-10

-15

-20
60
Phase (deg)

30

0
10 -2 10 -1 10 0 10 1 10 2 10 3
Frequency (rad/s)

s+1
Figure 9.7: G(s) = s+10
: single zero followed by a single pole.

10

8
Magnitude (dB)

0
0
Phase (deg)

-10

-20

-30
10 -2 10 -1 10 0 10 1 10 2
Frequency (rad/s)

s+3
Figure 9.8: G(s) = s+1
: single pole followed by a single zero.
9.8. Bode plots of common transfer functions 131

40
Magnitude (dB)

20

-20
180
Phase (deg)

135

90

45

0
10 -1 10 0 10 1
Frequency (rad/s)

Figure 9.9: G(s) = s2 + 0.1s + 1: complex-conjugate pair of

zeros.

80

60
Magnitude (dB)

40

20

-20
0
Phase (deg)

-45

-90

-135

-180
10 -1 10 0 10 1 10 2
Frequency (rad/s)

2
Figure 9.10: G(s) = ss2 +0.02s+1
+s+100
: complex-conjugate pair of zeros
and complex-conjugate pair of poles.
9.A. Derivation of the frequency response 132

9.A Derivation of the frequency response

Here we derive equation (9.2) for the frequency response. (You
are not expected to be able to derive this.) As described in § 9.2
we set the input to be
u(t) = ejωt .
1
Then U (s) = L{u(t)} = s−jω
and

Y (s) = G(s)U (s)

Πm (s − zi ) 1 Πm
i=1 (s − zi )
= k ni=1 = k n+1 .
Πi=1 (s − pi ) (s − jω) Πi=1 (s − pi )

(We have absorbed the 1/(s − jω) term into the denominator so
that we have an extra pole, pn+1 = jω.)
Now write Y (s) as a sum of partial fractions:
n+1
X B̃i
Y (s) =
i=1
s − pi

Now we need expressions for the B̃i terms:

B̃i = lim (s − pi )G̃(s)

s→pi
1
B̃i = lim (s − pi )G(s)
s→pi (s − jω)
| {z }
Bi

for i = 1, . . . , n, we have:
Bi
B̃i =
pi − jω
9.B. The decibel 133

while for i = n + 1, we have:

B̃n+1 = lim (s − jω)G̃(s)
s→jω
1
= lim (s − jω)G(s)
s→jω (s − jω)
= G(jω).
The output is then given by
Y (s) = G̃(s) = G(s)U (s)
n+1
X B̃i
=
i=1
s − pi
n
X Bi /(pi − jω) G(jω)
= + .
i=1
s − pi s − jω
Taking inverse Laplace transforms, we get
n n
X Bi X
y(t) = epi t + G(jω) ejωt
i=1
pi − jω i=1

which is equation (9.2).

9.B The decibel

The decibel is based on the measurement of power originally used
in telephony in the Bell Telephone Company. One decibel is one
tenth (deci-) of one bel:
 
Pout
gain in decibels (dB) = 10 log10
Pin
 2
Vout
= 10 log10
V
 in 
Vout
= 20 log10 .
Vin
9.B. The decibel 134

So the factor of 10 comes from the fact that we are measuring

decibels rather than bels (named in honour of Alexander Graham
Bell, above); and the factor of 2 comes from the fact that the peo-
ple at Bell chose to measure the ratio of power (Pout /Pin ), which
involves the ratio of voltages squared (Vout /Vin )2 . Despite starting
its life as a unit of measurement for telephony, the decibel is still
used today in a broad range of fields such as acoustics, electronics
and in control system design. (Once something becomes standard
it is unlikely to change. Everyone uses it because everyone uses
it!)
10 Nyquist plots

This chapter is all about Nyquist plots, which also go by the name
of polar plots.

10.1 Nyquist plots

A Nyquist plot (or polar plot) shows the same information as a
Bode plot but in a different way. We saw in the previous chapter
that the frequency response is obtained by setting s = jω in the
transfer function:

G(jω) = frequency response.

This frequency response G(jω) is in general complex and remem-

ber that there are two ways to plot a complex number:

G(jω) = |G(jω)|ej∠G(jω) (10.1)

G(jω) = Re G(jω) + Im G(jω) (10.2)

In the Bode plot we opt for (10.1) and plot the frequency response
in terms of its gain |G(jω)| and phase ∠G(jω). In the Nyquist
plot we opt for (10.2) and plot the frequency response in terms of
its real and imaginary parts.
One important difference between Bode plots and Nyquist plots is
that, in the Nyquist plot, we generally plot the frequency response

135
10.2. Nyquist plots of common transfer functions 136

for both positive and negative frequencies, i.e. ω ∈ [−∞, +∞].

For this it is useful to know that

G(−jω) = [G(jω)]∗ = Re G(jω) − Im G(jω).

So once we have worked out the Nyquist plot for positive ω, we

can get the response for negative ω quite easily by simply flipping
the sign of the imaginary part. In other words, the Nyquist plot of
G(−jω) is the mirror image of G(+jω) about the real axis.
Why bother learning about Nyquist plots? There are two main
reasons:
• It is useful to be able to ‘picture’ frequency responses from both
points of view (i.e. Bode and Nyquist views).
• Nyquist plots are essential for understanding closed-loop stabil-
ity and closed-loop control performance. You will study them
further in ELEN90055 so consider this a first introduction.
The Nyquist plots of a number of transfer functions are shown in
the next few pages. We will also go through lots of examples in
lectures.
I will finish by saying that there are a couple of things that we
are not covering on Nyquist plots:
• the Nyquist D-contour
• Nyquist plots for systems with imaginary-axis poles.
These will both be covered in ELEN90055 when you need them.

10.2 Nyquist plots of common transfer func-

tions
10.2. Nyquist plots of common transfer functions 137

10

2
Im(G)

-2

-4

-6

-8

-10
-10 -8 -6 -4 -2 0 2 4 6 8 10
Re(G)

Figure 10.1: G(s) = s: differentiator.

10

2
Im(G)

-2

-4

-6

-8

-10
-10 -8 -6 -4 -2 0 2 4 6 8 10
Re(G)

Figure 10.2: G(s) = s + 1: single LHP zero.

10.2. Nyquist plots of common transfer functions 138

0.5

0.4

0.3

0.2

0.1
Im(G)

-0.1

-0.2

-0.3

-0.4

-0.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Re(G)

1
Figure 10.3: G(s) = s+1
: single LHP pole.

60

40

20
Im(G)

-20

-40

-60
-30 -20 -10 0 10 20 30
Re(G)

1
Figure 10.4: G(s) = s2 +0.02s+1
: pair of complex LHP poles.
10.2. Nyquist plots of common transfer functions 139

0.8

0.6

0.4

0.2
Im(G)

-0.2

-0.4

-0.6

-0.8
-0.2 0 0.2 0.4 0.6 0.8 1
Re(G)

32
Figure 10.5: G(s) = (s+3)2
: pair of real LHP poles.

0.5

0.4

0.3

0.2

0.1
Im(G)

-0.1

-0.2

-0.3

-0.4

-0.5
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Re(G)

s+1
Figure 10.6: G(s) = s+10
: single LHP zero followed by a single
LHP pole.
10.2. Nyquist plots of common transfer functions 140

1.5

0.5
Im(G)

-0.5

-1
1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
Re(G)

s+3
Figure 10.7: G(s) = s+1
: single LHP pole followed by a single
LHP zero.

0.8

0.6

0.4

0.2
Im(G)

-0.2

-0.4

-0.6

-0.8

-1
-100 -80 -60 -40 -20 0 20
Re(G)

Figure 10.8: G(s) = s2 + 0.1s + 1: complex-conjugate pair of

zeros.
10.2. Nyquist plots of common transfer functions 141

5000

4000

3000

2000

1000
Im(G)

-1000

-2000

-3000

-4000

-5000
-3000 -2000 -1000 0 1000 2000 3000
Re(G)

2
Figure 10.9: G(s) = ss2 +0.02s+1
+s+100
: complex-conjugate pair of zeros
and complex-conjugate pair of poles.
11 Fourier analysis

Fourier series come up in all kinds of areas of science and engi-

neering. Some examples include mp3 audio compression; image
processing; audio and image CODECs; noise cancellation such
as used in some headphones; magnetic resonance imaging (MRI);
DNA sequencing; and, last but not least, the technology behind
shazam. This list is by no means exhaustive and there are many
more examples.

11.1 Fourier series

So clearly Fourier series are very useful. But what are they?
The basic idea, first formulated by Joseph Fourier, is this: any
periodic signal, no matter how complicated, can be broken down
into smaller periodic building blocks. These building blocks are
composed of sine terms like sin ωi t and cosine terms like cos ωi t.
With this quite simple idea we can now write down the expression
for the Fourier series. The Fourier series of some periodic function
f (t) with period T is
∞   ∞  
X 2πkt X 2πkt
f (t) = a0 + ak cos + bk sin . (11.1)
k=1
T k=1
T

So what does eq. (11.1) actually say? It says that some periodic
function, f (t), can be represented as a sum of building blocks,

142
11.1. Fourier series 143

and that those building blocks are sines and cosines of different
frequencies (and a constant). The amplitudes of the cosine build-
ing blocks (i.e. how much of each we have) are given by the ak
terms; and the amplitudes of the sine building blocks are given
by the bk terms. Finally, the amplitude of the constant term is a0
(which could of course be zero).
Let’s take note of a few important things here:
1. We have specified (in the sentence before eq. (11.1)) that f (t)
should be periodic with period T .
2. There is a constant term, a0 . This accounts for the fact that,
as well as being made up of sines and cosines, f (t) may also
have a contribution that is constant with t.
3. The two sums go from i = 1 to i = ∞. We can generally
truncate the series at some finite number, though, and still
achieve good results.
4. We have chosen f to be a function of t, f (t), because we had
to choose something, but we could equally have written, say,
f (x) or f (y) instead.
Equation (11.1) tells us how to get f (t) from knowledge of ak and
bk terms. But we are usually interested in the problem the other
way around: given f (t), how do we find the ak and bk terms? We
will go through this in more detail in lectures. Here we will just
write down the final result. The constants a0 , ak and bk are given
11.2. Complex-exponential form 144

by

1 T
Z
a0 = f (t) dt (11.2a)
T 0
2 T
Z  
2πkt
ak = f (t) cos dt (11.2b)
T 0 T
2 T
Z  
2πkt
bk = f (t) sin dt. (11.2c)
T 0 T

11.2 Complex-exponential form

It turns out that equation (11.1) is not the only way to represent
a Fourier series. Before we go on, a little reminder about Euler’s
identity, which says that

ejωt = cos ωt + j sin ωt. (11.3)

The complex exponential is therefore made up of a cosine and sine

which, remember, are the building blocks of the Fourier trans-
form. So linking (11.1) to (11.3), we might wonder whether the
Fourier series could instead be represented as a sum of complex
exponentials—whether, that is, we could use ejωt as the building
block. The answer is yes, we can—the periodic function f (t) can
also be represented as
∞
X
f (t) = ck ejωk t (11.4)
k=−∞

where
ωk = 2πk/T.
In this case the coefficients ck are given by

1 T
Z
ck = f (t)e−jωk t dt. (11.5)
T 0
11.3. Relationship between them 145

11.3 Relationship between them

What is the relationship between the constants ak , bk in (11.1)
and the constants ck in (11.4)? Can we use one set of constants
to find the other set? To get ak , bk from knowledge of ck we can
use
a0 = c 0 (11.6a)
ak = ck + c−k (11.6b)
bk = j(ck − c−k ) (11.6c)
and to get ck from knowledge of ak , bk we can use
c 0 = a0 (11.7a)
1
ck = (ak − jbk ) (11.7b)
2
1
c−k = (ak + jbk ). (11.7c)
2
You can find these equations by writing out the two versions of
the Fourier series and equating each term in the series.

11.4 Example: Square wave

Let’s now look at the specific example of a square wave of period
T. The square wave is shown in figure 11.1. We will find its Fourier
series in terms of sine and cosine building blocks, i.e. (11.1).
First let’s find a0 :
1 T
Z
a0 = f (t) dt
T 0
1 T /2 1 T
Z Z
= 1 · dt + 0 · dt
T 0 T T /2
1 T
= ·
T 2
11.4. Example: Square wave 146

nf t

O t
T 72 Tiz T

Figure 11.1: Square wave.

so
1
a0 =
2

Now for ak :

2 T
Z  
2πkt
ak = f (t) cos dt
T 0 T
2 T /2 2 T
Z   Z  
2πkt 2πkt
= 1 · cos dt + 0 · cos dt
T 0 T T T /2 T
2 T /2
Z  
2πkt
= cos dt
T 0 T
  T /2
2 T 2πkt
= · sin
T 2πk T 0
1
= [sin(πk) − sin(0)]
πk
= 0.

so
ak = 0 for all k = 1, 2, 3, . . .
11.4. Example: Square wave 147

And finally bk :

2 T
Z  
2πkt
bk = f (t) sin dt
T 0 T
2 T /2 2 T
Z   Z  
2πkt 2πkt
= 1 · sin dt + 0 · sin dt
T 0 T T T /2 T
2 T /2
Z  
2πkt
= sin dt
T 0 T
  T /2
2 T 2πkt
=− · cos
T 2πk T 0
1
= (cos(0) − cos(πk))
πk
1
= (1 − (−1)k )
πk
and therefore

bk = 0 if k is even
2
bk = if k is odd.
kπ
Putting this all together, we can therefore write the Fourier series
of the square wave as
∞  
1 X 2 2πkt
f (t) = + sin (11.8)
2 k=1,3,5,... kπ T

The Fourier series (11.8) is compared to the true square wave for
different numbers of Fourier modes in figure 11.2.
It makes sense that the ak terms came out as zero: the original
square wave is an odd function (f (−t) = −f (t)). Meanwhile
cosine waves are even functions (f (t) = f (−t)) and sine waves are
odd functions. All of the cosine waves come out to zero because
11.4. Example: Square wave 148

Figure 11.2: Square wave of figure 11.1 (with T = 1): original

square wave and its Fourier series when including (from top to
bottom) 1 Fourier mode; 3 modes; 5 modes; 11 modes; 101 modes;
and 1001 modes.
11.5. Application: Pulse Width Modulation (PWM) 149

you cannot build an odd function out of even building blocks!

(Try shifting the square wave along the time axis so that it is
symmetric about t = 0 and therefore even and try repeating what
we did. You should find that the sine terms are zero this time
around.)

11.5 Application: Pulse Width Modulation (PWM)

Pulse width modulation (PWM) varies the average power pro-
vided by an electrical signal by chopping it into smaller parts.
This is achieved by varying the pulse’s duty cycle, which is the
time the signal spends high versus the time it spends low. PWM is
used in many applications such as AC motors, photovoltaic solar
battery chargers and voltage regulators.
We will now look at the application of PWM for a DC motor.
We will see that this is as an application for which both Fourier
series (for the pulse wave modulation) and Bode plots (for the DC
motor) are needed.
So, what is the idea behind PWM? Here is a step input of ampli-
tude 1:
Ault
1

What if weAult
want the input to be (say)Ault
half as large?
1

t t
X
Tiz Tiz
11.5. Application: Pulse Width Modulation (PWM) 150
t

Ault Ault
1

t t
X
Tiz Tiz

One option (on the left above) would be to reduce the size of the
step. Another option (on the right above) would be to reduce the
energy of the input by placing ‘gaps’ in it. Pulse width modulation
uses the second option (on the right above) and we will now look
at this in more detail by considering the Fourier series of the
input; and the Bode plot of a DC motor; and combining them
together.

11.5.1 Complex Fourier of the PWM signal

The PWM input signal is shown below:

Notice that the period of the cycle is T ; that the signal is high
for a length of time of 2T0 ; and that it is low for the rest of the
cycle (i.e. for an amount of time (T − 2T0 )).
Our first task is to find the Fourier series of the pulse and we
will find it in complex-exponential form (11.4) because this will
make application of the Bode plot easier. This means finding the
11.5. Application: Pulse Width Modulation (PWM) 151

coefficients ck . For k = 0 we have

1 T /2
Z
c0 = u(t)e−j·0·t dt
T −T /2
1 T0
Z
= 1 · e0 dt
T −T0
T0
=2 ,
T
which is the average value of the signal. For k 6= 0 we have
T0
1
Z
1 · e−jk( T )t dt
2π
ck = (u(t) = 0 outside this range)
T −T0
" #T0
1 1
e−jk( T )t
2π
=
T −jk 2πT −T0
j h −jk( 2π )T0 i
− e−jk( T )(−T0 )
2π
= e T
2πk         
j T0 T0 T0 T0
= cos 2πk − j sin 2πk − cos 2πk − j sin 2πk
2πk T T T T
  
j T0
= −2j sin 2πk
2πk T
 
1 T0
= sin 2πk .
πk T

Now we can use these to form the Fourier series:

∞
X
u(t) = ck ejωk t
k=−∞

1 1 1
first noting that x
sin(cx) = −x
sin(−cx) (i.e. x
sin(cx) is an even
11.5. Application: Pulse Width Modulation (PWM) 152

function)
∞   −1  
T0 j0 X 1 T0 jk 2π
X 1 T0 2π
u(t) = 2 e + sin 2πk e T
t
+ sin 2πk e−jk T t
T πk T πk T
k=1 k=−∞
∞  
T0 X 1 T0  jk 2π t 2π

=2 + sin 2πk e T + e−j(−k) T t
T πk T
k=1
∞  
T0 X 2 T0 2π
=2 + sin 2πk ejk T t
T πk T
k=1

therefore the coefficients are

 
T0 2 T0
c0 = 2 , ck = sin 2πk
T πk T

and we see that

c−k = ck (u(t) is even so c−k = ck )

ck = 0 for k = ±4, ±8, ±12, . . .

The coefficients ck are shown below:

0.2

0
-20 -10 0 10 20

and the truncated approximation is

k=K  
T0 X 2 T0 jk 2π t
uK (t) = 2 + sin 2πk e T
T k=1
πk T

which is shown below for 1, 3, 5, 20, 50 and 200 Fourier modes.

11.6. Application: PWM of a DC motor 153

1 1

0.5 0.5

0 0
-10 -5 0 5 10 -10 -5 0 5 10

1 1

0.5 0.5

0 0
-10 -5 0 5 10 -10 -5 0 5 10

1 1

0.5 0.5

0 0
-10 -5 0 5 10 -10 -5 0 5 10

11.6 Application: PWM of a DC motor

A DC motor, which we have considered before, is shown be-
low:

The question is: How do we drive the DC motor with a pulse-

width-modulated input such that its angular velocity, ω = θ̇, is
the same as what would be achieved by applying a constant step
input of amplitude α? We want

θ̇ = ω = ωd = constant

Can we achieve this with an input that is pulse-width modulated,

i.e. constantly begin switched on and off? We will see that we
can, provided that we switch it on and off quickly enough.
11.6. Application: PWM of a DC motor 154

Dynamics of the motor:

d
V (t) = Ra ia (t) + La ia (t) + Ve (t)
dt
The torque is τ (t) = Kt ia (t). The back-EMF is Ve (t) = Ke ω(t) =
Ke θ̇(t).
Dynamics of loaded gear:
J θ̈(t) = −bθ̇(t) + τ (t)
The input to the system is the voltage V (t). The output is the
angular rotation of the gear ω = θ̇.
Now combining all equations gives
Kt V (t) = JLa ω̈(t) + (JRa + La b)ω̇(t) + (Ra b + Ke Kt )ω(t)
Ω(s) Kt
=⇒ G(s) = = 2
V (s) JLa s + (JRa + La b)s + Ra b + Ke Kt
Now we choose the parameters J = 1, La = 1, b = 5, Ra = 20,
Ke = 0, Kt = 100 to give
Ω(s) 100
G(s) = = .
V (s) (s + 5)(s + 20)

First, how would we achieve ω = constant = ωd using a step in-

put? For this we can use the final value theorem (chapter 7):
lim ω(t) = lim sΩ(s) (11.9)
t→∞ s→0

for a step of amplitude α

α
lim ω(t) = lim sG(s)V (s) = lim s G(s)
t→∞ s→0 s→0 s
= lim αG(s) = ωd
s→0
Kt
=α = ωd
Ra b + Ke Kt
11.6. Application: PWM of a DC motor 155

so
Ra b + Ke Kt
α= ωd
Kt
100
=
100
= ωd .

Now for the pulse: our approximation of the pulse is

K
X
uK (t) = ck ejωk t
k=−K

with
 
T0 2 T0 2π
c0 = 2 , ck = sin 2πk , ωk = k.
T πk T T

Using the frequency response G(jω), the output will be

K
X
yK (t) = G(jωk )ck ejωk t
k=−K

and what we want is that

• the zero-frequency component of y(t) equals ωd
• the other Fourier components of y(t) are small (i.e. the contri-
butions from terms with ωk 6= 0).
The zero-frequency component is

y0 = c0 ej0 G(j0)
T0 100 T0
=2 =2
T 5 × 20 T
11.6. Application: PWM of a DC motor 156

and we want this to be equal to ωd :

T0 T0 1
2 = ωd =⇒ = ωd
T T 2
Where do the other Fourier components of u(t) lie on the Bode
plot of G(jω)?
2π
ωk = k
T
the spacing in ωk is then
2π
∆ω = .
T
If we can make this spacing large enough then the other contri-
butions to the Fourier series will be filtered out by the frequency
response of the DC motor so that all that is left—to a reasonable
approximation—is the zero-frequency response

yK (t) ≈ G(j0)c0 = ωd .

Now let’s at what happens for a few different values of T (and

therefore ∆ω).
11.6. Application: PWM of a DC motor 157

Case 1: T = 4 =⇒ ω1 = 2π/T = 1.57 Hz

0.5

0
-10 -5 0 5 10

0.5

0
-10 -5 0 5 10

-100

-200
10 -2 10 0 10 2 10 4 10 6

The spacing between frequencies ∆ω is too large and higher fre-

quencies are still very noticeable in the output.
11.6. Application: PWM of a DC motor 158

Case 2: T = 0.4 =⇒ ω1 = 2π/T = 15.7 Hz

0.5

0
-10 -5 0 5 10

0.5

0
-10 -5 0 5 10

-100

-200
10 -2 10 0 10 2 10 4 10 6

The spacing ∆ω is now ten times larger. Although the higher

frequencies are much smaller, there is still a reasonable amount
of ‘wiggle’ in the output.
11.6. Application: PWM of a DC motor 159

Case 3: T = 0.04 =⇒ ω1 = 2π/T = 157 Hz

0.5

0
-10 -5 0 5 10

0.5

0
-10 -5 0 5 10

-100

-200
10 -2 10 0 10 2 10 4 10 6

The spacing ∆ω is now large enough that, to a reasonable ap-

proximation, the output is constant. This is because the non-zero
frequencies are ‘far enough down’ the Bode gain plot that they
are filtered out by the DC motor and are therefore not seen in the
output.
11.7. Links and further reading 160

11.6.1 Summary of PWM for the DC motor

We have seen that
• T0 /T sets the amplitude of the zero-frequency component of
u(t) and therefore the zero-frequency component of y(t).
• T sets the frequencies of the non-zero frequency components of
u(t) and therefore the non-zero frequency components of y(t).
• By making T small enough, the non-zero frequency components
of u(t) are at frequencies where the gain of the DC motor,
|G(jω)|, is small:

|G(jωk )| ≈ 0 for k 6= 0.

• The physical explanation is that the pulse is so fast (relative to

the dynamics of the DC motor) that the DC motor isn’t ever
able to ’catch up’ with the rapid duty cycle. The DC motor
therefore ‘sees’ the pulse as a constant input with a reduced
amplitude.

11.7 Links and further reading

• How does Shazam work
11.A. The discrete Fourier transform (DFT) 161

11.A The discrete Fourier transform (DFT)

Note: This material is not examinable but it may be helpful for
the final part of assignment 3.
How would we find the Fourier transform using a computer? A
simple example might be that I give you a square wave or triangle
wave and ask you to find its Fourier transform; more complicated
examples might be that I give you a music file or a time series
of measured wind speed from a possible wind-farm site and ask
you to find the dominant frequencies in the signal. For this we
need a discrete version of the Fourier transform since the data will
be discrete—for example the wind speed data might be measured
every 0.01 seconds (which is a sampling frequency of 1/0.01 =
100 Hz); and music is often sampled at a sampling frequency of
44.1 kHz.
The discrete Fourier transform pair are defined as
N −1
2π
X
fˆn = fk e−j N nk (11.10)
k=0

and
N −1
1 X ˆ +j 2π nk
fk = fn e N (11.11)
N n=0

where fk represents the time-domain signal and fˆn represents its

discrete Fourier transform.
11.B. The fast Fourier transform (FFT) 162

We can alternatively write (11.10) as a matrix equation:

    
 fˆ0  1 1 1 ··· 1   f0 
    
    
 fˆ1  1
    
   e−j2π/N e−j4π/N ··· e−j2π(N −1)/N  
 f1 

    
    
 ˆ =
 f2  1 e−j4π/N e−j8π/N ··· e−j2π2(N −1)/N  
  f2 

    
    
 ..   .. .. .. ..
...   ... 
  
 .  . . . .
    
    
    
fˆN −1 1 e−j2π(N −1)/N e−j2π2(N −1)/N ··· e −j2π(N −1)(N −1)/N
fN −1
(11.12)
(and similarly we could write (11.11) as a matrix equation).
The matrix in (11.12) is of dimension N × N and solving (11.12)
requires performing N × N = N 2 operations. This can become
very costly when N is large. For example, the song used in lectures
is about 256 seconds in duration and the sampling frequency is
44,100 samples per second. So the total number of samples in the
song is N = 11, 289, 600 and if we square this we get over 127
trillion operations! (N 2 = 127.455 × 1012 .)

11.B The fast Fourier transform (FFT)

Efficient computer algorithms for calculating discrete Fourier trans-
forms have been developed since the mid-1960s. These are known
as fast Fourier transform (FFT) algorithms and they exploit the
fact that the standard DFT involves a lot of redundant calcula-
tions. Looking again at (11.12), we can already see from the few
entries of the matrix that are shown that e−j4π/N shows up twice
and that (e−j0π/N = 1) shows up several times.
The fast Fourier transform exploits this redundancy to signifi-
cantly reduce the number of operations from N 2 operations to
to N2 log2 N operations. Table 11.1 lists the savings provided by
11.C. Sampling frequency 163

the fast Fourier transform (FFT) when compared to the discrete

Fourier transform (DFT) for different values of N .

N
N N 2 (DFT) 2
log2 N (FFT) saving

32 1,024 80 92%

256 65,536 1,024 98%

1,024 1,048,576 5,120 99.5%

Table 11.1: Computational savings provided by the fast Fourier

transform (FFT) when compared to the discrete Fourier trans-
form (DFT) for different values of N .

11.C Sampling frequency

For the DFT and the FFT the signal will always be discrete: that
is, defined at discrete points in time. We call the time between
neighbouring points the sampling period Ts . Then the sampling
frequency Fs is the reciprocal of this:
1
Fs = .
Ts

11.D Examples
Perhaps the best thing to do now is to look at the FFT in action
by looking at some examples. Some of these examples are taken
almost verbatim from Matlab’s documentation for its fast Fourier
transform fft() (see links at the end of this chapter).
11.D. Examples 164

We will also look at simple examples of phenomena known as i)

spectral leakage and ii) aliasing.
11.D. Examples 165

Example 1: two sinusoids

In this first example we form a signal that contains two dis-
tinct frequencies (50 Hz and 120 Hz). We will see that the fast
Fourier transform gives us back these frequencies exactly in its
spectrum.

Fs = 1000; % Sampling frequency

T = 1/Fs; % Sampling period
L = 1500; % Length of signal
t = (0:L-1)*T; % Time vector

% signal composed of 50 Hz and 120Hz sinusoids

y = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

% plot first 100 values in time

figure, set(gcf,'Position', [10 10 400 150])
plot(t(1:100), y(1:100), 'k')
xlabel('t (s)'), ylabel('y(t)')

% take Fourier transform

f = Fs*(0:(L/2))/L;
Y = fft(y);
P2 = abs(Y/L); % two-sided spectrum
P1 = P2(1:L/2+1); % one-sided spectrum
P1(2:end-1) = 2*P1(2:end-1);

% plot one-sided spectrum

figure, set(gcf,'Position', [10 10 400 150])
stem(f, P1, 'k')
xlabel('f (Hz)'), ylabel(' | P 1(f) | ')
11.D. Examples 166

1
y(t)

-1

-2
0 0.02 0.04 0.06 0.08 0.1
t (s)
1
|P1(f)|

0.5

0
0 100 200 300 400 500
f (Hz)

Figure 11.3: Example 1. Top: the first 0.1 s of the signal in the
time domain. Bottom: fast Fourier transform.
11.D. Examples 167

Example 2: two sinusoids and random noise

We now take the same signal from example 1 and add random
noise to it. The fast Fourier transform still comes back with the
two frequencies (50 Hz and 120 Hz) as the most dominant, but
with other frequencies in the signal as well.

Fs = 1000; % Sampling frequency

T = 1/Fs; % Sampling period
L = 1500; % Length of signal
t = (0:L-1)*T; % Time vector

% signal composed of 50 Hz and 120Hz sinusoids

y = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

% add zero-mean white noise to y

y = y + 1*randn(size(t));

% plot first 100 values in time

figure, set(gcf,'Position', [10 10 400 150])
plot(t(1:100), y(1:100), 'k')
xlabel('t (s)'), ylabel('y(t)')

% take Fourier transform

f = Fs*(0:(L/2))/L;
Y = fft(y);
P2 = abs(Y/L); % two-sided spectrum
P1 = P2(1:L/2+1); % one-sided spectrum
P1(2:end-1) = 2*P1(2:end-1);

% plot one-sided spectrum

figure, set(gcf,'Position', [10 10 400 150])
stem(f, P1, 'k')
xlabel('f (Hz)'), ylabel(' | P 1(f) | ')
11.E. Spectral leakage 168

2
y(t)

-2

-4
0 0.02 0.04 0.06 0.08 0.1
t (s)
1
|P1(f)|

0.5

0
0 100 200 300 400 500
f (Hz)

Figure 11.4: Example 2. Top: the first 0.1 s of the signal in the
time domain. Bottom: fast Fourier transform.

11.E Spectral leakage

We will now look very briefly at spectral leakage. Note that
you do not need this to complete the last part of assignment
3, but nonetheless it is useful to know roughly what it is and
why it occurs because it does appear in the signals you have been
given.
Spectral leakage can occur if the original time-domain signal is
not perfectly sinusoidal over the length of that signal. This is
best seen with an example.
11.E. Spectral leakage 169

Fs = 20; % Sampling frequency

T = 1/Fs; % Sampling period
L = 20; % Length of signal
t = (0:L-1)*T; % Time vector

% signal x is periodic over its length

x = sin(2*pi*1*t);

% signal y is not periodic over its length

y = sin(2*pi*1.05*t);

% plot x and y in time

figure, set(gcf,'Position', [10 10 400 150])
plot(t, x, 'k-o'), hold on, plot(t, y, 'r-o')
xlabel('t (s)'), ylabel('x(t), y(t)')
legend('x(t)', 'y(t)'), set(gca,'Box','On')

% find fft of x and y

f = Fs*(0:(L/2))/L;
% for x
X = fft(x);
P2 = abs(X/L); % two-sided spectrum
P1 = P2(1:L/2+1); % one-sided spectrum
P1(2:end-1) = 2*P1(2:end-1);
PX = P1; % set spectrum of x to one-sided spectrum
% and for y
Y = fft(y);
P2 = abs(Y/L); % two-sided spectrum
P1 = P2(1:L/2+1); % one-sided spectrum
P1(2:end-1) = 2*P1(2:end-1);
PY = P1; % set spectrum of y to one-sided spectrum

% plot one-sided spectrum

figure, set(gcf,'Position', [10 10 400 150])
stem(f, PX, 'k'), hold on, stem(f, PY, 'r')
xlabel('f (Hz)'), ylabel(' | P x(f)| , | P y(f) | ')
legend(' | P x(f) | ', ' | P y(f) | '), set(gca,'Box','On')
11.E. Spectral leakage 170

1
x(t)
0.5 y(t)
x(t), y(t)

-0.5

-1
0 0.2 0.4 0.6 0.8 1
t (s)
1
|Px (f)|
|Px (f)|, |Py (f)|

|Py (f)|

0.5

0
0 2 4 6 8 10
f (Hz)

Figure 11.5: Example of spectral leakage. Top: signals x(t)

and y(t) in the time domain. Bottom: fast Fourier transforms of
x(t) and y(t).

The signal x(t) is periodic over the length of the signal:

x(t) = sin(2π × 1t)

and its FFT is non-zero only at f = 1 Hz, as expected (see bottom

part of figure 11.5). The signal y(t), meanwhile, is not periodic
over the length of the signal:

y(t) = sin(2π × 1.05t)

and we can note two things. First, that the largest component
of the fast Fourier transform of y(t) is at f = 1 Hz. (Which
is the nearest frequency to the true frequency of f = 1.05 Hz.)
11.F. Aliasing 171

Second, that the fast Fourier transform also has energy at other
frequencies—energy has ’leaked’ into other frequencies.

11.F Aliasing
Aliasing does not occur for the data you are given for the last
part of assignment 3 so you don’t need to worry about it. But I
feel I should still tell you a little about it.
It can be shown that, for a sampling frequency of Fs , the highest
frequency that can be resolved is
1
Fnyq = Fs
2
and this is known as the Nyquist frequency. Aliasing can occur
when the frequency at which a signal is sampled is too low when
compared to the frequencies that are present in the signal. In
particular, aliasing will occur if the highest frequency present in
the signal is higher than the Nyquist frequency, Fnyq .
This is best seen with an example, the code and figures for which
follow on the next two pages. The frequency (f0) contained in
the signal is 15 Hz. The signal x1 (t) has a sampling frequency
of 100 Hz, which is high enough to resolve this frequency (since
FNyq = Fs /2 = 100/2 = 50 Hz > 15 Hz). The spectrum therefore
(correctly) shows a peak at f = 15 Hz. The signal x2 (t) has a
sampling frequency of only 20 Hz and so the frequency 15 Hz is
not resolved (since in this case FNyq = 20/2 = 10 Hz < 15 Hz).
Instead it is aliased to f = 5 Hz—it is as if the true frequency has
been ‘folded over’ to a lower frequency.
11.F. Aliasing 172

f0 = 15; % choose frequency of the signal

% sampling frequency of 100 Hz

Fs = 100; T = 1/Fs; L = 100;
t1 = (0:L-1)*T; % time vector
x1 = sin(2*pi*f0*t1); % signal
f1 = Fs*(0:(L/2))/L; % frequency vector
X1 = fft(x1);
P2 = abs(X1/L); P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
PX1 = P1; % one-sided spectrum

% sampling frequency of 20 Hz
Fs = 20; T = 1/Fs; L = 20;
t2 = (0:L-1)*T; % time vector
x2 = sin(2*pi*f0*t2); % signal
f2 = Fs*(0:(L/2))/L; % frequency vector
X2 = fft(x2);
P2 = abs(X2/L); P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
PX2 = P1; % one-sided spectrum

% plot x1 and x2 in time

figure, set(gcf,'Position', [10 10 400 150])
plot(t1, x1, 'k'), hold on, plot(t2, x2, 'ro-')
xlabel('t (s)'), ylabel('x(t), y(t)')
set(gca,'Box','On')

% plot one-sided spectrum

figure, set(gcf,'Position', [10 10 400 150])
stem(f1, PX1, 'k'), hold on, stem(f2, PX2, 'r')
xlabel('f (Hz)'), ylabel(' | P x(f)| , | P y(f) | ')
legend(' | P 1(f) | ', ' | P 2(f) | ')
axis([0 20 0 1]), set(gca,'Box','On')
11.G. Links and further reading 173

0.5
x(t), y(t)

-0.5

-1
0 0.2 0.4 0.6 0.8 1
t (s)
1
|P1 (f)|
|Px (f)|, |Py (f)|

|P2 (f)|

0.5

0
0 5 10 15 20
f (Hz)

Figure 11.6: Example of aliasing. Top: signals in the time do-

main. Bottom: spectra obtained using the fast Fourier transform.

11.G Links and further reading

• Fast Fourier transform
• 44,100 Hz
• Camera shutter speed synchronized with helicopter blade fre-
quency
• Guitar oscillations captured with iPhone 4
• rotor panorama
12 Discrete-time systems

So far all of the dynamical systems we have considered have been

continuous in time. In this chapter we look at how to model dy-
namical systems that are discrete in time. By discrete we mean
that, rather than the input and output varying continuously in
time, they are instead defined only at discrete instances in time.
(With the interval between them determined by the sampling pe-
riod.)
In this chapter we will look at i) how to analyse systems that are
discrete-time from the beginning; and ii) how to take a continuous-
time system and convert it into a discrete-time system.

12.1 Motivation: Line-follower robot

Let’s motivate our study of discrete-time systems by thinking
about a specific example: the line-following robot shown below.

174
12.2. Difference equations 175

What is continuous about this system?

• the dynamic motion of the vehicle and its motors
What is discrete about this system?
• discrete-time controller: on-board microprocessor samples at
its clock speed
• discrete sensor signal: the sensor operates in discrete time, mea-
suring the sensed signal at discrete points in time determined
by its sampling period
• signals sent to the motor: although the motor is continuous,
the signals it receives are discrete in time

12.2 Difference equations

So far we have considered systems governed by differential equa-
tions like

an y (n) (t) + · · · + a1 ẏ(t) + a0 y(t) = bm u(m) (t) + · · · + b1 u̇(t) + b0 u(t).

Using Laplace transforms, this leads to an equation of the form

(an sn + · · · + a1 s + a0 )Y (s) = (bm sm + · · · + b1 s + b0 )U (s)

and from there the transfer function:

Y (s) bm s m + · · · + b1 s + b 0
G(s) = = .
U (s) an s n + · · · + a1 s + a0

If we want to implement systems like this digitally on a computer,

say, then we need to represent the system discretely (in time)
using a difference equation:

an y[k−n]+· · ·+a1 y[k−1]+a0 y[k] = bm u[k−m]+· · ·+b1 u[k−1]+b0 u[k].

12.3. The z transform 176

We will see that, using z transforms, this leads to an equation of

the form

(an z −n + · · · + a1 z −1 + a0 )Y (z) = (bm z −m + · · · + b1 z −1 + b0 )U (z)

and a discrete-time transfer function, G(z):

Y (z) bm z −m + · · · + b1 z −1 + b0
G(z) = = .
U (z) an z −n + · · · + a1 z −1 + a0

12.3 The z transform

We call the relation between y[n] and Y (z) the z transform and
it is given by
∞
X
Y (z) = y[k]z −k .
k=0

The z transform maps a function of discrete time-step k to a

function of z. The z transform can be applied both to signals
and to systems.

12.3.1 Examples of z transform of signals

From now on we will often use the notation

y[k] = yk

just because it is a little bit neater.

z transform of a discrete impulse

The discrete impulse is 1 at the zeroth time-step and zero every-
where else:
12.3. The z transform 177

Ukr
1

Ukr
O 1
1 2 3 4 5 K

Its z transform isUkr

1
∞
X
U (z) = uk z −k
O
k=0 1 2 3 4 5 K
0 −1 −2
Ukr
= u0 z + u1 zUkr+ u2 z + ...
1 O 1
= 1 1i 2 3 4 5 K

Ukr
z transform of a1 delayed discrete impulse
On the left below we delay the impulse by one time-step; on the
O O i
right we 1delay2 it by
3
n time-steps:
4 5 K 1 2 3 4 5 K

Ukr Ukr
1 1
O i
n 2 n I n n 11 nt2 K

O i O i
1 2 3 4 5 K n 2 n I n n 11 nt2 K

IfUkr
we delay the impulse by one time-step (on the left above) then
1
we get a z transform of

U (z) = 0 + 1 · z −1 + 0 + · · ·
= z −1 .
O i
n 2 n I n n 11 nt2 K
12.3. The z transform 178

More generally, if we delay the impulse by n time-steps (on the

right above) then we get a z transform of

U (z) = 0 + 0 + · · · + 1 · z −n + 0 + · · ·
= z −n .

z transform of a discrete unit step

The discrete step is 1 for all non-negative time-steps:
Ukr
1 o

O e i e i l
1 2 3 4 5 K

Its z transform is

U (z) = 1 + z −1 + z −2 + z −3 + · · ·

We can simplify this by using

∞
X 1
rk = (for |r| < 1)
k=0
1−r

and then setting r = z −1 to give

1 z
U (z) = −1
= .
1−z z−1

12.3.2 z −1 as the unit delay

What is z? We can think of z −1 as representing the unit de-
lay:
 o i z z 4 I
a K
i
12.3. Theo z transform
z z 4 I K 179

o i z z 4 I K
un i
un i

o i z z 4 I K un i
−n
and we can think of z as representing a delay of n time-steps:
un n
un n
un i

un n

uan for which the unit

We can compare this to the Laplace transform,
uan
delay (with sampling time T ) would be represented by
n un

uan

h
h
y(t) = u(t − T ) −→
uan
Y (s) = e−sT U (s)

h
and for a delay of n time-steps we would have

h
y(t) = u(t − nT ) −→ Y (s) = e−snT U (s)

We might therefore expect that the relationship between s and z

would be
z −1 ←→ e−sT
and therefore
z ←→ esT
and this is indeed true. (We will use this later to convert between
the s- and z-domains.)
12.4. Mapping the s-plane to the z-plane 180

12.3.3 Properties of the z transform

There are two important properties of the z transform that we
will use:
1. linearity: If we have

ax[k] + by[k]

then when we take z transforms we get simply

aX(z) + bY (z).

2. delay: a delay in the discrete-time domain is like multiplica-

tion by z −1 in the z domain:
z
y[k] →
− Y (z)
z
− z −1 Y (z)
y[k − 1] →
z
− z −n Y (z).
y[k − n] →

In lectures we will look at some examples where we use these

properties to take z transforms of a difference equation.

12.4 Mapping the s-plane to the z-plane

How do points in the Laplace (s) domain map to points in the
z domain? For this we can use the fact that the way to convert
between them is using
z = esT .
Let’s now look at how points in the s-plane are mapped to points
in the z-domain. This is shown below in figure 12.1, the code for
which follows.
12.4. Mapping the s-plane to the z-plane 181

4 2

2 1
Im(s)

Im(z)
0 0

-2 -1

-4 -2
-1 -0.5 0 0.5 1 -2 -1 0 1 2
Re(s) Re(z)

Figure 12.1: Mapping of the s-plane (left) to the z-plane (right).

Red lines map to red lines; green to green; and blue to blue.

T = 1; % sampling period
% choose points in s-plane
s1 = -0.5 + 1j*pi*[-1:0.05:1]; % stable (LHP)
s2 = 0.0 + 1j*pi*[-1:0.05:1]; % marginally stable
s3 = +0.5 + 1j*pi*[-1:0.05:1]; % unstable (RHP)
% map each line to points in z
z1 = exp(s1*T); z2 = exp(s2*T); z3 = exp(s3*T);
% plot points in s plane
figure, set(gcf,'Position', [10 10 200 200])
plot(real(s1), imag(s1), 'bo-'), hold on
plot(real(s2), imag(s2), 'go-')
plot(real(s3), imag(s3), 'ro-')
xlabel('Re(s)'), ylabel('Im(s)')
set(gca, 'Box', 'On'), pbaspect([1 1 1])
axis([-1. 1. -4. 4.])
% plot points in z plane
figure, set(gcf,'Position', [10 10 200 200])
plot(real(z1), imag(z1), 'bo-'), hold on
plot(real(z2), imag(z2), 'go-')
plot(real(z3), imag(z3), 'ro-')
xlabel('Re(z)'), ylabel('Im(z)')
set(gca, 'Box', 'On'), pbaspect([1 1 1])

12.4.1 Going the other way

What if we want to map from the z-plane back to the s-plane?
For this we can take the mapping from s to z:

z = esT
12.4. Mapping the s-plane to the z-plane 182

and rearrange for s to get

ln z
s= .
T
Remember that z is in general complex. How do we take the
logarithm of a complex number? It is not too difficult provided
we use
z = reiθ .
Plugging this in:

ln z = ln(reiθ )
= ln r + ln(eiθ )
= ln r + iθ.

So finally we have
ln r + iθ
s=
T
where r and θ are defined such that z = reiθ .
12.4. Mapping the s-plane to the z-plane 183

12.4.2 The z transform and stability

First a reminder: for Laplace transforms (the s domain), a system
is stable if all of its poles lie inside the left half-plane:
Im

Im

For the z transform (the z domain), a system G(z) is stable if all

of its poles lie inside the unit circle:

We can see this from figure 12.1 in the previous section, where
i) points in the left half-plane in s mapped to points inside the
unit circle; ii) points in the right half-plane in s mapped to points
outside the unit circle; and iii) points on the imaginary axis in s
mapped to points on the unit circle.
12.5. Converting from s to z 184

12.5 Converting from s to z

What if we have a transfer function in continuous time (the trans-
fer function of a controller, for example) and want to convert it to
discrete time to implement it digitally? Clearly we need a method
to convert from continuous-time (s-domain) to discrete time (z-
domain). A number of methods exist to do this but we will focus
on only one: the bilinear transform.

12.5.1 The bilinear transform (Tustin’s method)

The bilinear transform is defined by the equation
2 z−1
s= .
T z+1
This tells us how to get from the s-domain to the z-domain—it
tells us that, every time we see an s, we should replace it with
2 z−1
T z+1
. It is also good to know how to go the other way—from the
z-domain to the s-domain. For this we have the inverse trans-
form:
1 + T2 s
z= ,
1 − T2 s
which tells us what we should replace z with to get an expression
in terms of s.
This method for converting between the s- and z-domains can go
by a number of different names:
• bilinear transform
• Tustin’s method
• trapezoidal integration rule
Where do these expressions come from? For more details on
this you can have a look at the additional material in section
12.A.
12.6. Frequency warping 185

12.6 Frequency warping

Tustin’s method (and the other methods mentioned) is an ap-
proximation. It will be a good approximation provided that the
frequencies of interest are low relative to the Nyquist frequency,
defined as
1
Fnyq = Fs
2
where Fs is the sampling frequency. If we denote by ωa the
continuous-time frequency and by ωd the discrete-time frequency
to which it maps following Tustin’s method then it turns out
that  
2 T
ωd = arctan ωa .
T 2
(For more details on where this expression comes from, see section
12.B.) This is plotted in figure 12.2. For small ωa the slope of the
curve is 1 giving ωd = ωa and there is no warping of frequencies.
However, as the frequency increases the slope of the curve flattens
and we have ωd < ωa . As ωa becomes larger, ωd asymptotes to
the value
π2 π
ωd = = = 2πfnyq = ωnyq .
2T T
This implies that higher frequencies all get ‘squashed’ into a rel-
atively tight band of discrete frequencies with these discrete fre-
quencies satisfying ωd ≤ π/T .
12.6. Frequency warping 186

2
d

0
0 5 10 15 20 25 30

Figure 12.2: Plot of ωd versus ωa for a sampling time of T = 1

(which corresponds to Fnyq = 1/2 and ωnyq = π).

T = 1; % sampling time

oma = [0:0.001:1] * 30;

omd = 2/T * atan(T/2*oma);

figure, set(gcf,'Position', [10 10 400 150])

plot(oma, omd, 'k')
xlabel('\omega a'), ylabel('\omega d')
axis([0 30 0 3.2])
12.A. Deriving the bilinear transform 187

12.A Deriving the bilinear transform

We will derive the bilinear transform by looking at a specific case.
Suppose that the output is related to the input by
Z
y(t) = u(τ ) dτ,

i.e. an integrator. This is a continuous function of time and taking

Laplace transforms gives
1
Y (s) = U (s).
s
What about in discrete time? For this we need to approximate
the integral—doing so will give us a difference equation which is
discrete in time. Below is plotted the input u(t) as a function of
time which, in discrete time, is sampled at the times t0 , t1 , t2 and
so on to give u0 , u1 , u2 and so on.

ru
Uz 43
44

the

t
to f Iz e ta

Using the trapezoidal rule to approximate the integration of u to

12.A. Deriving the bilinear transform 188

give y, assuming y0 = 0, gives

1
y1 = (u0 + u1 )T
2
1
y2 = y1 + (u1 + u2 )T
2
..
.
1
yk = yk−1 + (uk−1 + uk )T.
2
Rearranging the final equation we have
1 1
yk − yk−1 = T uk + T uk−1 .
2 2
Taking z transforms:
1
(1 − z −1 )Y (z) = T (1 + z −1 )U (z)
2
and rearranging to get the discrete-time transfer function, G(z):
1
2
+ z −1 )
T (1
G(z) =
1 − z −1
or
T z+1
G(z) =
2 z−1
Now remember that the transfer function in continuous time was
1
G(s) =
s
and therefore the correspondence between s and z must be
1 T z+1
=
s 2 z−1
or
2 z−1
s=
T z+1
which is the bilinear transform.
12.B. Relationship between ωa and ωd 189

12.B Relationship between ωa and ωd

The discrete-time transfer function Gd (z), comes by setting s =
2 z−1
T z+1
(bilinear transform) in the continuous-time transfer function
G(s):

Gd (z) = G(s)|s= 2 z−1

T z+1

2 z−1
 
=G .
T z+1

We want to know how this (approximation) comes out if we set

s = esT = ejωd T (which is the exact expression for z):
 jωd T
−1

jωd T 2e
Gd (z) = Gd (e )=G
T ejωd T + 1

!
2 ejωd T /2 ejωd T /2 − e−jωd T /2
=G
T ejωd T /2 (ejωd T /2 + e−jωd T /2 )

!
2 ejωd T /2 − e−jωd T /2
=G
T (ejωd T /2 + e−jωd T /2 )
!
2 ejωd T /2 − e−jωd T /2 /(2j)


=G j
T (ejωd T /2 + e−jωd T /2 ) /2
 
2 sin(ωd T /2)
=G j
T cos(ωd T /2)
 
2
= G j tan (ωd T /2)
T

This is like evaluating the continuous-time G(s) at the frequency

 
2 T
ωa = tan ωd .
T 2
12.B. Relationship between ωa and ωd 190

Rearranging, the discrete-time frequency ωd is warped (relative

to ωa ) by an amount
 
2 T
ωd = arctan ωa .
T 2
13 State-space and Laplace connections

In this chapter we will look at some connections between the state-

space and Laplace-transform descriptions of a linear system.

13.1 Introduction
We have seen that there are two different ways to represent a
linear dynamical system:
1. using a state-space model in chapter 2;
2. using Laplace transforms in chapter 4 (which leads to a transfer
function covered in chapters 5–7).
As we might expect, there are close connections between the trans-
fer function of a system and its state-space representation. We
will now look at these connections: first for autonomous systems
(in which there is no input); and second for non-autonomous sys-
tems (in which the effect of some input(s) is included).

13.2 Autonomous system ẋ = Ax

An autonomous system is one in which there is no input:

ẋ(t) = Ax(t)

191
13.2. Autonomous system ẋ = Ax 192

with initial condition x(0) = x0 . Remember (chapter 2) that we

can solve this using the matrix exponential:
x(t) = eAt x0 .
We can also solve using Laplace transforms:
L {ẋ(t) = Ax(t)}
sX(s) − x0 = AX(s)
so
sX(s) − AX(s) = x0
(sI − A)X(s) = x0 .
Then the solution in time is
x(t) = L−1 (sI − A)−1 x0


= L−1 (sI − A)−1 x0 .



See §13.A for some help on finding the inverse of a matrix.

13.2.1 Connection between eAt and (sI − A)

To summarize, the two solutions of
ẋ(t) = Ax(t) with x(0) = x0
are
x(t) = eAt x0
and
x(t) = L−1 (sI − A)−1 x0 .


The two solutions must of course give the same answer and so

eAt = L−1 (sI − A)−1 = Φ(t)



where Φ(t) is called the state-transition matrix.

13.3. Non-autonomous system ẋ = Ax + Bu 193

13.3 Non-autonomous system ẋ = Ax + Bu

This time we
1. include the effect of inputs u(t) (non-autonomous system)
2. have a second equation for the output(s) of interest y(t).
So this time we have

ẋ(t) = Ax(t) + Bu(t)

y(t) = Cx(t) + Du(t)

we assume x(0) = 0 because our interest is now in the system’s

response to inputs. Let us now list the different terms and their
dimensions:
• um×1 is the input; yp×1 is the output; xn×1 is the state
• An×n is the state matrix
• Bn×m is the input matrix
• Cp×n is the output matrix
• Dp×m is called the direct transmission matrix.
Now we take Laplace transforms to get a matrix of transfer func-
tions:

L {ẋ(t) = Ax(t) + Bu(t)}

sX(s) = AX(s) + BU (s)

so

sX(s) − AX(s) = BU (s)

(sI − A)X(s) = BU (s)
X(s) = (sI − A)−1 BU (s)
13.3. Non-autonomous system ẋ = Ax + Bu 194

and we do the same for the output equation:

L {Y = CX + DU }
Y (s) = CX(s) + DU (s)
putting them together we get
Y (s) = CX(s) + DU (s)
= C(sI − A)−1 BU (s) + DU (s)
= (C(sI − A)−1 B + D)U (s)
= G(s)U (s)
and so
G(s) = C(sI − A)−1 B + D
G(s) is a matrix of transfer functions of dimensions p × m, i.e.
(number of outputs) × (number of inputs).
If we were to write the whole thing out then we would get
      
y1 (s)  u1 (s)  G11 (s) · · · G1m (s)
  u1 (s) 
 
    
      
 u (s)   ... .. 
      
y (s) .   u (s) 
 2   2    2 
  = G(s)p×m  =  
 ..   ..     .. 
      
 .   .    . 
      
      
      
yp (s) um (s) Gp1 (s) Gpm (s) um (s)

Therefore the transfer function describing the dynamics from in-

put uj (t) to output yi (t) is Gij (s). In other words,
yi (s) = Gij (s)uj (s)
or
yi (s)
= Gij (s).
uj (s)
13.A. Calculating B −1 for a matrix Bn×n 195

13.A Calculating B −1 for a matrix Bn×n

The inverse of some general B is given by
1
B −1 = Adj(B)
det(B)

where
 
 (−1)1+1 M1,1 · · · (−1)j+1 Mj,1 · · · (−1)n+1 Mn,1 
 
 
 .. ... 

 . 

 

Adj (B) =  .. 
 (−1)1+i M1,i (−1)j+i Mj,i .


 
 
 .. .. 

 . . 

 
 
(−1)1+n M1,n ··· (−1)n+n Mn,n

where Mi,j is the “i, j” minor of B, which is the determinant of the

matrix (denoted B(i),(j) ) obtained by deleting row i and column j
of B.

Application to some matrix B

 
1 2 3
 
 
B = 4 5 6
 
 
 
 
7 8 9
13.A. Calculating B −1 for a matrix Bn×n 196

What is M1,2 ?
 
 
 
 
  4 6
B(1),(2) = 4 =
   
 6 





 7 9
7 9

and
M1,2 = detB(1),(2) = 4 · 9 − 7 · 6 = −6.

Application to (sI − A)
Applying this to B = (sI − A), we get

1 N (s)
(sI − A)−1 = · Adj(sI − A) = .
det(sI − A) D(s)

The scalar det(sI − A) is a polynomial in s and each element in

the matrix Adj(sI − A) is a polynomial in s.
   
1 −2 s − 1 2
A= =⇒ sI − A = 
   
 
   
0 0 0 s
13.A. Calculating B −1 for a matrix Bn×n 197

What is (sI − A)−1 ?

det(sI − A)−1 = (s − 1)s

 
s −2 
Adj(sI − A) = 
 

 
0 s−1
   
1
s −2   (s−1) −2 
−1 1 (s−1)s 
(sI − A) = =
  
(s − 1)s 
 
  
1
0 s−1 0 s
14 MDoF vibration

In this final chapter we will look at vibrating mechanical systems

when there are two or more degrees of freedom. We will only look
at the case where there is i) no damping and ii) no forcing. For
this case we will find the governing equations and then look at
how to find i) the natural frequencies; and ii) the mode shapes of
the vibrating system.

14.1 The two-mass case

We can write the two-mass case in matrix form as

M q̈(t) + Kq(t) = 0

with
   
m1 0  k1 + k2 −k2 
M = ; K= .
   
   
0 m2 −k2 k2 + k3

we will derive these equations in lectures.

Notice that both matrices are symmetric:

M = MT ; K = KT .

198
14.2. The three-mass case 199

Physically, this symmetry makes sense. If the first mass moves,

it generates a force on the second mass because of the spring
stiffness k2 ; and if the second mass moves, it generates a force on
the first mass because of that same spring stiffness.

14.2 The three-mass case

When there are three masses we can still write the governing
equations in matrix form as
M q̈(t) + Kq(t) = 0
where this time
   
m1 0 0 k1 + k2 −k2 0 
   
   
M =  0 m2 0  ; K =  −k2 .
   
   k2 + k3 −k 3 
   
   
0 0 m3 0 −k3 k3 + k4

14.3 The n-mass case

You might be able to see where this is going. When there are n
masses we can still write the governing equations in matrix form
as
M q̈(t) + Kq(t) = 0
where this time
   
m1 0 0  k1 + k2 −k2 0 
   
   
M =  0 ... 0  ; K =  −k2 .. .
   
   . −kn  
   
   
0 0 mn 0 −kn kn + kn+1
14.4. Natural frequencies and mode shapes 200

In this most general case we see that

• the mass matrix M is diagonal
• the stiffness matrix is tridiagonal, which means that all entries
are zero except for the diagonal, the line above the diagonal,
and the line below the diagonal.

14.4 Natural frequencies and mode shapes

Our equation of motion is

M q̈(t) + Kq(t) = 0. (14.1)

Now setting1 q(t) = q̃est :

(s2 M + K)q̃ = 0, (14.2)

and letting s2 = −λ:

K q̃ = λM q̃. (14.3)
We can make this look like an eigenvalue problem by multiplying
both sides by M −1 :
M −1 K q̃ = λq̃.
Then the eigenvalues can be found from

det(M −1 K − λI) = 0.

From this we can find the n values of the undamped natural fre-
quencies, ω12 , ω22 , . . . , ωn2 . Substituting any one of these back into
(14.2), we get a corresponding set of relative values for the vector
q̃, which contain the eigenvectors or mode shapes, which we will
1
By doing this we are saying that we expect the vibrating system to
oscillate at some frequency given by est ; and that the oscillation will have
some shape given by the vector q̃. In other words, the masses will oscillate
together at some common frequency and with some overall shape.
14.4. Natural frequencies and mode shapes 201

label ψ from now on. There are n of these too, ψ1 , ψ2 , . . . , ψn .

The mode shapes ψi are not unique: we can double them, say,
and they are still an eigenvector. (This is always true of eigen-
vectors.)
We will look at some examples of finding natural frequencies and
mode shapes in lectures.