Chapter 3

Continuous Models Using Ordinary

Differential Equations

3.1 Introduction to Continuous Models

Continuous models are systems whose inputs and outputs are capable
of changing at any instant of time. A continuous model consists of
a dependent continuous variable, varying with some other independent
continuous variables. We use a first-order ordinary differential equation (or
a system of first-order ordinary differential equations) to model a continuous
system if we have some information or assumption about the rate of change
of the dependent variable(s) with respect to the independent variable.

Let us model a situation of rumor propagation or spreading of rumors.

Suppose, a company has 2000 employees. On Monday, a rumor begins to
spread among them that the CEO of the company has been hospitalized after
being infected with COVID-19, the infectious disease caused by the newly
discovered coronavirus.

We assume that the spread of rumors is proportional to the number of

possible meetings between employees who have heard the rumor and those
who have not. Let x(t) be the number of employees who have heard the rumor
after t days; then the scenario is modeled with a first-order ordinary differential
equation as

dx
= kx(2000 − x), (3.1)
dt

where k is a constant of proportionality.

We also assume that initially (when t = 0), only 40 employees heard the
rumor, as they all attended the same meeting. After 2 days, 200 employees
heard the rumor. Mathematically, this means y(0) = 40 and y(2) = 200.

DOI: 10.1201/9781351022941-3 123

ISTUDY
 124 Mathematical Modeling: Models, Analysis and Applications

Equation(3.1) can be easily solved using the separation of variables as follows:

Z   Z
dx 1 1 1
= k dt ⇒ + dx = k dt,
x(2000 − x) 2000 x 2000 − x
 
1 x
⇒ ln = kt + C,
2000 2000 − x
x
⇒ = A e2000kt , where A = e2000C .
2000 − x
2000 A e2000kt 2000A 1
⇒ x(t) = 2000kt
. Given x(0) = 40 ⇒ 40 = ⇒ A= .
1+A e 1+A 49
2000 e2000kt
⇒ x(t) = . After 2 days, 200 employees heard the rumor.
49 + e2000kt
 
2000 e4000k 1 49
⇒ 200 = ⇒ k= ln .
49 + e4000k 4000 9

Therefore, the number of employees who have heard the rumor is given by

2000 e 2 ln( 9 )
t 49
2000
x(t) = ⇒ x(t) = . (3.2)
t
ln( 49 ) − 2t ln( 49
9 ) + 1
49 + e 2 9 49 e

We now pose a question that how much time will the rumor take to reach
1000 employees? Substituting x(t) = 1000, we obtain

2000 2 ln 49
1000 = ⇒ t= = 4.6 days.
− 2t ln( 49
9 ) ln(49/9)
49 e +1

Therefore, the rumor reaches 1000 employees in 4 days, 14 hours and 24

minutes. So, if the rumor starts on Monday 9:30 am, it will reach 1000
employees by Friday 11:54 pm.

We observe from (3.2) that x(t) → 2000 as t → ∞. The obvious question

is why x(t) approaches a constant value 2000 as time t becomes large? Why
is it not approaching any other value? What is the significance of the value
2000? Answer to these questions with possible reasons will be discussed in the
next section.

ISTUDY
 Continuous Models Using Ordinary Differential Equations 125

3.2 Steady-State Solution

We consider a system of n non-linear, autonomous (does not explicitly
depend on time, that is, time-invariant) differential equations
dx̃
= f (x̃),
dt
where x̃ = (x1 , x2 , .............., xn )T ,
and f (x̃) = (f1 (x̃), f2 (x̃), .........., fn (x̃))T .

A steady-state solution or equilibrium solution or critical point is a

constant solution, that is, where the value of x̃ does not change over time
and is obtained by putting
dx̃
= 0.
dt
Please note that in order for the values of x̃ to be the same (constant) over
time, there must not be any change in x̃, implying dx̃ dt = 0. Therefore, the only
value(s) of x̃ for which this can happen is f (x̃) = 0, and so f (x̃) = 0 gives a
steady-state solution or an equilibrium solution or a fixed point.

In short, if we consider a nonlinear time-invariant system

dx̃
= f (x̃), x̃(t0 ) = x˜0 , where f : ℜ → ℜ,
dt
a point x˜e is an equilibrium point or a steady-state solution of the system if
f (x˜e ) = 0.

A steady-state solution or equilibrium point or critical point is said to

be isolated if there is a neighborhood to the equilibrium point that does not
contain any other equilibrium points. In what follows, we shall assume that
every critical point is isolated. Consider a simple growth model of a population:
dx
= ax − bx2 (a, b > 0).
dt
Here, the population grows linearly and there is a crowding effect or
intra-specific competition (−bx2 ), which depresses the rate of growth of the
population. The steady-state solution is given by
dx a
= 0 ⇒ ax − bx2 = 0 ⇒ x = 0 and .
dt b
a
Therefore, 0 and are two steady states or equilibrium solutions of this
b
population.

ISTUDY
 126 Mathematical Modeling: Models, Analysis and Applications

Again, consider a two species predator–prey system

dx  x
= rx 1 − − αxy,
dt k
dx
= −βy + γxy.
dt
dx dy
The steady-state solution is given by dt = 0, = 0.
dt  Solving, we get
 the
β r(kγ − β)
steady-state solutions of the system as (0, 0), (k, 0) and , . The
γ kαγ
last solution is feasible if kγ − β > 0 because populations are positive in real
life.

3.3 Stability
In layman’s language, we say that an equilibrium point or a steady-state
solution x˜e is locally stable, if all solutions that start near x˜e (that is, the
initial conditions are in the neighborhood of x˜e ) remain near x˜e for all the time.
Furthermore, if all the solutions starting near x˜e approach x˜e as t → ∞, we say
that the equilibrium point or steady-state solution x˜e is locally asymptotically
stable.

