Mathematical Modelling, I Notes - 2

Applications of first-order ODEs

An ordinary differential equation (ODE) is an equation that involves the derivatives

of a dependent variable y with respect to a single independent variable x, e.g. 𝑦 / =
𝑑𝑦 𝑑2 𝑦
, 𝑦 // = . The equation may also involve y itself and some given functions of x.
𝑑𝑥 𝑑𝑥 2
Sometimes, there may also be more than one dependent variable. Consider the
following ordinary differential equations.

a) 𝑦 / = 4𝑥7
b) 𝑦 // − 𝑦 / + 12𝑦 = 5 cos 𝑥
2
c) 𝑦 // + 2(𝑦 / ) = 75𝑥
d) 𝑦𝑦// + 1 = 0
2
e) (𝑦/ ) + 𝑦 = 7𝑥
f) 4

Classification
In order to know whether or not we are in a case where an explicit solution is
possible, we need to be able to classify the ODE.

Order:
The order of an ODE is the largest number of times that the dependent variable is
differentiated in the ODE. Orders of (a)–(e) above are respectively: 1, 2, 2, 3, 1,

Linear:
An ODE is linear if it contains no products or powers (other than one) of the
dependent variable or its derivatives. A linear ODE cannot contain terms like 𝑦 2 , 𝑦𝑦 / ,
√𝑦, 𝑐𝑜𝑠𝑦 etc. Powers of the independent variable are allowed. So the examples (a)
and (b) above are linear and the others are nonlinear.

Why differential equations?

Typically, mathematical models can be classified in a number of ways, e.g., empirical
vs. mechanistic. Empirical models are often data-driven and use statistical methods
to infer relations between different variables in a system; Mechanistic models use
prior knowledge and understanding of the mechanism that causes change in a system.

1
In many applications, we are given the value of a quantity at the present time (for
example, the temperature of coffee in a cup, the number of people infected with a
virus, the concentration of carbon dioxide in the atmosphere) and we wish to predict
its value in the future. To do this, we must know how quickly the quantity is changing.
Mathematically, the rate of change of this quantity is its derivative.
If we can write down an equation relating this derivative to some knowledge of the
mechanism that causes the change, we can calculate how the temperature changes,
how the number of infected people changes, how the concentration of carbon dioxide
changes. This gives rise to a differential equation—which is a mechanistic model of
the quantity of interest.

Solving a problem using mathematical models often involves three steps:

(1) building the model,

(2) solving the equations and
(3) interpreting the solutions.

In this chapter, we apply the theory and solution methods of first-order ODEs to a
number of application areas. For each modelling problem, the goal is to formulate
the ODE that must be solved, along with its initial condition, and use an appropriate
method to find the particular solution.

The ‘maths’ in this chapter is no more difficult than what we have already seen in
the study of the course ODE. However, before getting to the maths, we have to
learn how to formulate the problem from a description (in words) of a certain
situation. So, while the maths might be no more difficult, the tricky part here is
coming up with the correct ODE in the first place.

Throughout this chapter, we consider the situation carefully and translate the word
describing this situation into an ODE—this is a key part of the process of
mathematical modelling.

Radioactivity and Carbon dating

(i) Decay and half-life
Atoms of a given radioactive isotope have a fixed percentage chance of decaying in
any given time period. Provided a sufficiently large number of atoms are considered
as a mass, the number of decays that occur in a fixed time period is simply
proportional to the number of atoms and the length of time.

2
Let N(t) be the number of atoms at time ‘t’ and consider this number a (short) time
later N(t + 𝛿t). Assume we know that there are N0 atoms at a particular time t0. In
the time interval between t and t + 𝛿t, the number of atoms decreases by an mount
proportional to

