MATH50004/50015/50019

Differential Equations

Lecture Notes

Spring Term 2024/25

Martin Rasmussen
Imperial College London
c Martin Rasmussen (2025).
These notes are provided for the personal study of students taking this
module. The distribution of copies in part or whole is not permitted.
Contents

Chapter 1. Introduction 5
§1. Ordinary differential equations and initial value problems 7
§2. Examples 9
§3. Visualisations 12
Chapter 2. Existence and uniqueness 17
§1. Picard iterates 17
§2. Lipschitz continuity 20
§3. Picard–Lindelöf theorem 25
§4. Maximal solutions 29
§5. General solutions and flows 32
Chapter 3. Linear systems 39
§1. Matrix exponential function 39
§2. Planar linear systems 43
§3. Jordan normal form 52
§4. Explicit representation of the matrix exponential function 54
§5. Exponential growth behaviour 56
§6. Variation of constants formula 58
Chapter 4. Nonlinear systems 61
§1. Stability 61
§2. Limit sets 75
§3. Lyapunov functions 79
§4. Poincaré–Bendixson theorem 87

3
 Chapter 1

Introduction

An algebraic equation over the real numbers is solved by real numbers. For
instance, x2 − 1 = 0 is solved by x = 1 and x = −1. In contrast, ordinary
differential equations have functions as their solutions. Let us look at an
example.

Example 1.1 (A first example). We consider the differential equation

ẋ = ax , (1.1)

where a ∈ R is a constant and ẋ means dx

dt . We say that a function λ : I → R,
where I ⊂ R is an interval, solves this differential equation if

dλ
λ̇(t) = (t) = aλ(t) for all t ∈ I .
dt
The argument t of the solution λ typically stands for time, and x in (1.1)
typically represents the state of a physical, ecological, or other system, so
the solution λ describes the evolution of the state x in time. It should be
noted that there are also lots of applications where t does not stand for time.
The differential equation (1.1) is a very simple but realistic model of appli-
cations in nature and society.
For instance, if a > 0, this models growth of capital with a interest rate
linked to a (note that normally interest rates are given as yearly rates, which
correspond to a discrete-time model; a here is a continuous-time interest
rate; try as an exercise to convert both rates). For positive capital x, the
right hand side of (1.1) is positive, so ẋ is increasing, and importantly for
the model of capital growth, the increase in capital is proportional to the
amount of capital available.

5
6 1. Introduction

In contrast, if a < 0, then, if x > 0, the right hand side of (1.1) is negative,
and ẋ is decreasing, proportional to x. This models, for instance, radioactive
decay, which describes the decay of certain atoms such as Uranium 238.
It is easy to see that for a given b ∈ R, the function λb : R → R,
λb (t) = beat for all t ∈ R ,
solves (1.1), see Figure 1.1. Are there more solutions to this differential
equation? Assume there is another solution µ : I → R, where I ⊂ R is an
interval. Then
d
(µ(t)e−at ) = µ̇(t)e−at − µ(t)ae−at = aµ(t)e−at − µ(t)ae−at = 0
dt
for all t ∈ I. Hence, µ(t)e−at ≡ b for some b ∈ R, so µ(t) = beat = λb (t) for
all t ∈ I, which is not a new solution, so all solutions to (1.1) are known to
us.

Figure 1.1. Solutions to the differential equation (1.1) for α = 1 > 0

(left) and α = −1 < 0 (right).

So in contrast to (simple) algebraic equations whose the solutions are in

finite-dimensional vector spaces (such as Rd ), differential equations are
solved by functions, and spaces of functions are typically infinite-dimensional
and studied in the mathematical discipline functional analysis. Infinite-
dimensional spaces are more difficult to grasp in general. However, for a
vast majority of the material covered in this course, a finite-dimensional
thinking and visualisation is enough to understand the material very well.
In places, however, we will need some material from functional analysis to
understand differential equations better.
You have encountered ordinary differential equations already last year. In
particular, you have learned how to solve certain types of ordinary differ-
ential equations. It should be noted that for the most interesting ordinary
differential equations, it is not possible to find solutions analytically. In
this course, we will learn techniques how to still understand the solutions
to these equations without knowing them explicitly. We also address the
1. Ordinary differential equations and initial value problems 7

important question when solutions to ordinary differential equations exist

and are unique.

1. Ordinary differential equations and initial value problems

In this section, we look at the definition of an ordinary differential equation
and an initial value problem, and we study basic examples.
In Example 1.1, we have studied a differential equation of the form
ẋ = f (x) (1.2)
with f : Rd → Rd , where d = 1 and f (x) = αx. Although such type of
equations (i.e. autonomous and first order ) will be mostly studied in this
course, we would also like to deal with nonautonomous differential equations,
i.e. where the right hand side of (1.2) depends on time t, for instance, ẋ = tx2
(see Example 1.8 below). We note that higher-order differential equations
also generalise the situation in (1.2), and you have studied such differential
equations already in Year 1. It is demonstrated in Repetition Material 1
that such differential equations can always be transformed to first-order
differential equations, so no separate treatment (with regard to the general
theory) is necessary. An example of a higher-order differential equation is
given by the harmonic oscillator ẍ = −x (see Example 1.10 below).
The setup for nonautonomous first-order differential equations is explained
in the next definition.
Definition 1.2 (Ordinary differential equation). Consider d ∈ N, an open
set D ⊂ R × Rd , and a function f : D → Rd . An equation of the form
ẋ = f (t, x) (1.3)
is called a d-dimensional (first-order) ordinary differential equation. A dif-
ferentiable function λ : I → Rd on an interval

I ⊂ R is called a solution to
the differential equation (1.3) if t, λ(t) ∈ D and


λ̇(t) = f t, λ(t) for all t ∈ I . (1.4)

An ordinary differential equation (1.3) is called autonomous if the right hand

side does not depend on t, i.e. (1.3) is of the form
ẋ = f (x) ,
where f : D → Rd for some open set D ⊂ Rd . In this case, we also use the
symbol D for the domain of the right hand side f , here as a subset of Rd
instead of R × Rd , but this should not cause confusion, as it will be clear
from the context. We note that any autonomous differential equation can be
interpreted as a nonautonomous differential equation (1.3), and the domain
8 1. Introduction

D ⊂ Rd then translates to the domain R × D, which is an open set if and

only if D is open.
We will only treat ordinary differential equations (ODEs) in this course. Of
great importance are also partial differential equations (PDEs), which are
solved by functions depending on more than one variable, so, in contrast
to ordinary differential equations, partial differentiation is needed to even
define a partial differential equation.
The easiest types of solutions are constant solutions, which are also called
equilibrium solutions. If the differential equation is autonomous, they are
found algebraically, by zeros of the right hand side.
Proposition 1.3 (Constant solutions to autonomous differential equations).
Consider an open set D ⊂ Rd and a function f : D → Rd , and consider the
autonomous differential equation
ẋ = f (x) .
Then there exists a constant solution λ : R → Rd of this differential equation
with a ∈ Rd , i.e. λ(t) = a for all t ∈ R, if and only if f (a) = 0.

Proof. (⇒) Suppose that λ : I → Rd is a constant solution, i.e. λ(t) = a

for all t ∈ I. The solution identity yields
λ̇(t) = f (λ(t)) for all t ∈ I , (1.5)
which implies f (a) = 0.
(⇐) Suppose that f (a) = 0 for some a ∈ Rd . Then for the constant function
λ : R → Rd , λ(t) = a, the solution identity (1.5) is clearly fulfilled, and thus,
the constant function λ is a solution to the given differential equation. 

This proposition says that constant solutions to autonomous ordinary differ-

ential equations are easy to find. For many differential equations, constant
solutions are the only solutions that can be given explicitly, which means
that there are no formulas for all other solutions. In fact, most of the in-
teresting differential equations cannot be solved analytically. There are two
approaches to overcome this deficit. Firstly, there are numerous schemes to
numerically approximate solutions of differential equations – this will not be
covered in this course. Secondly, the so-called qualitative theory of ordinary
differential equations provides insights into how solutions behave without
knowing them explicitly – you will learn some basic elements of this theory
in this course.
We are interested now to understand in solutions for a given pair of initial
time and initial conditions. For the model discussed in Example 1.1, this
would mean that we are interested in the time evolution of capital, given
that we have x0 capital at time t0 .
2. Examples 9

Definition 1.4 (Initial value problem). Consider d ∈ N, an open set D ⊂

R × Rd , and a function f : D → Rd . The combination of the ordinary
differential equation
ẋ = f (t, x)
with an initial condition of the form
x(t0 ) = x0 , (1.6)
where (t0 , x0 ) ∈ D, is called an initial value problem, and (1.6) is called
initial condition. A solution to the above initial value problem is a solution
λ : I → Rd to the differential equation such that t0 is in the interior of I
and
λ(t0 ) = x0 .

We now solve an initial value problem for the simple differential equation
ẋ = ax.

Example 1.5 (A first example revisited). Consider the ordinary differential

equation (1.1) from Example 1.1 with the solutions λb for b ∈ R. For fixed
t0 , x0 ∈ R, we show that there exists a unique solution µ : R → R to (1.1)
solving the initial condition x(t0 ) = x0 . This follows from
λb (t0 ) = x0 ⇔ beat0 = x0 ⇔ b = x0 e−at0 .
Hence, the solution to this initial value problem is given by µ(t) = x0 ea(t−t0 )
for all t ∈ R.

We will later find conditions for ordinary differential equations that guaran-
tee that all initial value problems have a unique solution. These conditions
are rather weak and apply to large classes of applications.

2. Examples
We have seen that the differential equation ẋ = x behaves as nicely as one
can imagine:
(i) a solution exists for every initial value problem,
(ii) the solution to each initial value problem is unique,
(iii) the solution to each initial value problem exists globally, i.e. can be
defined on I = R.
In this section, we look at examples for which not all of this properties
are satisfied. The first example demonstrates that solutions to initial value
problems do not need to exist.
10 1. Introduction

Example 1.6 (No solution to an initial value problem). Consider the one-
dimensional initial value problem

1 : x<0
ẋ = f (x) = , x(0) = 0 ,
−1 : x ≥ 0
which has a discontinuous right hand side. Show as an exercise that this
initial value problem does not have any solutions.
There may exist more than one solution to an initial value problem.
Example 1.7 (Many solutions to an initial value problem). Consider the
one-dimensional initial value problem
p
ẋ = f (x) := |x| , x(0) = 0 . (1.7)
Since f (0) = 0, Proposition 1.3 implies that there exists a constant solution
with value 0. In addition, for any b ≥ 0, the function λb : R → R,

0 : t≤b
λb (t) = 1 2 : t>b ,
4 (t − b)
is a solution to this initial value problem. To check this, we have to verify
the solution identity for t < b, t = b and t > b. Clearly, λ̇(t) = 0 = f (λ(t))
d 1 2 1
for all t < b. Note that dt 4 (t − b) = 2 (t − b), which is 0 at t = b, so the
1
p
identity also holds at t = b. For t > b, we have λ̇(t) = 2 (t − b) = |λ(t)|,
which finishes the proof.
Question: can you find even more solutions to this initial value problem?

Figure 1.2. Three different solutions to the initial value problem (1.7)
(b = 1, 2, 3).