Consider the dynamical system satisfying

dx̃
= f (x̃), x̃(t0 ) = x˜0 , where f : ℜ → ℜ, (3.3)
dt
In the sense of Lyapunov, the equilibrium point or the steady-state solution
x˜e of (3.3) is stable at t = t0 if for any ǫ > 0, there exists a δ(t0 , ǫ) > 0 such
that

k x̃(t0 ) − x˜e k< δ ⇒ k x̃(t) − x˜e k< ǫ, ∀ t ≥ t0 .

The equilibrium point or the steady-state solution x˜e of (3.3) is

asymptotically stable (in the sense of Lyapunov) at t = t0 if
(i) x˜e is stable,
(ii) there exists δ(t0 ) > 0 such that

k x̃(t0 ) − x˜e k< δ ⇒ lim x̃(t) = 0.

t→∞

An equilibrium point or a steady-state solution is unstable if it is not stable.

Please note that the definitions of stability and asymptotic stability are local
definitions as they describe the behavior of a system near an equilibrium point.

ISTUDY
 Continuous Models Using Ordinary Differential Equations 127

3.3.1 Linearization and Local Stability Analysis

This method is known as Lyapunov’s first method or reduced method,
where the stability analysis of a steady-state or equilibrium point is done
by studying the stability of the corresponding linearized system in the
neighborhood (vicinity) of the steady state.

We consider the model of the form

dx
= f (x),
dt
whose local stability analysis we want to perform about the equilibrium point
x∗ (obtained by putting f(x) = 0). We give a small perturbation to the system
about the equilibrium point x∗ . Mathematically, this means we put x = X +x∗
into the above equation and get

dX
= f (x∗ + X) = f (x∗ ) + Xf ′ (x∗ ) + ... (higher-order terms),
dt
dX
≈ f ′ (x∗ )X, since f (x∗ ) = 0 and neglecting higher-order terms.
dt
Therefore, we conclude that the system is stable if f ′ (x∗ ) < 0 (decreasing
function) and unstable if f ′ (x∗ ) > 0 (increasing function). If f ′ (x∗ ) = 0, no
definite conclusion can be drawn from linear stability analysis.

Let us now consider the model given by the system of differential equations
of the form
dx dy
= f (x, y), = g(x, y). (3.4)
dt dt
Let (x∗ , y ∗ ) be the steady-state solution of (3.4), then f (x∗ , y ∗ ) = 0 and
g(x∗ , y ∗ ) = 0. We now give a small perturbation to the system about the
steady state, and mathematically this means we put x = X + x∗ and y =
Y + y ∗ . This implies

dX
= f (x∗ + X, y ∗ + Y )
dt
= f (x∗ , y ∗ ) + Xfx (x∗ , y ∗ ) + Y fy (x∗ , y ∗ ) + ... higher-order terms
(by Taylor series expansion of two variables).
Similarly,
dY
= g(x∗ , y ∗ ) + Xgx (x∗ , y ∗ ) + Y gy (x∗ , y ∗ ) + ... higher-order terms,
dt

where fx (x∗ , y ∗ ) is ∂f ∗ ∗
∂x evaluated at the steady state (x , y ). Since by
definition, f (x∗ , y ∗ ) = 0, g(x∗ , y ∗ ) = 0, by neglecting second and higher-order

ISTUDY
 128 Mathematical Modeling: Models, Analysis and Applications

terms, we obtain
dX
= fx (x∗ , y ∗ )X + fy (x∗ , y ∗ )Y,
dt
dY
= gx (x∗ , y ∗ )X + gy (x∗ , y ∗ )Y,
dt
which can be put in matrix form as
   
dx̃ x fx fy
= Ax̃, where x̃ = and A = . (3.5)
dt y gx gy

Let x̃ = vbeλt be a trial solution of (3.5), where vb (6= 0) is some fixed vector
which needs to be determined. Then
dx̃
= vbλeλt = Ab
v eλt ⇒ Ab
v = λb
v.
dt
From linear algebra, it can be easily concluded that λ is the eigenvalue of the
matrix A, whose eigenvector is νb. The eigenvalues are obtained by solving

fx − λ fy
det(A − λI) = 0 ⇒ = 0,
gx gx − λ
⇒ λ2 − (fx + gy )λ + fx gy − fy gx = 0 ⇒ λ2 − trace(A) λ + det(A) = 0.

Let λ1 and λ2 be the two eigenvalues of the matrix A. The necessary and
sufficient condition that λ1 and λ2 will be negative (if real) or have negative
real parts (if complex) is

trace(A) = fx + gy < 0,
det(A) = fx gy − fy gx > 0.

Routh-Hurwitz Stability Criteria:

(i) Quadratic Equation: λ2 + a1 λ + a2 = 0,

Stability Criteria: a1 > 0, a2 > 0.

(ii) Cubic Equation: λ3 + a1 λ2 + a2 λ + a3 = 0,

Stability Criteria: a1 > 0, a2 > 0, a3 > 0, a1 a2 − a3 > 0.

(iii) Fourth-Order Equation: λ4 + a1 λ3 + a2 λ2 + a3 λ + a4 = 0,

Stability Criteria: a1 > 0, a2 > 0, a3 > 0, a4 > 0, a1 a2 − a3 > 0,
a1 a2 a3 − a23 − a21 a4 > 0.

Example 3.3.1 The fish growth model by Von Bertalanffy [42] is given by

dF (t)
= αF 3/2 (t) − βF (t),
dt

ISTUDY
 Continuous Models Using Ordinary Differential Equations 129

where F (t) denotes the weight of the fish, and α and β are positive constants.

Solution: The equilibrium solution of the model is given by

