C HA P T E R 4

Systems of ODEs. Phase Plane.

Qualitative Methods

Tying in with Chap. 3, we present another method of solving higher order ODEs in
Sec. 4.1. This converts any nth-order ODE into a system of n first-order ODEs. We also
show some applications. Moreover, in the same section we solve systems of first-order
ODEs that occur directly in applications, that is, not derived from an nth-order ODE but
dictated by the application such as two tanks in mixing problems and two circuits in
electrical networks. (The elementary aspects of vectors and matrices needed in this chapter
are reviewed in Sec. 4.0 and are probably familiar to most students.)
In Sec. 4.3 we introduce a totally different way of looking at systems of ODEs. The
method consists of examining the general behavior of whole families of solutions of ODEs
in the phase plane, and aptly is called the phase plane method. It gives information on the
stability of solutions. (Stability of a physical system is desirable and means roughly that a
small change at some instant causes only a small change in the behavior of the system at
later times.) This approach to systems of ODEs is a qualitative method because it depends
only on the nature of the ODEs and does not require the actual solutions. This can be very
useful because it is often difficult or even impossible to solve systems of ODEs. In contrast,
the approach of actually solving a system is known as a quantitative method.
The phase plane method has many applications in control theory, circuit theory,
population dynamics and so on. Its use in linear systems is discussed in Secs. 4.3, 4.4,
and 4.6 and its even more important use in nonlinear systems is discussed in Sec. 4.5 with
applications to the pendulum equation and the Lokta–Volterra population model. The
chapter closes with a discussion of nonhomogeneous linear systems of ODEs.
N OT AT I O N. We continue to denote unknown functions by y; thus, y1(t), y2(t)—
analogous to Chaps. 1–3. (Note that some authors use x for functions, x 1(t), x 2(t) when
dealing with systems of ODEs.)
Prerequisite: Chap. 2.
References and Answers to Problems: App. 1 Part A, and App. 2.

4.0 For Reference:

Basics of Matrices and Vectors
For clarity and simplicity of notation, we use matrices and vectors in our discussion
of linear systems of ODEs. We need only a few elementary facts (and not the bulk of
the material of Chaps. 7 and 8). Most students will very likely be already familiar

124
130 CHAP. 4 Systems of ODEs. Phase Plane. Qualitative Methods

This quadratic equation in l is called the characteristic equation of A. Its solutions are
the eigenvalues l1 and l2 of A. First determine these. Then use (14*) with l 5 l1 to
determine an eigenvector x (1) of A corresponding to l1. Finally use (14*) with l 5 l2
to find an eigenvector x (2) of A corresponding to l2. Note that if x is an eigenvector of
A, so is kx with any k 0.

EXAMPLE 1 Eigenvalue Problem

Find the eigenvalues and eigenvectors of the matrix

24.0 4.0
(16) A5 c d
21.6 1.2

Solution. The characteristic equation is the quadratic equation

24 2 l 4
det ƒ A 2 lI ƒ 5 5 l2 1 2.8l 1 1.6 5 0.
21.6 1.2 2 l

It has the solutions l1 5 22 and l2 5 20.8. These are the eigenvalues of A.

Eigenvectors are obtained from (14*). For l 5 l1 5 22 we have from (14*)

(24.0 1 2.0)x 1 1 4.0x 2 50

21.6x 1 1 (1.2 1 2.0)x 2 5 0.

A solution of the first equation is x 1 5 2, x 2 5 1. This also satisfies the second equation. (Why?) Hence an
eigenvector of A corresponding to l1 5 22.0 is

2 1
(17) x (1) 5 c d. Similarly, x (2) 5 c d
1 0.8

is an eigenvector of A corresponding to l2 5 20.8, as obtained from (14*) with l 5 l2. Verify this. j

4.1 Systems of ODEs as Models

in Engineering Applications
We show how systems of ODEs are of practical importance as follows. We first illustrate
how systems of ODEs can serve as models in various applications. Then we show how a
higher order ODE (with the highest derivative standing alone on one side) can be reduced
to a first-order system.

EXAMPLE 1 Mixing Problem Involving Two Tanks

