Malawi University of Science and Technology

Malawi Institute of Technology

Introduction to Mathematical Modeling

MATH 2107

Module Writer

Precious Juwawo

Module Reviewer

Date
January 2024

1
 TABLE OF CONTENTS

TABLE OF CONTENTS . . . . . . . . . . . . . . . . . . . . . 4

1 INTRODUCTION 1
1.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Mathematical Modeling . . . . . . . . . . . . . . . . . . 3
1.2.1 Examples of Mathematical models . . . . . . . . 5
1.2.2 Why Mathematical Modeling . . . . . . . . . . 5
1.2.3 The Modeling process . . . . . . . . . . . . . . 6
1.2.4 Classification of Mathematical models . . . . . 7
1.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4 End of unit test . . . . . . . . . . . . . . . . . . . . . . 9
1.5 Further Reading . . . . . . . . . . . . . . . . . . . . . . 9

2 DIFFERENCE EQUATIONS 11
2.1 Examples of difference equations . . . . . . . . . . . . . 13
2.2 Solving Difference equations by iteration . . . . . . . . 14
2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 Further Reading . . . . . . . . . . . . . . . . . . . . . . 19

3 MODELING WITH DIFFERENCE EQUATIONS (DISCRETE-

TIME MODELS) 20
3.1 Discrete-time models . . . . . . . . . . . . . . . . . . . 22

2
 3.1.1 Examples of discrete time models . . . . . . . . 23
3.1.2 Solutions and qualitative analysis of discrete-
time models . . . . . . . . . . . . . . . . . . . . 25
3.1.3 Fixed points/ steady states of discrete-time models 26
3.1.4 Stability criterion . . . . . . . . . . . . . . . . . 27
3.1.5 Model analysis . . . . . . . . . . . . . . . . . . 30
3.1.6 Application of discrete-time models to finance . 31
3.2 Systems of difference equation models . . . . . . . . . . 33
3.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.4 End of unit test . . . . . . . . . . . . . . . . . . . . . . 38
3.5 Further Reading . . . . . . . . . . . . . . . . . . . . . . 39

4 CONTINUOUS-TIME MODELS 40
4.1 Differential equations in Mathematical modeling . . . . 42
4.1.1 Solving differential equations by separation of
variables . . . . . . . . . . . . . . . . . . . . . . 42
4.2 Continuous-time models: Motivation and derivation . . 45
4.2.1 Examples of Continuous models . . . . . . . . . 46
4.2.2 Solution to the logistic population growth model 48
4.2.3 Stability Criterion of continuous models . . . . 49
4.3 Differential Equation models . . . . . . . . . . . . . . . 51
4.4 Systems of differential equation models . . . . . . . . . 52

3
 4.4.1 Analytical methods for systems of differential
equation models. . . . . . . . . . . . . . . . . . 56
4.4.2 Equilibrium points (or fixed points or steady states) 57
4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.6 Further Reading . . . . . . . . . . . . . . . . . . . . . . 64

4
unit 1

INTRODUCTION

Introduction
This chapter provides definitions of some terms used in Mathematical
modeling. We also define the term Mathematical modeling and pro-
vide some examples of Mathematical models. The modeling process
and types of Mathematical models are also studied.

Intended learning outcomes

By the end of this unit you should be able to;

1
 a . Define Mathematical modeling

b . Give examples of Mathematical models

c . Describe the process of Mathematical modeling

d . Formulate Mathematical models

Key terms
Ensure that you are familiar with the following terms;
Modeling, discrete, continuous, static, phenomena, empirical.

1.1 Preliminaries
The difficulty of problems treated in science and engineering is due
to the complex nature of the systems under consideration. Models
provide an adequate tool to breakup this complexity and make the
problem easy to deal with. A model is something that mimics relevant
features of the situation being studied.

The usefulness of model is how well it performs when applied to the

problems it was designed to handle. When a model is used, it may
lead to incorrect predictions. At that stage, the model may either
be modified, discarded, or used anyway. The best model is the sim-
plest model that still serves its purpose, that is, which is still complex

2
enough to help use understand a system and to solve problems.

Modelling is an iterative process to obtain a better understanding of

some observation from the real world. The Real world is abstracted
into a symbolic idea, which is expressed as mathematical equations
or a mathematical model. Empirical data are collected about the sys-
tem of interest, and these data are compared to the output from the
mathematical model. No model is THE right one, but only A right
one.

1.2 Mathematical Modeling

Now that you know the meaning of the word model-

ing, think of how this can help you to define the word Mathematical
modeling.

A mathematical model is a representation of a real system using math-

ematical language. Examples include Algebraic Equations, Differen-
tial Equations, Integral equations, algorithms, Formulae, and theo-
rems. All mathematical models are wrong, but some models are use-
ful to help us understand real-world systems.

The essence of good mathematical model is that it is simple in design

and exhibits the basic properties of the real system. That is, a good

3
mathematical model;

(a) Should be testable against empirical data.

(b) Should iteratively lead to improved mathematical models.

Modelling is an evolutionary process.

Figure 1.1: The modeling diagram

Definition 1.1. (Mathematical Modelling). Mathematical modelling is

a process that uses mathematics to represent, analyse, make predictions
or otherwise provide insight into real-world phenomena

4
1.2.1 Examples of Mathematical models

(a) Newton’s law F = ma, where F is force, m is mass, a is accelera-

tion

(b) Ohm’s law V = IR, where V is voltage, I is current, R is resis-

tance.

(c) Kepler’s third law T = cR3/2 , where T is the orbital period of the
planet, R is the mean distance from the planet to the sun.

(d) Einstein’s relativity theory E = M c2 , where E is energy, M is

mass, c is light speed.

(e) Metabolic theory of ecology B = B0 M 3/4 , where B is organism

metabolic rate, B0 is a mass-independent normalization constant,
M is organism mass.

dN N


(f) Logistic population model dt = rN 1 − K , where N is popula-
tion size, r is the maximum per capita growth rate, K is carrying
capacity.

1.2.2 Why Mathematical Modeling

(a) To gain understanding

(i) Through analysis of the model.

5
 (ii) By identifying factors that are most important in the system

(b) To predict or simulate Very often we wish to know how the real-
world system will do in the future, but it is expensive, impractical
and impossible to experiment directly with the system.

(c) To optimise some performance.

1.2.3 The Modeling process

Figure 1.2: The Modelling process

(a) Step 1. Identify the problem with specific goals and questions.

(b) Step 2. Formulate/Construct the model.

(i) Make assumptions

(ii) Define variables

6
 (iii) Determine the relationship among the variables defined.

(c) Step 3. Solve the Model

(i) Analyse the model

(ii) Simulate the model.

(d) Step 4. Validate the model with real phenomena or empirical

data.

(i) Does it address the problem?

(ii) Does it make sense?

(iii) Test it with real data.

(e) Step 5. Implement the model - Apply the model to make predic-
tions.