dF (t)
= 0 ⇒ αF 3/2 − βF = 0,
dt
β2
⇒ F (αF 1/2 − β) = 0 ⇒ F ∗ = 0 and F ∗ = 2 .
α
3/2 ′
∗ 3 ∗1/2
Let W (F ) = αF − βF, then W (F ) = αF − β.
2
′
Now, W (0) = −β < 0 implies that the equilibrium point F ∗ = 0 is stable
′
and W (β/α) = 12 β > 0 implies that the equilibrium point F ∗ = β 2 /α2 is
unstable.
Example 3.3.2 Consider the autonomous non-linear system
dx
= β(2 − x) + x2 y,
dt
dy
= x − x2 y (β > 1).
dt
Find the equilibrium point(s) and comment on the stability (local) of the
equilibrium point(s).
Solution: The equilibrium points of the model is given by
β(2 − x) + x2 y = 0, (3.6)
2
x − x y = 0. (3.7)
(3.7) gives x = 0 and xy = 1. But, x = 0 does not satisfy (3.6). Hence,
2β
substituting y = x1 in (3.6), we get x = β−1 . Therefore, the only equilibrium
 
2β β−1
point is (x∗ , y ∗ ) = β−1 , 2β . Linearizing about the steady state (x∗ , y ∗ ),
we obtain the system as
 dx  4β 2
! 
dt
2 − β (β−1) 2 x
dy = −4β 2
.
dt −1 y
(β−1)2

The characteristic equation is given by

4β 2
2−β−λ (β−1)2
−4β 2
= 0,
−1 (β−1)2 − λ

β 3 + 5β − 2 4β 2
⇒ λ2 + + = 0.
(β − 1)2 β−1
Clearly, β 3 + 5β − 2 > 0, for all β > 1. Therefore, by Routh Hurwitz’s
 criteria,

∗ ∗ 2β β−1
the system is always stable about the equilibrium point (x , y ) = β−1 , 2β .

ISTUDY
 130 Mathematical Modeling: Models, Analysis and Applications

3.3.2 Lyapunov’s Direct Method

Lyapunov’s direct method, also known as Lyapunov’s second method,
determines the stability of a system without explicitly integrating the
differential equation
dx̃
= f (x̃), x̃(t0 ) = x˜0 , where f : ℜ → ℜ.
dt
The method is a generalization of the idea that if the potential energy
has a relative minimum at the equilibrium point, then the equilibrium
point is stable; otherwise, it is unstable. Russian mathematician Aleksandr
Mikhailovich Lyapunov generalized this principle to obtain a method for
studying the stability of the general autonomous system.

[Link] Lyapunov’s Condition for Local Stability

Consider the autonomous system
dx̃
= f (x̃), x̃(t0 ) = x˜0 , where f : ℜ → ℜ, (3.8)
dt
having isolated critical point at the origin x̃ = 0 and f (x̃) has continuous
partial derivatives for all x̃. Let V (x̃) be a positive definite function in a
neighborhood S of the origin x̃ = 0 and the derivative V̇ (x̃) of V (x̃) with
respect to the system (3.8) is negative semi-definite in the neighborhood S
of the origin x̃ = 0, then V (x̃) is called a Lyapunov function for the system
(3.8). Mathematically, V (x̃) is called a Lyapunov function if
(i) V (x̃) > 0 in the nbd. of the origin x̃ = 0,

(ii) V (0) = 0 for all x̃ = 0,

(iii) V̇ (x̃) ≤ 0 in the nbd. of the origin x̃ = 0,

(iv) V̇ (0) = 0 for all x̃ = 0.

If there exists a Lyapunov function V (x̃) for the system (3.8) in the nbd.
of the origin x̃ = 0, then the steady-state solution or the equilibrium point
x̃ = 0 is stable (locally).

If the derivative V̇ (x̃) of V (x̃) is negative semi-definite in the neighborhood

S of the origin x̃ = 0, then the steady state solution or the equilibrium point
x̃ = 0 is locally asymptotically stable (LAS).

Please note that there is no general method for constructing a Lyapunov

function, but if one can construct a Lyapunov function for the system (3.8),
then one can directly obtain information about the steady state of the
equilibrium point x̃ = 0 and hence the name Lyapunov’s direct method.

ISTUDY
 Continuous Models Using Ordinary Differential Equations 131

[Link] Lyapunov’s Condition for Global Stability

The equilibrium point x̃ = 0 is globally asymptotically stable if it is a
basin of attraction in the entire state space or entire domain. Thus, x̃ = 0 is
globally asymptotically stable (GAS) for the system (3.8) if

(i) x̃ = 0 is locally asymptotically stable,

(ii) The Lyapunov function V (x̃) is radially unbounded, that is,

| V (x̃) |→ ∞ as norm(x̃) =k x̃ k→ ∞.

Example 3.3.3 Consider the system

dx dy
ẋ = = −x − y, ẏ = = x − y3.
dt dt
Check the local and global stability of the stable.

Solution: We first show that the system is locally asymptotically stable

(LAS).

Lyapunov’s Indirect Method

Clearly, (0, 0) is the only steady-state solution (equilibrium point) of the

given system. Linearizing about the steady state (0, 0), we obtain the system
as
   dx    
ẋ dt −1 −1 x
= dy = .
ẏ dt
1 0 y
  √
−1 −1
The variational matrix has eigen values −1
2 ± i 23 , whose
1 0
real part is negative. Hence, by Lyapunov’s indirect method, the system is
locally asymptotically stable (LAS) about (0, 0).

Lyapunov’s Direct Method

We consider a function V (x, y) = ax2 +by 2 , where a and b are some positive
constants to be determined. Clearly, the function V (x, y) is positive definite
as V (x, y) > 0, ∀ (x, y) 6= (0, 0) and V (0, 0) = 0. Also,