A mixing problem involving a single tank is modeled by a single ODE, and you may first review the
corresponding Example 3 in Sec. 1.3 because the principle of modeling will be the same for two tanks. The
model will be a system of two first-order ODEs.
Tank T1 and T2 in Fig. 78 contain initially 100 gal of water each. In T1 the water is pure, whereas 150 lb of
fertilizer are dissolved in T2. By circulating liquid at a rate of 2 gal>min and stirring (to keep the mixture uniform)
the amounts of fertilizer y1(t) in T1 and y2(t) in T2 change with time t. How long should we let the liquid circulate
so that T1 will contain at least half as much fertilizer as there will be left in T2?
SEC. 4.2 Basic Theory of Systems of ODEs. Wronskian 137

4.2 Basic Theory of Systems of ODEs.

Wronskian
In this section we discuss some basic concepts and facts about system of ODEs that are
quite similar to those for single ODEs.
The first-order systems in the last section were special cases of the more general system

y r1 5 f1(t, y1, Á , yn)

y r2 5 f2(t, y1, Á , yn)
(1)
Á

y rn 5 fn(t, y1, Á , yn).

We can write the system (1) as a vector equation by introducing the column vectors
y 5 [ y1 Á yn]T and f 5 [ f1 Á fn]T (where T means transposition and saves us
the space that would be needed for writing y and f as columns). This gives

(1) y r 5 f(t, y).

This system (1) includes almost all cases of practical interest. For n 5 1 it becomes
y r1 5 f1(t, y1) or, simply, y r 5 f (t, y), well known to us from Chap. 1.
A solution of (1) on some interval a , t , b is a set of n differentiable functions

y1 5 h 1(t), Á, yn 5 h n(t)

on a , t , b that satisfy (1) throughout this interval. In vector from, introducing the
“solution vector” h 5 [h 1 Á h n]T (a column vector!) we can write

y 5 h(t).

An initial value problem for (1) consists of (1) and n given initial conditions

(2) y1(t 0) 5 K 1, y2(t 0) 5 K 2, Á, yn(t 0) 5 K n,

in vector form, y(t 0) 5 K, where t 0 is a specified value of t in the interval considered and
the components of K 5 [K 1 Á K n]T are given numbers. Sufficient conditions for the
existence and uniqueness of a solution of an initial value problem (1), (2) are stated in
the following theorem, which extends the theorems in Sec. 1.7 for a single equation. (For
a proof, see Ref. [A7].)

THEOREM 1 Existence and Uniqueness Theorem

Let f1, Á , fn in (1) be continuous functions having continuous partial derivatives
0f1 >0y1, Á , 0f1 > 0yn, Á , 0fn > 0yn in some domain R of ty1 y2 Á yn-space
containing the point (t 0, K 1, Á , K n). Then (1) has a solution on some interval
t 0 2 a , t , t 0 1 a satisfying (2), and this solution is unique.
140 CHAP. 4 Systems of ODEs. Phase Plane. Qualitative Methods

4.3 Constant-Coefficient Systems.

Phase Plane Method
Continuing, we now assume that our homogeneous linear system

(1) y9 5 Ay

under discussion has constant coefficients, so that the n 3 n matrix A 5 [ajk] has entries
not depending on t. We want to solve (1). Now a single ODE y r 5 ky has the solution
y 5 Cekt. So let us try

(2) y 5 xelt.

Substitution into (1) gives y r 5 lxelt 5 Ay 5 Axelt. Dividing by elt, we obtain the
eigenvalue problem

(3) Ax 5 lx.

Thus the nontrivial solutions of (1) (solutions that are not zero vectors) are of the form
(2), where l is an eigenvalue of A and x is a corresponding eigenvector.
We assume that A has a linearly independent set of n eigenvectors. This holds in most
applications, in particular if A is symmetric (akj 5 ajk) or skew-symmetric (akj 5 2ajk)
or has n different eigenvalues.
Let those eigenvectors be x (1), Á , x (n) and let them correspond to eigenvalues
l1, Á , ln (which may be all different, or some––or even all––may be equal). Then the
corresponding solutions (2) are