(f) Step 6. Calibrate and extend the model

1.2.4 Classification of Mathematical models

Mathematical models are classified according to several criteria. We

look at the following models;

(a) Static vs Dynamic Dynamic models: Models where state vari-

7
 ables and relationships describing the system change with time.
Thus, time-dependent.
Static modelS: models where the time factor is not taken into
account. Hence the variables and relationships describing the
system are time-independent.

(b) Stochastic vs Deterministic

Stochastic models are models where values assumed by the vari-
ables or changes to the variables are not predictable with cer-
tainty. Thus, random effects are involved. For example, rolling a
dice, share prices, and so on.
Deterministic models: these does not depend on random ef-
fects. They predict outcomes from a given starting point. In other
words, there is certainty in prediction of values of variables or
changes to the variables.

(c) Discrete-time vs Continuous-time Discrete-time models: time

changes in incremental steps.
Continuous-time models: time changes continuously.

(d) Linear Model vs Non-linear Models Linear models involve linear

relationships amongst the variables of the system, hence has lin-
ear equations.
Non-linear models involve non-linear relationships among the
variables of the system, hence has non-linear equations.

8
(e) Simulation models are models that are solved on a computer due
to the complex nature which can not be solved analytically.

1.3 Summary
In this unit we have laid a foundation to the study of Mathematical
modeling. We have defined useful terms in this topic and given a good
classification of Mathematical modeling.

1.4 End of unit test

a . Come up with 3 different phenomena in real life which can be

modeled using Mathematically.

b . Give five more different examples of Mathematical models.

1.5 Further Reading

1. Wang H. (2012) Mathematical Modeling I - preliminary. ebook.

9
2. Giordano F. A first course in Mathematical modeling. 5th edition.
Cengage Learning.

10
unit 2

DIFFERENCE
EQUATIONS

Introduction
In this unit we are going to study difference equations and methods of
finding solutions to difference equations. This will help us to under-
stand Mathematical models that use difference equations in the next
unit.

11
 Learning outcomes
By the end of this unit you should be able to:

a . Define difference equations

b . Solve difference equations by iteration

Key terms
Make sure you familiarize yourself withe the following terms in this
unit;
Sequence, relations, iterations.

Definition 2.1. A difference equation for the sequence {an } is an equa-

tion that expresses an interms of one or more of the previous terms
a0 , a1 , . . . an−1 for integers n with n ≥ n0 where n0 is a nonnegative
integer.

Difference equations are also known as recurrence relations. A se-

quence is called a solution of recurrence relation if its terms satisfy
the recurrence relation. The Initial condition is required to precede
the first terms where relations take effect.

12
2.1 Examples of difference equations

(a) Let an = an−1 +6 for n ≥ 1 given that a0 = 3. Find a1 , a2 , a3 . . . a100

Solution:
a0 = 3
a1 = a0 + 6 = 3 + 6 = 9
a2 = a1 + 6 = 9 + 6 = 15
a3 = a2 + 6 = 15 + 6 = 21

We can an equation a+nd. Since

a4 = 3 + 4 · 6 = 27. Therefore
a100 = a + n · d = 3 + (100 × 6)
= 603

(b) If u1 = 4 and un = 2un−1 + 3n − 1, for n ≥ 2, find the values of u2

and u3 .

Solution
u2 = 2u1 + 3 × 2 − 1
=2×4+6−1
=8+6−1
= 13

13
 and
u3 = 2u2 + 3 × 3 − 1
= 2 × 13 + 9 − 1
= 26 + 9 − 1
= 34

2.2 Solving Difference equations by iteration

When solving recurrence relation we aim at finding a non-recursive
formula to calculate an where the solution is called a closed formula.
We can also use Iteration method involving substitution and back
tracking methods to solve these recurrence relation. In the follow-
ing example we are going to apply these methods.

Example 1: Let an be the sequence that satisfies the recurrence re-

lation an = an−1 + 3 for n ≥ 1 with initial condition a0 = 2. Use
Iteration method involving substitution or back tracking method to
solve for a1 , a2 , a3 , a4 . . . an

14
Solution: By Iteration method involving substitution we have

a0 = 2
a1 = 2 + 3
a2 = a1 + 3 = (2 + 3) + 3 = 2(3 + 3) = 2 + 2(3)
a3 = a2 + 3 = 2 + 2(3) + 3 = 2 + 3(3)
a4 = 2 + 4(3)
..
.
an = 2 + n(3)
an = 2 + 3n

The use of back tracking is left for you as an exercise.

Example 2. Solve
un = 2un−1 + 1, n≥2

For example,
u6 = 2u5 + 1
= 2 × 15 + 1
= 31.

In a similar way, un−1 can be written in terms of un−2 as

un−1 = 2un−2 + 1

15
Then equation (2) can be substituted into equation (1) to give

un = 2 (2un−2 + 1) + 1
= 4un−2 + 2 + 1

Repeating this process using un−2 = 2un−3 + 1 gives

un = 4 (2un−3 + 1) + 2 + 1
= 8un−3 + 4 + 2 + 1
⇒ un = 23 un−3 + 22 + 21 + 20 .

You can see a pattern developing. Continuing until un is expressed in

terms of u1 , gives

un = 2n−1 u1 + 2n−2 + 2n−3 + . . . + 22 + 21 + 20

= 2n−1 + 2n−2 + . . . + 22 + 21 + 20 , ( as u1 = 1)

We recognise the solution as a geometric progression (GP) with first

term 1 and common ratio 2 . Using the formula for the sum to n terms
of a GP,
1 (2n − 1)
un = = 2n − 1, n ≥ 1.
2−1

This is the solution to the differential equation un = 2un−1 + 1. This

process involved the repeated use of a formula and is known as itera-
tion.

16
You can now see how easy it is to calculate a value for un . For exam-
ple,
u100 = 2100 − 1
≈ 1.27 × 1030 .

Activity 1
Practice Questions
Find the next two elements of the following recurrence relations with
their initial conditions.

(a) an = 2an−1 an −2 + n2 ; a1 = 1, a2 = 2

(b) an = 2an−1 + 5 an −2 − 6nn−3 ; a1 = 1, a2 = 2 and a3 = 1

(c) an = nan−1 + 3an − 2; a1 = 1, a2 = 2

2.3 Summary
In this unit we have studied difference equations and methods of solv-
ing them. This will help us in studying modeling with difference equa-
tions in the next unit.

17
End of unit test

1 . For each equation you are given the first term of a sequence.
Find the 4 th term in each case;

a .u1 = 2 and un = un−1 + 3, n≥2

b . u1 = 1 and un = 3un−1 + n, n≥2

c . u1 = 0 and un − un−1 = n + 1, n ≥ 2.

2 . For each sequence write down a difference equation which

describes it:

a . 3 5 7 9 11

b . 251123 47

c . 125 1441.

3 . Solve by iteration, giving un in terms of u1 .

a . un = un−1 + 2, n≥2