dV ∂V ∂V
V̇ (x, y) = = ẋ + ẏ
dt ∂x ∂y
= 2ax(−x − y) + 2by(x − y 3 )
= −2ax2 − 2by 4 + 2xy(b − a).

ISTUDY
 132 Mathematical Modeling: Models, Analysis and Applications

If we choose b = a = 1, then
dV
V̇ (x, y) = = −2x2 − 2y 4 < 0, ∀ (x, y) 6= (0, 0),
dt
= 0, ∀ (x, y) = (0, 0).
This implies V̇ (x, y) is negative definite. Hence, by Lyapunov’s direct method,
(0,0) is locally asymptotically stable (LAS).

Global Stability

We now use Lyapunov’s conditions to prove that the system is globally

asymptotically stable (GAS). Already, we have shown that the function is
locally asymptotically stable. Now, the function V (x, y) = x2 +y 2p(a=1, b=1)
is also radially unbounded as V (x, y) → ∞ as norm(x, y) = x2 + y 2 →
∞. Therefore, by Lyapunov’s stability conditions, the system is globally
asymptotically stable (GAS).

3.4 Phase Plane Diagrams of Linear Systems

We consider a two-dimensional autonomous linear system of the form
dx
= λ1 x + λ2 y,
dt
dy
= λ3 x + λ4 y, (3.9)
dt
which can be written in matrix form as
   
dx̃ λ1 λ2 x
= A x̃, where A = and x̃ = . (3.10)
dt λ3 λ4 y

Clearly, the linear system has a unique steady state (0, 0), provided
det(A) = λ1 λ4 − λ2 λ3 6= 0. The solution of (3.10) can be visualized as
trajectories moving in the xy-plane and can be sketched, which are known
as phase portraits.
For a better understanding of the system, we consider a much more
simplified linear system of the form
dx dy
= λ1 x and = λ4 y, (3.11)
dt dt
where (0, 0) is the unique equilibrium solution of (3.11). This can be put in
matrix form as     
d x λ1 0 x
= .
dt y 0 λ4 y

ISTUDY
 Continuous Models Using Ordinary Differential Equations 133
 
λ1 0
The characteristic equation of the matrix B = is
0 λ4

λ1 − k 0
= 0 ⇒ (λ1 − k)(λ4 − k) = 0 ⇒ k = λ1 , λ4 ,
0 λ4 − k

which are the roots of the characteristic equation, also called the eigenvalues
of the matrix B.

The solution to (3.11) is

x = x(0)eλ1 t , y = y(0)eλ4 t ,

where x(0) and y(0) are the initial values.

Case I: If both the eigenvalues λ1 and λ4 are negative, all the trajectories
approach (0, 0), that is, all the solutions of the system converge to the
equilibrium solution (0, 0), no matter what the initial conditions may be (fig.
3.1(a)). The steady state (0, 0) is called a stable node (fig. 3.1(b)).

x(t), (t)
5 4

4 2

-2

1
-4

t
1 4 5 -4 -2 0 2 4

(a) Two negative eigenvalues. (b) Stable node.

FIGURE 3.1: Phase portrait showing all the trajectories approach (0, 0), a
stable node.

Case II: If both the eigenvalues λ1 and λ4 are positive, all the trajectories
move away from (0, 0), that is, all the solutions of the system diverge from the
equilibrium solution (0, 0), irrespective of the initial conditions (fig. 3.2(a)).
In this case, the steady state (0, 0) is called an unstable node (fig. 3.2(b)).

Case III: If the eigenvalues are opposite in sign, say, λ1 < 0 and λ4 > 0,
then x(t) decreases whereas y(t) increases exponentially (fig. 3.3(a)). All the
solutions, in this case, approach the steady state (0,0) for some time until
they come close to it, and then they start moving away from the steady state,

ISTUDY
 134 Mathematical Modeling: Models, Analysis and Applications

x(t)()

 4

2
0.15

0.10
-2

0.05
-4

t
  0.6  1.0 -4 -2 0 2 4

(a) Two positive eigenvalues. (b) Unstable node.

FIGURE 3.2: Phase portrait showing all the trajectories move away from
(0, 0), an unstable node.

irrespective of the initial solutions. The steady state (0, 0) is called a saddle
point (fig. 3.3(b)).

Case IV: If the eigenvalues are complex conjugates, say, λ1 = a + ib and

λ4 = a − ib, then the solutions are of the form

x(t) = x(0)eat (cos b + i sin b),

y(t) = y(0)eat (cos b − i sin b).

(i) If a < 0, then both the eigenvalues have negative real parts, the term
eat decays for increasing t, and hence the solution decreases (fig. 3.4(a)). In
this case, all the trajectories spiral towards the steady state (0, 0), irrespective
of the initial conditions and the steady state is known as a stable spiral or
stable focus (fig. 3.4(b)).

(ii) If a > 0, both the eigenvalues have positive real parts, the term eat
grows exponentially and so is the solution, for t > 0 (fig. 3.5(a)). In this case,
all the trajectories spiral away from the steady state (0, 0), irrespective of the
initial conditions and the steady state is called an unstable spiral or unstable
focus (fig. 3.5(b)).

(iii) If a = 0, both the eigenvalues are purely imaginary. In this case, all
the trajectories are closed orbits about the steady state (0, 0). The solutions
are periodic (fig. 3.6(a)) and the steady state is called a center (fig. 3.6(b)).

Please note that for simplicity, we have considered the linear system of the
form
dx dy
= λ1 x and = λ4 y. (3.12)
dt dt

ISTUDY
 Continuous Models Using Ordinary Differential Equations 135

x(t)()
4


2

0
150

100 -2