(i) the number itself at that time (i.e., 𝑁(𝑡)and

(ii) the length of the time interval (i.e.,𝛿𝑡). Mathematically, this translates to:
𝑁(𝑡 + 𝛿𝑡) = 𝑁(𝑡) − 𝑘𝑁(𝑡)𝛿𝑡

where k is a positive constant and we take the negative sign since the number is
decreasing. Rearranging and making –kN(t) or just -kN the subject, we have:

𝑁(𝑡+𝛿𝑡)−𝑁(𝑡)
= −𝑘𝑁(t)
𝛿𝑡

and taking the limit of small 𝛿𝑡 from the method of first principles with 𝑓 (𝑥) = 𝑁(𝑡),
we have

= −𝑘𝑁, with initial condition: (3.1)

𝑑𝑁
𝛿𝑡

Solving 3.1 by the method of separation of variables, we obtain

𝑑𝑁
𝑁
= −𝑘𝛿𝑡

Integrating both sides,

𝑑𝑁
∫ = ∫ −𝑘𝛿𝑡
𝑁

𝑙𝑛𝑁 = −𝑘𝑡 + 𝑐

exponentiating,
𝑒 𝑙𝑛𝑁 = 𝑒 −𝑘𝑡+𝑐

or 𝑁 = 𝑒 −𝑘𝑡 ∙ 𝑒 𝑐

or 𝑁 = 𝐴𝑒 −𝑘𝑡

using the initial condition, N(t0) = N0 to find A,

𝑁0 = 𝐴𝑒 −𝑘𝑡0 ,

3
making A the subject, we have

𝐴 = 𝑁0 𝑒 𝑘𝑡0 ,

replacing A with 𝑁0 𝑒 𝑘𝑡0 , we have

𝑁 = 𝑁0 𝑒 𝑘𝑡0 ∙ 𝑒 −𝑘𝑡

or 𝑵 = 𝑵𝟎 𝒆−𝒌(𝒕−𝒕𝟎 ) (3.2)

which results in the particular solution (3.2) for N(t):

Note that the symbol N is the dependent variable is our choice. We could also choose
to denote the number of atoms by, for example, x and the problem we must solve is
the same. For this problem (and many others), the independent variable is time. For
modelling problems, we try to keep a sensible naming convention, so if we are
discussing rates of change with respect to time, we almost always name the
independent variable as t. However, there is nothing wrong mathematically with
choosing another symbol for time.

4
Plotting a graph of N against t for equation (3.2) we obtain

The blue line in this figure plots N(t) for a general initial condition N(t0) = N0; it is
an example of exponential decay where k is the rate of decay.

Half-life is the time required for a quantity to reduce to half of its initial value. In
radioactivity, the half-life, which we will denote 𝑇ℎ , is the time after which half of
the isotope has decayed. Initially, at time t0, there are N0 atoms. So, we are
interested in the time 𝑇ℎ in the future at which there are 2 𝑁0 atoms.
1

Mathematically, this means N(𝑡0 + 𝑇ℎ ) = 2 𝑁0 . Using the solution (3.2) at time 𝑡0 + 𝑇ℎ :

1

i.e N(𝑡0 + 𝑇ℎ ) =
1
2
𝑁0 atoms at half-life,

from (3.2), or 𝑁 = 𝑁0 𝑒 −𝑘(𝑡−𝑡0) , replacing t with (𝑡 + 𝑇ℎ ), we have

5
 1
𝑁 = 𝑁0 𝑒 −𝑘(𝑡0+𝑇ℎ −𝑡0) = 2 𝑁0
or 𝑁0 𝑒 −𝑘𝑇ℎ = 2 𝑁0
1
, from which,

1
𝑒 −𝑘𝑇ℎ = 2
Taking reciprocals,
𝑒 𝑘𝑇ℎ = 2

Introducing Napierian logarithms and rearranging, to find an expression for the

half-life in terms of the rate of decay k, we obtain

𝑘𝑇ℎ
𝑙𝑛𝑒 = 𝑙𝑛2

𝑘𝑇ℎ = 𝑙𝑛2

or 𝑇ℎ = 𝑘 𝑙𝑛2
1

(ii) Carbon dating

There are three isotopes of Carbon: 12C, 13C and 14C. Almost all Carbon is made up of
the first two (12C and 13C) because 14C is radioactive: it decays with a half-life of
5730 years to form 14N, a Nitrogen isotope. Although it decays quite quickly, it is
constantly being produced in the upper atmosphere by the action of cosmic rays. The
equilibrium level of 14C is about 1 part per trillion (1012).

When an organism dies, it ceases to absorb Carbon from the environment, so the
amount of 14C it contains will decrease as this decays radioactively. The time since
the death of the organism can be estimated by measuring how much 14C remains.

14
For C, the half-life is 𝑇ℎ 5730 years, from the relation 𝑇ℎ = 𝑘 𝑙𝑛2
1

so
1
𝑘= 𝑙𝑛2
𝑇ℎ

0.000121/year

Example
14
A fossilized bone is found to contain 0.1% of its original C. Find the age of the
fossil.

6
Solution:
Let ‘𝑡’ be today’s date, and suppose the fossil is T years old, so it was fossilized at
time 𝑡 = 𝑡0 − 𝑇. At this time it had 𝑥0 14C, and now it has 0.001𝑥0 . So the number of
atoms ‘x’ at time ‘𝑡’ is

From (3.2), or 𝑁 = 𝑁0 𝑒 −𝑘(𝑡−𝑡0) , where 𝑁0 = 𝑥0 (or replacing 𝑁0 with 𝑥0 and N(t) with
x(t))

𝑥 (𝑡) = 𝑥0 𝑒 −𝑘(𝑡−𝑡0 ) = 0.001𝑥0

𝑥0 𝑒 −𝑘(𝑡−𝑡0) 0.001𝑥0 , replacing (𝑡 − 𝑡0 with T, from 𝑡 = 𝑡0 − 𝑇

𝑥0 𝑒 −𝑘𝑇 = 0.001𝑥0

𝑒 −𝑘𝑇 = 0.001

which on solving for T, we obtain

−1
𝑇= 𝑙𝑛0.001
𝑘

−1 1
= 𝑙𝑛 ( )
𝑘 1000

1 1 −1
= 𝑙𝑛 ( ) , from the third law of logarithms
𝑘 1000

= 𝑙𝑛1000, but k = 0.000121/year

1
𝑘

= 57,088.87

≈ 𝟓𝟕, 𝟏𝟎𝟎

Therefore, the fossil is approximately 57100 years old.

7
(iii) Newton’s law of cooling
(Sir Isaac Newton FRS: 1643–1727, English physicist, mathematician, astronomer,
natural philosopher, alchemist and theologian).

Newton’s law of cooling states that if an object is hotter than the ambient
temperature, then the rate of change of the object’s temperature is proportional to
the temperature difference between the object and its surrounding.
Mathematically, we write:

= −𝑘 (∅ − 𝐴), with (3.3)

𝑑∅
𝑑𝑡
∅(𝑡0 ) = ∅0

where ∅(𝑡) is the object’s temperature, A is the ambient temperature (a constant),

′𝑡’ is time and k is a positive constant. This is a first-order linear (separable) ODE:
(Note that, in more complicated examples, the ambient temperature may be time-dependent A = A(t).)
Equation (3.3) is rearranged to become

𝑑∅
𝑑𝑡
+ 𝑘∅ = 𝑘𝐴
Which on solving using the integrating factor 𝑅(𝑡) = 𝑒 ∫ 𝑘(𝑑𝑡) 𝑒 𝑘𝑡 , when multiplied by
(3.3) yields

𝑑∅
( + 𝑘∅) 𝑒 𝑘𝑡 = 𝑘𝐴𝑒 𝑘𝑡
𝑑𝑡

𝑑∅
𝑒 𝑘𝑡 + 𝑘∅𝑒 𝑘𝑡 𝑘𝐴𝑒 𝑘𝑡
𝑑𝑡

where RHS is the derivative of ∅𝑒 𝑘𝑡 by the product rule, with both ∅ and 𝑒 𝑘𝑡 ,
functions of ‘t’.

thus,
𝑑
(∅𝑒 𝑘𝑡 ) = 𝑘𝐴𝑒 𝑘𝑡
𝑑𝑡
rearranging,
𝑑(∅𝑒 𝑘𝑡 ) = 𝑘𝐴𝑒 𝑘𝑡 𝑑𝑡
integrating

∫ 𝑑(∅𝑒 𝑘𝑡 ) = ∫ 𝑘𝐴𝑒 𝑘𝑡 𝑑𝑡,

yielding,
𝑘𝐴𝑒 𝑘𝑡
∅𝑒 𝑘𝑡 = 𝑘
+𝐶

or ∅𝑒 𝑘𝑡 = 𝐴𝑒 𝑘𝑡 + 𝐶 …… which is the general solution.

8
applying the initial conditions, ∅(𝑡0 ) = ∅0 , to determine C, we have

𝐶 = ∅𝑒 𝑘𝑡 − 𝐴𝑒 𝑘𝑡 , but ∅ = ∅0 when 𝑡 = 𝑡0
Hence
𝐶 = ∅0 𝑒 𝑘𝑡0 − 𝐴𝑒 𝑘𝑡0

𝐶 = (∅0 − 𝐴)𝑒 𝑘𝑡0

Replacing C in the general equation with (∅0 − 𝐴)𝑒 𝑘𝑡0 , we obtain,

∅𝑒 𝑘𝑡 = 𝐴𝑒 𝑘𝑡 + (∅0 − 𝐴)𝑒 𝑘𝑡0

dividing all through by 𝑒 𝑘𝑡 ,

∅ = 𝐴 + (∅0 − 𝐴)𝑒 𝑘𝑡0 ∙ 𝑒 −𝑘𝑡

or ∅(𝑡) = 𝐴 + (∅ − 𝐴)𝑒 −𝑘(𝑡−𝑡0)

Exercise:
Verify that this expression is correct.
Note that as 𝑡 → ∞, we have ∅(𝑡) → 𝐴, i.e., the temperature decays exponentially to
the ambient temperature. (the exponent of ∞is zero)

Example 3.2. A horse shoe is heated to 1000 C and then placed in a room to cool. The
temperature of the room is 100 C. After 20 minutes, the temperature of the bar is
500 C. It is safe to handle at 300 C; how long must we wait to handle it?

Solution:
We consider time in minutes with 𝑡0 = 0 the time at which the metal bar is placed in
the room. We use the particular solution to the cooling problem with ∅(𝑡0 ) = 0 = ∅0 =
100, 𝐴 = 100.
∅(𝑡) = 10 + 90𝑒 −𝑘𝑡

where k is the (unknown) decay constant. Next we “translate” the remainder of the
question into equation form:

_ “After 20 minutes, the temperature of the bar is 500 C” translates to ∅(20) = 50

9
_ “It is safe to handle at 300 C; how long must we wait to handle it?” This is the crux
of the problem: find T such that ∅(T) = 30.
We use the first point (∅(𝑡) = 10 + 90𝑒 −𝑘𝑡 to determine k

i.e.∅(20) = 10 + 90𝑒 −𝑘(20) = 50

50−10
𝑒 −20𝑘 = 90

Introducing natural (Napierian) logarithms, we have

ln 𝑒 −20𝑘 = 𝑙𝑛 49

−20𝑘 = −0.810930216

𝑘 = 0.0405

Using this value of k, we can find T, such that ∅(T) = 30

i.e. ∅(𝑇) = 10 + 90𝑒 −𝑘𝑇 = 30

30−10
𝑒 −𝑘𝑇 =
90

Introducing Napierian logarithms again, we have

30−10
𝑙𝑛𝑒 −𝑘𝑇 = 𝑙𝑛 ( 90
)

2
−𝑘𝑇 = 𝑙𝑛 (9)

−0.0405𝑇 = −1.504077397

𝑇 = 37.1377135

𝑇 ≈ 37 minutes

Thus, we must wait approximately 37 minutes to handle the metal bar.

10
(iv) Population growth models
Let p(t) be the population of a country at time t. The rate of change of population in
a country is equal to the rate at which people enter the country (e.g., births,
immigration) minus the rate at which people leave the country (e.g, deaths,
emigration). Hence we can write

𝑑𝑃
𝑑𝑡
= 𝐵(𝑃, 𝑡) − 𝐷(𝑃, 𝑡) + 𝑀(𝑃, 𝑡)

Where,
𝐵(𝑃, 𝑡) represents births,
𝐷(𝑃, 𝑡) represents deaths,
𝑀(𝑃, 𝑡) represents net migration into the country (hence the plus sign).
The processes B, D and M may depend on the population p itself and the time t.

(a) The Malthusian model

(Thomas Robert Malthus FRS: 1766–1834, English clergyman, political economist
and demographer)
.
Malthus (1798) suggested a simple model that has no migration, so M = 0. For the
birth and death rates, he assumed they are proportional to the population p:

𝐵(𝑃, 𝑡) = 𝑏𝑝(𝑡) and 𝐷(𝑃, 𝑡) = 𝑑𝑝(𝑡)

where b and d are positive constants, so

(3.4)
𝑑𝑃
= (𝑏 − 𝑑)𝑃 = 𝛾𝑃
𝑑𝑡

where (𝑏 − 𝑑)𝑃 = 𝛾𝑃 is a constant called the growth rate. Note that depending on
the values of b and d , 𝛾 can be positive, zero or negative. Equation (3.4) is a separable
ODE which we can solve with the initial condition 𝑃(𝑡0 ) = 𝑃0

Given
𝑑𝑃
𝑑𝑡
= (𝑏 − 𝑑)𝑃 = 𝛾𝑃

or = 𝛾𝑃, separating the variables and integrating wrt t

𝑑𝑃
𝑑𝑡

𝑙𝑛𝑃 = 𝛾𝑡 + 𝐶

exponentiating,

11
 𝑒 𝑙𝑛𝑝 = 𝑒 𝛾𝑡+𝐶

𝑃 = 𝑒 𝛾𝑡 ∙ 𝑒 𝐶 = 𝐴𝑒 𝛾𝑡

Solving for A using the initial condition 𝑃(𝑡0 ) = 𝑃0 , we have

𝐴 = 𝑃 (𝑡0 )𝑒 −𝛾(𝑡0)

Substituting for A in 𝑃 = 𝐴𝑒 𝛾𝑡 , we obtain

𝑃 (𝑡) = 𝑃(𝑡0 )𝑒 −𝛾(𝑡0) . 𝑒 𝛾𝑡

𝑃(𝑡) = 𝑃 (𝑡0 )𝑒 𝛾(𝑡−𝑡0)

This solution suggests that it clearly depends on 𝛾; the plot below shows three
possible solutions 𝑃(𝑡) depending on the sign of 𝛾 . Its sign determines whether the
population will increase (more births than deaths) or decrease (more deaths than
births). If the birth and death rate are the same ( 𝛾 = 0), then the population is
steady, i.e., it does not change in time. As a model it is quite limited: it predicts that
the population will increase without bound if 𝛾 > 0 and the population will die out if
𝛾 < 0.

The three solutions for 𝛾 > 0 , 𝛾 = 0and 𝛾 < 0 are depicted in the graph below

12
Example:
In 1770, the population of Great Britain was estimated to be 6.4 million. By 1790,
the population had grown to 8 million. Estimate, and predict the population in the
year 2012.

Solution
Take 𝑡0 = 1770, and 𝑃(𝑡0 ) = 6.4 × 106 . We know 𝑃(1790) = 8 × 106
So
𝑃(1790) = 8 × 106 = 𝑃(𝑡0 )𝑒 (1790−𝑡0)𝛾

6.4 × 106 × 𝑒 20𝛾

Therefore
8×106
𝛾 = 20 𝑙𝑛 (6.4×106 ) = 0.0112 per year
1

Now we can get 𝑃(2012)

𝑃 (2012) = 𝑃 (𝑡0 )𝑒 (2012−𝑡0)𝛾

6.4 × 106 × 𝑒 242×0.0112

96 × 106

13
Simulating a graph of P(t) against time ‘t’ using a computer software like Matlab, we
get;

In fact, the population of Great Britain in 2012 was about 61 × 106 so the estimate
is approximately 50% too large. This is not surprising, since we have not included
effects like immigration, birth and death rates that change with time, changes in
agriculture that allow more food production etc.

14