b . un = 4un−1 − 1, n≥2

18
 c . un = 3un−1 + 2, n≥2

2.4 Further Reading

1. Wang H. (2012) Mathematical Modeling I - preliminary. ebook.

2. Giordano F. A first course in Mathematical modeling. 5th edition.

Cengage Learning.

19
unit 3

MODELING WITH
DIFFERENCE EQUATIONS
(DISCRETE-TIME
MODELS)

Introduction
In this unit we study Mathematical models that use difference equa-
tions. Methods of solving such models are going to be studied as well
as mathematical analysis of them.

20
 Learning outcomes
By the end of this unit you should be able to;

a . Formulate discrete-time models

b . Solve discrete-time models

c . Apply discrete-time models to solve predict population growth

d . Apply discrete-time models in finance

e . Analyze discrete-time models mathematically

f . Solve systems of discrete-time models

g . Analyse systems of discrete-time models mathematically

Key words
The following terms are useful in this unit;
Equilibrium, fixed points, stability, recurrence relation, eigenvalue.

21
3.1 Discrete-time models
Discrete-time models describe phenomena in terms of fixed-time steps.
The underlying assumption will always be that the population size or
the amount of quantity being modeled at each stage is determined by
the population sizes in the past generations, but that the intermediate
population sizes between generations are not needed.

A fundamental concept when working with discrete-time models is

called recurrence relation. Consider a sequence of quantities, x0 , x1 , x2 , . . .,
where xi denotes the quantity after i time steps. If xn+1 depends only
on xn , a discrete-time model is expressed by

xn+1 = f (xn ) n = 0, 1, 2, 3, . . .

with some initial condition x0 . This discrete-time model is called a

difference equation. This discrete model gives

x1 = f (x0 )
x2 = f (x1 ) = f (f (x0 )) = f [2] (x0 )
x3 = f (x2 ) = f (f (x1 )) = f (f (f (x0 ))) = f [3] (x0 )
..
.
xn = f (xn−1 ) = f (f (xn−2 )) = · · · = f (f (f (· · · f (x0 ) · · · ))) = f [n] (x0 )

The resulting sequence x0 , x1 , x3 , · · · , is called an orbit of the map.

22
3.1.1 Examples of discrete time models

An example - bacterial reproduction

Bacterial cells divide into more cells after one sampling time.

The number of bacterial cells in the next measurement will be some

multiple of the current number. We assume that this multiple is a
constant over several sample times. Note that this assumption is ob-
viously invalid for many sampling times due to resource and space
limitations. From the scenario above we can develop a dicrete-time
model and solve it. Denote Bn as the cell number observed at the nth
sampling time, then the model can be written as;

Bn+1 = rBn (Malthus model, 1798)

where the constant r is called the growth rate.

Solution: Given the initial cell number B0 ,

B1 = rB0
B2 = rB1 = r2 B0
B3 = rB2 = r3 B0
..
.
Bn = rn B0

which is the solution of this discrete model.

23
In reality, the growth rate r usually depends on the cell number be-
cause of competition for resource and space. r = r (Bn ) is a decreas-
ing function of Bn .
The more general model than Malthus model is of the form

Bn+1 = r (Bn ) Bn .

The solution is a simple linear difference equation that involves an ex-

pression of the form (some constant) n , where n is a generation num-
[Link] is true in general for linear difference [Link] magni-
tude of r determines whether the population grows or dwindles with
time.
-When |r| > 0, Mn increases over successive generations.
-When |r| < 0 Mn decreases over successive generations.
- When r = 1 Mn is constant.
The above model is very simple. It assumes unlimited space and re-
sources for the culture.

Suppose this models bacterial cells dividing. We expect competition

for space and [Link] reality, a depends on the cell number be-
cause competition for space and [Link], r = r (Mn ) is a de-
creasing function of Mn . The more general model than Malthus model
is of the form
Mn+1 = r (Mn ) Mn

2
For instance, Verhulst model assumes that r (Mn ) = 1+ MKn
For instance,

24
Verhulst model assumes

2
r (Mn ) =
1 + MKn

and thus the model becomes

2Mn
Mn+1 = .
1 + MKn

The function r (Mn ) is a decreasing function of Mn with the maximum

2 occurring at Mn = 0. When Mn = K, r (Mn ) = 1, one-half of its
maximum.

3.1.2 Solutions and qualitative analysis of discrete-

time models

We know that the solution of a linear discrete-time model xn+1 =

rxn , r ∈ R is xn = rn x0 , n = 0, 1, 2, · · · , as shown [Link] us investi-
gate the long time behavior as n → ∞.
- If r > 1, then limn→∞ xn = limn→∞ rn x0 = ∞.
- If r = 1, then limn→∞ xn = limn→∞ x0 = x0 .

- If −1 < r < 1, then limn→∞ xn = limn→∞ rn x0 = 0.

- If r = −1, then

x0 n even

lim =
n→∞
−x0 ,

n odd

25
and thus there is no convergence.

- If r < −1, then limn→∞ |xn | = limn→∞ |r|n |x0 | = ∞, thus no conver-
gence.
In addition rn is positive for n even, and negative for n odd.

3.1.3 Fixed points/ steady states of discrete-time mod-

els

Definition 3.1. ; A number x∗ is called an equilibrium point or fixed

point of xn+1 = f (xn ), if xn = x∗ for all n = 1, 2, 3, · · · when x0 =
x∗ .That is, xn = x∗ is a constant solution to the discrete model.

For our first model,

Mn+1 = rMn

an equilibrium point satisfies M ∗ = rM ∗ .

- if r ̸= 1, M ∗ = 0 is the only equilibrium point.

- If r = 1, every number is an equilibrium point.

Now, let us consider the following complicated discrete model

xn+1 = rxn + b

26
where r and b are constants. An equilibrium x∗ satisfies the

x∗ = rx∗ + b

- If r ̸= 0, x∗ = b
1−r

- If r = 1, then x∗ = x∗ + b.

-Every point is an equilibrium if b = 0, and no equilibrium point

exists if b ̸= 0 As a summary, results for equilibria are listed below:

- x∗ = b
1−r is the only equilibrium if r ̸= 1.

- Every number is an equilibrium if r = 1 and b = 0.

- no equilibrium exists if r = 1 and b ̸= 0.

3.1.4 Stability criterion

Theorem 1. Theorem (Stability Criterion):Let x∗ be an equilibrium of

xn+1 = f (xn ), then we have the results:

(a) If |f ′ (x∗ )| < 1, x∗ is stable.

(b) If |f ′ (x∗ )| > 1, x∗ is unstable.

(c) If |f ′ (x∗ )| = 1, there is no conclusion about the stability of x∗ .

Higher order terms need to be examined to determine stability.

The constant f ′ (x∗ ) is called the eigenvalue of the map f at x∗ .