50
-4

t
  0.6  1.0  -4 -2 0 2 4

(a) Positive and negative eigenvalues. (b) A saddle point.

FIGURE 3.3: Phase portrait shows all the trajectories approach the steady
state (0,0) for some time until they come close to it, and then they start
moving away from it, a saddle.

A similar approach with (3.9) would have given the same 

result. In that case,
λ1 λ2
the characteristic equation of the matrix A = is
λ3 λ4

λ1 − k λ2
= 0 ⇒ k 2 − (λ1 + λ4 )k + (λ1 λ4 − λ2 λ3 ) = 0,
λ3 λ4 − k
1h √ i
⇒ k 2 − Trace(A) k + Det(A) = 0 ⇒⇒ k = Trace(A) ± Discriminant ,
2
where Discriminant = (Trace(A))2 − 4Det(A).

Case I: Let Det(A) = −D (<0, D being a positive quantity), then

Discriminanth = (Trace(A))2 − 4Det(A)i = (Trace(A))
h
2
+ 4D > Trace(A).i
1
√ 1
√
Therefore, 2 Trace(A) + Discriminant and 2 Trace(A) − Discriminant
will have opposite signs and the phase portrait will be a saddle. The sign of
Trace(A) is irrelevant in this case.

Case II: Let Det(A) > 0 and Discriminant = (Trace(A))2 − 4Det(A) > 0,
that is,
s
√ Det(A)
Discriminant = Trace(A) 1 − < Trace(A).
(Trace(A))2
h √ i h √ i
1 1
Therefore, 2 Trace(A) + Discriminant and 2 Trace(A) − Discriminant
will have the same signs. The eigen values will be positive if Trace(A) >0 and
negative if Trace(A) < 0.

ISTUDY
 136 Mathematical Modeling: Models, Analysis and Applications

x(t)#$(%)
4

1.0
2

0.5 0

-2

t
1 ! " 5 6
-4

-0.5 -4 -2 0 2 4

(a) Eigenvalues with negative real parts. (b) Stable spiral.

FIGURE 3.4: Phase portrait showing all the trajectories spiral towards (0, 0),
a stable spiral or a stable focus.

Conclusions:
(i) Discriminant > 0, Trace(A) < 0 ⇒ phase portrait is a stable node.
(ii) Discriminant > 0, Trace(A) > 0 ⇒ phase portrait is an unstable node.
(iii) If Discriminant = 0, then the phase portrait is stable or unstable node
according as Trace(A) < 0 or Trace(A) > 0.

Case III: Let Discriminant = (Trace(A))2 − 4Det(A) < 0, then the values
of k are complex conjugates, namely,

1h √ i 1h √ i
Trace(A) + i Discriminant and Trace(A) − i Discriminant .
2 2

The phase portrait will be a stable spiral if Trace(A) <0 and an unstable
spiral if Trace(A) > 0. If Trace(A) = 0, the phase portrait will be center.

Let λ and µ be the roots of the characteristic equation or the eigenvalues

of the matrix A. The possible cases and the results are now shown in tabular
form.

ISTUDY
 Continuous Models Using Ordinary Differential Equations 137

x(')()(') 4

& × 106
2

2 × 106
0

1 × 106
-2

t
5 6 7 8 9 10 -4

6
-1 × 10
-4 -2 0 2 4

(a) Eigenvalues with positive real parts. (b) Unstable spiral.

FIGURE 3.5: Phase portrait showing all the trajectories spiral away from
(0, 0), an unstable spiral or an unstable focus.

Nature of the Nature of the Stability of the

eigenvalues λ, µ steady state (0,0) steady state (0,0)
Real, unequal, same sign:
(i) λ < 0, µ < 0 Node Asymptotically stable
(ii) λ > 0, µ > 0 Node Unstable
Real, unequal, opposite sign:
(iii) λ < 0, µ > 0 Saddle point Unstable
Real, equal, same sign:
(iv) λ = µ = λ∗ < 0 Node Asymptotically stable
(v) λ = µ = λ∗ > 0 Node Unstable
Complex conjugates (a+ib, a-ib):
(vi) a < 0 Spiral point Asymptotically stable
(vii) a > 0 Spiral point Unstable
(viii) a=0 Center Stable but not
asymptotically stable

Thus, a steady state or an equilibrium point can be stable, asymptotically

stable, or unstable. A point is stable if the orbit of the system is inside a
bounded neighborhood to the point for all time t after some t0 . A point is
asymptotically stable if it is stable and the orbit approaches the critical point
as t → ∞. If the equilibrium point is not stable, then it is unstable.

Example 3.4.1 Determine the type and stability of steady state at (0, 0) of
the following systems:

dx dy dx dy
(i) = 2x + 7y, = −5x − 10y. (ii) = 3x − 4y, = 2x − y.
dt dt dt dt

ISTUDY
 138 Mathematical Modeling: Models, Analysis and Applications

x(t)*+(t)
4

6
2
4

2 0

t
5 10 15 20
-2
-2

-4 -4

-6
-4 -2 0 2 4

(a) Eigenvalues with zero real parts. (b) Center.

FIGURE 3.6: Phase portrait showing all the trajectories are closed orbits
about (0, 0), a center.

Solution: (i) The characteristic equation is given by

2−λ 7
= 0 ⇒ λ2 + 8λ + 15 = 0 ⇒ λ = −3, −5.
−5 −10 − λ

The eigenvalues are real, unequal and negative, implying that the phase
portrait is a stable node (asymptotically stable).

(ii) The characteristic equation is given by

2−λ 7
= 0 ⇒ λ2 − 2λ + 5 = 0 ⇒ λ = 1 + 2i, 1 − 2i.
−5 −10 − λ

