CHAPTER 8

Solving Differential Equations of the First Order

First order differential equations arise in almost all scientific fields ranging from biology to engineering. We will
be considering those with analytic (exact) solutions and we will explore the physical laws leading to their
construction and source solutions from the integration techniques covered in previous Chapters.

Definitions
Differential equation: an equation containing a variable and its derivatives. In a first order differential equation
the first derivative is the highest derivative present in the equation, for example, (or ).

An example of a second order differential equation is (or ) as the

second derivative is the highest derivative present in the equation. In these examples is the dependent
variable and time the independent variable which is symbolised as

Solution of a differential equation: That that satisfies the differential equation (when you substitute back
into the differential equation, left hand side=right hand side).

Rate of change: the gradient is the rate of change in as a function of time . When this rate is

negative (negative gradient) it indicates a decrease in y as time increases and vice versa for a positive rate
(gradient).

Proportionality: when the rate of change in is proportional to it is indicated as . This means that

there exists a constant such that we can forge the differential equation .

Steady state: This is the state the solution settles into after a long period of time. It is calculated by applying
. The in this Chapter refers to the dependent variable in the model. Before the steady state can

be determined, the dependent variable must be made subject of the solution equation after solving the
differential equation.

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 1

8.1 Three solution methods
8.1.1 Separation of variables
If a differential equation of the first order can be split across the sign into a left hand side containing only the
dependent variable and the right hand side only the independent variable, it is referred to as a separable
differential equation.

To determine if the differential equation is separable, place the derivative term (for example ) on the left

hand side, all other terms on the right hand side and 'multiply' the equation right through by . This will
establish the side (left) and the side (right). You must then:

1. Separate the variables by multiplying and/or dividing terms until all the terms are on the right hand
side and all the terms on the left hand side.
2. Integrate the side with regards to and the side with regards to . Place the constant of integration
on only one of the two sides (usually the side on the right).

For example (or ). Establish sides by

1. Separate variables:

2. Integrate:

This is termed the general solution as we don't know the value of . To calculate a specific solution you must
know one detail about the process being modelled, for instance, suppose that at time also indicated
as . You will then substitute into the general solution to calculate a specific value for :

When you substitute the specific in the general solution, you have the specific solution.

It will be required in many cases (and if possible) that you express the solution in terms of :

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 2

8.1.2 Direct integration
Differential equations only requiring and no further intervention to separate variables can be integrated
directly. The example in 8.1.1 can't be directly integrated as we first had to move to the left hand side to
complete the separation process. An example of a differential equation that can be directly integrated is
(or ) as it only requires for the separation process to be completed

Note that you can easily spot a separable equation by the fact that the dependent variable (in this case ) is only
located in the derivative term and nowhere else. This type can also be seen as the direct reversal of

differentiation:

Paced Exercise (8A)

Calculate the solutions to the following first order differential equations:

1. The streamlines of fluid flow are given by (or ). It is given that .

• You should find that this equation requires separation of variables which should lead to

• After integrating this equation, you should find that the general solution is
• After substituting into the general solution you should find that
• The specific solution is then

2. The acceleration of a particle is given by . It is given that the initial displacement was
and the initial velocity was . Calculate expressions for the velocity and the
displacement of this particle as a function of time. It is known that

• We were given and that . This leads to the differential equation

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 3

 • This equation can be directly integrated as is only contained in the term and nowhere else.

• After substituting into the general expression, the specific is

• To solve the displacement differential equation , we use the specific calculated above

• This can also be integrated directly as is only contained in the term

• You should find that the general solution is

• Using , you should find that the specific solution is

3. Calculate the solutions to . State if you are applying separation of variables or if you can

integrate directly.

• Separation of variables (after you also had to move the to their

respective sides)
• Integrate (the right hand side will require DI integration)

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 4

4. Calculate the solutions to . State if you are applying separation of variables or if you can integrate

directly.

• Separation of variables (after you also had to move the to the side)

• Integrate

• You should find that

5. Calculate the solutions to . State if you are applying separation of variables or if you can

integrate directly.