27
Examples
Example 1: Consider the discrete logistic equation yn+1 = ryn (1 − yn /K),
where the parameters r ≥ 0, K > 0. Find all equilibria and determine
their stability.
1. Solution: Let xn = yn /K, then xn+1 = rxn (1 − xn ).
Note that xn ≥ 0 represents population size, and r ≥ 0 is the maxi-

mum growth rate. Equilibria satisfy

x∗ = rx∗ (1 − x∗ )

which leads to x∗ = 0 or x∗ = r−1

r .

The nontrivial equilibrium x∗ = r−1

r > 0 if and only if r > 1.

To check stability, we compute the eigenvalue:

f ′ (x∗ ) = r (1 − 2x∗ )

For the trivial equilibrium x∗ = 0, the eigenvalue f ′ (0) = r, thus

x∗ = 0 is stable if 0 ≤ r < 1 and unstable if r > 1.

For the nontrivial equilibrium x∗ = r−1

the eigenvalue f ′ r−1


r , r =
2 − r, thus x∗ = r−1
r is stable if 1 < r < 3 (from |2 − r| < 1 ) and

28
unstable if r > 3 (from |2 − r| > 1 and r > 1 ). Example 2: Consider
rxn
the Beverton-Holt model xn+1 = 1+ r−1
with r > 0 and K > 0. Find
K xn

all equilibria and determine their stability.

[Link]: Equilibria x∗ satisfy

∗ rx∗
x =
1 + r−1
K x
∗

which leads to x∗ = 0 or x∗ = K (if r ̸= 1 ). Note that for r = 1 the

model becomes xn+1 = xn , which is not interesting.

To check stability, we compute the eigenvalue

r
f ′ (x∗ ) =
.
r−1 ∗ 2
1+ K x

For the trivial equilibrium x∗ = 0, the eigenvalue f ′ (0) = r, thus

x∗ = 0 is stable if 0 < r < 1 and unstable if r > 1.

For the nontrivial equilibrium x∗ = K, the eigenvalue f ′ (K) = 1r , thus

x∗ = K is stable if r > 1 and unstable if 0 < r < 1.

Activity 2
Practice Questions
For the following models, find all equilibria and determine their sta-

29
bility
(1)Bn+1 = rBn
2Bn
(2)Bn+1 =
1 + BKn
rBn2
(3)Bn+1 =
2
1 + BKn

3.1.5 Model analysis

Analytical Solution to xn+1 = rxn + b

xn = rxn−1 + b
= r (xn−2 + b) + b
= r2 xn−2 + b(1 + r)
= r2 (rxn−3 + b) + b(1 + r)
= r3 xn−3 + b 1 + r + r2

..
.
= rn x0 + b 1 + r + r2 + r3 + · · · + rn−1

n (1 − rn )
= r x0 + b ( if r ̸= 1)
1−r
1 b n
= rn x0 + − r
1−r 1−r

If r = 1, xn+1 = xn + b. Thus,

xn = xn−1 + b = xn−2 + 2b = xn−3 + 3b = · · · = x0 + nb

30
Long term behaviour of model can be discussed for b ̸= 0.
b 1


(1) If r ̸= 1, then xn = rn x0 − 1−r + 1−r , there are two cases

b
- When |r| < 1, xn → 1−r as n → ∞.

- When |r| > 1, xn is not convergent as n → ∞

(2) If r = 1, xn = x0 + nb.

- If b > 0, xn → ∞ as n → ∞.

- If b < 0, xn → −∞ as n → ∞

3.1.6 Application of discrete-time models to finance

Example :Consider the annuity for retirement with 0.5% as the monthly
interest rate and a monthly withdrawal of $2000. Develop a discrete-
time model to describe the annuity problem. Determine equilibria
and stability. How much of an initial deposit is needed to deplete the
annuity in 30 years?
Solution: Let xn be the amount in the account after n months, then

the discrete-time model is provided by

xn+1 = xn (1 + 0.5%) − 2000

which can be simplified as

xn+1 = 1.005xn − 2000

31
with the initial deposit x0 .

Equilibria x∗ satisfy x∗ = 1.005x∗ − 2000 which leads to x∗ = 400000.

The eigenvalue is f ′ (x∗ ) = 1.005 > 1, thus the only equilibrium x∗ =
400000 is unstable. Recall that the solution of xn+1 = rxn + b is
 
b b
xn = r n x0 − + if r ̸= 1
1−r 1−r
xn = x0 + nb if r = 1.

In this example, r = 1.005 ̸= 1 and b = −2000, thus we use the first

formula to obtain
 
−2000 −2000
0 = a360 = 1.005360 x0 − +
1 − 1.005 1 − 1.005

from which we solve for the initial deposit: x0 = 333580.

As a conclusion, an initial deposit of $333580 allows the withdrawal

of $2000 per month from the account from 30 years. The total with-
drawal is $720000, and at the end of 30 years the account is de-
pleted.

32
3.2 Systems of difference equation models
We consider a two-dimensional discrete-time system:

xn+1 = f (xn , yn ) ,
yn+1 = g (xn , yn ) .

The qualitative analysis of a higher dimensional discrete-time system

is, although more complicated, similar to a two-dimensional discrete-
time system.

Same logic as the scalar case, equilibria (x∗ , y ∗ ) satisfy

x∗ = f (x∗ , y ∗ ) ,
y ∗ = g (x∗ , y ∗ ) .

Stability of an equilibrium (x∗ , y ∗ ) can be determined by the following

theorem:

Theorem 2. (Linear stability analysis): Let (x∗ , y ∗ ) be an equilibrium

of
xn+1 = f (xn , yn ) ,
yn+1 = g (xn , yn ) ,

and f, g are at least twice continuously differentiable. Let

 
∂f ∗ ∗ ∂f ∗ ∗
(x , y ) (x , y )
J (x∗ , y ∗ ) =  ∂xn ∂yn 
∂g ∗ ∗ ∂g ∗ ∗
∂xn (x , y ) ∂yn (x , y )

33
  
f
be the Jacobian matrix of  , evaluated at (x∗ , y ∗ ). Then (x∗ , y ∗ )
g
is stable if all eigenvalues of J have magnitude less than one; (x∗ , y ∗ ) is
unstable if at least one of the eigenvalues has magnitude greater than
one.

Notice the following;

(a) Eigenvalues λ of J are obtained from the characteristic equation

det(J − λI) = 0 where I denotes the identity matrix.

(b) Magnitude of a real eigenvalue is absolute value, while magni-

√
tude of a complex eigenvalue |a + bi| = a2 + b2 .

(c) This stability theorem can be easily extended to a higher dimen-

sional system. For a system of m difference equations, the Ja-
cobian matrix will be m × m, and there will be m eigenvalues
(counting multiplicity).

(d) For a two-dimensional system, the characteristic equation det(J −

λI) = 0 is equivalent to λ2 − (tr J)λ + det J = 0. We can show
that | tr J| < 1 + det J < 2 (Jury conditions) are sufficient and
necessary conditions for all eigenvalues of J to have magnitude
less than one, then the equilibrium (x∗ , y ∗ ) is stable.