The eigenvalues are imaginary with positive real part, implying that the phase
portrait is an unstable focus.

3.5 Continuous Models

Different models are formulated and discussed in this section.

3.5.1 Carbon Dating

Carbon dating (Carbon 14 dating) is a method, developed by W.F. Libby
at the University of Chicago in 1947 [97], that can be used to accurately date
archaeological samples to determine the ages of the plant (wood fossil) or any

ISTUDY
 Continuous Models Using Ordinary Differential Equations 139

material which got its carbon from the air. Carbon 14 (C 14 ), a radioactive
isotope of carbon, is a result of constant bombardment by radiation from the
sun in the atmosphere. During this bombardment, neutrons hit nitrogen 14
atoms and transmute them to carbon.

In a living organism, the absorption rate of C 14 balances the disintegration

rate of C 14 . When the organism dies and the body is preserved, it does
not absorb C 14 but disintegration continues. As mentioned before, C 14 is
radioactive in nature and has a half-life (the time taken by a substance
undergoing decay to decrease to half) of 5730 years. Scientists use this
information. The method of carbon dating involves measuring the strength of
C 14 archaeological samples or fossils and then comparing it with the expected
strength of C 14 in the atmosphere, to calculate the accurate age.

Suppose an archaeological sample was found whose age needs to be

determined. Let A(t) be the amount of C 14 present in the sample at any
time t, then
dA
= −λA (following radioactive decay law), (3.13)
dt
where λ is the decay constant of the sample. Integrating, we get,
A(t) = A0 e−λt ,
where A0 = A(0) is the amount of C 14 present in the sample when it
was discovered. From equation (3.13), we can obtain the present ratio of
disintegration of C 14 in the archaeological sample, given by
dA M (t) 1 M (0)
M (t) = − = λA0 e−λt ⇒ = e−λt ⇒ t = log , (3.14)
dt M (0) λ M (t)
M (0) (= λA0 ), being the original rate of disintegration. From (3.14), we can
determine the age, provided we can measure M (t) and M (0). The rate of
disintegration of C 14 in a comparable amount of archaeological sample or
fossil (living wood) is equal to M (0).
The half-life of the sample under investigation can be obtained from
equation (3.14). Since the half-life is the amount of time required by the
decreasing substance to reduce to half, we get
A0 1 0.6931
= A0 e−λτ ⇒ τ = loge 2 = , where τ is the half-life.
2 λ λ
Example 3.5.1 A fossil is found that has 20% C 14 compared to the living
sample. How old is the fossil, knowing that the C 14 half-life is 5730 years?

Solution: Let A(t) be the amount of C 14 present in the sample at any

time t, then
A(t) = A0 e−λt (following exponential decay),

ISTUDY
 140 Mathematical Modeling: Models, Analysis and Applications

where A0 is the initial amount of C 14 present in the sample and λ is the decay
constant of the sample. The fossil contains 20% of the original amount at time
t, which implies

A(t) 20 20 1
= ⇒ e−λt = ⇒ t = ln(5). (3.15)
A0 100 100 λ

Since the half-life is the amount of time required by the decaying substance
to reduce to half, we get,

A0
= A0 e−λτ , where τ is the half-life.
2  
1 1 τ
⇒ − λτ = ln ⇒ = . (3.16)
2 λ ln(2)

From (3.15) and (3.16), we obtain,

ln(5) 1.6094
t=τ = 5730 × ≃ 13305 years.
ln(2) 0.6931

3.5.2 Drug Distribution in the Body

The study of the movement of drugs in the body is called pharmacokinetics.
The science of pharmacokinetics uses mathematical equations and utilizes
them to describe the movement of the drug through the body [65].

We now study a simple problem in pharmacology, where we will be dealing

with the dose-response relationship of a drug. In this problem, the drug present
in the system follows certain laws. Let us assume that the rate of decrease
of the concentration of the drug is directly proportional to the square of its
amount present in the body and C0 be the initial dose of the drug given to the
patient at time t = 0. The mathematical model that captures this dynamics
is given by
dC(t)
= −kC 2 , (3.17)
dt
where k is a constant depending on the drug used and its value can be obtained
from the experiment. Solving (3.17), we get

C0
C(t) = , where C(0) = C0 .
1 + C0 kt

Let an equal dose of drug C0 be given to the body at equal time intervals, T .
Then, immediately after the second dose, the concentration of the drug inside
the body is
C0
C1 = C0 + .
1 + C0 kT

ISTUDY
 Continuous Models Using Ordinary Differential Equations 141

Immediately after the third dose, the concentration of the drug inside the
body is
C1
C2 = C0 + .
1 + C1 kT
In a similar manner, we can conclude that
Cn−1
Cn = C0 + , (3.18)
1 + Cn−1 kT

which is a non-linear difference equation. Now,

Cn − Cn−1
Cn+1 − Cn = . (3.19)
(1 + kT Cn )(1 + kT Cn−1 )

From (3.18), we conclude that Cn > C0 , which implies Cn+1 − Cn and Cn −

Cn−1 have the same sign. Noting that Cn is an increasing function of n, we
attempt to find the limiting value of the concentration by taking limits on
both sides of (3.18), that is,
 
Cn−1
lim Cn = lim C0 + ,
t→∞ t→∞ 1 + Cn−1 kT

lim Cn−1
t→∞
⇒ lim Cn = C0 + ,
t→∞ 1 + kT lim Cn−1
t→∞

C∞
⇒ C∞ = C0 + where C∞ = lim Cn = lim Cn−1 ,
1 + C∞ kT t→∞ t→∞