• Direct integration (after all terms were on the side they belong)

• Integrate (the right hand side will require partial fractions)

• You should find that

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 5

8.1.3 Integrating Factor (IF)
If the differential equation can't be directly integrated and also can't be separated, you must attempt to write it as

Note that the coefficient in front of must always be 1. and may be constants or functions containing

but they may not contain . The only 's are the two you see pointed out: one in the rate expression and

the other raised to the power one. If this format is possible, an integrating factor can be forged that will lead
to the solution. An example of a differential equation that can be expressed in this form and which is neither
separable nor can it be directly integrated is :

Step 1: Direct integration test= fail as is also in another term and not just in .

Step 2: Separability test= fail

. It is impossible to separate the from the term

Step 3: test= pass

If the format is possible, the solution is determined as follows:

• Write as and identify

• Calculate
• Put the solution together as

For those whom appreciate the art of mathematics:

Theorem: The solution to can be calculated by solving where

Proof:

First of all, take the time derivative of the product :

Note that in the first term on the right is the derivative of an integral - an inverse

action resulting in just P remaining (like e and ln, sin and arcsin etc.) The above simplifies to:

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 6

For all differential equations of type

use (*) to replace the first term on the left:

Integrate both sides (the left will be an inverse combination - an integral of a derivative):

Use this method when separation of variables and direct integration is not possible and the equation fits the
format

Example

Solve the equation

Step 1: The coefficient in front of must be 1 so first divide the equation by

Move the term to the left

where you identify

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 7

Step 2:

Step 3: The solution is then put together as

Step 4: general solution

Paced Exercise (8B)

[Link] is the level of alcohol in the blood after a heavy night of drinking (or substance abuse), is the rate of
intake (per unit volume of the person's blood) and is the rate of decay (broken down by the liver), the
differential equation modelling the state of alcohol in the blood stream is

Solve this equation in terms of and . What will happen if the drinking continues at this steady pace?

Note that this equation does pass the separability test as it can be separated as

This equation also passes the test

•
• The
• The general solution is calculated from

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 8

 • From substituting you should find that and the specific solution is

or rather

Determine the steady state: Note that therefore, if drinking carries on at a

steady pace.

2. Calculate the solution to , given that

• This equation won't pass the separability or direct test. You will however find that it passes the test:

• The
• The solution formula is

• the integral will require DI integration:

general solution:

• Substituting into the general solution, the value of

• The specific solution is then

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 9

3. Calculate the general solution to

• First divide the equation by as the coefficient of must be 1

• Identify that

•
• The solution formula is

• You should find that (using DI integration)

• The solution is

4. Calculate the general solution to

• Identify

•
• Solution formula

• You should find that (using DI integration)

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 10

 • The solution is

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 11

Models resulting in first order Differential Equations
8.2 Terminal velocity and steady states
Terminal velocity is the constant velocity an object "falls" at when the forces acting on the object are in balance
(steady state). Before an object reaches terminal velocity, its velocity changes with time and this part will be
modelled by a first order differential equation relating velocity and time.

The modelling process is as follows: A body with mass falling through the air experiences gravitational
acceleration. If at the body is released (initial velocity is zero or ), it will fall downwards with
acceleration and experience a downwards force . The body will also experience an upwards force
referred to as drag (friction) which is either modelled as

• Laminar flow: the drag force is proportional to the velocity of the falling body ( is the
proportionality constant). Laminar flow or streamline flow in pipes (or tubes) occurs when a fluid flows in
parallel layers, with no disruption between the layers. At low velocities, the fluid tends to flow without
lateral mixing, and adjacent layers slide past one another like playing cards. There are no cross-currents
perpendicular to the direction of flow, nor eddies or swirls of fluids. In laminar flow, the motion of the
particles of the fluid is very orderly with all particles moving in straight lines parallel to the pipe walls.
• Turbulent flow: the drag force is proportional to the square of the velocity ( is the
proportionality constant). Turbulent flow is a flow regime characterized by chaotic property changes. This
includes rapid variation of pressure and flow velocity in space and time. In contrast to laminar flow the
fluid no longer travels in layers and mixing across the tube is highly efficient. Flows at Reynolds
numbers larger than 4000 are typically (but not necessarily) turbulent, while those at low Reynolds
numbers below 2300 usually remain laminar. Flow in the range of Reynolds numbers 2300 to 4000 and
known as transition.

The proportionality constant is calculated as where

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 12

a) b) c)