If additionally (tr J)2 − 4 det J > 0, then (x∗ , y ∗ ) is a stable node (real
eigenvalues); and if additionally (tr J)2 − 4 det J < 0, then (x∗ , y ∗ ) is

34
a stable spiral (complex eigenvalues).

Example
We consider love affairs between Romeo and Juliet. Let xn be Romeo’s
love/hate for Juliet on day n, and let yn be Juliet’s love/hate for
Romeo on day n. xn > 0 implies Romeo loves Juliet, xn < 0 implies
Romeo hates Juliet, and xn = 0 implies Romeo is neutral to Juliet.
The larger |xn |, the stronger feeling of love/hate. Parallel assump-
tions hold for yn . A simple linear model is provided by

xn+1 = axn + byn ,

yn+1 = cxn + dyn ,

where a and d are response rates to their own feelings, b and c are re-
sponse rates to the feelings of the other. The sign of b or c determines
a particular romantic style. Analyze this linear model mathemati-
cally.

Solution: Equilibria (x∗ , y ∗ ) satisfy

x∗ = ax∗ + by ∗
y ∗ = cx∗ + dy ∗

(a − 1)x∗ + by ∗ = 0
cx∗ + (d − 1)y ∗ = 0

35
     
∗
a−1 b x 0
  = 
∗
c d−1 y 0
   
a−1 b a−1 b
Two cases: i) det   ̸= 0; ii) det   = 0.
c d−1 c d−1
 
a−1 b
We discuss these two cases one by one. i) If det   ̸= 0,
c d−1
then (x∗ , y ∗ ) = (0, 0) is the only equilibrium. Stability ofthis unique

a b
equilibrium is determined by the Jacobian matrix J =   ⇒
c d
tr J = a + d, det J = ad − bc. According to Jury conditions, (0, 0) is
stable if |a + d| < 1 + ad − bc < 2.
 
a−1 b a−1
If det   = 0, then (a − 1)(d − 1) − bc = 0, that is =
c
c d−1
b
d−1 , thus the two algebraic equations are identical, we only need to
∗ 1−a ∗
solve from one equation, e.g. the first equation: = b x , which
is the only condition for equilibria to satisfy. Hence, all points of the
form (x∗ , y ∗ ) = x∗ , 1−a ∗


b x are equilibria. All equilibria have the same
Jacobian matrix because the system is linear. Therefore, any equilib-
rium x∗ , 1−a ∗


b x is stable if |a+d|
 < |1+ad−bc 