2
⇒ kT C∞ − kT C0 C∞ − C0 = 0,
p
kT C0 ± k 2 T 2 C02 + 4C0 kT
⇒ C∞ = ,
2kT
r
C0 C0 4
⇒ C∞ = + 1+ (taking positive sign only).
2 2 C0 kT
This implies C0 < Cn < C∞ (since Cn is an increasing function), that is, the
concentration is bounded.

3.5.3 Growth and Decay of Current in an L-R Circuit

We consider an L-R circuit where L is the inductance of the coil and R is
the resistance. The coil is connected to a battery of voltage V through a key
K (fig. 3.7).
In the ON position, the current flows through the coil. When the current
i(t) starts to flow, the negative lines of force move outward from the coil and
an electromotive force (e.m.f.) will induce across L. According to the law of
electromagnetic induction, this e.m.f. will oppose the voltage, as a result of

ISTUDY
 142 Mathematical Modeling: Models, Analysis and Applications

R L

V
K

FIGURE 3.7: The inductance-resistance (L-R) circuit, connected to a battery

of voltage V through a key K.

which there will be a voltage drop across R, which will also oppose the applied
voltage. Let, at any time t, i be the current in the circuit increasing from 0 to
a maximum value at a rate of increase di dt . Now, the potential difference across
di
the inductor is V1 = L dt and across the resistor is V2 = iR. The differential
equation modeling of this scenario is given by

di
V = L + iR (since V = V1 + V2 ),
dt
Z i Z  
di R t i R
⇒ V
= − dt ⇒ log e i − = − t,
0 i − R
L 0 V /R L
V  R

⇒ i(t) = 1 − e− L t ,
R
V
which shows that the current grows exponentially. As t → ∞, i → R = I
(say), a steady value (fig. 3.8(a)).

Current Current
60
50
50
40
40
30
30
20 20
10 10

Time 0 Time
100 200 300 400 500 0 100 200 300 400 500

(a) Current grows and reaches a (b) Current decays to zero.

steady value.

FIGURE 3.8: Graphs showing the (a) growth and (b) decay of current, with
L = 50, R = 1 and V = 50.

We now put the key in the OFF position. Initially, when the key was in
the ON position, a steady current I = VR was flowing. With no current flowing
in the circuit, the flux will reduce gradually, resulting in a voltage drop iR
di
across the resistance R and the induced e.m.f. L dt across the inductance L.

ISTUDY
 Continuous Models Using Ordinary Differential Equations 143

Now, since the key is OFF, the current becomes open, implying that the
impressed voltage is zero.
The differential equation showing this scenario is given by

di
0 = L + iR (since V = V1 + V2 = 0),
dt
Z i Z t
di R
⇒ =− dt (since at t = 0, i = I),
i 0 L
I
i R V −Rt
ln = − t ⇒ i(t) = e L .
I L R

Thus, the current decays exponentially as time increases and ultimately goes
to zero (fig. 3.8(b)).

3.5.4 Rectilinear Motion under Variable Force

x
P

O F a A

FIGURE 3.9: A particle moving in a straight line towards the origin O

(fixed) and acted upon by a force F directed towards O.

Let a particle move in a straight line and be acted upon by a force F =

mµ
, µ (> 0) being the constant of proportionality, which is always directed
x
towards a fixed point O. Here, m is the mass of the particle and P is the
position of the particle at any time t such that OP = x (fig. 3.9). The equation
of motion modeling the given scenario is given by

dv µ
mv = −m ,
dx x
dv
where v dx is the acceleration of the particle of mass m. Since the force is
attractive, the sign of right-hand side is negative. Integrating, we get,

v2
= −µ ln x + constant.
2

Initially, let the particle start from rest (from A) at a distance a from the
fixed point O (origin), then at time t = 0, x = a, v = 0,
a
⇒ constant = µ ln a ⇒ v 2 = 2µ ln .
x

ISTUDY
 144 Mathematical Modeling: Models, Analysis and Applications
r a
dx
⇒v= = − 2µ ln (negative sign as distance decreases with time).
dt x
p Z T Z 0
dx
⇒ 2µ dt = − q
, where T is the time taken by the particle
0 a ln xa
to reach the origin. Therefore,
Z ∞ −y2
p ae (−2y)dy  a 
2µT = − Substitute ln = y2
0 y x
Z ∞ Z ∞ −z
2 e
= 2a e−y dy = a √ dz (Put z = y 2 )
0 0 z
Z ∞   r
−z 12 −1 1 √ π
=a e z dz = a Γ =a π ⇒T =a .
0 2 2µ

3.5.5 Mechanical Oscillations

In this section, we discuss three different oscillations, namely, horizontal
oscillations, vertical oscillations and damped forced oscillations.

[Link] Horizontal Oscillations

d1 d2

O1 A a P O2

FIGURE 3.10: A particle at rest is attracted by two forces towards two fixed
centers O1 and O2 , respectively.

We consider a particle of mass m resting in equilibrium at a point A,

being attracted by two forces equal to mµn1 ×(distance)n and mµn2 ×(distance)n
respectively, towards two fixed centers O1 and O2 (fig. 3.10).
In this equilibrium position of the particle at A, the forces acting on it are
equal and opposite, and hence O1 , A and O2 are collinear.
Let O1 O2 = a, O1 A = d1 and AO2 = d2 . Then d1 + d2 = a and µn1 dn1 =
n n
µ2 d2 , since the forces are equal and opposite at A.

µn1 dn1 = µn2 dn2 ⇒ µ1 d1 = µ2 d2 ,

µ1 µ2 µ1 + µ2 µ1 + µ2
⇒ = = = ,
d2 d1 d1 + d2 a
aµ2 aµ1
⇒ d1 = and d2 = .
µ1 + µ2 µ1 + µ2
Now, let the particle be slightly displaced towards the fixed center O2 and P
be the new position of the particle at any time t, such that AP = x (small).