(4) y (4) 5 x (1)el1t, Á , y (n) 5 x (n)elnt.

Their Wronskian W 5 W(y (1), Á , y (n)) [(7) in Sec. 4.2] is given by

x (1)
1 e
l1t Á x (n)
1 e
lnt
x (1)
1
Á x (n)
1

x (1)
2 e
l1t Á x (n)
2 e
lnt
x (1)
2
Á x (n)
2
l1t1 Á 1lnt
W 5 (y , Á , y (n)) 5
(1)
5e .
# Á # # Á #
x (1)
n e
l1t Á x (n)
n e
lnt
x (1)
n
Á x (n)
n

On the right, the exponential function is never zero, and the determinant is not zero either
because its columns are the n linearly independent eigenvectors. This proves the following
theorem, whose assumption is true if the matrix A is symmetric or skew-symmetric, or if
the n eigenvalues of A are all different.
148 CHAP. 4 Systems of ODEs. Phase Plane. Qualitative Methods

4.4 Criteria for Critical Points. Stability

We continue our discussion of homogeneous linear systems with constant coefficients (1).
Let us review where we are. From Sec. 4.3 we have

a11 a12 y r1 5 a11 y1 1 a12 y2

(1) y r 5 Ay 5 c d y, in components,
a21 a22 y r2 5 a21 y1 1 a22 y2.

From the examples in the last section, we have seen that we can obtain an overview of
families of solution curves if we represent them parametrically as y(t) 5 [ y1(t) y2(t)] T
and graph them as curves in the y1 y2-plane, called the phase plane. Such a curve is called
a trajectory of (1), and their totality is known as the phase portrait of (1).
Now we have seen that solutions are of the form

y(t) 5 xelt. Substitution into (1) gives y r (t) 5 lxelt 5 Ay 5 Axelt.

Dropping the common factor elt, we have

(2) Ax 5 lx.

Hence y(t) is a (nonzero) solution of (1) if l is an eigenvalue of A and x a corresponding

eigenvector.
Our examples in the last section show that the general form of the phase portrait is
determined to a large extent by the type of critical point of the system (1) defined as a
point at which dy2 >dy1 becomes undetermined, 0>0; here [see (9) in Sec. 4.3]

dy2 y r2 dt a21 y1 1 a22 y2

(3) 5 5 .
dy1 y 1r dt a11 y1 1 a12 y2

We also recall from Sec. 4.3 that there are various types of critical points.
What is now new, is that we shall see how these types of critical points are related
to the eigenvalues. The latter are solutions l 5 l1 and l2 of the characteristic equation

a11 2 l a12
(4) det (A 2 lI) 5 5 l 2 2 (a11 1 a22)l 1 det A 5 0.
a21 a22 2 l

This is a quadratic equation l2 2 pl 1 q 5 0 with coefficients p, q and discriminant ¢

given by

(5) p 5 a11 1 a22, q 5 det A 5 a11a22 2 a12a21, ¢ 5 p 2 2 4q.

From algebra we know that the solutions of this equation are

(6) l1 5 12 ( p 1 1¢), l2 5 12 ( p 2 1¢).

152 CHAP. 4 Systems of ODEs. Phase Plane. Qualitative Methods

4.5 Qualitative Methods for Nonlinear Systems

Qualitative methods are methods of obtaining qualitative information on solutions
without actually solving a system. These methods are particularly valuable for systems
whose solution by analytic methods is difficult or impossible. This is the case for many
practically important nonlinear systems

y r1 5 f1( y1, y2)

(1) y r 5 f(y) , thus
y 2r 5 f2( y1, y2).

In this section we extend phase plane methods, as just discussed, from linear systems
to nonlinear systems (1). We assume that (1) is autonomous, that is, the independent
variable t does not occur explicitly. (All examples in the last section are autonomous.)
We shall again exhibit entire families of solutions. This is an advantage over numeric
methods, which give only one (approximate) solution at a time.
Concepts needed from the last section are the phase plane (the y1 y2-plane), trajectories
(solution curves of (1) in the phase plane), the phase portrait of (1) (the totality of these
trajectories), and critical points of (1) (points (y1, y2) at which both f1( y1, y2) and f2( y1, y2)
are zero).
Now (1) may have several critical points. Our approach shall be to discuss one critical
point after another. If a critical point P0 is not at the origin, then, for technical
convenience, we shall move this point to the origin before analyzing the point. More
formally, if P0: (a, b) is a critical point with (a, b) not at the origin (0, 0), then we apply
the translation