< 2| which is equivalent 
a−1 b
to |a+d| < a+d < 2 (since det   = 0 ⇒ ad − bc + 1 = a + d.
c d−1
However, the first Jury condition is never satisfied, then at least one
eigenvalue has magnitude greater than or equal to one. To determine
stability, We need to examine more details.

36
Actually when a + d > 0, |a + d| = a + d, at least one eigenvalue has
magnitude one with the proof:

λ2 − (tr J)λ + det J = 0 ⇒λ2 − (a + d)λ + ad − bc = 0 ⇒

p
a + d ± (a + d)2 − 4(ad − bc)
λ=
p 2
a + d ± (a + d)2 − 4(a + d − 1)
=
p 2
a + d ± (a + d − 2)2
=
2
a + d ± (a + d − 2)
=
2

which leads to eigenvalues λ1 = a + d − 1, λ2 = 1. In general, we may

need to look at higher order terms of the right hand side of the system
when one or a few eigenvalues have magnitude one and all other
eigenvalues have magnitudes less than one. Fortunately the system
here is linear (no high order terms), then the stability is determined
by the eigenvalue whose magnitude is not one. Hence, all equilibria
x∗ , 1−a ∗


b x are stable if |a + d − 1| < 1 ⇔ 0 < a + d < 2.

When a + d < 0, |a + d| > a + d, the eigenvalues λ1 = a + d − 1 <

−1, λ2 = 1. Since |λ1 | = |a + d − 1| > 1, then the stability theorem
implies that all equilibria x∗ , 1−a ∗


b x are unstable.

Activity 3
Consider host-parasitoid interactions (see Figure 3.1). Parasitoids lay
eggs on hosts, and thus hosts are separated into two groups. Par-
asitized hosts give rise to the next generation of parasitoids, while

37
non-parasitized hosts give rise to the next generation of hosts. Write
up a simple model for host-parasitoid interactions and perform math-
ematical analysis.

Figure 3.1: host-parasite model

3.3 Summary
In this unit we have discussed discrete-time models and solutions of
solving such models. We have also studied analysis of such models
and stability of equilibrium points. Systems of difference equation
models have also been studied.

3.4 End of unit test

1 . Consider the following model

xn+1 = xn + k (N − xn ) xn

38
 where N is the maximum population that can be sustained by
the environment. N often is referred to as the carrying capacity
of the population.

a . Solve the model

b . Analyze the model

3.5 Further Reading

1. Wang H. (2012) Mathematical Modeling I - preliminary. ebook.

2. Giordano F. A first course in Mathematical modeling. 5th edition.

Cengage Learning.

39
unit 4

CONTINUOUS-TIME
MODELS

Introduction
In this unit we are going to study continuous-time models. We first
study differential equations and how they are used in mathematical
modeling. Analysis of continuous models is also done which is fol-
lowed by studying systems of differential model.

dN (t)
= (b − d)N (t)
dt
These are the models that treat time as a continuous variable.

40
 Learning outcomes
By the end of this unit you should be able to;

a . Solve differential equations by separation of variables

b . Solve continuous-time models

c . Analyze continuous time models

d . Perform qualitative analysis of solutions to continuous-time

models

c . Solve systems of differential equation models

Key words
The following terms are useful in this unit;
Differential equations, continuous, qualitative analysis, eigenvalue,
stability, Jacobian matrix, determinant.

41
4.1 Differential equations in Mathematical
modeling

Theorem 3. Separation of variables

Consider the nonautonomous differential equation

dx
= f (t, x)
dt

which is separable if
f (t, x) = p(t)q(x)

Note that p(t) or q(x) may be a constant function. For a separable

equation, we can express it as
Z Z
dx dx
= p(t)q(x) ⇒ = p(t)dt
dt q(x)

from which we can obtain the general solution to the differential equa-
tion.

4.1.1 Solving differential equations by separation of

variables

Examples:

42
 dx
(a) Solve the differential equation dt = t2 + (xt)2 .
Solution: This equation is separable because we can rewrite it as

dx
= t2 1 + x2


dt

then Z Z
dx
= t2 dt
1 + x2
tan−1 (x) = t3 /3 + C
x = tan t3 /3 + C

where C is a constant determined by the initial condition.

dx t
(b) Solve dt = et+x .

Solution: This equation is separable because it is equivalent to

dx t 1
= t x
dt e e

Thus, Z Z
x t
e dx = dt
et
te−t dt = integration by parts = − te−t − e−t dt +C =


⇒ ex =
R R

− (te−t + e−t )+C = −(t+1)e−t +C. Hence, the general solution is

x(t) = ln [−(t + 1)e−t + C], where the constant C depends on the
initial condition. For instance, given the initial condition x(0) =
0, then
x(0) = ln[−1 + C] = 0,

43
 which leads to C = 2. Then the solution of the IVP

 dx = t ,

dt et+x
(0) = 0,


is x(t) = ln [−(t + 1)e−t + 2].

x

 dx = e

(c) Solve the IVP dt t2 x
(1) = 0


Solution: We apply separation of variables to obtain

Z Z
xe−x dx = t−2 dt,
 Z 
− xe−x − e−x dx = −t−1 + C,

− xe−x + e−x = −t−1 + C,

−(x + 1)e−x = −t−1 + C.

Then use the initial condition x(1) = 0 to obtain −(0 + 1)e−0 =

−1−1 + C, which gives C = 0. Hence, the solution to IVP is

1 1
−(x + 1)e−x = − − ,
t 2

Activity 4

1 . Solve the following differential equations by separation of vari-

44
 ables

dy
a x dx =y

b 4 − y 2 y ′ = x2

c x2 + y 2 y ′ = y 2

4.2 Continuous-time models: Motivation and

derivation
We derive a continuous model from a discrete model for population
prediction. Let N (t) be the population size at time t. In a small time
period ∆t, a percentage b of the population is born, and a percent-
age d of the population dies. Thus the change of the population size
during the time period ∆t is

N (t + ∆t) − N (t) = bN (t)∆t − dN (t)∆t

Note that b and d have the unit per time, and they are called growth
and death rates. Divide both sides by ∆t to obtain

N (t + ∆t) − N (t)
= (b − d)N (t)
∆t

45
Let ∆t → 0, we have

dN (t)
= (b − d)N (t)
dt

using the definition of derivative. Let α = b − d, then the model

becomes
dN
= αN
dt

Give the initial condition N (t0 ) = N0 , the continuous-time model

dN
dt = αN is defined for t ≥ t0 . Solution: We apply separation of
variables to obtain Z Z
dN
= αdt
N
which leads to
ln N = αt + C

for some constant C. We then use the initial condition N (t0 ) = N0 to

obtain ln N0 = αt0 + C, thus C = ln N0 − αt0 . Substitute it back into
the solution, then ln N = αt + ln N0 − αt0 , which leads to ln (N/N0 ) =
α (t − t0 ), and thus
N (t) = N0 eα(t−t0 )

which is the solution of the continuous model.

4.2.1 Examples of Continuous models

Suppose the 2000 census for the population of China was 1262600000
and in 1980 it was 981235000 . Substitute these values into the so-

46
lution by letting t0 = 1980 and N0 = 981235000 :

1262600000 = 981235000eα(2000−1980) ,

from which we solve for α : α = 0.0126. Hence the model be-

comes
N (t) = 981235000e0.0126(t−1980) ,

which can be used to predict future population. For example, in 2010

the Chinese population size should be N (2010) = 981235000e0.0126(2010−1980) =
1432000000, overestimate the realistic number 1338300000. How
about the year 2100?N (2100) = 981235000e0.0126(2100−1980) = 4450700000,
obviously unsustainable in China. Clearly the model is oversimplified.
Improved Continuous Model
We improve the model with limited growth. The constant α should
be a decreasing function of the population size N and becomes zero
when N reaches the sustainable maximum populations size M . The
simplest to incorporate the population ceiling is

α = r(M − N )

which is
- positive when N < M ;
- zero when N = M ;
- negative when N > M .

47
Therefore, the improved model is provided by

dN
= r(M − N )N
dt

called logistic population growth model.

4.2.2 Solution to the logistic population growth model

Solution: We apply separation of variables to obtain

Z Z
dN
= rdt
(M − N )N
Z   Z
1 1 1
+ dN = rdt
M M −N N
1
(− ln |M − N | + ln N ) = rt + C
M
ln N − ln |M − N | = M (rt + C)

Use the initial condition N (t0 ) = N0 < M ⇒ ln N0 − ln (M − N0 ) =

M (rt0 + C), from which we solve for C : C = 1
M ln MN−N
0
0
− rt0 .
Hence,
 
1 N0
ln N − ln |M − N | = M rt + ln − rt0
M M − N0

Consider the case N < M which is usually valid according to the

48
definition of M , then

N (M − N0 )
ln = rM (t − t0 )
N0 (M − N )
N (M − N0 ) = N0 (M − N )erM (t−t0 )
N (M − N0 ) + N0 N erM (t−t0 ) = N0 M erM (t−t0 )
N0 M erM (t−t0 )
N (t) =
M − N0 + N0 erM (t−t0 )
M N0
N (t) =
N0 + (M − N0 ) e−rM (t−t0 )

This is the solution of the improved model with limited [Link]

figure 4.1 We can easily observe that

N (t) → M as t → ∞.

4.2.3 Stability Criterion of continuous models

Theorem 4. (Stability Criterion) Let x∗ be an equilibrium of an au-

dx
tonomous equation dt = f (x). Then
- x∗ is stable when f ′ (x∗ ) < 0;
- x∗ is unstable when f ′ (x∗ ) > 0;
- there is no conclusion about the stability of x∗ when f ′ (x∗ ) = 0.

f ′ (x∗ ) is called the eigenvalue of x∗ .

49
 dN
Figure 4.1: The logistic curve: the solution of the model dt = r(M −
N )N with the initial condition N (t0 ) = N0 .

Example 1. Determine equilibrium values and their stability for the

differential equation

dx
= (x + 1)(x − 2).
dt

dN
Example 2. Revisit the logistic population growth model dt =

r(M − N )N to determine equilibrium values and their stability. [Link]-

lution: Equilibrium values are x∗ = −1 and x∗ = 2.
f (x) = (x + 1)(x − 2) = x2 − x − 2, then f ′ (x) = 2x − 1, thus

50
f ′ (x∗ ) = 2x∗ − 1. f ′ (−1) = −3 < 0

Hence the equilibrium solution x = −1 is stable. f ′ (2) = 3 > 0

Hence the equilibrium solution x = 2 is unstable. [Link]: Equi-

librium values are N ∗ = M and N ∗ = 0.
f (N ) = r(M − N )N = rM N − rN 2

Then f ′ (N ) = rM − 2rN , thus f ′ (N ∗ ) = rM − 2rN ∗ .

f ′ (M ) = −rM < 0, hence the equilibrium solution N = M is stable.

f ′ (0) = rM > 0, hence the equilibrium solution N = 0 is unsta-

ble.

4.3 Differential Equation models

Consider a general first-order differential equation:

dx
= f (t, x)
dt

dx
If f (t, x) is independent of t, i.e. dt = f (x), then the differential equa-
tion is called autonomous.
dx
The values of x such that dt = 0 are called equilibrium values or

steady states or fixed points.

51
 A phase line is a plot on the x axis that shows all fixed points to-

dx
gether with the intervals where we can determine the signs of dt and
d2 x
dt2 , from which we know the monotonicity and concavity of solution
curves. Examples
dx
1: For the autonomous differential equation dt = x(x − 1), determine
equilibrium values and perform the phase line analysis.
Solution Equilibrium values are x∗ = 0 and x∗ = 1, i.e. the equation
has two constant solutions x = 0 for all t and x = 1 for all t.
dx
It is easy to determine that dt = x(x − 1) is positive if x < 0 or x > 1;

negative if 0 < x < 1. We can compute

d2 x d dx d dx


dt2 = dt dt = dt (x(x − 1)) = (2x − 1) dt = (2x − 1)x(x − 1)

which is positive if x > 1 or 0 < x < 1/2; negative if 1/2 < x < 1 or

x < 0.

dx
Figure 4.2: The phase line analysis of dt = x(x − 1)

4.4 Systems of differential equation models

Some motivation examples

52
(a) Example 1 (A competition model). Consider a pond that sup-
ports two fish types: trout and bass. They compete for food. Let
x(t) be the population of trout at time t and y(t) be the population
of bass at time t. The model is provided by

dx
dt = ax − bxy, ( growth-competition)
dy
dt = my − nxy, (growth-competition)

where a, b, m, n > 0, x, y ≥ 0.
The parameters a and m represent the per capita growth rates of

trout and bass, respectively.

The parameters b and n represent the competition strengths

between trout and bass.

(b) Example 2 (The Lotka-Volterra competition model).

Let y1 (t) be the population of species 1 and y2 (t) be the popula-

tion of species 2. The model incorporating both between-species

and within-species competitions can be provided by

dy1
= r1 y1 (1 − y1 − a1 y2 )
dt
dy2
= r2 y2 (1 − y2 − a2 y1 )
dt

where r1 , r2 , a1 , a2 > 0, y1 , y2 ≥ 0. The right hand side of either

equation is formulated as growth within species competition -
between species competition. This model is usually called the

53
 Lotka-Volterra competition model.

(c) Example 3 A basic epidemic model. Epidemic models are con-

structed for understanding the spread of an infectious disease in
a host population.

We define the variables and the parameters as follows:

S - the number of susceptible individuals.

I - the number of infected-infectious individuals.

R - the number of recovered individuals.

β > 0 - the transmission coefficient.

α > 0 - the recovery rate.

γ > 0 - the rate for the loss of immunity Infected individuals can

Figure 4.3: A flowchart of SIR epidemic model.

54
infect susceptible individuals and can be recovered. It is possible
for recovered individuals to lose their immunity after some time.
The classical SIR (susceptible-infected-recovered) model describ-
ing the above process is provided by

S ′ = −βSI + γR,
I ′ = βSI − αI,
R′ = αI − γR,

where βSI is the rate of new incidences via direct contact be-
tween S and I, using the mass-action form. (S + I + R)′ = 0 ⇒
S + I + R = constant, say N , then

S ′ = −βSI + γ(N − S − I),

I ′ = βSI − αI,
R = N − S − I,

If we assume permanent immunity, then γ = 0, then

S ′ = −βSI,
I ′ = βSI − αI,
R′ = αI,

or equivalently
S ′ = −βSI
I ′ = βSI − αI
R=N −S−I

55
4.4.1 Analytical methods for systems of differential
equation models.

Consider the following system of differential equation models.

x′1 = f1 (x1 , x2 ) ,
x′2 = f2 (x1 , x2 ) .

     
x1 f1 (x) f1 (x1 , x2 )
Let x =   , f (x) =   =  , we have the
x2 f2 (x) f2 (x1 , x2 )
vector form:
x′ = f (x).

Actually, for any number of autonomous differential equations, we

can always write the system as this vector form.

Definition 4.1. Nullclines:

The x1 -nullcline is the set of points (x1 , x2 ) such that x′1 = f1 (x1 , x2 ) = 0.
The x2 -nullcline is the set of points (x1 , x2 ) such that x′2 = f2 (x1 , x2 ) =

0.

56
4.4.2 Equilibrium points (or fixed points or steady
states)

All equilibrium points (x∗1 , x∗2 ) satisfy

f1 (x∗1 , x∗2 ) = 0
f2 (x∗1 , x∗2 ) = 0

the intersection(s) of x1 - and x2 -nullclines.

Stability of equilibrium points.  
∂f1 ∂f1
∂x1 (x∗1 , x∗2 ) ∂x2 (x∗1 , x∗2 )
Jacobian matrix at (x∗1 , x∗2 ) : J (x∗1 , x∗2 ) =  
∂f2 ∂f2
∂x1 (x∗1 , x∗2 ) ∂x2 (x∗1 , x∗2 )
Characteristic equation: det (J (x∗1 , x∗2 ) − λI) = 0 whose roots λ1 , λ2

are eigenvalues.

Theorem 5. Stability criterion:

The equilibrium point (x∗1 , x∗2 ) is stable if all eigenvalues have negative
real parts. The equilibrium pint (x∗1 , x∗2 ) is unstable if at least one of the
eigenvalues has positive real part.

For a system of two differential equations, the condition ”all eigenval-

ues have negative real parts” is equivalent to det J (x∗1 , x∗2 ) > 0 and
tr J (x∗1 , x∗2 ) < 0.
Example

57
(a) Find all equilibria and determine their stability for the basic com-
petition model below.
A pond supports two fish types: trout and bass. They compete

for food. Let x(t) be the population of trout at time t and y(t) be
the population of bass at time t. The model is provided by

dx
dt = ax − bxy, ( growth-competition)
dy
dt = my − nxy, (growth-competition)

where a, b, m, n > 0, x, y ≥ 0.
The parameters a and m represent the per capita growth rates of

trout and bass, respectively.

The parameters b and n represent the competition strengths

between trout and bass. Equilibrium points (x∗ , y ∗ ) satisfy

ax∗ − bx∗ y ∗ = 0
my ∗ − nx∗ y ∗ = 0

There are two possibilities: (0, 0), (m/n, a/b).

We can determine their stability by checking the Jacobian matrix

 
a − by −bx
J(x, y) =  .
−ny m − nx

58
  
∗ ∗
a − by −bx
At each equilibrium point (x∗ , y ∗ ) , J (x∗ , y ∗ ) =  .
∗ ∗
−ny m − nx
For the equilibrium (0, 0),
 
a 0
J(0, 0) =  
0 m

det(J(0, 0) − λI) = 0
 
a−λ 0
det  
0 m−λ

(a − λ)(m − λ) = 0

Eigenvalues: λ1 = a > 0, λ2 = m > 0. Thus, the equilibrium (0, 0)

is an unstable node. For the equilibrium (m/n, a/b),
 
0 −bm/n
J(m/n, a/b) =  
−na/b 0

det(J(m/n, a/b) − λI) = 0

 
−λ −bm/n
det  =0
−na/b −λ

λ2 − am = 0
√
λ = ± am

√
Since am > 0, we obtain two real eigenvalues: λ1 = am >
√
0, λ2 = − am < 0.

59
 Thus the equilibrium (m/n, a/b) is a unstable.

(b) Consider the Lotka-Voltera competition model in example 2

dy1
= r1 y1 (1 − y1 − a1 y2 )
dt
dy2
= r2 y2 (1 − y2 − a2 y1 )
dt

where r1 , r2 , a1 , a2 > 0, y1 , y2 ≥ 0. Solve the model and perform

stability analysis of all equilibrium values. Equilibrium points
(ȳ1 , ȳ2 ) satisfy
r1 ȳ1 (1 − ȳ1 − a1 ȳ2 ) = 0
r2 ȳ2 (1 − ȳ2 − a2 ȳ1 ) = 0

There are four possibilities: (0, 0), (1, 0), (0, 1), (y1∗ , y2∗ ) where y1∗ , y2∗
satisfy
1 − y1∗ − a1 y2∗ = 0
1 − y2∗ − a2 y1∗ = 0
1 − a1 1 − a2
⇒ y1∗ = , y2∗ = .
1 − a1 a2 1 − a1 a2

60
Jacobian matrix
 
∂f1 ∂f1
∂y1 ∂y2
J (ȳ1 , ȳ2 ) =  
∂f2 ∂f2
∂y1 ∂y2 (ȳ1 ,ȳ2 )
 
r1 (1 − 2y1 − a1 y2 ) −r1 a1 y1
= 
−r2 a2 y2 r2 (1 − 2y2 − a2 y1 )
 
r1 (1 − 2ȳ1 − a1 ȳ2 ) −r1 a1 ȳ1
= .
−r2 a2 ȳ2 r2 (1 − 2ȳ2 − a2 ȳ1 )

For the equilibrium

  point (0, 0), the Jacobian matrix J(0, 0) =
r 0
 1  has eigenvalues λ1 = r1 > 0, λ2 = r2 > 0, thus the
0 r2
equilibrium (0, 0) is an unstable node.

For the equilibrium point (1, 0), the Jacobian matrix J(1, 0) =

−r1 −r1 a1
  has eigenvalues λ1 = −r1 < 0, λ2 = r2 (1 − a2 ).
0 r2 (1 − a2 )
Thus the equilibrium (1, 0) is asymptotically stable (a stable node)
if a2 > 1, and it is unstable (a saddle) if a2 < 1.

For the equilibrium point (0, 1), the Jacobian matrix J(0, 1) =

r (1 − a1 ) 0
 1  has eigenvalues λ1 = r1 (1 − a1 ) , λ2 = −r2 <
−r2 a2 −r2
0. Thus the equilibrium (0, 1) is asymptotically stable (a stable
node) if a1 > 1, and it is unstable (a saddle) if a1 < 1. For the co-

61
existence/internal equilibrium point (y1∗ , y2∗ ), the Jacobian matrix
 
r1 (1 − 2y1∗ − a1 y2∗ ) −r1 a1 y1∗
J (y1∗ , y2∗ ) =  
−r2 a2 y2∗ r2 (1 − 2y2∗ − a2 y1∗ )
   
2−2a1 a1 −a1 a2 a1 −a21
r1 1 − 1−a1 a2 − 1−a1 a2 −r1 1−a 1 a2
=  
a2 −a22 2−2a2 a2 −a1 a2
−r2 1−a 1 a2
r2 1 − 1−a1 a2 − 1−a1 a2
 
r1 (1−a1 a2 −2+2a1 −a1 +a1 a2 ) r1 (a21 −a1 )
1−a1 a2 1−a1 a2
= r2 (a22 −a2 )

r2 (1−a1 a2 −2+2a2 −a2 +a1 a2 )
1−a1 a2 1−a1 a2
 
1 r (a − 1) r1 a1 (a1 − 1)
=  1 1 .
1 − a1 a2 r2 a2 (a2 − 1) r2 (a2 − 1)

Now we check the conditions tr J (y1∗ , y2∗ ) < 0 and det J (y1∗ , y2∗ ) >
0. All eigenvalues have negative real parts if and only if

r1 (a1 − 1) + r2 (a2 − 1)
tr J = <0
1 − a1 a2

and

r1 r2 (a1 − 1) (a2 − 1) a1 a2 r1 r2 (a1 − 1) (a2 − 1)

det J = −
(1 − a1 a2 )2 (1 − a1 a2 )2

r1 r2 (a1 − 1) (a2 − 1)
= > 0.
1 − a1 a2
r1 r2 (a1 −1)(a2 −1)
We consider four cases: (i) a1 > 1 and a2 > 1 ⇒ 1−a1 a2 <
0, the second inequality is violated, thus the coexistence equilib-
rium (y1∗ , y2∗ ) is unstable. (ii) a1 < 1 and a2 < 1 ⇒ both in-

62
 equalities are valid, thus the coexistence equilibrium (y1∗ , y2∗ ) is
asymptotically stable. (iii) a1 > 1 and a2 < 1 ⇒ it is difficult
to determine whether any of the two inequalities are valid. (iv)
a1 < 1 and a2 > 1 ⇒ it is difficult to determine whether any of
the two inequalities are valid.

Since populations y1 , y2 are nonnegative, we only consider the

first quadrant of the phase plane ( y1 − y2 plane). We need to
examine whether y1∗ , y2∗ are nonnegative. For case (iii) and case
(1−a1 )(1−a2 )
(iv), y1∗ y2∗ = (1−a1 a2 )
2 < 0, thus either y1∗ or y2∗ is negative. This
means that the coexistence equilibrium (y1∗ , y2∗ ) is unfeasible.

For case (i) and case (ii), the coexistence equilibrium (y1∗ , y2∗ ) is
obviously feasible in the first quadrant.

As a conclusion, there are four cases under different values of a1

and a2 . y1 -nullcline: f1 (y1 , y2 ) = 0 ⇔ y1 = 0 or 1 − y1 − a1 y2 =
0. y2 -nullcline: f2 (y1 , y2 ) = 0 ⇔ y2 = 0 or 1 − y2 − a2 y1 = 0.
Intersections of different nullclines provide equilibrium points.

4.5 Summary
In this unit we have looked at continuous models and methods of
solving them and their analysis. We have also studied systems of dif-

63
ferential equation models and analysis of stability of equilibria using
method of Jacobian matrix.

End of unit test

1 . Think of a real life scenario which can be modeled by a

continuous-time model

a . Develop the model

b . Solve the model

c . Analyse the model

4.6 Further Reading

1. Wang H. (2012) Mathematical Modeling I - preliminary. ebook.

2. Giordano F. A first course in Mathematical modeling. 5th edition.

Cengage Learning.

64