ISTUDY
 Continuous Models Using Ordinary Differential Equations 145

The differential equation modeling of this physical scenario is given by

d2 x
m = −mµn1 O1 P n + mµn2 O2 P n . (3.20)
dt2
When the particle is slightly displaced from the equilibrium position A to
the position P , it tends to come back to its original position. Doing so, O1
becomes the origin of attraction and O2 , the origin of repulsion. This explains
the reason for the first term on the right-hand side of (3.20) to be negative
and the second term to be positive.

d2 x
Therefore, = −µn1 (d1 + x)n + µn2 (d2 − x)n ,
dt2  n  n
x x
= −µn1 dn1 1 + + µn2 dn2 1 − ,
d1 d2
 
nx n(n − 1) x2
= −µn1 dn1 1 + + + ...
d1 2! d21
 
nx n(n − 1) x2
+ µn2 dn2 1 − + + ... ,
d2 2! d22
d2 x
or, = −µn1 dn1 − nxµn1 d1n−1 + µn2 dn2 − nxµn2 d2n−1 ,
dt2
(Since x is small, neglecting higher powers of x)
2

d x
⇒ = −n µn1 d1n−1 + µn2 d2n−1 x,
dt2 (  n−1  n−1 )
n aµ2 n aµ1
= −n µ1 + µ2 x,
µ1 + µ2 µ1 + µ2
nµ1n−1 µ2n−1 an−1
= − (µ1 + µ2 )x,
(µ1 + µ2 )n−1
d2 x (µ1 µ2 a)n−1
⇒ = −n x.
dt2 (µ1 + µ2 )n−2

Therefore, we conclude that the motion of the particle is simple harmonic

about A and the period of oscillation is
s
(µ1 + µ2 )n−2
2π .
(µ1 µ2 a)n−1

[Link] Vertical Oscillations

Consider an elastic string of upstretched length AB (= a), fixed at point
A. A particle of mass m is attached to the end of the string so that the string is
extended to the length AO (= b) when the mass is at rest. In the equilibrium
position, the tension in the string is balanced by the weight of the particle,

ISTUDY
 146 Mathematical Modeling: Models, Analysis and Applications

that is,  
AO − AB
mg = T0 = λ (by Hooke’s law)
AB
 
b−a
⇒ mg = λ , λ is the modulus of elasticity. (3.21)
a

Now, the particle is pulled to the point D and A

then released. Let P be the position of the particle

at any time t such that OP = x (fig. 3.11). The
a
differential equation that models the scenario is b

given by
d2 x B
m = mg − T, where T is the tension of the string.
dt2 O
d2 x b+x−a b − a λx
⇒ m 2 = mg − λ = mg − λ − , -
dt a a a
d2 x λx P
⇒ =− (using (3.21)).
dt2 am D

Hence, the motion is simple harmonic aboutr the FIGURE 3.11:

am
center O, the period of oscillation being 2π . Vertical oscillation.
λ

[Link] Damped and Forced Oscillations

2kv

O
v P

FIGURE 3.12: A disturbing force 2kv acting on a particle moving in a

straight line towards a fixed point O.

We consider a particle moving in a straight line with an acceleration µ2

(distance) towards a fixed point O in the line. A disturbing acceleration of
2k (velocity)(µ > k) also acts on the particle along with a periodic additional
acceleration F cos(bt) [54] (fig. 3.12).

The equation of motion that models the scenario is given by

d2 x dx
= −µ2 x − 2k + F cos(bt).
dt2 dt
Let x = Aemt (A 6= 0) be a trial solution of
d2 x dx
2
+ 2k + µ2 x = 0,
dt dt

ISTUDY
 Continuous Models Using Ordinary Differential Equations 147

then the required auxiliary equation is

p p
m2 + 2mk + µ2 = 0 ⇒ m = −k ± k 2 − µ2 = −k ± i µ2 − k 2 (µ > k).
p
Complimentary function is e−kt A cos( µ2 − k 2 t + ε1 ) and the particular
integral is given by
1 D2 − 2Dk + µ2
F cos(bt) = F cos(bt)
D2 + 2Dk + µ2 (D2 + µ2 )2 − 4D2 k 2
(µ2 − b2 ) cos(bt) + 2kb sin(bt)
=F = B cos(bt − ε2 ),
(µ2 − b2 )2 + 4k 2 b2
F 2kb
where B = p and tan ε2 = 2 .
2 2 2
(µ − b ) + 4k b 2 2 µ − b2

Therefore, the general solution is

p
x = Ae−kt cos( µ2 − k 2 t + ε1 ) + B cos(bt − ε2 ). (3.22)

From (3.22), it is concluded that motion is the resultant of two oscillations,

2.0
200

/:5
100
1.0

.:5
;< 100 =;< 200 >;<
5. 100 /5. 200 65.
-100
-.:5
-1.0 -200

(a) k = 0.1, F = 4, b = 0.1, µ = 2.0. (b) k = 0.1, F = 4, b = µ = 0.1.

FIGURE 3.13: A motion which is the resultant of two oscillations, namely,

a free oscillation and a forced one.

namely, the free oscillation (first part) and forced oscillation (second part).
The arbitrary constants A and B can be obtained from the initial conditions.
From the expression (3.22), it is clear that the amplitude of free oscillation
decreases with time t because of the factor e−kt and ultimately vanishes for
large t. However, the amplitude of the forced oscillation persists as there is no
diminishing factor, whose period of oscillation is 2π
b . This is also evident from
fig. 3.13(a).

Special Case: If the period of forced oscillation is equal to the period of

free oscillation, that is, 2π 2π
b = µ ⇒ b = µ, then the amplitude of the forced
F
oscillation is B = . If k is small, then the amplitude of the forced oscillation
2kb

ISTUDY