~
y 1 5 y1 2 a, ~
y 2 5 y2 2 b

which moves P0 to (0, 0) as desired. Thus we can assume P0 to be the origin (0, 0), and
for simplicity we continue to write y1, y2 (instead of ~ y 1, ~
y 2). We also assume that P0 is
isolated, that is, it is the only critical point of (1) within a (sufficiently small) disk with
center at the origin. If (1) has only finitely many critical points, that is automatically
true. (Explain!)

Linearization of Nonlinear Systems

How can we determine the kind and stability property of a critical point P0: (0, 0) of
(1)? In most cases this can be done by linearization of (1) near P0, writing (1) as
y r 5 f( y) 5 Ay 1 h( y) and dropping h(y), as follows.
Since P0 is critical, f1(0, 0) 5 0, f2(0, 0) 5 0, so that f1 and f2 have no constant terms
and we can write

y 1r 5 a11 y1 1 a12 y2 1 h 1( y1, y2)

(2) y r 5 Ay 1 h(y), thus
y2r 5 a21 y1 1 a22 y2 1 h 2( y1, y2).

A is constant (independent of t) since (1) is autonomous. One can prove the following
(proof in Ref. [A7], pp. 375–388, listed in App. 1).
160 CHAP. 4 Systems of ODEs. Phase Plane. Qualitative Methods

13. y s 1 sin y 5 0 15. Trajectories. Write the ODE y s 2 4y 1 y 3 5 0 as a

14. TEAM PROJECT. Self-sustained oscillations. system, solve it for y2 as a function of y1, and sketch
(a) Van der Pol equation. Determine the type of the or graph some of the trajectories in the phase plane.
critical point at (0, 0) when m . 0, m 5 0, m , 0. y2
c=5
(b) Rayleigh equation. Show that the Rayleigh
equation5
c=4
Y s 2 m(1 2 13Y r 2)Y r 1 Y 5 0 (m . 0)
also describes self-sustained oscillations and that by
differentiating it and setting y 5 Y r one obtains the van
der Pol equation. –2 c=3 2
(c) Duffing equation. The Duffing equation is y1

ys 1 v20y 3
1 by 5 0
where usually ƒ b ƒ is small, thus characterizing a small
deviation of the restoring force from linearity. b . 0
and b , 0 are called the cases of a hard spring and a
soft spring, respectively. Find the equation of the
trajectories in the phase plane. (Note that for b . 0 all
these curves are closed.) Fig. 98. Trajectories in Problem 15

4.6 Nonhomogeneous Linear Systems of ODEs

In this section, the last one of Chap. 4, we discuss methods for solving nonhomogeneous
linear systems of ODEs

(1) y9 5 Ay 1 g (see Sec. 4.2)

where the vector g(t) is not identically zero. We assume g(t) and the entries of the
n 3 n matrix A(t) to be continuous on some interval J of the t-axis. From a general
solution y (h)(t) of the homogeneous system y r 5 Ay on J and a particular solution
y (p)(t) of (1) on J [i.e., a solution of (1) containing no arbitrary constants], we get a
solution of (1),

(2) y 5 y (h) 1 y (p).

y is called a general solution of (1) on J because it includes every solution of (1) on J.

This follows from Theorem 2 in Sec. 4.2 (see Prob. 1 of this section).
Having studied homogeneous linear systems in Secs. 4.1–4.4, our present task will be
to explain methods for obtaining particular solutions of (1). We discuss the method of

5
LORD RAYLEIGH (JOHN WILLIAM STRUTT) (1842–1919), English physicist and mathematician,
professor at Cambridge and London, known by his important contributions to the theory of waves, elasticity
theory, hydrodynamics, and various other branches of applied mathematics and theoretical physics. In 1904 he
was awarded the Nobel Prize in physics.