pic a) shows values for C for different type of body shapes, pic b) shows values for for different temperatures
if the medium is air and pic c) shows how you calculate the area A.

In all the previous models the proportionality constants can be calculated through similar processes (the
ingredients will be situation specific). It is however simpler to calculate these constants through "model fitting"
which is possible if both the initial condition of the process is known as well as one other condition later on in the
process.

Newton's second law: For a falling body, the frictional force is in an upwards direction (opposite to ). From
Newton's second law, the resultant force (where is the body's acceleration), the body is
experiencing is given by

When these two opposite forces are in balance, the resultant force is zero and which can be
used to calculate the terminal velocity by solving from

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 13

a) Laminar flow

To calculate the velocity of the object as a function of time before it reaches terminal velocity, one must solve
the force equation:

which turns into the first degree differential equation:

Equation (3) is separable and suited to .

Separation:

Step 1:

Step 2: Integrate

Step 3: General solution

Step 4: Calculate using

To calculate one more condition is required.

Step 5: The updated general solution is

The specific solution (with both known) can be used to calculate the object's velocity for different values
of time t.

Step 6: In making subject of the equation, the limit will determine the body's terminal velocity (steady

state)

as

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 14

Integrating factor:

Step 1: write as

Step 2:

Step 3: which simplifies to

Step 4: initial condition

Step 5: updated general solution

syms v(t) g k
eq=diff(v)==g-k*v;
y=dsolve(eq,v(0)==0)
y =

Example

Given , solve the differential equation modelling a falling object experiencing laminar flow
(the body falls from rest). Calculate the terminal velocity for this body.

The differential equation is

Step 1:

separate variables

Step 2: Integrate

Step 3: General solution

Step 4: Calculate by using

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 15

Step 5: Specific solution

Step 6: The terminal velocity is as the negative exponential

Confirm calculations using Matlab:

syms v(t)
eq=diff(v)==9.8-0.5*v;
v(t)=dsolve(eq,v(0)==0)

v(t) =

fplot(v,[0,10])

title('Laminar flow')
xlabel('time')
ylabel('velocity')

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 16

b) Turbulent flow

For this case, , the resultant force experienced by the falling body is calculated as

This differential equation is not suited to but can be separated:

Step 2:

The left hand side requires hyperbolic substitution

with

substitution:

Example

Given , calculate the velocity of a falling body experiencing turbulent flow (the body falls
from rest). Calculate the terminal velocity (steady state) of this object.

Step 1:

Separate variables

Step 2: Integrate

hyperbolic substitution:

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 17

from

Step 3: The general solution is then

Step 4: Calculate using

Step 5: The specific solution is

Step 6: terminal velocity

from the definition

limits of negative exponentials:

therefore

syms v(t)
dv=diff(v);
eq=dv==9.8-0.5*v^2;
v=dsolve(eq,v(0)==0)

v =

v=vpa(v,3) %simplify solution to show 3 digits

v =
fplot(v,[0,3.5])

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 18

title('Turbulent flow')
xlabel('time')
ylabel('velocity')

Note that we use eval to simplify an expression not containing unknown values but vpa when the expression

does contain unknown values, for example, the in the expression is unknown and

the eval command won't be able to simplify this expression.

Paced Exercise (8C)

1. A 60 kg sky diver jumps out of a plane with no initial velocity and proportionality constant =0.8. Determine his
velocity after 10 seconds if and laminar flow is presumed. What will his terminal velocity be?

• You must first calculate

• The general solution is

• The value of is

• The specific solution is

• After 10 seconds,
DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 19
 • Terminal velocity