We now study a differential equation for which there are solutions that do
not exist for all times, i.e. they can only be defined on a proper subset
I ( R of the real numbers. To solve this differential equation, we need the
separation of variables technique that you have learned in your first year.
This technique is useful to compute initial value problems of the form
ẋ = g(t)h(x) , x(t0 ) = x0 , (1.8)
2. Examples 11

where g : I → R and h : J → R are two continuous functions, with two

open intervals I, J ⊂ R, and we assume that h(x0 ) 6= 0 (otherwise the initial
value problem has the constant solution).
The formal procedure is as follows:
dx dx
= g(t)h(x) , x(t0 ) = x0 =⇒ = g(t)dt , x(t0 ) = x0
dt h(x)
Z x Z t
1
=⇒ dy = g(s)ds ,
x0 h(y) t0

and solving this equation with respect to x will give a solution to the differ-
ential equation (1.8).
We study the separation of variables procedure for the following example.

Example 1.8 (Solutions do not need to exist for all times). Consider the
initial value problem
ẋ = tx2 , x(t0 ) = x0 ,
where x0 6= 0. Using the above procedure, we get
dx 1 1 t2 t20
= tdt , x(t0 ) = x0 =⇒ − = −
x2 x0 x 2 2
2x0
=⇒ x = .
2 + x0 (t20 − t2 )
It is easy to see that all solutions with x0 > 0 are defined only on a bounded
subinterval of R: there always exist two times t ∈ R such that the denomi-
nator 2 + x0 (t20 − t2 ) is equal to 0, and the solution converges to ∞ if these
times are approached, see Figure 1.3 for an illustration of √ the √
solution

with
t0 = 0 and x0 = 1. This solution exists on the interval − 2, 2 .

Figure 1.3. Solution to the initial value problem ẋ = tx2 , x(0) = 1.

Question: what happens for x0 < 0 and x0 = 0?

12 1. Introduction

3. Visualisations
Two different ways to visualise the solutions to ordinary differential equa-
tions are discussed in this section. The first one concerns solutions curves
of nonautonomous differential equations in the (t, x)-space, while the sec-
ond one concerns projections of solution curves of autonomous differential
equations in the x-space.

3.1. Solution portrait. We consider a function f : D ⊂ R × Rd → Rd

and the corresponding ordinary differential equation
ẋ = f (t, x) .
A solution to this equation is a differentiable function λ : I → Rd fulfilling
λ̇(t) = f (t, λ(t)) on the interval I. Then the graph of this solution, given by
the so-called solution curve
G(λ) = {(t, λ(t)) : t ∈ I} ⊂ R × Rd ,
is a differentiable curve. The derivative of this curve in the point t0 ∈ I is
d
given by dt (t, λ(t)) t=t0 = (1, λ̇(t0 )) = (1, f (t0 , λ(t0 ))).
This implies in particular that the vector field
(t, x) 7→ (1, f (t, x)) , (1.9)
defined on D, is crucial for the shape of the solution curves, in the sense
that the solution curves are tangential to the vector field (1.9).
A solution portrait is given by a visualisation of several solution curves in
the (t, x)-space, the so-called extended state space. The x-space is normally
called state space, and it is extended by the time axis.
See Figure 1.4 for a solution portrait of the differential equation ẋ = tx2
(note that we have studied this differential equation in Example 1.8).

3.2. Phase portrait. A different visualisation is possible for autonomous

differential equation, i.e. for differential equations that do not explicitly de-
pend on time t.
We first demonstrate that in this context, solutions stay solutions when we
shift them in time, which leads to some kind of redundancy when visual-
ising them via a solution portrait. This implies that a visualisation in the
state space (which is one dimension lower) is meaningful to illustrate certain
properties of solutions.
Proposition 1.9 (Translation invariance). Let λ : I → Rd be a solution to
the autonomous differential equation
ẋ = f (x) .
3. Visualisations 13

Figure 1.4. Solution portrait with several solutions to the differential

equation ẋ = tx2 (in blue), and the vector field (1.9) (in red).

Then for all τ ∈ R, the function µ : I˜ → Rd , where I˜ := {t ∈ R : t + τ ∈ I}

and
µ(t) := λ(t + τ ) for all t ∈ I˜ ,
is also a solution to this differential equation.

Proof. Since λ is a solution, we have λ̇(t) = f (λ(t)) for t ∈ I. The chain

˜ and we get
rule implies that µ̇(t) = λ̇(t + τ ) for all t ∈ I,
µ̇(t) = λ̇(t + τ ) = f (λ(t + τ )) = f (µ(t)) for all t ∈ I˜ .
This finishes the proof. 

We demonstrate visualisation using phase portraits by means of the har-

monic oscillator.
Example 1.10 (Harmonic oscillator). Consider the differential equation of
the harmonic oscillator
ẍ = −x ,
which, according to Repetition Material 1, can be re-written as
    
ẋ 0 1 x
= . (1.10)
ẏ −1 0 y
Note that this is a so-called linear differential equation, since the right hand
side of (1.10) is a matrix multiplied by the state space vector. (We have
encountered another linear differential equation already in Example 1.1.)
This differential equation is also autonomous (i.e. the matrix does not de-
pend on time t), and we will see later in Chapter 3 that all such systems
14 1. Introduction

are solvable explicitly. The solution to the initial value problem (1.10),
(x, y)(0) = (x0 , y0 ), is given by
  
cos(t) sin(t) x0
λ(t) = .
− sin(t) cos(t) y0
A visualisation of one solution (in the extended state space, and for the
initial condition (x, y)(0) = (1, 1)) is given in Figure 1.5.

Figure 1.5. Solution to the initial value problem (1.10), (x, y)(0) = (1, 1).

Note that due to the translation invariance, any translation of the above
solution is also a solution to the differential equation. Therefore, a visuali-
sation in the extended state space contains redundant information, and it is
common to visualise in the state space itself, and not in the extended state
space. This is done by a projection of all solutions to the x-space, and we
obtain the so-called phase portrait, see Figure 1.6 for the differential equa-
tion (1.10). Note that a projected solution is called an orbit or trajectory,
see Definition 2.25 below for a precise definition.
Question: What information is lost if we do such a projection?
Let us consider now an arbitrary autonomous differential equation
ẋ = f (x) ,
with a right-hand side f : D ⊂ Rd → Rd . Similarly to the vector field (1.9)
for the solution portrait, also a vector field is tangential to the trajectories
in the phase portrait. It is simply given by the projected version of vector
field (1.9), given by
x 7→ f (x) ,
defined on D, indicated by the red arrows in Figure 1.6.
3. Visualisations 15

Figure 1.6. Phase portrait of the differential equation (1.10).

 Chapter 2

Existence and
uniqueness

We have seen in Example 1.6 that it is not guaranteed to have solutions to

an initial value problem, and we have seen in Example 1.7 that solutions
to an initial value problem do not need to be unique. In this chapter, we
present a theory that guarantees existence and uniqueness for solutions to
initial value problems. Luckily, the conditions we need to impose to obtain
existence and uniqueness are rather weak and fulfilled by the vast majority
of differential equation that come from applications.

1. Picard iterates
We first want to establish a procedure to show that solutions to specific ini-
tial value problems exist. The following proposition is (maybe surprisingly)
extremely helpful for this purpose, although it looks a bit like a triviality: we
just reformulate the differential equation equivalently as an integral equation
by integrating it.
Proposition 2.1 (Reformulation as integral equation). Consider the initial
value problem
ẋ = f (t, x) , x(t0 ) = x0 , (2.1)
where f : D ⊂ R × Rd → Rd is continuous and (t0 , x0) ∈ D, and let λ : I →
Rd be a function on an interval I such that t0 ∈ I and (t, λ(t)) : t ∈ I ⊂ D.
Then the following two statements are equivalent.
(i) λ solves the initial value problem (2.1), i.e.
λ̇(t) = f (t, λ(t)) for all t ∈ I , and λ(t0 ) = x0 . (2.2)

17
18 2. Existence and uniqueness

(ii) λ is continuous, and we have

Z t
λ(t) = x0 + f (s, λ(s)) ds for all t ∈ I . (2.3)
t0

Proof. We integrate (2.2) from t0 to t and obtain (2.3), and differentiating

(2.3) with respect to t yields (2.2). 

Note that the Riemann integration in (2.3) is higher-dimensional in general,

and higher-dimensional integration of a function g : R → Rd is defined
componentwise by
R t 
Z t t0 g1 (s) ds
g(s) ds =  ..
,
 
.
t0 Rt
t0 gd (s) ds
provided that all integrals on the right hand side exist.
Why is (2.3) so helpful to obtain the solution λ to the initial value problem
(2.1)? The reason is that λ appears both on the left and right hand side of
(2.3), and we can exploit an iterative scheme to approximate it.
We first study this in a simpler setting, i.e. for an algebraic equation. Con-
sider for a > 0 the algebraic equation
a 1
a= + . (2.4)
2 a
This equation is obviously equivalent to a2 = 2, and as above for λ, the
unknown quantity a appears on the left and the right hand side of (2.4).
We define the iterative scheme
an 1
an+1 = + for all n ∈ N0 , (2.5)
2 an
with an arbitrary starting value a0 > 0. This defines a sequence {an }n∈N0 .
If we can show that this sequence converges with limit a∞ , i.e. limn→∞ an =
a∞ , then we have
a 1 a∞ 1
n
lim an+1 = lim + =⇒ a∞ = + ,
n→∞ n→∞ 2 an 2 a∞
so a∞ solves (2.4).
It is easy to see that, indeed, the sequence {an }n∈N0 converges, since it is
both monotone (for n ≥ 1) and bounded.
Question: Can you establish this rigorously? Figure 2.1 will help you with
that.
Motivated by Proposition 2.1 and the above iterative scheme for the al-
gebraic equation (2.4), we define the Picard iterates for the initial value
problem (2.1).
1. Picard iterates 19

Figure 2.1. In blue: the function g(a) = a2 + a1 , so that (2.5) reads as

an+1√= g(an ). In red: the identity function. Both graphs intersect at
a = 2.

Definition 2.2 (Picard iterates). Consider the initial value problem (2.1),
and choose an interval J that contains t0 . We define a initial function

λ0 (t) ≡ x0 for all t ∈ J ,

and inductively, the Picard iterates

Z t
λn+1 (t) := x0 + f (s, λn (s)) ds for all t ∈ J and n ∈ N0 . (2.6)
t0

Note that the interval J has to be chosen appropriately, but we do not worry
about this now, since we are interested in the general principle in the first
instance. It will follow from a question on the second problem sheet that
if this sequence is uniformly convergent with the limiting function λ∞ , we
obtain
Z t
λ∞ (t) = x0 + f (s, λ∞ (s)) ds for all t ∈ J .
t0
Rt
Note that uniform convergence is needed for limn→∞ t0 f (s, λn (s)) ds =
Rt
t0 limn→∞ f (s, λn (s)) ds. Thus, Proposition 2.1 yields that λ∞ is solution
to the integral equation (2.3) and therefore solves the initial value problem
ẋ = f (t, x), x(t0 ) = x0 .
Please remind yourself what uniform convergence means! Can you give an
example of a sequence of functions that does not converge uniformly?
We study the procedure of Picard iteration for a simple example.

Example 2.3 (Picard iterates for ẋ = ax). We would like to compute the
Picard iterates for the initial value problem

ẋ = ax , x(t0 ) = x0 ,
20 2. Existence and uniqueness

where a ∈ R is fixed. The first three iterates are

λ0 (t) = x0 ,
Z t

λ1 (t) = x0 + ax0 ds = x0 1 + a(t − t0 ) ,
t0
Zt
ax0 1 + a(s − t0 ) ds = x0 1 + a(t − t0 ) + 21 a2 (t − t0 )2 .



λ2 (t) = x0 +
t0
It is easy to prove by induction that
n
X ai (t − t0 )i
λn (t) = x0 for all n ∈ N0 and t ∈ R ,
i!
i=0
and this sequence of functions converges to the limit function
λ∞ (t) = x0 ea(t−t0 ) for all t ∈ R .
This coincides with the solution to this initial value problem we have iden-
tified in Example 1.5.

Note that we obtain global convergence of the Picard iterates for this exam-
ple (i.e. convergence for all t ∈ R, although this convergence is not uniform
on R, but uniform on any compact interval). We will see later that, in
general, under suitable but weak conditions, we obtain local convergence,
i.e. we choose a small enough compact interval J around t0 . This will es-
tablish existence and uniqueness of solutions to initial value problems on a
theoretical level. From a practical perspective, however, it is not a standard
approach to look at Picard iterates to find solutions for specific systems.

2. Lipschitz continuity
We aim at uniform convergence of the Picard iterates {λn : J → Rd }n∈N0
to a limit function λ∞ : J → Rd , where J is compact interval. How can we
express this type of uniform convergence? We will do so by considering the
space of continuous functions on a compact interval J, denoted by C 0 (J, Rd ).
As you have seen in the analysis course, this space is a complete normed
vector space when equipped with the supremum norm, and convergence in
this norm corresponds to uniform convergence (see Repetition Material 2 ).
We will see later in Section 3 that Lipschitz continuity plays an important
role in establishing uniform convergence of the Picard iterates. Firstly, the
above required uniform convergence for the Picard iterates will follow from
an application of Banach’s fixed point theorem (see Repetition Material 3 ),
and we need a Lipschitz constant less than 1 for this. Secondly, we will
see that a Lipschitz condition for Banach’s fixed point theorem holds due
to a certain Lipschitz condition on the right hand side f of the differential
2. Lipschitz continuity 21

equation under consideration (and we note that this Lipschitz constant does
not necessarily need to be less than 1).
Recall that a vector space V over the reals is an abelian group (V, +) with
an additional scalar multiplication (V, ·). In particular, for x, y ∈ V , we have
ax + by ∈ V for all a, b ∈ R.
Definition 2.4 (Normed vector space). A norm on a vector space V over
the reals is a map k · k : V → R+
0 such that

(i) kxk = 0 ⇔ x = 0 (positive definiteness),

(ii) kaxk = |a|kxk for all a ∈ R and x ∈ V (absolute homogeneity),
(iii) kx + yk ≤ kxk + kyk for all x, y ∈ V (triangle inequality).
A vector space with norm is called a normed vector space.

Examples for normed vector spaces are the finite-dimensional Euclidean

spaces Rd , where d ∈ N. In Euclidean spaces Rd , we will normally use
√
the Euclidean norm kxk := x1 + · · · + xd for x ∈ Rd .
As motivated above, we will be also interested in the infinite-dimensional
normed vector space C 0 (J, Rd ), the space of continuous functions on a com-
pact interval J.
The norm k · k of a normed vector space V naturally describes the distance
between two vectors x, y ∈ V . This distance is given by kx−yk. In fact, every
normed vector space is a metric space, where the metric d : V × V → R+ 0 is
given by
d(x, y) := kx − yk for all x, y ∈ V .
The normed vector space (V, k·k) is called complete if every Cauchy sequence
converges in V . A complete normed vector space is called a Banach space.
As mentioned above, Lipschitz continuity is crucial for rigorously establish-
ing convergence of Picard iterates.
Definition 2.5 (Continuous and Lipschitz continuous functions). Let X be
a subset of a normed vector space (V, k · kV ) and Y be a subset of a normed
vector space (W, k · kW ). Then a function f : X → Y is called
(i) continuous if for all x ∈ X and ε > 0, there exists a δ > 0 such
that
kx − x̄kV < δ =⇒ kf (x) − f (x̄)kW < ε .
(ii) Lipschitz continuous if there exists a constant K > 0 such that
kf (x) − f (x̄)kW ≤ Kkx − x̄kV for all x, x̄ ∈ X .
The constant K is called a Lipschitz constant.
22 2. Existence and uniqueness

It is easy to see that Lipschitz continuous functions are continuous (for a

given ε > 0, choose δ := Kε ), but the reverse is not true, as Example 2.6
below shows.

2.1. Lipschitz continuity and the mean value theorem. In this sub-
section, we explore in dimension one how Lipschitz continuity is related to
the mean value theorem. Consider a differentiable function f : I → R, where
I ⊂ R is an interval. Recall that the mean value theorem says that for any
x, y ∈ I, there exists an ξ between x and y such that
f (x) − f (y) = f 0 (ξ)(x − y) .
This implies
|f (x) − f (y)| = |f 0 (ξ)||x − y| , (2.7)
and it is clear if the derivative f0
is bounded on the interval I, then f is
Lipschitz continuous. In particular, this holds when I is compact and f is
continuously differentiable. We now look at some examples.
Example 2.6 (Lipschitz continuity in dimension one). We consider several
real-valued functions defined on intervals.
√
(i) The function x 7→ x, where x ∈ [0, 1], is continuous, but not
Lipschitz continuous. Note that the function is differentiable in the
open interval (0, 1) with unbounded derivative, and one can prove
this rigorously using this fact and the mean value theorem.
(ii) The function x 7→ x2 , where x ∈ R, is not Lipschitz continuous.
Note that the derivative x 7→ 2x of this function is unbounded, and
one can argue as outlined in (i).
(iii) The function x 7→ x2 , where x ∈ [0, 1], is Lipschitz continuous with
Lipschitz constant 2, since the derivative x 7→ 2x of this function
is bounded by 2 on the interval [0, 1]; this follows from the mean
value theorem as outlined above.

2.2. Lipschitz continuity and the mean value inequality. The above
example shows that a Lipschitz condition is closely connected to derivatives
(although, in general, Lipschitz continuous functions do not need to be dif-
ferentiable). In the one-dimensional context of this example, this followed
from the mean value theorem. We explore now in what sense this result can
be generalised to higher dimensions. It turns out that we only get a mean
value inequality (in contrast to an equality that holds in dimension one),
but that is good enough to obtain Lipschitz continuity.
Although our main interest are nonlinear mappings, we first look at linear
mappings. It is clear that for any matrix A ∈ Rm×n , the linear mapping
x 7→ Ax is continuous. To see that this mapping is Lipschitz continuous, we
2. Lipschitz continuity 23

introduce the so-called operator norm, which is a norm on the vector space
of all matrices in Rm×n .

Definition 2.7 (Operator norm of a matrix). For a given matrix A ∈ Rm×n ,

the operator norm of A is defined by

kAxk
kAk := sup . (2.8)
x∈Rn \{0} kxk

Note that the three norms used in (2.8) are different (unless n = m). In
addition to the operator norm (which we define), we also use the Euclidean
norm k · k on Rm and Rn here. Show as an exercise that A 7→ kAk is indeed
a norm.
Note that due to linearity of A, we have

kAxk x
kAk = sup = sup A kxk = sup kAxk .
n
x∈R \{0} kxk n
x∈R \{0} x∈Rn ,kxk=1

In particular, since the linear mapping x 7→ Ax is continuous, the above

supremum is a maximum (on the compact set Sn−1 := {x ∈ Rn : kxk = 1}),
so kAk is a finite real number. We have kAx − Ayk = kA(x − y)k ≤
kAkkx − yk, so the mapping x 7→ Ax is even Lipschitz continuous with
Lipschitz constant kAk.
The following result is the appropriate analogue of the mean value theorem
in higher dimensions. As mentioned above, we do not have an equality in
general. Can you find an example? Note that an equality can be established
in the case of m = 1.

Theorem 2.8 (Mean value inequality). Consider an open set D ⊂ Rn , and

let f : D → Rm be continuously differentiable. Then for all x, y ∈ D with
[x, y] ⊂ D, there exists a ξ ∈ [x, y] such that

kf (x) − f (y)k ≤ kf 0 (ξ)kkx − yk .

n
 x, y ∈ R , the closed
Here, for any line segment connecting x and y is given
by [x, y] := αx + (1 − α)y ∈ Rn : α ∈ [0, 1] .

Proof. Consider the function g : [0, 1] → Rm , g(α) := f (αx + (1 − α)y .

By the fundamental theorem of calculus, we have

Z 1 Z 1
0
f 0 αx + (1 − α)y (x − y) dα .


f (x) − f (y) = g(1) − g(0) = g (α) dα =
0 0
24 2. Existence and uniqueness

This implies
R1
f 0 αx + (1 − α)y (x − y) dα


kf (x) − f (y)k = 0
Lemma 2.9 R1
f 0 αx + (1 − α)y (x − y) dα


≤ 0
(2.8) R1
f 0 αx + (1 − α)y


≤ 0 dαkx − yk
0


≤ max f αx + (1 − α)y kx − yk
α∈[0,1]
0
= kf (ξ)kkx − yk

for some ξ ∈ [x, y]. Note that continuous differentiability of f was used in
the last step of this proof. 

It remains to prove the triangle-like inequality we used in the above proof.

Lemma 2.9 (Triangle-like inequality for integrals). Let I ⊂ R be an interval

and f : I → Rm be a continuous function. Then
Z t Z t
f (s) ds ≤ kf (s)k ds for all t, t0 ∈ I .
t0 t0

Proof. We first show the case t
0 < t. For n ∈ N, we look at the Riemann
t−t0 Pn−1 i
sum n
Rt i=0 f t0 + n (t − t0 ) , which in the limit n → ∞ converges to

t0 f (s) ds. The triangle inequality implies

n−1   n−1  
t − t0 X i t − t0 X i
f t0 + (t − t0 ) ≤ f t0 + (t − t0 ) ,
n n n n
i=0 i=0
Rt
and since the right hand side converges to t0 kf (s)k ds, the statement fol-
Rt Rt
lows. The case t0 ≥ t follows due to t0 f (s) ds = − t 0 f (s) ds. 

The following statement is an immediate corollary from the mean value

inequality. It follows from the fact that continuous functions attain their
(finite) maximum on compact sets.

Corollary 2.10 (Lipschitz continuity and the mean value inequality). Let
U ⊂ Rn be open and f : U → Rm be continuously differentiable. Then given
a compact and convex set C ⊂ U , the restricted function f |C : C → Rm is
Lipschitz continuous.

Note here that convexity of C means that for any two points x, y ∈ C, the
closed line segment lies in C, i.e. [x, y] ⊂ C.
3. Picard–Lindelöf theorem 25

3. Picard–Lindelöf theorem
We aim at an easily verifiable condition that the Picard iterates correspond-
ing to an initial value problem converge (at least locally in a neighbourhood
of the initial time). It turns out that a Lipschitz condition in the state space
variable x for the right hand side of a differential equation is an appropriate
condition. We first study the easiest situation where a system is globally
defined and has a global Lipschitz constant in x. After understanding the
global case, we then show that a local Lipschitz condition is sufficient for
local existence and uniqueness of solutions.
We note that the proof of the next theorem crucially makes use of Banach’s
fixed point theorem, which is applied to a complete normed vector space,
given by the space of continuous functions on a compact interval, see Repe-
tition Material 2 and 3.
Theorem 2.11 (Picard–Lindelöf theorem, global version). Consider an or-
dinary differential equation
ẋ = f (t, x) (2.9)
such that the function f : R × Rd → Rd is continuous and satisfies a global
Lipschitz condition of the form
kf (t, x) − f (t, y)k ≤ Kkx − yk for all t ∈ R and x, y ∈ Rd , (2.10)
1
where K > 0 is a constant. Define h := 2K . Then every initial value
problem (2.9), x(t0 ) = x0 , admits a unique solution λ : [t0 − h, t0 + h] → Rd .

Proof. The proof is divided in three steps and relies on the construction of
a contraction P : X → X on the Banach space X := C 0 [t0 − h, t0 + h], Rd .

It turns out that a fixed point of P , which is obtained by Banach’s fixed

point theorem, solves the above initial value problem.
Step 1. Definition of the function P : X → X.
Due to Proposition 2.1, it follows that solving the initial value problem with
a solution λ : [t0 − h, t0 + h] → Rd is equivalent to the finding a continuous
function λ : [t0 − h, t0 + h] → Rd solving the integral equation
Z t
λ(t) = x0 + f (s, λ(s)) ds for all t ∈ [t0 − h, t0 + h] .
t0
This in turn follows from finding a fixed point of the function P : X → X,
defined by
Z t
P (u)(t) := x0 + f (s, u(s)) ds for all t ∈ [t0 − h, t0 + h] .
t0
Note that P assigns to a function u ∈ X another function which we denote
by P (u). To define the function P (u), we need to evaluate its value at all
26 2. Existence and uniqueness

t ∈ [t0 − h, t0 + h], and this is what P (u)(t) means. Note that X is a Banach
space with the supremum norm
kuk∞ = sup ku(t)k for all u ∈ X (2.11)
t∈[t0 −h,t0 +h]

(see also Repetition Material 2 ). Note that the operator P is well-defined,

since the continuity of f guarantees that the integral exists. Note that P
is the function for iteratively constructing the Picard iterates; in fact, the
sequence {λn }n∈N0 from Definition 2.2 satisfies λn+1 = P (λn ) for all n ∈ N0 .
Step 2. P is a contraction.
We will prove that kP (u1 ) − P (u2 )k∞ ≤ 21 ku1 − u2 k∞ for all u1 , u2 ∈ X. To
do so, let u1 , u2 ∈ X. Then for all t ∈ [t0 − h, t0 + h], we have
Z t


kP (u1 )(t) − P (u2 )(t)k = f (s, u1 (s)) − f (s, u2 (s)) ds
t0
Lemma 2.9
Z t
≤ f (s, u1 (s)) − f (s, u2 (s)) ds
t0
(2.10)
Z t (2.11)
Z t
≤ K ku1 (s) − u2 (s)k ds ≤ K ku1 − u2 k∞ ds
t0 t0
1
≤ Khku1 − u2 k∞ = ku1 − u2 k∞ .
2
This implies by taking the supremum over all t ∈ [t0 − h, t0 + h] that
1
kP (u1 ) − P (u2 )k∞ = sup kP (u1 )(t) − P (u2 )(t)k ≤ ku1 − u2 k∞ ,
t∈[t0 −h,t0 +h] 2
which shows that P is a contraction on X.
Step 3. Application of the Banach fixed point theorem.
Since P is a contraction on a Banach space (which naturally is a complete
metric space), the Banach fixed point theorem (see Repetition Material 3 )
implies that there exists a unique fixed point λ : [t0 − h, t0 + h] → Rd . As
outlined in Step 1, with the help of Proposition 2.1, such a fixed point solves
the initial value problem under consideration, and since it is the unique fixed
point, the solution to this initial value problem is unique on the interval
[t0 − h, t0 + h]. 

Note that the Lipschitz condition (2.10) is a strong assumption, and al-
though some very important differential equations fulfill it (such as au-
tonomous linear systems, studied later in Chapter 3), it is not true for most
interesting differential equations. As an example, consider, for instance, the
differential equation ẋ = tx2 from Example 1.8. However, it turns out that
for this differential equation, solutions to any initial value problem exist lo-
cally and are unique. The reason is that the existence of a local Lipschitz
3. Picard–Lindelöf theorem 27

condition is enough. In the following definition, we distinguish between the

global Lipschitz condition (2.10) and its local version.
Definition 2.12 (Global and local Lipschitz continuity). Let D ⊂ R × Rd
be open, and consider a function f : D → Rd .
(i) f is said to be globally Lipschitz continuous (with respect to x) if
there exists a constant K > 0 such that
kf (t, x) − f (t, y)k ≤ Kkx − yk for all (t, x), (t, y) ∈ D .
(ii) f is said to be locally Lipschitz continuous (with respect to x) if
for all (t0 , x0 ) ∈ D, there exists a neighbourhood U ⊂ D of (t0 , x0 )
and a constant K > 0 such that
kf (t, x) − f (t, y)k ≤ Kkx − yk for all (t, x), (t, y) ∈ U .

The global version of the Picard–Lindelöf theorem says that, under global
Lipschitz continuity, solutions to all initial value problems exist locally and
are unique. Such a statement is true even under the weaker assumption
of local Lipschitz continuity. However, in contrast to Theorem 2.11, the
interval length 2h on which the solution exists will depend on the specific
initial value.
Theorem 2.13 (Picard–Lindelöf theorem, local version). Let D ⊂ R × Rd
be open, and consider a function f : D → Rd that is continuous and locally
Lipschitz continuous with respect to x. Consider for a fixed (t0 , x0 ) ∈ D the
initial value problem
ẋ = f (t, x) , x(t0 ) = x0 . (2.12)
Then the following two statements hold:
(i) Qualitative version. The initial value problem (2.12) has locally a
uniquely determined solution, i.e. there exists a h = h(t0 , x0 ) > 0
such that (2.12) has exactly one solution on [t0 − h, t0 + h].
(ii) Quantitative version. Consider for some τ, δ > 0 the set
τ,δ
W (t0 , x0 ) := [t0 − τ, t0 + τ ] × Bδ (x0 ), where Bδ (x0 ) := x ∈
Rd : kx − x0 k ≤ δ is the closed δ-neighbourhood of x0 . We assume
that W τ,δ (t0 , x0 ) ⊂ D, and we suppose that there exist K, M > 0
such that
kf (t, x) − f (t, y)k ≤ Kkx − yk for all (t, x), (t, y) ∈ W τ,δ (t0 , x0 ) (2.13)
and
kf (t, x)k ≤ M for all (t, x) ∈ W τ,δ (t0 , x0 ) . (2.14)
Then (2.12) has exactly one solution on [t0 − h, t0 + h], where h =
1 δ
h(t0 , x0 ) := min{τ, 2K ,M }.
28 2. Existence and uniqueness

The proof of the local version of the Picard–Lindelöf theorem is similar to

the global version, and it will be skipped here for this reason. The full proof
is given in Extra Material 1.
The following proposition shows that if the right hand side of a differential
equation is continuously differentiable, then it is locally Lipschitz continuous,
and local existence and uniqueness of solution holds.
Proposition 2.14 (Continuous differentiability and Lipschitz continuity).
Consider an open set D ⊂ R × Rd and a continuously differentiable function
f : D → Rd . Then f is locally Lipschitz continuous with respect to x, and
thus, every initial value problem involving a differential equation with right
hand side f can be solved locally uniquely.

Proof. Since D is open, for each fixed (t0 , x0 ) ∈ D, there exists a compact
and convex neighbourhood U of (t0 , x0 ). Since f is continuously differen-
tiable, and U is compact, there exists a K > 0 such that
∂f
(t, ξ) ≤ K for all (t, ξ) ∈ U . (2.15)
∂x
Due to the mean value inequality (Theorem 2.8), for any (t, x), (t, y) ∈ U ,
there exists a ξ ∈ [x, y] with
∂f (2.15)
kf (t, x) − f (t, y)k ≤ (t, ξ) kx − yk ≤ Kkx − yk .
∂x
Hence, f is locally Lipschitz continuous (with respect to x), and thus, The-
orem 2.13 implies the assertion. 

The following lemma shows that two solutions cannot cross.

Lemma 2.15 (Solutions cannot cross). Let D ⊂ R × Rd be open, and
consider a function f : D → Rd that is continuous and locally Lipschitz
continuous with respect to x. Consider two solutions of
ẋ = f (t, x) ,
given by λ : I → Rd and µ : J → Rd , where I and J are intervals. Then
either
λ(t) = µ(t) for all t ∈ I ∩ J
or
λ(t) 6= µ(t) for all t ∈ I ∩ J .

Proof. Assume to the contrary that λ(t0 ) = µ(t0 ) for some t0 ∈ I ∩ J, and
λ(t̃) 6= µ(t̃) for some t̃ ∈ I ∩ J. Without loss of generality, we assume that
t̃ > t0 . Define
t∗ := sup t > t0 : λ(t0 ) = µ(t0 ) for all t0 ∈ [t0 , t] .

4. Maximal solutions 29

Due to Theorem 2.13, we get t∗ > t0 , and due to continuity of λ and µ,

we have λ(t∗ ) = µ(t∗ ). Now both λ and µ solve the differential equation
with the initial condition x(t∗ ) = λ(t∗ ). Thus, Theorem 2.13 implies unique
solvability of this initial value problem around t∗ . This contradicts the
definition of t∗ and finishes the proof of this lemma. 

4. Maximal solutions
Let D ⊂ R × Rd be open, and consider a function f : D → Rd that is
continuous and locally Lipschitz continuous with respect to x. For a given
initial pair (t0 , x0 ) ∈ D, consider the initial value problem
ẋ = f (t, x) , x(t0 ) = x0 . (2.16)

In the last section, we have shown that we obtain a local solution to (2.16),
given by an interval of length 2h around t0 . In this section, we prove that
there exists a maximal time interval around t0 (containing [t0 − h, t0 + h])
on which the solution to (2.16) exists. We also discuss the behaviour of this
solution when time approaches the left and right points of this maximal time
interval.
Definition 2.16 (Maximal existence interval). Consider the initial value
problem (2.16). We define

I+ (t0 , x0 ) := sup t+ ≥ t0 : there exists a solution to (2.16) on [t0 , t+ ] ,

I− (t0 , x0 ) := inf t− ≤ t0 : there exists a solution to (2.16) on [t− , t0 ] ,
and the interval Imax (t0 , x0 ) := (I− (t0 , x0 ), I+ (t0 , x0 )) is called the maximal
existence interval for the initial value problem (2.16).

In the following theorem, we clarify the question of existence of a solution on

the time interval Imax (t0 , x0 ) and the boundary behaviour of this solution.
Theorem 2.17 (Existence of the maximal solution and boundary be-
haviour). There exists a maximal solution λmax : Imax (t0 , x0 ) → Rd to the
initial value problem (2.16), i.e. any other solution to this initial value prob-
lem is defined on an interval that is a subset of Imax (t0 , x0 ). The maximal
solution has the following two properties:
(i) If I+ (t0 , x0 ) is finite, then either
lim kλmax (t)k = ∞ , (2.17)
t%I+ (t0 ,x0 )

or the boundary ∂D of D is nonempty, and we have

lim dist (t, λmax (t)), ∂D = 0 . (2.18)
t%I+ (t0 ,x0 )
30 2. Existence and uniqueness

(ii) If I− (t0 , x0 ) is finite, then either

lim kλmax (t)k = ∞ ,
t&I− (t0 ,x0 )

or the boundary ∂D of D is nonempty, and we have

lim dist (t, λmax (t)), ∂D = 0 .
t&I− (t0 ,x0 )

Here, for a given set A ⊂ Rn , the function dist(·, A) : Rn → R+ 0 is defined

by
for all y ∈ Rn .

dist(y, A) := inf ky − ak : a ∈ A
Furthermore, if a solution λ : I → Rd to the initial value problem (2.16),
defined on an open interval I ⊂ R, satisfies


lim kλ(t)k = ∞ or lim dist (t, λ(t)), ∂D = 0 , (2.19)
t%sup I t%sup I

then sup I = I+ (t0 , x0 ). Furthermore, if one has

lim kλ(t)k = ∞ or lim dist (t, λ(t)), ∂D = 0 ,
t&inf I t&inf I

then inf I = I− (t0 , x0 ).

λmax (t)
x0 D

I− (t0 , x0 ) t0 I+ (t0 , x0 ) t

Figure 2.2. Illustration of possible boundary behaviour of the maximal solution.

Proof. Step 1. The existence of the maximal solution.

Choose t̄ ∈ Imax (t0 , x0 ). Due to Definition 2.16, there exists a solution
µ : I → Rd of the initial value problem (2.16) such that t̄ ∈ I and t0 ∈ I.
We note that due to Lemma 2.15, all solutions to this initial value problem
having t̄ in their domain must coincide at the time t̄, and we define λmax (t̄) =
µ(t̄), and clearly also λ̇max (t̄) = f (t̄, λmax (t̄)). Since t̄ was chosen arbitrarily,
this defines the maximal solution in the open interval Imax (t0 , x0 ). We also
4. Maximal solutions 31

note that the maximal solution cannot be defined on the endpoints of the
(open) interval Imax (t0 , x0 ), since via the local version of the Picard–Lindelöf
theorem, we would be able to extend this solution beyond either I− (t0 , x0 )
or I+ (t0 , x0 ), and thus contradicting Definition 2.16.
Step 2. Proof of (i).
Assume that both (2.17) and (2.18) do not hold. Hence, there exists an
M > 0 and a sequence {tn }n∈N with limn→∞ tn = I+ (t0 , x0 ) and

1
kλmax (tn )k ≤ M and dist (tn , λmax (tn )), ∂D ≥ for all n ∈ N .
M
(2.20)
The sequence {(tn , λmax (tn ))}n∈N is bounded, and thus, there exists a con-
vergent subsequence (tnk , λmax (tnk )) → (t∗ , x∗ ) ∈ D as k → ∞ (we have
t∗ = I+ (t0 , x0 )). Note that (t∗ , x∗ ) ∈ D, because of the second part of (2.20).
Due to an exercise on the current problem sheet, in a neighbourhood W of
(t∗ , x∗ ), there exists a h = h(W ) such that all initial value problems with
initial values (t0 , x0 ) ∈ W

, there exists a solution on the interval [t0 −h, t0 +h].
Due to tnk , λmax (tnk ) → (t , x∗ ) ∈ D, there exists an N ∈ N such that
∗

(tnk , λmax (tnk )) ∈ W for all k ≥ N . This implies that tnk + h > t∗ for large
k ∈ N, and thus, the solution can be extended beyond t∗ , which contradicts
the maximality of the solution λmax . The proof of (ii) is analogous.
Step 3. The converse of (i) and (ii).
Assume that a solution λ : I → Rd satisfies (2.19), and assume that the
existence interval can be extended beyond sup I. This implies that the
solution is defined at the time sup I, which means that limt→sup I λ(t) =
λ(sup I) and (sup I, λ(sup I)) ∈ D. This is a contradiction to (2.19) (note
that D is open). Hence the solution can not be extended beyond sup I. A
similar argument can be made for extension below inf I. 

Example 2.18. For a fixed parameter α > 0, we consider the autonomous

differential equation

ẋ = xα ,

the right hand side of which is defined for all x > 0 (and formally for
all t ∈ R), so the domain is given by D = R × R+ . We consider the
initial condition x(0) = 1. Using separation of variables (see description
before Example 1.8), we can compute the maximal solution depending on
the parameter α > 0. It is clear that for the linear case α = 1, we have
λ1,max (t) = et , and for α 6= 1, we get

1
λα,max (t) = 1 + (1 − α)t 1−α .
32 2. Existence and uniqueness

Note that the maximal existence intervals, on which these solutions exist,
are given by
1


 α−1 , ∞ : α ∈ (0, 1) ,
Iα,max (0, 1) = (−∞, ∞)
: α = 1 ,
1
− ∞, α−1 : α ∈ (1, ∞) ,


and depend on α (see also Figure 2.3).

Figure 2.3. λα,max for α ∈ (0, 1) (left), and λα,max for α ∈ (1, ∞) (right).

Both situations described in Theorem 2.17 (explosion or convergence against

boundary of D) can occur. For α ∈ (0, 1), Iα,max is bounded below, and the
solution convergence to the boundary of D (which is the t-axis R × {0}). On
the other hand, for α ∈ (1, ∞), Iα,max is bounded above, and we see that
the solution converges to infinity when approaching the upper boundary
of Iα,max .

5. General solutions and flows

After studying solutions for specific initial value problems, we introduce
notion of a general solution and a flow that comprise all solutions of a
differential equation (provided conditions for local existence and uniqueness
are satisfied).

5.1. General solutions. We first will deal with nonautonomous differ-

ential equations before providing a simpler approach for the autonomous
special case.
Consider an open subset D ⊂ R × Rd and a continuous and locally Lipschitz
continuous right hand side f : D → Rd of the differential equation
ẋ = f (t, x) . (2.21)

The notion of a general solution is explained in the following definition.

5. General solutions and flows 33

Definition 2.19 (General solution to a nonautonomous differential equa-

tion). Consider the nonautonomous differential equation (2.21), and define

Ω := (t, t0 , x0 ) ∈ R1+1+d : (t0 , x0 ) ∈ D and t ∈ Imax (t0 , x0 ) .



Then the function λ : Ω → Rd , defined by

λ(t, t0 , x0 ) := λmax (t, t0 , x0 ) ,

where λmax is defined as in Theorem 2.17, is called the general solution of

(2.21).

Note that in the setting of Theorem 2.17, the initial pair (t0 , x0 ) was fixed,
but here we vary it to combine the maximal solutions to all initial value
problems in one notion, so we indicate the dependence of λmax (t, t0 , x0 ) on
this initial pair here.
This means that a general solution is a function that brings together all
maximal solutions of initial value problems. Using the notion of a general
solution, the solution identity then reads as
∂λ

(t, t0 , x0 ) = f t, λ(t, t0 , x0 ) for all (t, t0 , x0 ) ∈ Ω , (2.22)
∂t
where a partial derivative has to be used, since the time t is not the only
argument in the general solution.
We first study a simple example.

Example 2.20. We consider the ordinary differential equation (1.1) from

Example 1.1, given by
ẋ = ax ,
where a ∈ R. We have seen already in Example 1.1 that each initial condition
x(t0 ) = x0 leads to a unique solution

λmax (t) = x0 ea(t−t0 ) for all t ∈ R ,

and we have Imax (t0 , x0 ) = R. Hence the general solution is given by

λ(t, t0 , x0 ) = x0 ea(t−t0 ) for all (t, t0 , x0 ) ∈ Ω ,

where the domain Ω is given by Ω = R × R × Rd .

The general solution has the following fundamental and important proper-
ties.

Proposition 2.21 (Properties of the general solution). Consider the nonau-

tonomous differential equation (2.21), and let (t0 , x0 ) ∈ D. Then for all
34 2. Existence and uniqueness

s ∈ Imax (t0 , x0 ), we have

Imax (s, λ(s, t0 , x0 )) = Imax (t0 , x0 ) , (2.23)
λ(t0 , t0 , x0 ) = x0 , (2.24)


λ t, s, λ(s, t0 , x0 ) = λ(t, t0 , x0 ) for all t ∈ Imax (t0 , x0 ) . (2.25)
The identity (2.24) is called initial value property, while (2.25) is called
cocycle property.

Question: Can you explain the identities (2.24) and (2.25) in words?

Proof. The identity (2.24) is clear. Due to s ∈ Imax (t0 , x0 ), we have

(s, λ(s, t0 , x0 )) ∈ D, and we consider this as initial pair. Then the func-
tions µ1 (t) := λ(t, t0 , x0 ) and µ2 (t) := λ(t, s, λ(s, t0 , x0 )) are maximal so-
lutions to the initial value pairs (t0 , x0 ) and (s, λ(s, t0 , x0 )), and obviously,
µ1 (s) = µ2 (s) = λ(s, t0 , x0 ). Due to uniqueness of solutions and the fact
that both µ1 and µ2 are maximal solutions, we get µ1 = µ2 . This implies
(2.23) and (2.25). 

5.2. Flows. A simpler situation is given when the differential equation un-
der consideration is autonomous. We consider an open subset D ⊂ Rd and
a locally Lipschitz continuous right hand side f : D → Rd of the differential
equation
ẋ = f (x) . (2.26)
The simplification is due to the fact that, because of the translation invari-
ance of autonomous differential equations, as proved in Proposition 1.9, the
general solution to an autonomous differential equation does not depend on
the actual time and initial time separately, but only on the elapsed time
(which is the difference between these two times). More precisely, let λ
denote the general solution to (2.26). Then we get
λ(t, t0 , x0 ) = λ(t − t0 , 0, x0 ) for all t ∈ Imax (t0 , x0 ) ,
as well as

Imax (t0 , x0 ) = Imax (0, x0 ) + t0 := t0 + t : t ∈ Imax (0, x0 ) , (2.27)
and this motivates the definition of a flow of the autonomous differential
equation (2.26).
Definition 2.22 (Flow of an autonomous differential equation). Consider
the autonomous differential equation (2.26), and define for any initial value
x0 ∈ D,
Jmax (x0 ) := Imax (0, x0 ) (2.28)
and
ϕ(t, x0 ) = λ(t, 0, x0 ) for all t ∈ Jmax (x0 ) .
5. General solutions and flows 35

The function (t, x0 ) 7→ ϕ(t, x0 ) is called the flow of the autonomous differ-
ential equation (2.26).

Analogously to (2.22) in the nonautonomous case and for the general so-
lution, using the flow of an autonomous differential equation, the solution
identity reads as
∂ϕ

(t, x0 ) = f ϕ(t, x0 ) for all x0 ∈ D and t ∈ Jmax (x0 ) .
∂t
We first consider a simple example.

Example 2.23. The differential equation

ẋ = ax
considered in the recent Example 2.20 to illustrate the notion of a general
solution is an autonomous differential equation, so we can consider this in
the setting of flows also. As explained above, the general solution of this
differential equation is given by

λ(t, t0 , x0 ) = x0 ea(t−t0 ) for all (t, t0 , x0 ) ∈ R1+1+d .

Hence, the flow of this differential equation is given by the function
ϕ(t, x0 ) = λ(t, 0, x0 ) = x0 eat .

We note that normally flows are denoted as ϕ(t, x) rather than ϕ(t, x0 ), and
we will use the more usual notation in the following.
Flows have the following properties, which correspond naturally to the prop-
erties of general solutions studied in Proposition 2.21.

Proposition 2.24 (Properties of the flow). Let ϕ be the flow of the au-
tonomous differential equation (2.26). Then for any x ∈ D, the following
statements hold.
Jmax (ϕ(t, x)) = Jmax (x) − t for all t ∈ Jmax (x) , (2.29)
ϕ(0, x) = x , (2.30)
ϕ(t, ϕ(s, x)) = ϕ(t + s, x) for all t, s with s, t + s ∈ Jmax (x) , (2.31)
ϕ(−t, ϕ(t, x)) = x for all t ∈ Jmax (x) . (2.32)

The identity (2.30) is called initial value condition, while (2.31) is called
group property.

Proof. Let λ denote the general solution to (2.26), and let Imax denote the
maximal existence intervals corresponding to λ.
36 2. Existence and uniqueness

(2.29): We have
(2.28)
Jmax (ϕ(t, x)) = Jmax (λ(t, 0, x)) = Imax (0, λ(t, 0, x))
(2.23)
= Imax (0, λ(0, −t, x)) = Imax (−t, x)
(2.27)
= Imax (0, x) − t = Jmax (x) − t ,

where we used translation invariance in the third equality.

(2.30): From (2.24), it follows that ϕ(0, x) = λ(0, 0, x) = x.
(2.31): We have

ϕ(t, ϕ(s, x)) = λ(t, 0, λ(s, 0, x)) = λ(t + s, s, λ(s, 0, x))

(2.25)
= λ(t + s, 0, x) = ϕ(t + s, x) .

(2.31) (2.30)
(2.32): This follows from ϕ(−t, ϕ(t, x)) = ϕ(−t + t, x) = x. 

As explained in Subsection 3.2 of Chapter 1, for visualising autonomous

systems, we project solution curves from the extended state space to objects
(called orbits or trajectories) in the state space. By doing so, we obtain the
phase portrait of the differential equation.

Definition 2.25 (Orbits or trajectories). Let ϕ be the flow of the au-

tonomous differential equation (2.26). For all x ∈ D, we call the set

O(x) := ϕ(t, x) ∈ D : t ∈ Jmax (x)

the orbit (or trajectory) through x. In addition, we call O+ (x) := ϕ(t, x) ∈


+
 : t ∈ Jmax (x) ∩ R0 the positive
D half-orbit through x, and O− (x) :=
ϕ(t, x) ∈ D : t ∈ Jmax (x) ∩ R−0 the negative half-orbit through x.

The geometric picture is that the domain D of the right hand side f is
partitioned into orbits of ϕ. There are essentially three different types of
orbits O(x) for x ∈ D.
(i) O(x) is a singleton. This implies f (x) = 0 (see also Proposi-
tion 1.3), and Jmax (x) = R. The point x is called equilibrium.
(ii) O(x) is a closed curve, i.e. there exists t > 0 such that ϕ(t, x) = x,
but f (x) 6= 0. This implies Jmax (x) = R. The point x is called
periodic, and O(x) is called periodic orbit.
(iii) O(x) is not a closed curve, i.e. the function t 7→ ϕ(t, x) is injective
on Jmax (x).
We identify all three types of points in the following example.
5. General solutions and flows 37

Example 2.26. Consider the autonomous two-dimensional differential

equation
ẋ = y + x(1 − x2 − y 2 ) ,
ẏ = −x + y(1 − x2 − y 2 ) ,
for which it is possible (via polar coordinates leading to the differential
equation ṙ = r(1 − r2 ), φ̇ = −1) to obtain the expression
 
1 x cos(t) + y sin(t)
ϕ(t, x, y) = q
y cos(t) − x sin(t)
x2 + y 2 + 1 − x2 − y 2 e−2t
for the flow of the system, see Figure 2.4 for the phase portrait.

Figure 2.4. Phase portrait of the differential equation from Example 2.26.

It is easy to see that all three trajectory types are present in this example.
The point (0, 0) is the only equilibrium, and there exists a periodic orbit
O(0, 1). All other orbits are not closed curves, and they converge in forward
time to the periodic orbit. It is clear that it makes sense to call the periodic
orbit stable, while the equilibrium (0, 0) is called unstable. We will formally
introduce these so-called notions of stability later in this course.
The situation is a bit simpler in the one-dimensional case, where periodic
orbits cannot occur. The proof of the following proposition is left as an
exercise.
Proposition 2.27 (Orbits of one-dimensional differential equations). Con-
sider the autonomous differential equation (2.26), where d = 1. Then all
solutions are monotone, and there do not exist periodic orbits. This means
that a trajectory is either a equilibrium or a non-closed curve.
38 2. Existence and uniqueness

The following remark addresses smoothness properties of general solutions

and flows.
Remark 2.28 (Continuity and differentiability of general solutions and
flows). Note that the general solution (as a function of three variables) and
the flow (as a function of two variables) is continuous. The proof of this fact
is lengthy and not very insightful, so we do not cover this in this course. In
addition, if we have more regularity of the right hand side, such as continuous
differentiability, one can even prove that the general solution and the flow
are also continuously differentiable. These are important results, since in
many problems, one is interested in variation of the initial conditions (and
also of parameters of the system, for which similar results hold).
 Chapter 3

Linear systems

Although the most interesting differential equations in applications are non-

linear, the class of linear systems is very important, because they allow to
describe a first-order approximation of the behaviour of solutions close to a
given reference solution. Consider such a solution µ : I → Rd of a differential
equation
ẋ = f (t, x) . (3.1)
Then the linearisation along the reference solution µ is given by
∂f
ẋ = (t, µ(t)) x
|∂x {z }
=:A(t)∈Rd×d

where the right hand side is a (time-dependent) linear function, the matrix
of which is given by the derivative of the (differentiable) function f evaluated
along the solution. In general, i.e. if the differential equation (3.1) is nonau-
tonomous, or the solution µ is not constant, then the matrix A depends on
time t, which makes the situation quite complicated. This is due to the fact
that nonautonomous linear systems ẋ = A(t)x are not solvable in general,
while there exists an explicit representation for the flow of an autonomous
linear system ẋ = Ax. For this reason, we focus on the autonomous case,
which is obtained when we linearise an autonomous differential equation
ẋ = f (x) in an equilibrium x∗ , leading to the linear system ẋ = f 0 (x∗ )x.

1. Matrix exponential function

We consider the linear differential equation
ẋ = Ax , (3.2)

39
40 3. Linear systems

where A ∈ Rd×d . Since kAx − Ayk = kA(x − y)k ≤ kAkkx − yk, this system
is globally Lipschitz continuous with Lipschitz constant kAk, and due to
the global version of the Picard–Lindelöf theorem (and an exercise from a
problem sheet), solutions to every initial value problem exist on R and are
unique, and this generates a globally defined flow ϕ : R × Rd → Rd .
Due to the Picard–Lindelöf theorem (Theorem 2.11), solutions can be ob-
tained locally by convergence of Picard iterates (λn )n∈N0 on some interval
J ⊂ R as defined in (2.6). We fix an initial value x0 ∈ Rd and are inter-
ested in how the solution corresponding to the initial condition x(0) = x0
looks in a neighbourhood J of 0. Note that the following analysis generalises
Example 2.3.
We define the initial function λ0 (t) := x0 for all t ∈ J, and the iterates as
defined in (2.6) read as
Z t
λn+1 (t) = P (λn )(t) = x0 + Aλn (s) ds .
0

Then for all t ∈ J, we have

Z t
λ1 (t) = x0 + Aλ0 (s) ds = x0 + tAx0 ,
0
Z t
t2
λ2 (t) = x0 + Aλ1 (s) ds = x0 + tAx0 + A2 x0 .
0 2
By induction, we obtain
n
t2 tn X tk Ak
λn (t) = x0 + tAx0 + A2 x0 + · · · + An x0 = x0 .
2 n! k!
k=0

So the solution to the initial value problem (3.2), x(0) = x0 , is given locally
around t = 0 by
∞ k k
X t A
λ∞ (t) = ϕ(t, x0 ) = eAt x0 , with eAt := with A0 = Idd . (3.3)
k!
k=0

The proof of the Picard–Lindelöf Theorem 2.11 shows that this infinite sum
exists (i.e. the series converges) whenever |t| ≤ h for some h > 0. We will
demonstrate later that it exists for all t ∈ R.
For a given matrix A ∈ Rd×d , the function t 7→ eAt is called the matrix
exponential function.
To prove that the matrix exponential function exists for all t ∈ R, we need
the following lemma (recall the definition of the operator norm of a matrix
from Definition 2.7).
1. Matrix exponential function 41

Lemma 3.1 (Sub-multiplicativity of the matrix norm). For two matrices

B, C ∈ Rd×d , we have
kBCk ≤ kBkkCk . (3.4)

Proof. For any 0 6= x ∈ Rd , we have

kBCxk ≤ kBkkCxk ≤ kBkkCkkxk ,
kBCxk
and thus, kBkkCk is an upper bound for kxk , which finishes the proof. 

Proposition 3.2 (Existence of the matrix exponential). Consider a matrix

B ∈ Rd×d . Then its matrix exponential
∞
X 1 k
eB := B
k!
k=0

exists and is a matrix in Rd×d .

Proof. Due to Lemma 3.1, we have kB k k ≤ kBkk for all k ∈ N. This

implies
∞ ∞
X kBkk X kB k k
ekBk = ≥
k! k!
k=0 k=0
kB k k
Hence, using the comparison test, we see that the series ∞
P
k=0 k! is con-
vergent. Define the sequences
n n
X kB k k X Bk
an := and bn := for all n ∈ N ,
k! k!
k=0 k=0

and note that due to the above observation, the sequence {an }n∈N converges.
For n > m, the triangle inequality implies that
n m
X 1 k X 1 k
kbn − bm k = B ≤ B = |an − am | .
k! k!
k=m+1 k=n+1

Since {an }n∈N is a Cauchy sequence, this inequality shows that {bn }n∈N is
a Cauchy sequence as well, and it is convergent in Rd×d . 

We have seen that in (3.3) that the flow of (3.2) (ẋ = Ax) is given locally
around t = 0 by ϕ(t, x) = eAt x. We show now that this holds for all t ∈ R.
Theorem 3.3 (The flow of an autonomous linear differential equation).
Consider the autonomous linear differential equation (3.2) with coefficient
matrix A ∈ Rd×d . Then the flow ϕ : R × Rd → Rd generated by this
differential equation is given by
ϕ(t, x) = eAt x for all t ∈ R .
42 3. Linear systems

Proof. Step 1. We show that eA(t+s) = eAt eAs for any t, s ∈ R.

We have
∞ k k X ∞ ` ` ∞ X ∞ k ` k+`
X t A sA X t sA
eAt eAs = = ,
k! `! k!`!
k=0 `=0 k=0 `=0

and setting n = k + `, so that k = n − `, we get

∞ X
∞ k n−k n ∞ n
X t s A X An X n!
eAt eAs = = tk sn−k
k!(n − k)! n! k!(n − k)!
k=0 n=k n=0 k=0
∞ n
X A (t + s)n
= = eA(t+s) .
n!
n=0

Step 2. We show that ϕ(t, x0 ) = eAt x0 for all t ∈ R.

We know already from (3.3) that there exists h > 0 such that

ϕ(t, x0 ) = eAt x0 for all t ∈ [−h, h] .

t
Let t ∈ R and choose N ∈ N with N ∈ [−h, h]. Then the group property
(2.31) of ϕ implies

ϕ(t, x0 ) = ϕ( Nt , ϕ( Nt , . . . , ϕ( Nt , x0 ) . . . ))
| {z }
N times
N
Y t Step 1 t
= e N A x0 = eAN N x0 = eAt x0 .
i=1

This finishes the proof. 

The matrix exponential has the following important properties. The proof
is left as an exercise.

Proposition 3.4 (Properties of the matrix exponential). Consider matrices

B, C, T ∈ Rd×d such that T is invertible. Then the following statements hold.
(i) If C = T −1 BT , then eC = T −1 eB T .

−1
(ii) e−B = eB .
(iii) If BC = CB, then eB+C = eB eC .
(iv) If B is a block diagonal matrix B = diag(B1 , . . . , Bp ) with matrices
B1 , . . . , Bp , then eB = diag(eB1 , . . . , eBp ).

Question: Why does (iii) of Proposition 3.4 not hold in general? Can you
find matrices B, C ∈ Rd×d such hat eB+C 6= eB eC ?
2. Planar linear systems 43

2. Planar linear systems

We look first at the two-dimensional case, and we would like to have an
explicit representation of the flow eAt for an autonomous linear differential
equation

ẋ = Ax ,

where A ∈ R2×2 . In addition, we would like to understand what different

types of phase portraits we can get.
We transform the matrix A in Jordan normal form J = T −1 AT , where the
transformation matrix T ∈ R2×2 is invertible. Using Proposition 3.4 (i), we
get eAt = T eJt T −1 , and it remains to understand eJt for two-dimensional
Jordan normal forms J.
There are four different cases:
a0


(C1) A has two different real eigenvalues a, b ∈ R: J = 0 b .
(C2) A has a double real eigenvalue

a ∈ R with two linearly independent
eigenvectors: J = a0 a0 .
(C3) A has a
double real eigenvalue a ∈ R with only one eigenvector:
J = a0 a1 .
a b .


(C4) A has complex pair a ± ib of eigenvalues with b 6= 0: J = −b a

We first study the situation for all cases when the matrix is not singular,
i.e. 0 is not an eigenvalue of A.

I. The matrix A is not singular.

I.(C1) J = a0 0b , where a, b ∈ R \ {0} and a 6= b.
We get
 at 
e 0
eJt = for all t ∈ R .
0 ebt

This means that for any(x0 , y0 ) ∈ R2
, the trajectory is given by O(x0 , y0 ) =
{eJt xy00 : t ∈ R} = x0 eat , y0 ebt : t ∈ R . We see that apart from

the equilibrium (0, 0), the four half axes are trajectories. Outside of these
trajectories, we obtain the representation
n
b  x o
O(x0 , y0 ) = x, y0 x
x0
a
∈ R2 : >0 .
x0

We obtain the following phase portraits, depending on the order and sign of
a and b.
44 3. Linear systems

a<b<0
stable knot
with two tangents

a<0<b
saddle

0<a<b
unstable knot
with two tangents
2. Planar linear systems 45

a0


I.(C2) J = 0a , where a ∈ R \ {0}.
We get

eat 0
 
Jt
e = for all t ∈ R .
0 eat

This means that for any(x0 , y0 ) ∈ R2
, the trajectory is given by O(x0 , y0 ) =
Jt x0

{e y0 : t ∈ R} = x0 eat , y0 eat : t ∈ R . We see that apart from
the equilibrium (0, 0), the four half axes are trajectories. Outside of these
trajectories, we obtain the representation

n x o
x, x xy00 ∈ R2 :


O(x0 , y0 ) = >0 .
x0

We obtain the following phase portraits, depending on the sign of a.

a<0
stable knot
with many tangents

a>0
unstable knot
with many tangents
46 3. Linear systems

a1


I.(C3) J = 0a , where a ∈ R \ {0}.
It follows from an exercise that
 at
teat

Jt e
e = for all t ∈ R .
0 eat

This means that for any (x0 , y0 ) ∈ R2 , the trajectory is given by O(x0 , y0 ) =
Jt x
 at at at


{e y0 : t ∈ R} = x0 e + y0 te , y0 e
0
: t ∈ R . We see that apart
from the equilibrium (0, 0), the two half x-axes are trajectories. Outside of
these trajectories, we obtain the representation

n
y
 y o
O(x0 , y0 ) = x0
y0 y + a ln yy0 , y ∈ R2 : >0 .
y0

We obtain the following phase portraits, depending on the sign of a.

a<0
stable knot
with one tangent

a>0
unstable knot
with one tangent
2. Planar linear systems 47

a b


I.(C4) J = −b a , where b ∈ R \ {0}.
It follows from an exercise that

 
cos(bt) sin(bt)
eJt = eat for all t ∈ R . (3.5)
− sin(bt) cos(bt)

This means that for any (x0 , y0 ) ∈ R2 , the trajectory is given by


  
x0 Jt
O(x0 , y0 ) = e :t∈R
y0
   
at x0 cos(bt) + y0 sin(bt)
= e :t∈R .
y0 cos(bt) − x0 sin(bt)

We obtain the following phase portraits, depending on the signs of a and b

a < 0, b > 0
stable focus

a < 0, b < 0
stable focus
48 3. Linear systems

a = 0, b > 0
centre

a = 0, b < 0
centre

a > 0, b > 0
unstable focus
2. Planar linear systems 49

a > 0, b < 0
unstable focus

II. The matrix A is singular.

II.(C1) J = a0 00 , where a ∈ R \ {0}.
We get

 at 
Jt e 0
e = for all t ∈ R .
0 1

This means that for any (x0 , y0 ) ∈ R2 , the trajectory is given by O(x0 , y0 ) =
x
{eJt y00 : t ∈ R} = (eat x0 , y0 ) : t ∈ R .



We obtain the following phase portraits, depending on the sign of a.

a<0
50 3. Linear systems

a>0

00


II.(C2) J = 00 .
In this trivial case, the whole state space consists of equilibria.


II.(C3) J = 00 10 .
We get

 
1 t
eJt = for all t ∈ R .
0 1

2
This means
Jt x0
that for any (x0 , y0 ) ∈ R , the trajectory is given by O(x0 , y0 ) =
e y0 = (x0 +y0 t, y0 ) : t ∈ R . Note that the x-axis consists of equilibria.
The phase portrait is given as follows.

II.(C4) This case does not exist, since a two-dimensional matrix without
real eigenvalues cannot have eigenvalue zero.
2. Planar linear systems 51

We note that while we know now everything about two-dimensional phase

portraits for differential equations in Jordan normal form, the question re-
mains how phase portraits look like when the system is not in Jordan nor-
mal form. Firstly, note that, given a general two-dimensional linear system
ẋ = Ax, we transform the matrix A ∈ R2×2 in Jordan normal form via
J = T −1 AT , where the transformation matrix T ∈ R2×2 is invertible. Using
Proposition 3.4 (i), we have eAt = T eJt T −1 , but what does this mean for
the phase portrait? In an exercise, you show that the phase portrait of the
original system ẋ = Ax is the result of applying the linear transformation T
to the phase portrait of system in Jordan normal form.
We close this section on two-dimensional linear systems with a remark con-
cerning the importance of the eigenvalues of A with regard to both expo-
nential growth and rotation.
Remark 3.5 (Meaning of real and imaginary part of the eigenvalues of A).
It turns out that in all two-dimensional examples, the eigenvalues are char-
acteristic for both the strength of exponential growth and rotation of the
solutions of the corresponding linear differential equation ẋ = Ax.
(i) Rate of exponential growth. The real part of the eigenvalues of
A determine rates the exponential growth behaviour of solutions
λ(t) = eAt (x0 , y0 )> , where t ∈ R. Note that if a function µ : R →
R \ {0} is growing exponentially, for instance, if µ(t) = eat for some
a ∈ R, then the exponential growth rate can be obtained as
ln µ(t) ln eat
lim = = a.
t→∞ t t
This is still true if the exponential growth is not purely expo-
nential, for instance, an easy calculation shows that the function
µ(t) = tn eat , where n ∈ N, will give the same exponential growth
rate a (note that we see such a function in Case I.(C3)).
This motivates the definition of a Lyapunov exponent correspond-
ing to the above solution λ, where we assume that the initial con-
dition (x0 , y0 ) 6= (0, 0).
ln kλ(t)k
σLyap (λ) := lim ,
t→∞ t
provided the limit exists. As can be checked easily, it turns out in
our all of our two-dimensional examples that the limit exists and
is equal the real part of one of the eigenvalues of the matrix A.
Note also that a solution exponentially decays if σLyap (λ) < 0 and
exponentially increases if σLyap (λ) > 0.
(ii) Rate of rotation. As can be seen in the examples, solutions only ro-
tate in Case I.(C4). This is the only case where the eigenvalues are
52 3. Linear systems

not real, and it can be seen that the absolute value of b determines
the speed of rotation, while the sign of b determines the orientation
of rotation: it goes in clockwise if b > 0, and anti-clockwise if b < 0.
It will turn out in a moment that these two-dimensional observations are
true in higher-dimensions as well, see the formulas for the matrix exponential
in Proposition 3.8 below.

3. Jordan normal form

We now aim at studying the higher-dimensional case, and we would like
to have an explicit representation of the flow eAt for an autonomous linear
differential equation
ẋ = Ax ,
where A ∈ Rd×d . As before, the Jordan normal form plays a crucial role,
since, as we will see later, it is easy to compute the matrix exponential
function corresponding to the Jordan normal form of A.
We first look at the complex Jordan normal form of A. The transformation
matrix (called T below) leading to this Jordan normal form may be complex
in this case when there are complex eigenvalues. When looking at the real
Jordan normal form in the next step, we will be able to remove the complex
entries in the transformation matrix, in order to arrive at the real Jordan
normal form.

Theorem 3.6 (Complex Jordan normal form). Consider a matrix A ∈

Rd×d . Then there exists a matrix T ∈ Cd×d so that under a basis transfor-
mation with the matrix T , we obtain the complex Jordan normal form
 
J1 0
J := T −1 AT = 
 .. ,

.
0 Jp

with the so-called Jordan blocks

 
ρj 1 0 0
 0 ρj 1 0
 
Jj = 
 .. .. 
for all j ∈ {1, . . . , p} , (3.6)
 . . 

0 ρj 1 
0 0 0 ρj

where the ρj , j ∈ {1, . . . , p}, are complex eigenvalues of the matrix A (some
of which may be the same).
3. Jordan normal form 53

Note that if Jj is a 1 × 1 matrix, then Jj = (ρj ), and if Jj is a 2 × 2 matrix,

then  
ρj 1
Jj = .
0 ρj
You have encountered the complex Jordan normal form already, and we will
not prove this theorem here, and we refer to Repetition Material 4 for a
more details how to compute the complex Jordan normal form.
We cover in more detail now the real Jordan normal form.

Theorem 3.7 (Real Jordan normal form). Consider a matrix A ∈ Rd×d .

Then there exists a matrix T ∈ Rd×d so under basis transformation with the
matrix T , we obtain the real Jordan normal form
 
J1 0
J := T −1 AT = 
 .. ,

.
0 Jp
where the Jordan blocks Jj are either as in Theorem 3.6, i.e. given by (3.6),
in case the eigenvalue ρj is real, or, in case the eigenvalue ρj is complex,
 
Cj Id2 0 0
 0 Cj Id2 0 
 
Jj = 
 . .
.. .. ,

(3.7)
 
0 Cj Id2 
0 0 0 Cj
   
aj bj 1 0
where Cj = with ρj = aj + ibj , and Id2 = .
−bj aj 0 1

Note that if Jj is a 2 × 2 matrix, then

 
aj bj
Jj = .
−bj aj

We also do not prove this theorem, but would like to understand how to
construct the matrix T to obtain the real Jordan form, given by the matrix J.
Question 1. Why does the real Jordan normal form consist of blocks of the
form Cj ? We demonstrate this here for the two-dimensional case. Assume
that we are in the situation of (C4) in Section 2, i.e. the matrix A ∈ R2×2
has the complex eigenvalues ρ = a±ib with b 6= 0. Then the complex Jordan
normal form reads as
 
a + ib 0
Jcomplex = .
0 a − ib
54 3. Linear systems

Assume that u + iv with u, v ∈ R2 is complex eigenvector for the eigenvalue

a + ib, i.e.
A(u + iv) = (a + ib)(u + iv) . (3.8)
We note that the vectors u and v are linearly independent in R2 . To see
this, assume they are linearly dependent, i.e. v = γu for some real value
γ 6= 0. Hence u + iv = (1 + iγ)u. Then (3.8) implies that (1 + iγ)Au =
(a + ib)(1 + iγ)u, which leads to
Au = (a + ib)u ,
and we get a contradiction, since the left hand side is vector in R2 , but not
the right hand side. We now look at action of the linear mapping on the
basis vectors u and v of R2 . Comparing real part and imaginary part of
(3.8) implies that
Au = au − bv
Av = bu + av ,
a b


which gives the desired Jordan normal form −b a corresponding to the
basis {u, v}. Note that the sign of b is not unique in the Jordan form. In
fact, taking −v instead of v as basis vector changes the sign of b. Note also
that taking complex conjugation on both sides of (3.8) implies that u − iv
is eigenvector for the eigenvalue a − ib.
Question 2. How do we compute the matrix T ∈ Rd×d ? We assume that
we know the complex Jordan normal form already, with the transformation
matrix. Note that the non-real eigenvalues appear in complex conjugate
pairs aj + ibj and aj − ibj with bj > 0. Now ignore all blocks with bj < 0
and replace the bases of complex generalised eigenvectors and eigenvectors
corresponding to aj + ibj with bj > 0, given by w1j , . . . , wdj j , by the column
vectors
Re w1j , Im w1j , . . . , Re wdj j , Im wdj j .
Then we will get a real matrix T ∈ Rd×d that guarantees transformation in
real Jordan normal form according to Theorem 3.7.

4. Explicit representation of the matrix exponential function

Consider a matrix A ∈ Rd×d . We explain in this section how to obtain an
explicit representation of eAt for all t ∈ R.
We assume that the invertible matrix T ∈ Rd×d transforms A into the real
Jordan normal form
 
J1 0
J = T −1 AT = 
 .. ,

.
0 Jp
4. Explicit representation of the matrix exponential function 55

as in Theorem 3.7, where the Jordan blocks Jj are given by either (3.6)
(corresponding to real eigenvalues ρi ) or (3.7) (corresponding to complex
eigenvalues ρj = aj + ibj ). Due to Proposition 3.4 (i) and (iv), we get
 Jt 
e 1 0
eAt = T eJt T −1 = T 
 ..  −1
T ,
.
0 J
e p t

so it remains to find an explicit representation for the matrix exponentials

t 7→ eJj t of each Jordan block Jj , where j ∈ {1, . . . , p}.

Proposition 3.8. Consider the matrix A ∈ Rd×d , and let Jj for j ∈

{1, . . . , p} be the Jordan blocks for the real Jordan normal form with eigen-
values ρj . The matrix exponentials eJj t for t ∈ R are then given as follows.
(i) If ρj is real, i.e. Jj ∈ Rdj ×dj is of the form (3.6), we obtain
 
ρ 1
 j .
0
..
 
. 2 tdj −1
 . .   1 t t2 · · · (dj −1)! 

. .. 1  t
..
  
0 1 t . . .
  
 
. 
e 0 ρj ρ t
=e j  .
 
.. .. t2
 . . 2

 
0 1 t 
0 0 0 1

(ii) If ρj = aj + ibj ∈ C is not real, i.e. Jj ∈ R2dj ×2dj is of the form

(3.7), we obtain
 
C Id
 j . 2 .
0   
. . t2 tdj −1
 . .  G(t) tG(t) 2 G(t) · · · (dj −1)! G(t)

.. t

. Id ..
  
 2   . . 
   0 G(t) tG(t) . . 
e 0 C j a t
=e j  ,
 
.. .. t2
 . . 2 G(t) 
 
 0 G(t) tG(t) 
0 0 0 G(t)
 
cos(bj t) sin(bj t)
where G(t) = for all t ∈ R.
− sin(bj t) cos(bj t)

Proof. (i) If ρj is real, then Jj ∈ Rdj ×dj has the form J = P + D, where
 
0 1 0 
ρj 0

. .
 .. ..   ...
P = ..  and D =  .

. 1 ..
  .
0 0 0 ρj
56 3. Linear systems

Obviously, we have P D = λP = DP , i.e. Proposition 3.4 (ii) is applicable,

and we get
eJj t = e(P +D)t = eP t eDt for all t ∈ R .
The matrix P is nilpotent, since we have P ` = 0 for all ` ≥ dj , and thus we
get
dj −1 `
X t
eP t = P ` for all t ∈ R ,
`!
`=0
and this implies the assertion, since eDt = eρj t Iddj for all t ∈ R.
(ii) If ρj = aj +ibj ∈ C is not real, then Jj ∈ R2dj ×2dj has the form J = Q+B
with
 
02×2 Id2 0 
Cj 0

. .
.. .. ..
Q=

..

 and B =  . .

. Id2 ..
.
  
0 02×2 0 Cj
Due to QB = BQ, we can apply Proposition 3.4 (ii), and we get
eJj t = e(Q+B)t = eQt eBt for all t ∈ R .
The matrix Q is nilpotent, since we have Q` = 0 for all ` ≥ dj , and thus we
get
dj −1 `
X t
Qt
e = Q` for all t ∈ R ,
`!
`=0
and this implies the assertion, since eCj t = eaj t G(t) for all t ∈ R due to (3.5).


5. Exponential growth behaviour

We are interested in the exponential growth of the flow induced by a linear
system ẋ = Ax, where A ∈ Rd×d , and motivated by Remark 3.5, we define
the spectrum of A as

Σ(A) := Re ρ : ρ is an eigenvalue of A = {s1 , . . . , sq } .
One can show that for the linear system
ẋ = Ax ,
one gets a decomposition
Rd = E1 ⊕ · · · ⊕ Eq
into linear spaces Ej , j ∈ {1, . . . , q}, which are invariant in the sense that
x ∈ Ej for some j ∈ {1, . . . , q} =⇒ ϕ(t, x) ∈ Ej for all t ∈ R .
5. Exponential growth behaviour 57

Moreover, the Lyapunov exponent associated to non-trivial solutions start-

ing in Ej is given by sj , where j ∈ {1, . . . , q}:
ln kϕ(t, x)k
x ∈ Ej \ {0} =⇒ σLyap (ϕ(·, x)) = lim = sj .
t→∞ t
We do not prove this here, but note that this does not say anything about
Lyapunov exponents for solutions starting in points x ∈ Rd that not in the
spaces Ej , j ∈ {1, . . . , q}, which are the vast majority of points (if we do
not consider the trivial case q = 1).
Question: What is the exponential growth behaviour for solutions starting
in a general element x ∈ Rd ?
Instead of analysing the above decomposition and the associated Lyapunov
exponents in detail, due to time constraints, we only aim at estimating the
exponential growth of the norm of the matrix exponential function t 7→ eAt ,
which will be important later when analysing nonlinear systems. To do this
in a sharp way, we need the notion of a semi-simple eigenvalue. Recall that
an eigenvalue is semi-simple if its algebraic multiplicity equals its geometric
multiplicity. Equivalently, an eigenvalue is semi-simple if all the Jordan
blocks associated to this eigenvalue in the real Jordan normal form from
Theorem 3.7 are one-dimensional for real eigenvalues and two-dimensional
for non-real eigenvalues.
Proposition 3.9 (Exponential estimate for the matrix exponential func-
tion). Consider a matrix A ∈ Rd×d , and choose γ ∈ R such that

γ > max Re ρ : ρ is an eigenvalue of A .

If all eigenvalues ρ with Re ρ = max Re ρ : ρ is an eigenvalue of A are
semi-simple, we can use a smaller γ, given by

γ := max Re ρ : ρ is an eigenvalue of A .
Then there exists a K > 0 such that
eAt ≤ Keγt for all t ≥ 0 .

Proof. Let J be the real Jordan normal form of the matrix A, i.e. there
exists a T ∈ Rd×d such that J = T −1 AT . We first consider exponential
bounds for eJt . In an exercise, you have proved that there exists a constant
C ≥ 1 with
1
kBk∞ ≤ kBk ≤ CkBk∞ for all B ∈ Rd×d .
C
Here, the Euclidean operator norm k · k, as introduced in (2.8), is compared
to the infinity norm k · k∞ , which is defined as
kBk∞ := max |bij | for all B = (bij )i,j∈{1,...,d} ∈ Rd×d .
i,j∈{1,...,d}
58 3. Linear systems

Note that a general result says that all norms on finite-dimensional spaces
are equivalent. It follows now from Proposition 3.8 that all entries of eJt are
of the form g(t)tn eρt , where ρ is an eigenvalue real part, g(t) is a bounded
function and n ∈ {0, . . . , d−1}. Note that in case of a semi-simple eigenvalue,
n is always equal to 0. All entries of eJt can thus be estimated by K̃eγt for
t ≥ 0, where K̃ > 0 is chosen appropriately. Note here that we need the
estimate tn eρt ≤ K 0 eγt for ρ < γ and t ≥ 0 here, where K 0 > 0 has to be
chosen appropriately. Finally, we get the inequality
(3.4)
eAt = T eJt T −1 ≤ kT k T −1 eJt ≤ CkT k T −1 eJt ∞
−1 γt
≤ CkT k T K̃ e for all t ≥ 0 ,
| {z }
=:K

which finishes the proof of this theorem. 

6. Variation of constants formula

We know now that the flow ϕ of a general autonomous (homogeneous) lin-
ear system ẋ = Ax, where A ∈ Rd×d , is given by the matrix exponential
function ϕ(t, x) = eAt x. We are now interested in the general solution to
the corresponding inhomogeneous equation
ẋ = Ax + g(t) , (3.9)
where g : I → Rd is a continuous function on an interval I ⊂ R.
Proposition 3.10 (Variation of constants formula). The general solution
to (3.9) is given by
Z t
λ(t, t0 , x0 ) = eA(t−t0 ) x0 + eA(t−s) g(s) ds for all t, t0 ∈ I and x0 ∈ Rd .
t0

Proof. We first show that for fixed t0 ∈ I, the function µg : I → Rd , given

by
Z t Z t
µg (t) = eA(t−s) g(s) ds = eAt e−As g(s) ds
t0 t0
is a solution to (3.9). This follows from
Z t
µ̇g (t) = AeAt e−As g(s) ds + eAt e−At g(t) = Aµg (t) + g(t) ,
t0

where we have used the product rule and the fundamental theorem of cal-
culus. Note that the general solution to the homogeneous system ẋ = Ax is
given by
λh (t, t0 , x0 ) = eA(t−t0 ) x0 ,
6. Variation of constants formula 59

which follows from the definition of a flow and Theorem 3.3. We show now
that for a fixed (t0 , x0 ), the function
νt0 ,x0 (t) := λh (t, t0 , x0 ) + µg (t) for all t ∈ I
is a solution to (3.9). This follows from
ν̇t0 ,x0 (t) = λ̇h (t, t0 , x0 ) + µ̇g (t) = Aλh (t, t0 , x0 ) + Aµg (t) + g(t)
= A(λh (t, t0 , x0 ) + µg (t)) + g(t) = Aνt0 ,x0 (t) + g(t) .
Obviously, the function νt0 ,x0 satisfies the initial condition x(t0 ) = x0 , which
proves that λ as given in the statement of the proposition is the general
solution to (3.9). 

Note that a nonautonomous (homogeneous) linear system of the form

ẋ = A(t)x ,
cannot be solved analytically in general, that means in the case when the
matrices A(t), t ∈ I, do not commute. Of course, one-dimensional matrices
always commute, and thus, we can compute the general solution to the
inhomogeneous linear differential equation
ẋ = a(t)x + g(t) ,
where a : I → R and g : I → R are continuous functions. Similarly to above,
one can show that its general solution is given by
Rt Z t R
a(s) ds t
λ(t, t0 , x0 ) = e 0
t x0 + e s a(τ ) dτ g(s) ds for all t, t0 ∈ I and x0 ∈ R .
t0
 Chapter 4

Nonlinear systems

This chapter deals with nonlinear autonomous differential equations. Such

systems are not solvable in general, in contrast to the linear autonomous
systems in the previous chapter, and we aim at understanding their be-
haviour in the spirit of the Russian mathematician and physicist Aleksandr
M. Lyapunov (1857–1918) and the French mathematician Henri Poincaré
(1854–1912), who are the so-called fathers of the Qualitative Theory of Dy-
namical Systems. This theory aims at understanding dynamical systems
(for instance, given in the form of ordinary differential equations) from a
qualitative point of view, i.e. without being able to solve them.
This chapter gives an introduction into the basic elements of this theory.
In particular, the basic elements of stability theory will be explained, and
the direct method of Lyapunov will be introduced. Finally, the asymptotic
behaviour of two-dimensional differential equations will be analysed using
the Poincaré–Bendixson theory.

1. Stability
We are interested in the dynamical behaviour of the flow of an autonomous
differential equation in the vicinity of equilibria. We will distinguish between
different types of stability close to equilibria, and we will learn about criteria
that indicate such behaviour.

1.1. Basic definitions. We introduce different notions of stability for an

autonomous differential equation

ẋ = f (x) , (4.1)

61
62 4. Nonlinear systems

where f : D → Rd is locally Lipschitz continuous and D ⊂ Rd is an open

set. We denote the flow of this differential equation by ϕ.
Definition 4.1 (Notions of stability). Let x∗ be an equilibrium of (4.1),
i.e. f (x∗ ) = 0.
(i) x∗ is called stable if for all ε > 0, there exists a δ > 0 such that
kϕ(t, x) − x∗ k < ε for all x ∈ Bδ (x∗ ) and t ≥ 0 .
(ii) x∗ is called unstable if x∗ is not stable.
(iii) x∗ is called attractive if there exists a δ > 0 such that
lim ϕ(t, x) = x∗ for all x ∈ Bδ (x∗ ) .
t→∞

(iv) x∗ is called asymptotically stable if x∗ is both stable and attractive.

(v) x∗ is called exponentially stable if there exist δ > 0, K ≥ 1 and
γ < 0 such that
kϕ(t, x) − x∗ k ≤ Keγt kx − x∗ k for all x ∈ Bδ (x∗ ) and t ≥ 0 .
(vi) x∗ is called repulsive if there exists a δ > 0 such that
lim ϕ(t, x) = x∗ for all x ∈ Bδ (x∗ ) .
t→−∞

Figure 4.1. Notions of stability, from left to right: stable, unstable,

attractive, asymptotically stable.

Example 4.2 (Stability of one-dimensional linear differential equations).

We study the trivial equilibrium x∗ = 0 of the linear differential equation
ẋ = αx ,
where α ∈ R. This differential equation has the flow ϕ : R × R → R,
ϕ(t, x) = xeαt . Depending on the parameter α, the equilibrium x∗ has
differential stability properties.
(i) x∗ is stable for α ≤ 0. Choose δ := ε for a given ε > 0. We have
|ϕ(t, x) − x∗ | = |x|eαt ≤ |x| < ε = δ for all t ≥ 0 and x ∈ (−δ, δ).
(ii) x∗ is unstable for α > 0. We fix ε = 1 and choose δ > 0 arbitrarily.
Then |ϕ(t, 2δ )| = 2δ eαt > ε for some t > 0, since eαt → ∞ as t → ∞.
1. Stability 63

(iii) x∗ is exponentially stable for α < 0, since |ϕ(t, x)| = |x|eαt =

Keγt |x| for t ≥ 0 and x ∈ Bδ (0) with γ := α < 0 and K := δ = 1.
(iv) x∗ is repulsive for α > 0. Choose δ := 1. Then for all x ∈ (−δ, δ),
we have |ϕ(t, x)| = |xeαt | = |x||eαt | < δ|eαt | → 0 as t → −∞.

We see in this example that the eigenvalues for this one-dimensional linear
system are crucial for the stability of the system. This applies to higher-
dimensional linear systems as well, and we will make this more precise later.
It is clear from the definitions that exponential stability implies asymptotic
stability. However, perhaps surprisingly, there is no relation between stabil-
ity and attractivity.
Example 4.3 (On the relation between stability and attractivity). We show
in this example that, in general, the notions of stability and attractivity are
not related.
(i) For the differential equation ẋ = 0, every point is an equilibrium,
and all equilibria are stable (choose δ := ε for a given ε > 0). It is
clear that no equilibrium is attractive in this example. Note that
stability does not imply attractivity can also been seen when look-
ing at the harmonic oscillator, which we studied in Example 1.10.
(ii) The two-dimensional differential equation
p
ẋ = x + xy − (x + y) x2 + y 2 ,
p
ẏ = y − x2 + (x − y) x2 + y 2
can be understood well when looking at the corresponding system
in polar coordinates
ṙ = r(1 − r) , φ̇ = r(1 − cos φ) .
The phase portrait is plotted in Figure 4.2. We note that this
system has exactly two equilibria: (0, 0), which is an unstable knot
with many tangents, and (1, 0), which admits a so-called homoclinic
orbit, given by the unit circle. Note that an orbit is called homo-
clinic if converges forward and backward in time to the same equi-
librium. It is clear, from looking at the polar coordinate system,
that the equilibrium (1, 0) is attractive. It can be even proved that
all orbits starting outside of the equilibrium (0, 0) converge to the
other equilibrium (1, 0) in forward time: limt→∞ ϕ(t, (x, y)) = (1, 0)
for all (x, y) 6= (0, 0). This is clearly seen in Figure 4.2, and can
be shown easily using the Poincaré–Bendixson theory, which will
be developed in Section 4 below. Although the equilibrium (1, 0) is
attractive, it is also clear that it is not stable, since orbits starting
in (x, y) on the unit circle very close to (1, 0), but with positive y
64 4. Nonlinear systems

take very long to complete the journey round the circle to come
close to (1, 0) from below, and for y → 0, this time converges to ∞.
Hence (−1, 0) ∈ ϕ(t, Bδ (1, 0)) for all δ > 0 and t sufficiently large.
This means that (1, 0) is unstable.

Figure 4.2. The attractive equilibrium (0, 1) is not stable.

In addition to the homoclinic orbit (connecting to (1, 0) in both time direc-

tions), there are also so-called heteroclinic orbits in this example (connect-
ing to (1, 0) in forward time and (0, 0) in backward time). In fact, all orbits
starting inside the unit circle have this property.
Definition 4.4 (Homoclinic and heteroclinic orbits). Consider the differ-
ential equation (4.1) with associated flow ϕ.
(i) An orbit O(x) for some x ∈ D is called a homoclinic orbit if there
exists an equilibrium x∗ ∈ D \ {x} such that
lim ϕ(t, x) = x∗ and lim ϕ(t, x) = x∗ .
t→∞ t→−∞

(ii) An orbit O(x) for some x ∈ D is called a heteroclinic orbit if there

exists two different equilibria x∗1 6= x∗2 ∈ D such that
lim ϕ(t, x) = x∗1 and lim ϕ(t, x) = x∗2 .
t→∞ t→−∞

1.2. Stability of linear systems. Before looking more into nonlinear sys-
tems, we study stability of the trivial equilibrium of autonomous linear sys-
tems, and we note that in the Example 4.2, we have already explored this
in the one-dimensional case. Recall also that in Remark 3.5, we identified
1. Stability 65

that the real part of eigenvalues play a role in the exponential growth be-
haviour in two-dimensional linear systems, and this was confirmed for higher
dimensions in Proposition 3.9.
Theorem 4.5 (Stability of linear systems). Consider the autonomous linear
system
ẋ = Ax , (4.2)
d×d ∗
where A ∈ R . Then the trivial equilibrium x = 0 of this system is
(i) stable if and only if the following two statements hold:
(a) the real part of all eigenvalues of A is non-positive, i.e. we
have Re ρ ≤ 0 for all eigenvalues ρ of A, and
(b) the eigenvalue ρ is semi-simple for all eigenvalues ρ of A with
Re ρ = 0.
(ii) exponentially stable if and only if Re ρ < 0 for all eigenvalues ρ
of A.

Proof. Let J be the real Jordan normal form of the matrix A, i.e. there
exists a T ∈ Rd×d such that J = T −1 AT . Using Proposition 3.4 (i), we get
ϕ(t, x) = eAt x = T eJt T −1 x (4.3)
for the flow of (4.2).
(i) (⇒) We show that if either (a) or (b) does not hold, then x∗ = 0 is not
stable. If either (a) or (b) is false, this means that either there exists an
eigenvalue ρ with positive real part, or there exists an eigenvalue ρ with real
part 0 that is not semi-simple. In both cases, it follows from Proposition 3.8
that t 7→ eJj t , for a Jordan block corresponding to the eigenvalue ρ is un-
bounded for t ∈ [0, ∞). Hence it follows from (4.3) that one element of the
matrix eAt is unbounded for t ∈ [0, ∞), say the element in the k-th row and
`-th column. Set ε := 1 and choose δ > 0. This implies that
ϕ t, 2δ e` = 2δ ϕ t, e` = 2δ eAt e` ∈



/ Bε (0) for some t ≥ 0 .
|{z}
∈Bδ (0)

Hence x∗ is not stable.

(⇐) It follows from Proposition 3.9 that (a) and (b) imply that there exists
a K > 0 such that
eAt ≤ K for all t ≥ 0 .
To prove that x∗ = 0 is stable, choose ε > 0 arbitrarily, and define δ := Kε .
Then for all x ∈ Bδ (0), we get
kϕ(t, x)k = eAt x ≤ eAt kxk < Kδ = ε for all t ≥ 0 .
Hence, x∗ = 0 is stable.
(ii) The strategy of the proof is similar to (i) and is left as an exercise. 
66 4. Nonlinear systems

It should be noted, although only stability and exponential stability are de-
scribed in the above theorem, the results immediately give a characterisation
for instability due the equivalence formulation in (i). Moreover, it can be
shown easily that attractivity of autonomous linear systems is equivalent to
exponential stability, and thus, all stability notions from Definition 4.1 can
clearly be understood in the context of linear systems.

1.3. Hyperbolicity. Since we understand the stability of linear systems

very well now, we focus our attention on nonlinear systems. As motivated
in the beginning of Chapter 3 on linear systems, we aim at using linear
systems to understand nonlinear systems locally in the neighbourhood of a
reference solution. In the autonomous context, this works well for reference
solutions that are constant, i.e. given by equilibria. To motivate the results
that will follow, we consider two different two-dimensional linear systems
from Section 2 in Chapter 3, to which we add some nonlinear perturbation.

Example 4.6 (Nonlinear perturbations of linear systems).

We first consider
−1 0
a perturbation of a saddle equilibrium. Let A1 := 0 1 , and consider the
linear system
   
ẋ x
= A1 , (4.4)
ẏ y
as well as the nonlinearly perturbed system
1 2
     
ẋ x 5 y
= A1 + 3 2 1 2 , (4.5)
ẏ y 10 x + 5 y

The phase portraits of both systems a neighbourhood [−2, 2] × [−2, 2] of

the trivial equilibrium (0, 0) are given in Figure 4.3, and we can see that
the nonlinear perturbation induces just a slight perturbation of the phase
portrait.

Figure 4.3. Phase portraits of (4.4) (on the left) and (4.5) (on the right).
1. Stability 67

The situation

is much different when we consider a centre equilibrium. Let
A2 := 01 −1
0 , and consider the linear system
   
ẋ x
= A2 , (4.6)
ẏ y
as well as the nonlinearly perturbed system
−5x(x2 + y 2 )
     
ẋ x
= A2 + , (4.7)
ẏ y −5y(x2 + y 2 )
The phase portraits of both systems in a neighbourhood [−2, 2] × [−2, 2] of
the trivial equilibrium (0, 0) are given in Figure 4.4, and in contrast to the
example above, the nonlinearly perturbed system behaves much differently
than the linear system, in the sense that all orbits converge forward to in
time to the trivial equilibrium (0, 0).

Figure 4.4. Phase portraits of (4.6) (on the left) and (4.7) (on the right).

The reason for this difference that if a linear system has only eigenvalues
with nonzero real parts, this implies exponential attractivity or repulsivity
in invariant linear subspaces, and this exponential behaviour cannot be de-
stroyed locally by a nonlinear perturbation. This is formulated precisely in
the so-called Hartman–Grobman theorem, which we will skip due to time
constraints. However, the boundary between attractivity and repulsivity is
given by an eigenvalue with real part 0, and different nonlinear perturba-
tions can make systems with zero real part eigenvalue attractive or unstable
(but not exponentially stable or exponentially unstable).
This observation leads to concept of hyperbolicity.
Definition 4.7 (Hyperbolicity). A matrix A ∈ Rd×d is called hyperbolic if
all eigenvalues λ of A have non-zero real part, i.e. Re λ 6= 0. An equilibrium
x∗ of a differential equation
ẋ = f (x) ,
d d
where f : D ⊂ R → R is continuously differentiable, is called hyperbolic
if the matrix f 0 (x∗ ) ∈ Rd×d is hyperbolic.
68 4. Nonlinear systems

We close this section by looking at the non-hyperbolic one-dimensional case.

Example 4.8 (One-dimensional non-hyperbolicity). We consider the linear
one-dimensional differential equation
ẋ = 0 .
Then the trivial equilibrium x∗= 0 is stable (see also Example 4.3 (i)), and
0 ∗
the derivative f (x ) of the right hand side (f (x) = 0 for all x ∈ R) is 0,
and thus this equilibrium (as well as all the other equilibria – every point
is an equilibrium) is non-hyperbolic. We consider the following nonlinear
perturbations and discuss its effect on the stability of x∗ = 0:
(i) ẋ = x2 : the equilibrium x∗ = 0 becomes unstable, although it
attracts all points starting in the negative half-line.
(ii) ẋ = −x2 : the equilibrium x∗ = 0 becomes unstable, although it
attracts all points starting in the positive half-line.
(iii) ẋ = x3 : the equilibrium x∗ = 0 becomes unstable, and all points
(except x∗ ) move away from x∗ forward in time.
(iv) ẋ = −x3 : the equilibrium x∗ = 0 becomes asymptotically stable,
but it is not exponentially stable.

1.4. Linearised stability. In this subsection, we consider a hyperbolic

equilibrium, for which all eigenvalues of the linearisation have negative real
parts. This situation is significantly easier to analyse than the case when
there are both positive and negative eigenvalue real parts (and due to time
constraints, we will not cover this case in detail, but we give some in insights
in the Subsection 1.5 below).
Under the assumption that the real parts of all eigenvalues of the hyperbolic
equilibrium are negative, it follows from Theorem 4.5 that the trivial equilib-
rium of the linearised system is exponentially stable. We show now that this
exponential stability is transferred to the nonlinear system. Importantly in
the proof, we make use of the following lemma, the Gronwall lemma. Via
the Gronwall lemma, we obtain an explicit exponential estimate from an
implicit inequality (and such implicit inequalities appear quite frequently in
the context of differential equations).
Lemma 4.9 (Gronwall lemma). We consider a continuous function u :
[a, b] → R defined on an interval [a, b], and let c, d ≥ 0. We assume that the
function u satisfies the implicit inequality
Z t
0 ≤ u(t) ≤ c + d u(s) ds for all t ∈ [a, b] . (4.8)
a
Then we have the explicit estimate
u(t) ≤ c ed(t−a) for all t ∈ [a, b] . (4.9)
1. Stability 69

Proof. Because the function u is continuous on the compact set [a, b], it is
bounded, i.e. there exists an M > 0 such that
u(t) ≤ M for all t ∈ [a, b] .
Using this in (4.8), we get
u(t) ≤ c + M d(t − a) for all t ∈ [a, b] .
Using this improved estimate in (4.8), we arrive at
1
u(t) ≤ c + cd(t − a) + M d2 (t − a)2 for all t ∈ [a, b] .
2
Inductively, we get after n steps that
n−1
X dk (t − a)k M dn (t − a)n
u(t) ≤ c + ,
k! n!
k=0 | {z }
| {z } →0
→ced(t−a)
which implies the claim in the limit n → ∞. 

We obtain exponential stability for an equilibrium of a nonlinear system, for

which the real parts of all eigenvalues of the linearisation are negative. This
is referred to as linearised stablity.
Theorem 4.10 (Linearised stability). Let D ⊂ Rd be open and f : D →
Rd be continuously differentiable, and consider the autonomous differential
equation
ẋ = f (x) . (4.10)
∗ ∗
Assume that x is an equilibrium of (4.10) (i.e. f (x ) = 0) such that for
all eigenvalues λ ∈ C of the linearisation f 0 (x∗ ) ∈ Rd×d , we have Re λ < 0.
Then the equilibrium x∗ of (4.10) is exponentially stable.

Proof. To simplify notation in the proof, we assume that the equilibrium

is given by x∗ = 0. To see this, one can make a change of variables using a
transformation y = x − x∗ , which yields a differential equation for y having
the zero equilibrium. It should be noted that this transformation (since it
is just a translation) does not change any stability properties.
Step 1. Some useful estimates.
We first write
ẋ = f (x) = f 0 (0) x + r(x) , (4.11)
| {z }
=:A
where the term
r(x) := f (x) − f 0 (0)x = o(kxk) for all x ∈ D
is a higher order term, defined on a neighbourhood of 0, and clearly satisfies
r(0) = 0 and r0 (0) = 0. We denote the flow of (4.11) by ϕ.
70 4. Nonlinear systems

Since the real parts of the eigenvalues of the matrix A are negative, due to
Proposition 3.9, there exist constants K > 0 and γ < 0 such that
eAt ≤ Keγt for all t ≥ 0 . (4.12)
γ
We choose a positive number M < − K . Since r is continuously differen-
0
tiable, there exists ρ > 0 such that kr (x)k ≤ M for all x ∈ Bρ (0). Due to
the mean value inequality (Theorem 2.8), we get then get
kr(x)k ≤ M kxk for all x ∈ Bρ (0) . (4.13)
Finally, we define for each initial value x ∈ Bρ (0) the escape time

Te (x) := sup T > 0 : kϕ(t, x)k ≤ ρ for all t ∈ [0, T ) .
Note that Te (x) can be ∞, which is desirable case, and, as we will see later,
this is true for small enough x.
Step 2. We show that for all x0 ∈ Bρ (0), we have
kϕ(t, x0 )k ≤ Ke(KM +γ)t kx0 k for all t ∈ [0, Te (x0 )) . (4.14)
We note that the solution t 7→ ϕ(t, x0 ) of the differential equation (4.10) is
also a solution to the nonautonomous linear differential equation
ẋ = Ax + r(ϕ(t, x0 )) ,
for which the variation of constants formula can be applied (see Proposi-
tion 3.10), and we obtain
Z t
ϕ(t, x0 ) = eAt x0 + eA(t−s) r(ϕ(s, x0 )) ds .
0
Hence, for all t ∈ [0, Te (x0 )), we have
Z t
At
kϕ(t, x0 )k ≤ e · kx0 k + eA(t−s) · kr(ϕ(s, x0 ))k ds
0
(4.12)
Z t
γt
≤ Ke kx0 k + Keγ(t−s) kr(ϕ(s, x0 ))k ds
0
(4.13)
Z t
γt
≤ Ke kx0 k + Keγ(t−s) M kϕ(s, x0 )k ds .
0

We multiply this inequality with e−γt and obtain the implicit estimate
Z t
−γt −γs
u(t) := e kϕ(t, x0 )k ≤ Kkx0 k + KM e|{z} kϕ(s, x0 )k ds .
0
u(s)

We use the Gronwall lemma (Lemma 4.9) for the above defined function u,
and get the explicit estimate
e−γt kϕ(t, x0 )k ≤ Kkx0 keKM t for all t ∈ [0, Te (x0 )) ,
1. Stability 71

which finishes the proof of this step.

Step 3. Finalisation of the proof.
From the fact that t 7→ e(KM +γ)t is monotonically decreasing, we get
(4.14)
kϕ(t, x0 )k ≤ Ke(KM +γ)t kx0 k
t=0
≤ Kkx0 k ≤ ρ for all t ≥ 0 and x0 ∈ Bρ/K (0) ,
and it follows that Te (x0 ) = ∞. Hence, (4.14) implies that (4.11) is expo-
nentially stable, since KM + γ < 0. 

Question: Does this result imply anything for hyperbolic equilibria, for
which the eigenvalues of the linearisation have only positive real parts?
Yes, a time-reversed version of Theorem 4.10 implies that such equilibria
are exponentially attractive backward in time, which means that they are
repulsive.
A more complicated situation is given when there eigenvalues with both pos-
itive and negative real parts, and we will discuss this in the next subsection,
but firstly, we now show that this result can be applied to the pendulum sys-
tem, which turns out to have an exponentially stable equilibrium whenever
friction is taken into account.
Example 4.11 (Pendulum, exponentially stable equilibrium). Consider a
pendulum moving along a circle of radius r > 0, with a mass m > 0 and
friction coefficient k > 0. Let x denote the angle from the vertical. The
force tangential to the circle depends on both the position x and the speed
ẋ of the pendulum, and is given by
Ftan (x, ẋ) = −mg sin(x) − krẋ ,
see Figure 4.5.

x
r

mg sin(x)
mg cos(x)
mg

Figure 4.5. The pendulum.

72 4. Nonlinear systems

Newton’s law reads as mrẍ = Ftan (x, ẋ), and we thus get the second-order
one-dimensional differential equation
g k
ẍ = − sin(x) − ẋ ,
r m
which we can transform into the first-order two-dimensional system
ẋ = y ,
g k
ẏ = − sin(x) − y .
r m
(see Repetition Material 1 ). The equilibria of this system are given by
(nπ, 0), where n ∈ Z, which corresponds to the pendulum being in vertical
position (pointing down for even n and pointing up for odd n. We linearise
this system in (0, 0) (or equivalently in (nπ, 0) for even n), and obtain the
linearisation  
0 1
,
− gr − mk
 q

k 2 g
which gives two eigenvalues λ± := 12 − m k
± m − 4 r . It follows that
the real parts of both eigenvalues λ± are negative, and hence, Theorem 4.10
implies that the equilibria (nπ, 0) for n even are exponentially stable. Note
) −4 gr < 0 then both eigenvalues are complex, and if ( m
k 2
that if ( m ) −4 gr ≥ 0,
k 2

then both eigenvalues are real and negative. The stability of the equilibria
(nπ, 0) for n odd will be discussed below in Example 4.16.

1.5. Stable and unstable sets, invariant sets. In the previous subsec-
tion, we have looked at the case where all eigenvalues of the linearisation in
an equilibrium are negative, which implies exponential stability of the equi-
librium. We now look at situations where some eigenvalues have positive
and some eigenvalues have negative real parts.
Definition 4.12 (Stable and unstable set). Consider the differential equa-
tion (4.1) with associated flow ϕ, and let x∗ be an equilibrium. We define
the stable set of x∗ as
n o
W s (x∗ ) := x ∈ D : lim ϕ(t, x) = x∗ ,
t→∞
and the unstable set of x∗ is defined as
n o
W u (x∗ ) := x ∈ D : lim ϕ(t, x) = x∗ ,
t→−∞

Note that if x∗ is an attractive equilibrium, then W s (x∗ ) is called domain

of attraction, and it follows from the definition of attractivity that this is
a neighbourhood of x∗ . Moreover, W s (x∗ ) is an open set in this case (see
exercises).
1. Stability 73

Example 4.13 (Stable and unstable set of a linear system). Given a matrix
A ∈ Rd×d , we consider the linear system
ẋ = Ax ,
which we assume to be hyperbolic, i.e there are no eigenvalues with zero
real part (which is equivalent to {s1 , . . . , sq } ∩ {0} = ∅, where Σ(A) =
{s1 , . . . , sq }). It was mentioned in Section 5 of Chapter 3 that one can show
that
Rd = E1 ⊕ · · · ⊕ Eq ,
with Lyapunov exponents s1 < · · · < sq associated to these spaces. Due to
hyperbolicity, there exists an k ∈ {1, . . . , q + 1} such that s` < 0 for all ` < k
and s` > 0 for all ` ≥ k. One can show that
k−1
M q
M
W s (0) = Ei and W u (0) = Ei ,
i=1 i=k
so that the above decomposition can be rewritten as
Rd = W s (0) ⊕ W u (0) .

We now look at nonlinear perturbations of these objects in a two-dimensional

example.
Example 4.14 (Stable and unstable set of a hyperbolic equilibrium). Recall
the first two systems discussed in Example 4.6. With A1 := −1 0


0 1 , we
consider the linear system
   
ẋ x
= A1 , (4.15)
ẏ y
as well as the nonlinearly perturbed system
1 2
     
ẋ x 5 y
= A1 + 3 2 1 2 , (4.16)
ẏ y 10 x + 5 y
From Example 4.13, it follows for the linear example (4.15) that R2 = E1 ⊕E2
with W s (0) = E1 being the x-axis and W u (0) = E2 being the y-axis, which
are both one-dimensional objects. If we look at the phase portrait of the
nonlinearly perturbed system (4.16), given in Figure 4.3, we notice that this
one-dimensional object survives – both the stable and unstable set seem to
be given by the union of two trajectories and the equilibrium point, and the
trajectories are tangential to the linear spaces E1 and E2 in the equilibrium
0 respectively.

We note that stable and unstable sets in these two examples are also called
stable and unstable manifolds of the equilibrium 0. We do not give a proper
definition of manifolds here, and only note that manifolds are objects that
locally look like Euclidean spaces in each point.
74 4. Nonlinear systems

Stable and unstable manifolds are special cases of so-called invariant man-
ifolds, and the notion of an invariant set is very important in the theory of
differential equations. We require from invariant sets that orbits starting in
it do not leave the set. Recall the definition of orbits O(x) and half-orbits
O+ (x), O− (x) from Definition 2.25.
Definition 4.15 (Invariance). Consider the differential equation (4.1).
Then a set M ⊂ D is called
(i) positively invariant if for all x ∈ M , we have O+ (x) ⊂ M ,
(ii) negatively invariant if for all x ∈ M , we have O− (x) ⊂ M ,
(iii) invariant if for all x ∈ M , we have O(x) ⊂ M .

Note that sets that consist of equilibria or periodic orbits are invariant. Sta-
ble and unstable sets are also invariant, and any union of orbits is invariant.
Accordingly, unions of half-orbits of the form O+ (x) or O− (x) are positively
invariant or negatively invariant, respectively.
We identify stable and unstable sets for the pendulum.
Example 4.16 (Pendulum, saddle equilibrium). We consider again the pen-
dulum from Example 4.11, given by the one-dimensional system of order two
g k
ẍ = − sin(x) − ẋ ,
r m
which we can transform due to Extra Material 1 into the first-order two-
dimensional system
g k
ẋ = y , ẏ = − sin(x) − y
r m
We have seen that the equilibria of this system are given by (nπ, 0), where
n ∈ Z, and the analysis in Example 4.11 showed that the equilibria (nπ, 0)
for even n ∈ Z are exponentially stable. We analyse now the equilibria
(nπ, 0) for odd n ∈ Z and obtain the the linearisation
 
0 1
g k ,
r −m
 q

k 2 g
which gives two eigenvalues λ± := 12 − m k
± m + 4 r . It follows that
we have λ− < 0 < λ+ , and one can see in the phase portrait in Figure 4.6
that the equilibria (nπ, 0) for n odd have stable and unstable sets/manifolds
that are given by one-dimensional curves.
Question: Can you give a physical interpretation of these two invariant
sets/manifolds?
We note that under the assumption that the right hand side of a differen-
tial equation is continuously differentiable, the so-called stable and unstable
2. Limit sets 75

Figure 4.6. Phase portrait of the pendulum.

manifold theorem says that for a hyperbolic equilibrium in dimension two

with one negative and one positive eigenvalue, the stable and unstable sets
are given by one-dimensional curves that intersect in the equilibrium and
are tangential to the eigenspaces of the linearisation. A similar statement
holds also in higher dimensions, where the curves need to be replaced by
higher-dimensional manifolds. We do not attempt to formulate the stable
and unstable manifold theorem here, but note that an application of this
theorem yields that instability of the linearisation (in form of an eigenvalue
with positive real part) carries over to the equilibrium of the nonlinear sys-
tem.

2. Limit sets
The asymptotic behaviour of differential equations (that is the limiting be-
haviour for t → ∞ and t → −∞) is determined by certain types of invariant
set, so-called limit sets. We study this for an autonomous differential equa-
tion of the form
ẋ = f (x) , (4.17)
where f : D → Rd is locally Lipschitz continuous and D ⊂ Rd is an open
set. We denote the flow of this differential equation by ϕ.
We now introduce two important classes of invariant sets, so-called omega
and alpha limits sets. Their importance is due to the fact that they describe
the asymptotic behaviour. Note that invariance is not part of the following
definition, but it will follow later from it.
76 4. Nonlinear systems

Definition 4.17 (Omega and alpha limit sets). Consider the flow ϕ of the
differential equation (4.17), and let x ∈ D.
(i) A point xω ∈ D is called omega limit point of x, if there exists a
sequence {tn }n∈N such that limn→∞ tn = ∞ and
xω = lim ϕ(tn , x) .
n→∞

We denote by ω(x) the set of all omega limit points of x. ω(x) is

called omega limit set of x.
(ii) A point xα ∈ D is called alpha limit point of x, if there exists a
sequence {tn }n∈N such that limn→∞ tn = −∞ and
xα = lim ϕ(tn , x) .
n→∞

We denote by α(x) the set of all alpha limit points of x. α(x) is

called alpha limit set of x.

Note that the omega limit set of a point x is empty if sup Jmax (x) < ∞, and
the alpha limit set to be nonempty requires inf Jmax (x) = −∞.
We look at omega limit sets for the differential equation from Example 2.26.
Example 4.18. Consider the autonomous two-dimensional differential
equation
ẋ = y + x(1 − x2 − y 2 ) ,
(4.18)
ẏ = −x + y(1 − x2 − y 2 ) .
We have seen already in Example 2.26 that this differential equation has
the trivial equilibrium (0, 0), which is unstable, and there exists a periodic
orbit, given by the unit circle S1 . All orbits of this system that do not start
in the unstable trivial equilibrium approach this periodic orbit in forward
time. It is possible to show that


{(0, 0)} : (x, y) = (0, 0) ,
ω (x, y) =
S1 : (x, y) 6= (0, 0) ,
and 

 {(0,10)} : k(x, y)k < 1 ,
α (x, y) = S : k(x, y)k = 1 ,
∅ : k(x, y)k > 1 .


Question: Can you establish this rigorously? You can use the explicit rep-
resentation of the flow from Example 2.26. Proving this will become easier
when we have established the Poincaré–Bendixson theorem, see Section 4
below.

We now derive an alternative characterisation of omega and alpha limit sets.

2. Limit sets 77

Proposition 4.19 (Alternative characterisation of limit sets). Consider the

flow ϕ of the differential equation (4.17), and let x ∈ D. Then we have
\
\

ω(x) = O+ ϕ(t, x) and α(x) = O− ϕ(t, x) . (4.19)
t≥0 t≤0

Before proving this proposition, we would like to understand its content

better using the following example.

Example 4.20. We again consider the differential equation (4.18). As de-

scribed above, for any point (x, y) 6= 0, the half-orbit O+ (x, y) approaches

the periodic orbit S1 , see Figure 4.7, so we have that

O+ (x, y) = O+ (x, y) ∪ S1 .

Hence,
\
\ + 
O ϕ(t, (x, y)) ∪ S1


O+ ϕ(t, (x, y)) =
t≥0 t≥0
\
= S1 ∪ O+ ϕ(t, (x, y)) = S1 .

t≥0

+ 6 1) or S1 (if
T

This follows, since t≥0 O ϕ(t, (x, y)) is either ∅ (if k(x, y)k =
k(x, y)k = 1).

Figure 4.7. Explanation of Proposition 4.19 using a half-orbit ap-

proaching a periodic orbit.

Proof of Proposition 4.19. We will only prove ω(x) = t≥0 O+ ϕ(t, x) ,
T
since the statement concerning alpha limit sets can be shown similarly.
78 4. Nonlinear systems

(⊂). Choose y ∈ ω(x). Then there exists a sequence {tn }n∈N with
limn→∞ tn → ∞ and y = limn→∞ ϕ(tn , x). This implies that for any t ≥ 0,
we have


y ∈ O+ ϕ(t, x) ,
since there exists an n0 = n0 (t) such that ϕ(tn , x) ∈ O+ ϕ(t, x) for all




n ≥ n0 . Hence y ∈ t≥0 O+ ϕ(t, x) .
T



(⊃). Choose y ∈ t≥0 O+ ϕ(t, x) . This implies that y ∈ O+ ϕ(t, x) for
T
all t ≥ 0, and in particular this means that
B1/n (y) ∩ O+ ϕ(n, x) 6= ∅ for all n ∈ N .

Hence, for all n ∈ N, there exists tn ≥ n such that ϕ(tn , x) ∈ B1/n (y). This
implies that limn→∞ tn → ∞ and y = limn→∞ ϕ(tn , x), and the proof is
finished, since y ∈ ω(x). 

Omega and alpha limit sets have the following important properties.
Proposition 4.21 (Properties of omega and alpha limit sets). Consider the
differential equation (4.17), and let x ∈ D. Then the following statements
hold.
(i) The omega limit set ω(x) is invariant. In addition, if O+ (x) is
bounded and O+ (x) ⊂ D, then ω(x) is non-empty and compact.
(ii) The alpha limit set α(x) is invariant. In addition, if O− (x) is
bounded and O− (x) ⊂ D, then α(x) is non-empty and compact.

Proof. We only prove (i), since (ii) can be shown similarly, and we denote
the flow of (4.17) by ϕ.
Step 1. ω(x) isnonempty.
The sequence ϕ(n, x) n∈N is bounded and thus has a convergent subse-

quence ϕ(nk , x) k∈N with limit in D. The limit of this subsequence is an
omega limit point of x, and thus, ω(x) is nonempty.
Step 2. ω(x) is compact.


This follows directly from (4.19), since each of the sets O+ ϕ(t, x) ⊂ O+ (x)
are compact, since they are bounded and closed, and the intersection ω(x)
over these sets is also bounded and closed, and thus compact.
Step 3. ω(x) is invariant.
To show that ω(x) is invariant, we need to show that for all x0 ∈ ω(x), we
have O(x0 ) ⊂ ω(x). Choose τ ∈ Jmax (x0 ). Since x0 ∈ ω(x), there exists a
sequence {tn }n∈N with limn→∞ tn = ∞ and
x0 = lim ϕ(tn , x) . (4.20)
n→∞
3. Lyapunov functions 79

Now consider the sequence sn := tn + τ , where n ∈ N, which also con-

verges to ∞ as n → ∞. Note that the flow ϕ is a continuous function (see
Remark 2.28). This implies that

lim ϕ(sn , x) = lim ϕ(tn + τ, x) = lim ϕ(τ, ϕ(tn , x))

n→∞ n→∞ n→∞
  (4.20)
= ϕ τ, lim ϕ(tn , x) = ϕ(τ, x0 ) ,
n→∞

and this means that ϕ(τ, x0 ) ∈ ω(x), which finishes the proof of this propo-
sition. 

3. Lyapunov functions
We have seen in Subsection 1.4 and Subsection 1.5 that stability can be de-
duced from the linearisation around hyperbolic equilibria. In particular, if
all eigenvalues of the linearisation f 0 (x∗ ) around an equilibrium x∗ of a dif-
ferential equation have negative real parts, then x∗ is asymptotically stable
(or even exponentially stable), and it follows that the domain of attraction
W s (x∗ ) contains a neighbourhood of the equilibrium x∗ .
In that sense, the method of linearisation provides local information, but
often it is useful to know more about global properties of the domain of
attraction W s (x∗ ). In addition, sometimes the method of linearisation can
not be used to determine stability of nonlinear systems, and this is the case
when the equilibrium is non-hyperbolic. It turns out that so-called Lyapunov
functions can be of help with regard to these restrictions to the methods
discussed so far. They are useful tools to prove stability (or instability)
of (not necessarily hyperbolic) equilibria and to determine their basin of
attraction.
In this section, we consider autonomous differential equations of the form

ẋ = f (x) , (4.21)

where f : D → Rd is locally Lipschitz continuous and D ⊂ Rd is an open

set. We denote the flow of this differential equation by ϕ.
Lyapunov functions are real-valued functions V : D → R and can be thought
of as energy functions. Their main property is that they decrease along
solutions of (4.21), for instance in systems with friction where energy is
lost. To model this, we will define the notion of an orbital derivative, i.e. the
derivative of V along solutions.

Definition 4.22 (Orbital derivative). Consider the differential equation

(4.21), and let V : D → R be a continuously differentiable function. Then
80 4. Nonlinear systems

the orbital derivative V̇ of the function V is defined by

d
X ∂V
V̇ (x) := V 0 (x) · f (x) = (x)fi (x) .
∂xi
i=1

Here, the row vector V 0 (x) ∈ R1×d is the gradient of V at x ∈ D.

Note that while the gradient V 0 does not depend on (4.21), the orbital
derivative V̇ does. In fact, V̇ describes the derivative of V along solutions
µ : I → D of (4.21). This follows from
d
V (µ(t)) = V 0 (µ(t)) · µ̇(t) = V̇ (µ(t)) for all t ∈ I , (4.22)
dt
where we have used the chain rule.
We consider the energy of the pendulum and study its orbital derivative.
Example 4.23 (Pendulum, orbital derivative of energy function). We con-
sider the pendulum
ẋ = y ,
g k
ẏ = − sin(x) − y .
r m
which we have studied first in Example 4.11 and then in Example 4.16.
Its kinetic energy is given by 21 m(ry)2 , and its potential energy is given by
mg(r − r cos(x)). The sum of these two parts constitutes the function
1
V (x, y) := m(ry)2 + mgr 1 − cos(x) .


2
V is equal to 0 in the asymptotically stable equilibria (nπ, 0) for even n, and
V is positive outside of these equilibria. We compute the orbital derivative
of V , given by
 
2 y
V̇ (x, y) = (mgr sin(x), mr y) = −kr2 y 2 .
− gr sin(x) − m
k
y
We considered so far only the case with positive friction k > 0, and in this
case V (x, y) < 0 whenever y 6= 0. Interesting is also the case without friction
(k = 0), since in this case, V̇ (x, y) = 0. This immediately implies that
solutions stay on the level sets of V , see Figure 4.8 for the phase portrait.
We concentrate only on functions V that do not increase along solutions,
which are so-called Lyapunov functions.
Definition 4.24 (Lyapunov function). Consider the differential equation
(4.21), and let V : D → R be a continuously differentiable function. Then
V is called a Lyapunov function if
V̇ (x) ≤ 0 for all x ∈ D .
3. Lyapunov functions 81

Figure 4.8. Phase portrait of the pendulum without friction.

Note that any Lyapunov function decreases along solutions, i.e.

V (ϕ(t, x)) ≤ V (x) for all t ∈ [0, sup Jmax (x)) , (4.23)
which can be seen as follows. Integrating (4.22) implies that for all t ∈
[0, sup Jmax (x)), we have
Z t
V (ϕ(t, x)) − V (ϕ(0, x)) = V̇ (ϕ(s, x)) ds ≤ 0 ,
| {z } 0 | {z }
=x ≤0

which shows (4.23). This property implies immediately that sublevel sets of
Lyapunov functions are positively invariant.
Proposition 4.25 (Sublevel sets of Lyapunov functions are positively in-
variant). Consider the differential equation (4.21) with a Lyapunov function
V : D → R. Then any sublevel set of the form

Sc := x ∈ D : V (x) ≤ c ,
where c ∈ R, is positively invariant.

Proof. Assume that Sc is not positively invariant. Then there exists an

x ∈ Sc and t > 0 such that ϕ(t, x) ∈/ Sc . Since x ∈ Sc implies V (x) ≤ c, and
ϕ(t, x) ∈
/ Sc implies that V (ϕ(t, x)) > c, this contradicts (4.23) and finishes
the proof. 

The following theorem says that if a Lyapunov function has a strict local
minimum in an equilibrium, then the equilibrium is stable. The result is
called a direct method, since it can get information on the stabilty behaviour
82 4. Nonlinear systems

of solutions directly from the right hand side of the differential equation,
and there is no need to solve the differential equation.
Theorem 4.26 (Lyapunov’s direct method for stability). Consider the dif-
ferential equation (4.21) with an equilibrium x∗ ∈ D, and let V : D → R be
a Lyapunov function such that
V (x∗ ) = 0 and V (x) > 0 for all x ∈ D \ {x∗ } .
Then the equilibrium x∗ is stable.

Proof. To prove stability of x∗ , we choose ε > 0. Since D is open, there

exists an ε̃ ∈ (0, ε] such that Bε̃ (x∗ ) ⊂ D, and we define
m := min V (x) : kx − x∗ k = ε̃ > 0 ,


which is positive due to V (x) > 0 outside of x = x∗ . Since V is a continuous

function, there exists a δ = δ(ε̃) ∈ (0, ε̃) such that
m
0 ≤ V (x) ≤ for all x ∈ Bδ (x∗ ) .
2
This implies with (4.23) that
m
V (ϕ(t, x)) ≤ V (x) ≤ for all x ∈ Bδ (x∗ ) and t ≥ 0 .
2
The orbit starting in x ∈ Bδ (x∗ ) can thus not leave the ε̃-neighbourhood of
x∗ , since the V is at least m on the boundary of this neighbourhood. Since
ε ≥ ε̃, this implies that x∗ is stable. 

Proving stability of an equilibrium is often checked by looking at the lineari-

sation, and exponential stability follows if the real parts of the eigenvalues of
the linearisation are negative (see Theorem 4.10). The above theorem using
Lyapunov functions can be helpful if the linearisation is non-hyperbolic and
thus not amenable to an analysis using Theorem 4.10.
Example 4.27 (Application of Lyapunov’s direct method for stability). We
consider the two-dimensional differential equation
ẋ = −y − xy 2 ,
ẏ = x − yx2 .
The only equilibrium of this system is the trivial equilibrium (0, 0). Indeed,
if ẋ = 0, then either y = 0 or 1 + xy = 0, but if y = 0, then x(1 − xy) = 0
implies x = 0, and if 1 + xy = 0, then 0 = x(1 − xy) = 2x implies x =
0
which contradicts 1 + xy = 0. The linearisation in (0, 0) is given by 10 −1
0 ,
and thus the system is non-hyperbolic and stability cannot be deduced from
Theorem 4.10.
3. Lyapunov functions 83

We first show that the trivial equilibrium is stable by considering the qua-
dratic function V (x, y) = x2 + y 2 for all (x, y) ∈ R2 . Then
−y − xy 2
 
V̇ (x, y) = (2x, 2y) = −4x2 y 2 ≤ 0
x − yx2
so (0, 0) is stable. We will see later in Example 4.32 that (0, 0) is even
asymptotically stable.

We now show that the existence of a Lyapunov function gives information

about the location of omega limit sets.
Theorem 4.28 (La Salle’s invariance principle). Consider the differential
equation (4.21) with a Lyapunov function V : D → R. Then

ω(x) ⊂ y ∈ D : V̇ (y) = 0 for all x ∈ D .

Proof. Assume to the contrary that

z ∈ ω(x) and V̇ (z) < 0 .
Then for some τ > 0, we have V (ϕ(τ, z)) < V (z). Since z ∈ ω(x), there
exists a sequence {tn }n∈N with limn→∞ tn = ∞ such that
lim ϕ(tn , x) = z .
n→∞

It is clear that we can choose the sequence {tn }n∈N such that tn+1 − tn > τ
for all n ∈ N. Due to V̇ (y) ≤ 0 for all y ∈ D, this implies that
V (ϕ(tn+1 , x)) ≤ V (ϕ(τ + tn , x)) = V (ϕ(τ, ϕ(tn , x))) for all n ∈ N .
Taking the limit n → ∞ on both sides of this inequality gives
V (z) ≤ V (ϕ(τ, z)) ,
which contradicts V (ϕ(τ, z)) < V (z). Note that here we have used that the
flow ϕ is continuous in x (see also Remark 2.28) 

We use La Salle’s principle to understand the asymptotic behaviour of the

pendulum better.
Example 4.29 (Pendulum, all orbits converge to equilibria). For a positive
friction coefficient k > 0, we consider the pendulum
ẋ = y ,
g k
ẏ = − sin(x) − y .
r m
which we have studied several times (in Examples 4.11, 4.16 and 4.23).
In particular, we have analysed the stability of all equilibria (nπ, 0), where
n ∈ Z, and we proved that the equilibria with n even are exponentially stable
84 4. Nonlinear systems

and the equilibria with n odd are unstable. We also found the Lyapunov
function
1
V (x, y) := m(ry)2 + mgr 1 − cos(x) ,


2
and we showed that
V̇ (x, y) = −kr2 y 2 .
Hence, V̇ (x, y) = 0 if and only if y = 0. Then La Salle’s principle implies
that
ω((x, y)) ⊂ R × {0} for all (x, y) ∈ R2 . (4.24)
It is possible to show, but this is left as a challenging exercise, that O+ ((x, y))
is bounded for all (x, y) ∈ R2 . Given this boundedness, Proposition 4.21
implies that ω(x) is nonempty. Proposition 4.21 also implies that ω((x, y))
is invariant. Suppose now that ω((x, y)) contains a point that is not an
equilibrium. Then (4.24) implies that there exists a (x̄, 0) ∈ ω((x, y)) with
x̄ 6= nπ for n ∈ Z. At this point (x̄, 0), we thus have ẏ = − gr sin(x̄) 6= 0.
This implies that the the flow starting in (x̄, 0) will leave the x-axis R × {0}
immediately (both forward and backward in time). Since the omega limit
set ω((x, y)) is invariant and we have (4.24), this cannot happen, and this
implies that such a point (x̄, 0) is not part of the the omega limit set ω((x, y)),
so
for all (x, y) ∈ R2 .

ω((x, y)) ⊂ (nπ, 0) : n ∈ Z
Since ω((x, y)) is connected (see exercises), this implies that ω((x, y)) is a
singleton. Since it is possible to show that when ω((x, y)) is a singleton,
then the flow starting in any (x, y) ∈ R2 converges in forward time to an
equilibrium (do this as an exercise, or look into the proof of Theorem 4.31
below).

These observations lead to a better formulation of La Salle’s principle.

Corollary 4.30 (Reformulation of La Salle’s invariance principle). Consider
the differential equation (4.21) with a Lyapunov function V : D → R. Then
for any x ∈ D, the omega limit set ω(x) is contained in the largest invariant
subset of y ∈ D : V̇ (y) = 0 . Here the largest invariant subset is given by

the union of invariant subsets of y ∈ D : V̇ (y) = 0

Proof. Let M be the largest invariant subset of y ∈ D : V̇ (y) = 0 .
Assume that for some x ∈ D, we have ω(x) \ M 6= ∅. Since ω(x) is invariant
due to Proposition 4.21, the set ω(x)∪M is invariant. Because of ω(x)∪M )
M , this contradicts the maximality of M . 

We apply La Salle’s principle to obtain asymptotic stability via Lyapunov

functions.
3. Lyapunov functions 85

Theorem 4.31 (Lyapunov’s direct method for asymptotic stability). Con-

sider the differential equation (4.21) with an equilibrium x∗ ∈ D, and let
V : D → R be a Lyapunov function such that

V (x∗ ) = 0 and V (x) > 0 for all x ∈ D \ {x∗ } ,

V̇ (x∗ ) = 0 and V̇ (x) < 0 for all x ∈ D \ {x∗ } .

Then the equilibrium x∗ is asymptotically stable.

Proof. The proof is divided in two steps.

Step 1. We show that there exists a δ > 0 such that ω(x) = {x∗ } for all
x ∈ Bδ (x∗ ).
Note first that x∗ is stable due to Theorem 4.26. Since D is open, there
exists an ε > 0 such that B2ε (x∗ ) ⊂ D, and since x∗ is stable, there exists a
δ > 0 such that for all δ > 0 such that for all x ∈ Bδ (x∗ ), we have

ϕ(t, x) ∈ Bε (x∗ ) for all t ≥ 0 .

We fix an x ∈ Bδ (x∗ ). Since ϕ(t, x) ∈ Bε (x∗ ) ⊂ D for all t ≥ 0, we get

O+ (x) ⊂ D, and this means that Proposition 4.21 implies that ω(x)  is
nonempty. La Salle’s principle (Theorem 4.28) implies that ω(x) ⊂ x̄ ∈
D : V̇ (x̄) = 0 = {x∗ }, and ω(x) being nonempty means that ω(x) = {x∗ }.
Step 2. We show that for all x ∈ Bδ (x∗ ), we have limt→∞ ϕ(t, x) = x∗ .
Assume to the contrary that we do not have limt→∞ ϕ(t, x) = x∗ . This
means that there exists an η > 0 and a sequence {tn }n∈N converging to ∞
such that
kϕ(tn , x) − x∗ k ≥ η for all n ∈ N ,
and since the sequence {ϕ(tn , x)}n∈N is bounded and bounded away from
the boundary of D, it has an accumulation point in D. This accumulation
point is an omega limit point of x, which contradicts the above observation
that the omega limit set ω(x) is the singleton {x∗ }. This finishes the proof
of this theorem. 

We use a slightly modified version of Theorem 4.31 (in the spirit of the
reformulation of La Salle’s principle in Corollary 4.30) to prove asymptotic
stability for the differential equation considered in Example 4.27.

Example 4.32 (Application of Lyapunov’s direct method for asymptotic

stability). We reconsider the differential equation

ẋ = −y − xy 2 ,
ẏ = x − yx2 ,
86 4. Nonlinear systems

from Example 4.27. We found a Lyapunov function for this differential

equation, given by V (x, y) = x2 + y 2 for all (x, y) ∈ R2 , and we showed that
V (0, 0) = 0 and V (x, y) > 0 for all (x, y) ∈ R \ {(0, 0)} .
In addition,
V̇ (x, y) = −4x2 y 2 for all (x, y) ∈ R2
shows that in the equilibrium (0, 0), we have V̇ (0, 0) = 0, but
(x, y) ∈ R2 : V̇ (x, y) = 0 = {0} × R ∪ R × {0} .


This means that Theorem 4.31 is not applicable directly. However, in Step 1
of the proof of Theorem 4.31, La Salle’s principle was applied to show that
the omega limit set is a singleton, and we will demonstrate now that instead,
we can apply the reformulation of this principle given by Corollary 4.30. This
is possible, since the largest invariant subset of
(x, y) ∈ R2 : V̇ (x, y) = 0 = {0} × R ∪ R × {0}


is given by the equilibrium {(0, 0)}. This follows from the fact that if we
start on the y-axis outside of trivial equilibrium (i.e. x = 0 and y 6= 0), then
ẋ = −y 6= 0, so we leave the y-axis immediately. The same holds for the
x-axis outside of trivial equilibrium (i.e. y = 0 and x 6= 0). In this case, we
get ẏ = x 6= 0, so we leave the x-axis immediately. It follows that a modified
version of Step 1 of Theorem 4.31 then implies that ω((x, y)) = {(0, 0)},
for (x, y) from a neighbourhood around (0, 0). This implies that (0, 0) is
asymptotically stable.

The following corollary to Theorem 4.31 shows that sublevel sets of Lya-
punov functions are part of the domain of attraction.
Corollary 4.33 (Sublevel sets of Lyapunov functions are subsets of the
domain of attraction). Under the assumptions of Theorem 4.31, we consider
the sublevel sets of the Lyapunov function V , which are of the form

Sc := x ∈ D : V (x) ≤ c ,
where c > 0. Then Sc is a subset of the domain of attraction W s (x∗ ) if
Sc ⊂ D is compact.

Proof. Let x ∈ Sc , and we need to show that limt→∞ ϕ(t, x) = x∗ . Proposi-

tion 4.25 implies that ϕ(t, x) ∈ Sc for all t ≥ 0, and since Sc is compact, we
get O+ (x) ⊂ D. This means due to Proposition 4.21, ω(x) is nonempty. La
Salle’s principle implies that ω(x) ⊂ x ∈ D : V̇ (x) = 0 = {x∗ }, and ω(x)


being nonempty means that ω(x) = {x∗ }. Then Step 2 of Theorem 4.31
finishes the proof. 
4. Poincaré–Bendixson theorem 87

4. Poincaré–Bendixson theorem
In this final section, we complete the set of tools that helps us to analyse two-
dimensional autonomous differential equations. This analysis decomposes
into a local and a global analysis.
In the local analysis, we
(i) locate the fixed points, and
(ii) determine their stability by linearising (which is possible if they are
hyperbolic).
In the global analysis, we
(i) look at nullclines in order to understand the global behaviour better
(see exercises),
(ii) try to find Lyapunov functions to understand stability of equilibria
(if non-hyperbolic) and domains of attractions,
(iii) locate periodic orbits (done in this section).
We consider the two-dimensional differential equation
ẋ = f (x) , (4.25)
where f : D → R2 is continuously differentiable function on an open set
D ⊂ R2 . We denote the flow of (4.25) by ϕ. Note that one can show (and
this is important for us) that the flow ϕ is continuously differentiable (see
Remark 2.28).

Theorem 4.34 (Poincaré–Bendixson theorem). Consider the differential

equation (4.25), and assume that for some x ∈ D, the positive half-orbit
O+ (x) lies in a compact subset K of D, which contains not more than finitely
many equilibria. Then one of the following three statements holds for the
omega-limit set ω(x).
(i) ω(x) is a singleton consisting of an equilibrium.
(ii) ω(x) is a periodic orbit.
(iii) ω(x) consists of equilibria and non-closed orbits. The non-closed
orbits in ω(x) converge forward and backward in time to equilibria
in ω(x), so they are either homoclinic or heteroclinic orbits.

We note that an analogous statement holds for alpha limit sets, when we
look at the negative half-orbit. The theorem of Poincaré–Bendixson shows in
particular that in two-dimensional differential equation, there is only very
regular behaviour and no chaos. Chaotic differential equations occur in
dimension three, for instance in the famous Lorenz system.
88 4. Nonlinear systems

The proof of the Poincaré–Bendixson theorem is quite involved, and we skip

it due to time constraints. The full proof is written down in Extra Material
2.
The Poincaré–Bendixson theorem is often applied in the following form to
prove the existence of a periodic orbit.
Corollary 4.35 (Existence of a periodic orbit). Consider the differential
equation (4.25), and assume that for some x ∈ D, the positive half-orbit
O+ (x) lies in a compact subset K of D that does not contain an equilibrium.
Then ω(x) is a periodic orbit.

We can apply the Poincaré–Bendixson theorem to prove the existence of a

periodic orbit.
Example 4.36 (Existence of a periodic orbit). We consider the two-
dimensional differential equation
ẋ = y ,
ẏ = −x + y(1 − x2 − 2y 2 ) .
We first show that M := (x, y) ∈ R2 : 31 ≤ x2 + y 2 ≤ 2 is positively


invariant. We show that the vector field of the right hand side points inwards
at the boundary of M . More precisely, we consider the scalar-valued function
V (x, y) = x2 +y 2 and show that the orbital derivative V̇ satisfies V̇ (x, y) < 0
for x2 +y 2 = 2 and V̇ (x, y) > 0 for x2 +y 2 = 31 . Firstly, the orbital derivative
reads as
V̇ (x, y) = 2xy + 2y − x + y(1 − x2 − 2y 2 ) = 2y 2 (1 − x2 − 2y 2 ) .

For x2 + y 2 = 2, we have V̇ (x, y) = 2y 2 (−1 − y 2 ) < 0 and for x2 + y 2 = 31 , we

have V̇ (x, y) = 2y 2 ( 32 − y 2 ) ≥ 0. This shows the positive invariance of M .
Note that the only equilibrium is clearly given by (0, 0) ∈ / M . We apply the
corollary to the Poincaré–Bendixson theorem (Corollary 4.35) and conclude
that the positively invariant set M contains a periodic orbit.