2. A 60 kg sky diver jumps out of a plane with no initial velocity and proportionality constant =0.8. Determine his
velocity after 10 seconds if and turbulent flow is presumed. What will his terminal velocity be?

• You must first calculate

• The substitution you will use is

• The general solution is
• You should find that results in
• The specific solution is
• After 10 seconds
• Terminal velocity (steady state)

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 20

8.3 Draining water from a container

In building the model for the rate of change in the water level, we must first review the principles of
conservation of mass and conservation of energy (we will assume a constant fluid density, constant container
surface area and draining tap area )

• Conservation of mass: the rate of change in the volume of the fluid in the container = rate of volume
loss through the tap opening (volume loss is the product of the velocity and opening area )

The volume of the fluid in the container is decreasing with time as the height is changing with time as
. Substituting in the above:

• Conservation of energy: The potential energy of a small mass of water at height is

and when the fluid flows out of the opening it has kinetic energy . From the conservation of

energy,

• Putting it together: substituting into the right hand side of relates

as or rather

This differential equation is only suited to separation. Since are all constants, we can combine
, then (5) is the more manageable:

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 21

Step 1: Separate the variables

Step 2: Integrate

Step 3: The general solution

Step 4: To calculate we use

Step 5: The specific solution is

Example

A cubical tank (sides 2 meters) is filled with water to height . The water is drained through a
circular hole with diameter .

Solve the differential equation modelling the draining process ( ). How long does it take to empty the
tank?

To start, the radius of the drainage opening is , therefore,

and

Step 1: Equation (5) simplifies to

Separate variables:

Step 2: Integrate

Step 3: General solution

Step 4: Calculate using

Step 5: Specific solution which simplifies to

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 22

Step 6: The container is empty when

Lets confirm this solution using Matlab:

syms h(t)
dh=diff(h); p=pi*(0.1)^2 ; A=2*2; g=9.8;
eq=dh==-(p/A)*sqrt(2*g)*h^0.5;
h(t)=dsolve(eq,h(0)==1.3);
h=vpa(h,3)
h(t) =

Note: the second entry of is the solution (why not the first one?)

The plot:

syms t
h(t)=( (-0.03*t+2.28)/2 )^2;
fplot(h,[0,77])

title('Draining Water')
xlabel('time')
ylabel('height')

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 23

Paced Exercise (8D)

1. A rectangular container with dimensions: and has initial water level .

How long will it take for this tank to empty through a circular drainage hole with radius . For which value of
will this tank be empty in 2 minutes? (Assume )

• You should find that

• The differential equation is which simplifies to

• General solution

• Specific solution
• The tank must be empty when . You should find that

2. A car fuel container has (approximate) rectangular dimensions: . If the initial fuel
level is , how long will it take for this tank to empty through a circular drainage hole with diameter
? (Assume )

• Since the formula contains the area fraction , you can keep the measurements in terms of (as long

as have the same units) as the units will cancel. Otherwise, you can convert to meters: You
should find that

• The differential equation is which simplifies to

• General solution
• Specific solution (when calculating , must be in meters )

• The tank will be empty when . You should find that

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 24

8.4 Heat transfer

Newton's law of cooling states that the rate at which a body cools from an initial temperature is proportional
to the difference in the body's (changing) temperature and the regulated constant temperature of it's
surrounding environment

The temperature can be in Kelvin (K) or in degrees. The order inside the bracket indicates the direction
of energy flow in the system. While , the energy flow will be away from the body and is negative (as

the temperature of the coffee decreases) and vice versa if . Both separation of variables and can be
used to solve the heat equation.

Separate variables:

Step 1: Separate variables

Step 2: Integrate

Step 3: The general solution is

Step 4: Calculate using

The updated general solution is

One more condition will be required to calculate .

Integrating factor:

Step 1:

Step 2:

Step 3: Put the solution together

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 25

Step 4: Calculate

Prove that the two solutions and are the same.

Example

A body is at temperature 400K and is submerged in an environment with constant temperature 300K. After 5
minutes the body is at 360K. How long will it take for the body to cool to 310K? What is the steady state
temperature?

Here and we will interpret time in terms of minutes and temperature in terms of
kelvin.

Note: If you choose to do this problem in terms of seconds, you must convert 5 minutes to 300 seconds, that is,
after 300 seconds the body is at 360 K. If you choose to work with degrees you must convert kelvin to degrees
by subtracting 273.15 from the kelvin amount, that is, .

Step 1:

Step 2: Integrate

Step 3: General solution

Step 4: Calculate using

Calculate using

Step 5: Specific solution

Step 6: Solving for T=310K:

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 26

Step 7: The steady state is determined by . This makes sense as eventually,

the coffee will have the same temperature as the environment.

syms T(t)
eq=diff(T)==-0.1022*(T-300);
T(t)=dsolve(eq,T(0)==400)
T(t) =
syms t
fplot(T,[0,100])

title('Heat Equation')
xlabel('time')
ylabel('Kelvin')

Paced Exercise (8E)

A cup of coffee is freshly made at 10 am Monday morning with temperature 99 degrees Celsius. It is put on the
breakfast table outside where the temperature is 16 degrees Celsius. After 2 minutes the coffee cools to 80
degrees Celsius. How long will it take for the coffee to cool to 60 degrees Celsius? What is the steady state
here?

• The heat equation is

• The general solution is

• The value of c is

• The value of k is
DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 27
 • The specific solution is

• When minutes

• The state is (make T subject of the solution) , therefore, the steady state is

2. A forensics team arrives at a crime scene with a dead body at 10 pm and measure the temperature of the
body to be . One hour later they measure the temperature again and found that it is . A healthy
person's temperature is on average . The room the body was found in is kept at a constant temperature
. At what time did this person die?

• Start the clock at 10 pm and assume with . When and .

• The general solution to is

• Using you should find that

• Using you should find that
• The specific solution is
• When did this person die: When , what is

. This is approximately 4 and a half hours before the forensic

team arrived at 10 pm. This person therefore died at 5:30 pm.

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 28

8.5 RL circuits
In this section, we will only consider circuits with an inductor , resistor and an energy source

DC: direct current (battery) AC: alternating current (plugs)

Kirchhoff's voltage law states that for a closed loop circuit:

sum of the voltage rises = sum of the voltage falls

The voltage fall across circuit components in terms of the current are modelled as follows:

• Across the resistor

• Across the inductor

The voltage rise is modelled as

• Across the power source

where

Voltage is measure in Volt ( )

Current is measured in Ampere ( )

Resistance is measured in ohm ( )

Inductance is measure in Henry ( )

Before the switch is opened, the current in the circuit is zero. Kirchhoff's law states that

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 29

a) When the power source E is constant (DC)

DC stands for direct current which results from a constant value of such as the voltage associated with
a battery. The differential equation is

it is both separable and suitable to :

Separable:

Step 1: Separate variables

Step 2: Integrate

Step 3: General solution

Step 4: Calculate when

Step 5: Specific solution

Integrating factor:

Step 1:

Step 2: IF

Step 3: Put the solution together which simplifies to

Step 4: Calculate when

Step 5: Specific solution

• Prove that the two solutions and are the same.

(Hint: Make the state subject of the equation in both cases)

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 30

Example

i) Calculate the current in the RL circuit below as a function of time, that is, solve the differential equation
modelling the voltage drops across the circuit components when the switch is turned on. Time zero is when the
switch is turned on, therefore . Calculate an expression for the current using separation of variables.
Calculate the steady state current for this circuit.

Step 1:

Separate variables:

Step 2: Integrate

Step 3: General solution

Step 4: Calculate when

Step 5: Specific solution

Step 4: the steady state in this circuit . The negative exponential therefore

Confirm the solution using Matlab:

syms i(t)
di=diff(i);
eq=di+5*i==20;
i(t)=dsolve(eq,i(0)==0) % no initial conditions were given
i(t) =

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 31

ii) Calculate the solution to i) using the method.

Step 1:

Step 2:

Step 3: Solution which simplifies to

Step 4: Calculate using

Step 5: Specific solution

b) With a time varying (AC) power source E(t)

AC stands for alternating current which results from an alternating (time varying) power sourse such as the
voltage associated with our wall sockets.

For the format, Q represents a power source. If the power source is constant, both the method and
separation of variables can be used to calculate the solution. If the power source is not constant, only the
method can be considered. These cases will usually also require DI integration when putting the solution
together.

In general, an alternating power source is represented as a sine function . The source is

specified by its frequency and root mean square voltage . From this information:

•
•

Examples

i) Given the AC circuit below. Calculate the steady state current by solving the differential equation modelling the
voltage drop in this circuit. The initial current is zero.

For a start we must convert from :

To express the power source as the sine wave

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 32

 • We must calculate the angular velocity using the frequency :
• The 10 V shown is the RMS value from which the amplitude can be calculated as:

Step 1:

write as:

Calculate using DI integration:

Step 4: the general solution

Step 5: initial conditions

Step 6: specific solution

syms i(t)
eq=0.01*diff(i)+5*i==14.14*sin(377*t);
DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 33
i(t)=dsolve(eq,i(0)==0)

i(t) =

i=vpa(i,3)
i(t) =
syms t
fplot(i,[0,0.5])
ylabel('i')
xlabel('t')

Note that the negative exponential (called the transient state) will go to zero as time progresses but the
sinusoidal part (called the steady state) will persist with time.

Paced Exercise (8F)

1. The power source and components of an RL circuit is given below. Model the differential equation modelling
the current in this circuit in and solve it.

The differential equation won't be separable with E not being constant and you will have to use the IF integrating
factor method. The integral will require DI.

• The differential equation is

• You should find that

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 34

 • You should find that the solution formula is
• You should find that (choosing )

• The general solution is

• The value of
• The specific solution is

The steady state is

which can be converted to a shifted sine or cosine:

2. The inductance of an RL circuit is , the resistance is and the power source is .

Calculate the steady state of the current in this circuit. It is given that .

• The equation is

• You should find that

• You should find that the solution formula is

• You should find that

• The general solution is

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 35

 • You should find that

• The specific solution is

• The steady state is

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 36

8.6 Population modelling
a) Unlimited population growth

Suppose that at time zero, a population has density and that the population increases at a rate
proportional to the difference between its birth rate and death rate as

Assuming and to be constant is a simplification that ignores many biological processes but this choice will
lead to an analytically solvable model. The above model can further be simplified to

This differential equation passes both the separable and IF tests :

Separation of variables:

Step 4: In calculating , you must substitute in the general solution

The specific solution is

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 37

Integrating factor:

Step 1: write as and identify

Step 2:

Step 3: put solution together which simplifies to

Step 5: In calculating , you must substitute in the general solution

The specific solution is

Example

A specific bacterial population modelled by exponential growth has density 0.022 units at time zero and doubles
its density after 10 seconds. What will the density be after 1 minute? Use separation of variables.

We were given that . These will be necessary to calculate

.

Step 1: The growth is unconstrained exponential

Step 4: Calculate by substituting

update the general solution:

Now substitute the additional condition

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 38

Step 5: The specific solution is

Step 6: After 60 seconds (1 minute)

syms k y(t)
dy=diff(y);
eq=dy==k*y;
y(t)=dsolve(eq,y(0)==0.022)
y(t) =

syms t
y(t)=0.022*exp(0.07*t);
fplot(y,[0,100])

title('Exponential growth')
xlabel('time')
ylabel('population density y')

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 39

b) Density dependent growth

Exponential growth as described in subsection a) is biologically unrealistic. The observation is that populations
exponentially grow until they reach a density limitation/carrying capacity (population cap) where if they exceed
this density, they will run out of resources forcing population densities to fall back to values below . When
resources recover, the population starts growing again until it again exceeds and the process repeats over
and over again. This process is modelled by the logistics model

is the intrinsic growth rate which is the maximum growth rate this population is capable of. Note from (2) that

• While , the growth rate is positive and increases

• When reaches saturation ( ), the growth rate is zero and the population remains constant
• When , the population depletes its resources and the growth rate is negative resulting in a
decrease in population density.
• When population densities are very small in comparison to , , the logistics model (2) can
be expressed as

which is the unlimited exponential population model (1) with

As the density increases, the density cut off regulation activates and the model restores to .

The growth rate adjusts to keep the population within their resource limitation with rapid growth when densities
are low and slowing down as densities come closer to .

Equation (2) is separable as it can be separated into a side and a side. The logistics equation is however not
friendly as it contains a and will not fit the format . For simplicity we will scale :

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 40

Calculating the constants A and B:

Calculating will require an initial condition and another condition associated with the growth. We were
only given the initial condition which will be used to determine :

The updated general solution is then

Example

A specific bacterial population modelled by density dependent growth has density 0.022
units at time zero and doubles its density after 10 seconds. What will the density be after 1 minute?

We were given that . We will need these values to calculate

The growth is constrained exponential so we use the logistics equation ( )

Step 3: The general solution is

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 41

Step 4: Calculate by using the initial condition

Update the general solution:

Calculate by using the additional condition

Step 5: The specific solution is

Step 6: After 60 seconds (1 minute)

syms r y(t)
dy=diff(y);
eq=dy==r*y*(1-y);
y(t)=dsolve(eq,y(0)==0.022)
y(t) =

syms t
y(t)=0.023*exp(0.07*t)/(1+0.023*exp(0.07*t));
fplot(y,[0,150])

title('Density regulated growth')

xlabel('time')
ylabel('population density y')

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 42

Note that this model won't exceed a density of 1 but the previous model will keep growing indefinitely.

Paced Exercise (8G)

A culture grows according to the density dependent growth model (logistics with ). The initial
density is units and the culture grows to 0.5 units after 30 seconds (0.5 minutes). How long will it take
this culture to reach a density of 0.9 units?

• For consistency, time will be measured in minutes

• This is the same differential equation used in the previous example with general

solution

• You should find that the initial condition determines that

• You should find that the additional condition minutes and determines that
• The specific solution is then

• When units, minute

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 43

8.7 Mixing problems
A dam initially contains litres of water containing kg/litre of pollutant (uniformly concentrated). A river
flows into this dam containing kg/litre of pollutant at a rate of litres/minute. The dam drains into another
river at a rate of litres/minute. If represents the kg pollutant content of the dam at time then the rate at
which changes is given by

rate of change in rate of entering the dam rate of is exiting the dam

where

which leads to the differential equation

Examples

i) A near empty dam initially contains litres of water containing no pollutants. A river flows into this dam
containing kg/litre of pollutant at a rate of 2 litres/minute. The dam drains into another river at a rate of
1 litre/minute. Calculate the concentration of pollutant in the dam after one hour.

Equation (7) simplifies to

This equation can only be solved with the IF method:

Step 1: write as

Step 2:

Step 3: solution:

DI integration for :

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 44

general solution

Step 4: initial condition

Step 5: specific solution

Step 6: After 1 hour (60 minutes)

syms x(t)
dx=diff(x);
eq=dx+(1/(t+100))*x==4*exp(-0.5*t);
x(t)=dsolve(eq,x(0)==0)

x(t) =

eval(x(60))
ans = 5.1000

syms t
fplot(x,[0,60])

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 45

ii) Change the draining rate in question i) to 2 litres per minute (same as the inflow). How does this change the
pollutant concentration in the dam after one hour?

Equation (7) simplifies to

This equation can only be solved with the IF method:

Step 1: write as

Step 2:

Step 3: solution which simplifies to

Step 4: initial condition

Step 5: specific solution

Step 6: After 1 hour (60 minutes)

syms x(t)
dx=diff(x); eq=dx+(1/(50))*x==4*exp(-0.5*t);
x(t)=dsolve(eq,x(0)==0)

x(t) =

eval(x(60))
ans = 2.5100

syms t
fplot(x,[0,60])

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 46

DR ERASMUS MEM12A, CMA12A, EMA125C CHAPTER 8 SECTION 1 47