A P P L I E D PH Y S I C S : M AT H E M AT I C A L M O D E L L I N G

PH Y S 1 0 3 4

Course Notes by
Warren Carlson and Jarrod Shipton

2025
Preface

These are the course notes for PHYS1034. These notes are intended to complement the lectures
and other course material, and are by no means intended to be complete. Students should consult
the references at the end of this document and those provided by the course instructor to find
additional material, and different views on the subject matter.
This material is under development. Please would you report any errors or problems (no
matter how big or small). Any suggestions would be gratefully appreciated.

School of Computer Science and Applied Mathematics,

University of the Witwatersrand,
Johannesburg,
South Africa

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivs

License. To view a copy of this license, visit,

• [Link]

Send a letter to Creative Commons, 543 Howard Street, 5th Floor, San Francisco, California, 94105,
USA. In summary, this means you are free to make copies provided it is not for commercial
purposes and you do not change the material.

i
Contributors

This course material was prepared with the help of the following people, Eric Mubai, Ezekiel
Lengaram, Kendall Born and Abdul Hamid Carrim. The authors of this manuscript are indebted
to these individuals for their help in preparing this material.

iii
Contents

Preface i

Contributors iii

Contents iv

1 Introduction 1
1.1 The Modelling Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Mathematical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Construction of Mathematical Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.1 Formal Model Construction Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.2 Conceptual and Mathematical Formulations . . . . . . . . . . . . . . . . . . . . . . 7
1.3.3 Dimensional Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2 Discrete Models 13
2.1 Modelling Discrete Change . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.1 Discrete Economic Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2 Formal Discrete Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3 General Discrete Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1 Discrete Population Growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3.2 Discrete Supply and Demand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.3 Critical Values and Equilibria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3 Continuous Models 31
3.1 Modelling Continuous Change . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.1.1 Euler’s Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Formal Continuous Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.1 Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.2 Euler’s Number and Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3 General Continuous Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.3.1 Continuous Population Growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

iv
 3.3.2 Continuous Supply and Demand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3.3 Newton’s Model of Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

Appendices 47

A Units and Dimensions 49

B Software 51
B.1 Getting Started . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
B.2 Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
B.3 Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
B.4 Scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Bibliography 61

v
Chapter 1

Introduction

Mathematical Modelling is a general procedure in the Mathematical Sciences wherein ideas about
observed processes are presented and interrogated using mathematical principles and techniques.
There may exist many different conceptual formulations of observed processes and so there may
be many equivalent models for those processes. In the sections that follow, we shall discuss some
of the conceptual and technical aspects of building models that describe physical phenomena.

1.1 The Modelling Process

To gain an understanding of the processes involved in mathematical modelling, consider two
worlds which can be depicted with the properties below,

• Real-world system:

– Observed behaviour

• Mathematical World:

– Models
– Mathematical operations and rules
– Mathematical conclusions

Suppose we want to understand some behaviour or phenomenon in the real world. We may
wish to make predictions about that behaviour in the future and analyse the effects that various
situations have on it. For example, when studying the populations of two interacting species,
we may wish to know if the species can coexist within their environment or if one species will
eventually dominate and drive the other to extinction. In the case of the administration of a drug
to a person, it is important to know the correct dosage and the time between doses to maintain a
safe and effective level of the drug in the bloodstream. In order to construct and use models in
the mathematical world to help us better understand real-world systems, it is important to gain
an understanding of how we link the two worlds together. To begin, we would first have to define
what we mean by a real-world system.

1
Definition 1 (System) A system is an assemblage of objects with some interaction or
interdependence.

Definition 2 (Model) A model is a conceptual representation of a given object or process that

captures some collection of aspects of that object or process.

The modeller is interested in understanding how a particular system works, which features
cause changes in the system, and how sensitive the system is to certain changes. They are also
interested in predicting what changes might occur and when they occur.

Remark 1 Physical systems can be very complicated, with subtle intricacies that are either not
well understood, or difficult to describe mathematically. Therefore, the impetus of mathematical
modelling is to produce an approximation that is sufficient to describe the important parts of the
system. This means that the domain of applicability of a model can be very narrow compared to
the conceptually perfect mathematical representation of such a system.

Suppose the goal is to draw conclusions about an observed phenomenon in the real-world.
One approach would be to conduct some experiments or real-world behaviour trials to observe the
effects thereof. This would be the process of drawing conclusions from the real-world behaviour.
However, this is not always a viable method. In many cases there may be prohibitive financial
or social costs which inhibit our ability to run such experiments. For example, determining the
concentration level at which a drug would be fatal, or studying the effects of a failure in a nuclear
power plant. It may even be the case where it is not possible to carry out such experiments, for
example, investigating a specific change in composition of the ionosphere and its corresponding
effect on the polar ice caps. We may also be interested in expanding our conclusions beyond that
of the very specific trial we have run. For example, we may have run a trial in Johannesburg with
a temperature of 32◦C, humidity 42%, and other very specific conditions. Even if we do succeed
in drawing a conclusion from the real-world behaviour under these very specific conditions, we
may not have necessarily explained why the particular behaviour has occurred. This gives us
motivation to derive indirect methods of studying real-world behaviour.
The alternative approach to draw conclusions about an observed phenomenon in the real-
world can be outlined as follows. The first step would be to make observations of the system we are
attempting to model. When making observations about real-world behaviour we will usually not
be able to identify or consider all the factors associated with the phenomenon. As such we would
make simplifying assumptions that eliminate some of the factors. Take for example a situation
such as the effects of a failure in a nuclear power plant near a city such as Cape Town. Initially it
might be found that we do not need information about the humidity levels and we would neglect
this factor in our model. Our next step would then be to conjecture tentative relationships among
the factors to create a model. After this, we would construct a model using these relationships.
After having constructed a model we would move to the next step of applying a comprehensive
mathematical analysis to the resultant model. We must be careful at this step to note that any
conclusions we obtain from this model are strictly to do with the model and not the real-world

2
behaviour, since we have made simplifications to the model, and our observations on which we
have based our model could contain errors. In summary, we have a rough outline of the procedure:

1. Through observation, identify the primary factors involved in the real-world behaviour.

2. Conjecture tentative relationships among the factors.

3. Apply mathematical analysis to the resultant model.

4. Interpret mathematical conclusions in terms of the real-world problem.

We shall consider several scenarios in detail using this approach.

1.2 Mathematical Models

Let us begin with a useful definition.

Definition 3 (Mathematical Model) A Mathematical model is a mathematical construct designed

to study a particular real-world system or phenomenon.

There are two approaches we can take when modelling a phenomenon mathematically. We
can either take pre-existing models of well studied phenomenon and build on these models, or
we can construct a completely new model for our special phenomenon.
When constructing a model it can often be the case that we reach a solution so complex and
intractable that there is little hope of analysing or even solving the model. This can be caused by
a number of conditions, and can often cause our model to have very little utility in understanding
the phenomenon. These complexities may arise when attempting to use a model given by a system
of equations. The problem may be so large that it is impossible to capture all the information in
a single mathematical model. An example of such a problem is the global effects of population
interaction, resource use and pollution. In such cases we can attempt to replicate the behaviour
directly by conducting various experiments in order to collect data and analyse the phenomenon
using statistical techniques or curve-fitting.
It is important to note that there is a trade-off among factors in selecting a specific type of
model. Three of these factors are listed below.

Fidelity: The preciseness of a model’s representation of reality.

Cost: The total cost of the modelling process.

Flexibility: The ability to change and control conditions affecting the model as required data are
gathered.

First let us consider fidelity. Observations made from real-world behaviour would be a direct
representation of reality and thus we would have the greatest fidelity from these types of models.
After this we would expect that experiments are the next best since they are scaled down versions

3
of the real phenomenon, followed by simulations which may not account for all the intricacies of
the system. Constructed models would be tailored for the specific problem and would probably
outperform pre-existing models.
The cost of real-world observations would generally be high since gathering data is an
expensive process. The cost of experiments and simulations would also generally be high since
the equipment required is usually expensive. Constructed models will generally have some cost
associated to the tailoring to the specific problem and pre-existing models would generally be the
cheapest.
The flexibility of constructed models would generally be the greatest since we can vary the
models easily. Pre-existing models can be used as a base and built upon and are generally available.
Experiments and simulations would only be as flexible as the technology, which is available at the
time, allows for. Finally, real-world observations would generally be the least flexible since we
cannot change the environment in which we make these.

1.3 Construction of Mathematical Models

Modelling is viewed as a process or translating real-world behaviour into understandable
conclusions, and there are various forms a model can take. In this section we will outline possible
steps that can be helpful in the construction of a mathematical model. It is a useful idea to
approach this with the mindset that this is an iterative process. We will often find that we missed
an important factor at a later step and will have to go back to a previous step to improve the
model. The concept of trial and error is central to the modelling process and model building.

1.3.1 Formal Model Construction Process

We shall follow an iterative process to build and verify each model. In each iteration we shall
implement an ordered sequence of tasks, namely

1. Problem Identification

2. Assumption Specification

3. Variable Classification

4. Solution and Interpretation

5. Verification

6. Implementation

7. Maintenance

Each of these exists for a specific reason, and at any point in the sequence it may be necessary to
review a previous step before continuing to the next.
Let us interrogate each of these tasks in some detail.

4
Problem Identification: First, identify the problem. This amounts to asking the question "What is
the problem you wish to explore?" This is usually not an easy step, since in real-life situations
there is usually no one who will hand you the problem that is to be solved mathematically.
This might require to sift through large amounts of data and identify some aspect of the
situation to study. It is important that the problem formulation be sufficiently precise to
allow for the translation of the verbal description of the problem into mathematical symbols.
This process is made clear in the subsequent steps of the modelling process.

Assumption Specification: Once the problem is identified, the next step is to make assumptions
about the problem. Usually we cannot identify all the factors affecting the problem in a
usable mathematical manner. So, we must simplify the task by reducing the number of
factors we need to consider to draw conclusions. Then, we must find which relationships
between the remaining factors could exist in the model. This both reduces the complexity of
the problem and gives us insight into which factors are necessary when drawing conclusions.

Variable Classification: This involves identifying all the variables, failing that, as many as
possible, and listing them. Any of the variables that we want to explain are considered to be
the dependent variables. There can be several dependent variables in a single model. The
remaining variables are considered to be the independent variables. All the variables
should be classified as dependent, independent or neither.

An independent variable may be neglected for one of two reasons. First, if the relative
influence of a given independent variable is small compared to another, then we might
choose to neglect it. The second reason could be that it affects the other alternatives in a
similar way. An example of such variables could be bullet colour and hat size
considerations when modelling wildlife hunting. While it may be an important factor in the
real-life situation, it would not have much impact on the model and may add unnecessary
complexity.

Next, we must determine the relationships among the remaining variables. Some problems
are sufficiently complicated that identifying direct relationships is difficult. In these cases,
it is productive to study some smaller model with only the variables that have obvious
direct relationships. After building relationships between each sub-collection of variables,
a complete and comprehensive model is built from the constituent sub-models.
Considering relationships among variables such as proportionality, will aid in
hypothesising the relationships among the different variables.

Solution and Interpretation: Once the relationships among the different variables are
determined, there remains to solve and/or interpret the model. At this step we take all the
sub-models we made in the previous step and ask what these models are telling us. These
sub-models can come in various forms and may consist of mathematical equations or
inequalities. The problem statement these relate to could be to find the behaviour of a
real-world situation given certain conditions, or to find the best solution, or optimal

5
 solutions to the model. In some cases we might find at this stage that we have a model that
is so complex, or vague that we cannot solve it. At this point it is often a good idea to revisit
our assumptions and possibly add or replace factors. It might be necessary to redefine the
problem statement entirely.

Verification: The next step is to verify the model. The model must be properly tested before it
can be used. Before embarking on this specific task, we consider whether the model does
indeed supply an answer to the problem statement, and whether it is usable. If a model
requires input data that is unavailable, then it is not usable. Finally, we consider whether
the model is sensible.

Once a model passes these checks, it can be tested against empirical observations. The
design and implementation of experiments to test a given model is outside of the scope of
this course, however it is useful to add some general commentary. Tests of models should
apply to the independent parametre ranges that are appropriate for the model in question.
As an example, Newton’s Second Law of motion is applicable where the momentum of an
object is approximately equal to the product of the mass of the object and velocity of the
object. This is violated when the velocity of the object approaches the speed of light.

Note that models are approximate descriptions of a system and not necessarily a governing
law. Unlike in pure Mathematics, the proof of a theorem is not merely demonstrated by a
single example, but rather a general statement which requires only formal mathematical
statements and is independent of physical measurement or observation. Similarly, a model
can only be demonstrated to be appropriate for a given system; we cannot in general
extrapolate any conclusions about a given model to broad generalisations from the evidence
we gather. A model does not become law just because it is verified repeatedly in certain
conditions. Instead, it is only the reasonableness of a model that is corroborated by the data
we collect.

Implementation: We want to explain our model in terms that the the people who use the model
to make decisions will understand. Furthermore we will want to make it as simple as possible
to use, otherwise these users will choose another model.

Maintenance: The model we have designed is often derived for a specific problem which is
identified in a problem statement or from the assumptions we have made. There may be a
number of reasons the model might need updating. For example, advances in measurement
or data collection and availability might allow previously unknown or immeasurable factors
to be included, or possibly the elimination of a previously important factor from a given
model.

There are two main types of Mathematical models that we shall like to build in the chapters to
follow. These are the so-called Discrete Mathematical Models and the Continuous Mathematical
Models. For each type of model, we shall follow the same development procedure, as outlined

6
above. We shall consider each in turn in the chapters to follow, and show the inter-relation of
discrete and continuous mathematical models in different scenarios.
A central part of checking the consistency of a mathematical model is the kind of output given
by the model. The output of a model usually describes some physical quantity. Physical quantities
carry units. The first and simplest consistency check that can be made is validity of the physical
output units. Since this is such an important part of the model building process, we take some
time to discuss it in detail.

1.3.2 Conceptual and Mathematical Formulations

A conceptual formulation invariably precedes any mathematical formulation of a situation. The

conceptual formulation can be expressed by means of a graph, diagram, physical model, or it
can reside "purely in the mind". A verbal language is used to describe these concepts which
are then translated into formal mathematical statements wherein the idea that motivates the
model is made explicit in mathematical language. We shall touch briefly on this process of verbal
language to mathematical language translation. Consider first the language used to describe the
relationship among variables. Table 1.1 presents these translations in detail. Similarly, we can
describe the change in the value of a function with respect a given parameter in that function.
Table 1.2 presents these relationships explicitly.

Table 1.1: Mathematical and verbal statements for the relationships among variables. Here, a , b ,
c and d are constants, k is a positive constant (0 < k ) and s , x , y and z are variables.

Mathematical Statement Verbal Statement

y =kx y is directly proportional to x
y = kxn y is directly proportional to the n -th power of x
y = kx y is inversely proportional to x
y = xkn y is inversely proportional to the n -th power of x
z =kx y z is proportional to x and y
s =ax +b y −c z +d s is linearly related to x , y and z

Table 1.2: Mathematical and verbal statements for the relationships among variables. Consider a
model of a physical situation where x dependends on t . We write x = x (t ).

Mathematical Statement Verbal Statement

x (t
) The value of the dependent variable, x , at any time t .
dx
€ d2t Š The rate of change of x at any time t .
d x
dt 2 The rate of change of the rate of change of x at any time t .

7
1.3.3 Dimensional Analysis
Consider the equation F⃗ = m a⃗ from Newtonian Mechanics. This physical law is taken to be true
as long as the units used in the measurement are consistent. If mass was measured in kilograms,
acceleration in meters per second squared and force in Newtons, then we could accept this as
true. For comparison, suppose we use the system of units where mass is measured in slugs and
acceleration in feet per second squared, then it would be incorrect to measure force in Newtons,
instead it should be measured in pounds to keep the units of measurement consistent. Similar
steps should be taken to ensure our units in our variables are all considered correctly.
The case of keeping with the example of F⃗ = m a⃗, let us use the three physical quantities related
to this, those being mass, length, and time, each with corresponding dimensions M , L , and T
respectively. There are other quantities which have dimensions not included in this example,
however many other quantities can be defined in these dimensions. Area, which is the product of
two lengths can be seen to have the dimension of L 2 . Velocity which is the ratio of the distance
travelled to the time taken has the dimension LT −1 . Area and velocity, being but two examples,
are simple products of these dimensions. Each product has dimensions M n L p T q , where n , p ,
and q are real numbers which can be positive, negative or zero. If a product is missing one of
these dimensions, the power of the dimension is zero. When all the powers are zero, the product
is said to be dimensionless.
One cannot add two values with different dimensions. For example, one cannot add a velocity,
−1
LT , to an acceleration, LT −2 . These are said to be dimensionally incompatible. If one wanted to
add a force to the product of mass and acceleration however, this would be possible, since these
are dimensionally compatible. If an equation exists that is true regardless of the units of measure,
the equation is said to have dimensional homogeneity.

Theorem 1 (Buckingham’s Theorem) An equation is said to be dimensionally homogeneous if

and only if it can be put into the form

f (Π1 , Π2 , . . . , Πn ) = 0. (1.1)

where f is some function of n arguments and {Π1 , Π2 , . . . , Πn } is a complete set of dimensionless

products.

When we have a function of the form of (1.1) we can find the product containing the dependent
variable and solve for it with a function in terms of the other products. That is, if our dependent
variable is in Πn we then have

Πn = H (Π1 , Π2 , . . . , Πn−1 ) = 0. (1.2)

We can then solve (1.2) for our dependent variable. The steps taken in dimensional analysis are
outlined below.

Step 1: Decide which variables are important to the problem being investigated.

8
Step 2: Determine the set of dimensionless products, {Π1 , Π2 , . . . , Πn }, among the variables.
Ensure that the dependent variable only appears in one of the dimensionless products.

Step 3: Check that all of the products found are indeed dimensionless and independent.

Step 4: Apply Theorem 1 to produce all the dimensionally homogeneous equations that are
possible among the variables.

Step 5: Solve the equation in Step 4 for the dependent variable.

Remark 2 The argument that motivates Theorem 1 is simply stated as follows. Any function that
takes as input a parametre that carries units must generate an output that also carries units.
Consider the function
X
∞
xn
f (x ) = .
n=0
n!
Clearly, any unit carrying input to this function will produce a sequence of terms carrying distinct
units in summation. The units of these terms are incompatible. Non-dimensionalising the input
parameter will remove the dimensional dependence. Therefore, dimensionless inputs to this
function will always produce dimensionless outputs. This argument generalises directly to any
function and so Theorem 1 applies to general mathematical formulations.

A collection of useful units and their corresponding dimensions are listed in Table A.1 in the
Appendix. The next example demonstrates how dimensional analysis is used in the modelling the
period of a simple pendulum.

Example 1.1 (The Period of a Simple Pendulum) Consider the motion of a simple pendulum.
How might we describe the period of its motion about the equilibrium position.

Problem Identification: We want to find a relationship among the important descriptive

parametres of the pendulum that will determine the time taken to complete a single
oscillation.

Assumption Specification: We do not consider any special arrangement of the pendulum. We have
specific knowledge of its motion, in particular, the simple pendulum should move with simple
harmonic motion and the position of the pendulum should be defined by a periodic function
such as sine or cosine. In addition, we should assume that the length of the pendulum ℓ, its
mass m and the effect of gravity g on the pendulum are all fixed. For a first iteration of the
model, we expect that additional complicating factors, such as resistance to motion, can be
neglected.

Variable Classification: We can collect all the variables that govern the process. For a simple
pendulum, we need to consider the mass m, its length ℓ, the initial angle, θ0 marks the start
of motion, and the period τ between swings under the effect of gravity g . The dimensions of
these variables are presented in Table 1.3.

9
 Table 1.3: The dimensions of the variables associated with the motion of a simple pendulum.

Variable Dimensions
m [M ]
ℓ [L ]
θ0 1
τ [T ]
L
g T2

Solution and Interpretation: We can analyse the dimensions of the problem using Buckingham’s
method to obtain a relation between the period and the other variables. First, note that the
input to the periodic function that will define the position of the pendulum should be unitless.
Collect the parametres into a dimensionless product,
 a b c d e
m ℓ θ0 τ g = 1 = M 0 L 0 T 0 .

Then, substitute in the dimensions of each parametre,

•  ‹e ˜
a b c d L
M L 1 T = 1 = M 0 L 0T 0.
T 2

Now, equate exponents,

M :0=a
L :0= b +e
T : 0 = d − 2e

It is now clear that

a =0 c =c b = −e
d = 2e
e = e.

Clearly, the motion of the pendulum is independent of mass since a = 0 means that mass
does not contribute to the dynamics of the problem. Now, we construct our two dimensionless
collections to be based on c and e . Setting c = k yields

Π1 = θ0k = constant.

Setting e = k we get,
 g k
Π2 = ℓ−k g k τ2k = τ2 .
ℓ
Make the assumption that f (Π1 ) = Π2 that is,
 v
g k tℓ
f (θ0 ) = τ2 or τ = h (θ0 ) ,
ℓ g

10
 for some function h . Note that, as stated in the Assumption Specification of this example, the
motion of the pendulum is periodic and determined by either a sine or cosine function. These
functions have a periodicity of 2π. Therefore the periodic angular dependence determines
that h = 2π, so v
tℓ
τ = 2π .
g

The formal construction of the equations of motion from Newton’s Second Law of motion are beyond
the scope of this course, so the verification, implementation and maintenance steps of the model
building process for this example are omitted. Students are encouraged to consult texts on mechanics
to verify that the output of this model is correct.

1.4 Exercises
Exercise 1.1 Translate into mathematics: "r is directly proportional to x and y 2 , and inversely
proportional to z 3 ".

Exercise 1.2 Consider the damped pendulum, where the damping force is given by F⃗ = −k v⃗ where
k is the damping coefficient, v⃗ is the velocity of the pendulum and F⃗ is the damping force applied
to the pendulum. Determine an expression for the period of the pendulum for its motion about the
equilibrium position as a function of the damping factor k .

Exercise 1.3 Consider a pipe of diameter D and length ℓ between pressure taps. The change in
pressure ∆p along the pipe has been experimentally studied and for laminar flow the change in
pressure is directly related to ℓ. You know that the change in pressure ∆p is a function of D , ℓ as
well as the fluids velocity,V , and the fluids viscosity, µ. Using dimensional analysis, deduce how
the change in pressure varies with the pipe’s diameter.

Exercise 1.4 Consider a cylinder of diameter D floating upright in a fluid. The cylinder is perturbed
in the vertical direction and begins to oscillate about its equilibrium position with some frequency,
ω. You are given that the frequency is a function of the cylinder’s diameter, mass, m, and the specific
weight, γ, of the fluid. Using dimensional analysis, deduce what would happen to the frequency if
the mass were increased.

11
Chapter 2

Discrete Models

The first models of interest to us are those where the state of the system changes in discrete steps,
such that at each state of the system has a single, well defined nature that is distinct from all others.
The set of all states of the system is then a collection of step and value pairs. If the system follows
some set of rules that dictates how it changes from one state to the next then the idea of a next
state of the system implies an ordering on the states that the system occupies. We can label these
states with symbols. Suppose that the current state of a system is labelled as b0 , the next state
labelled b1 , the state after that labelled b2 and so on. After n steps, the system will have traversed
n discrete relabelling, following that given sequence. If there is a well-defined set of rules that
dictates how any one state is related to the next state in the sequence, then it should be possible
to determine mathematically the state of the system after some given number of steps m ahead
of the current state. This Chapter discusses some ideas and examples of systems that change in
discrete steps.

2.1 Modelling Discrete Change

A powerful paradigm to use in modelling change is

future value = present value + change.

Often, we wish to predict the future on the basis of what we know now, in the present, and add
the change that has been carefully observed. In such cases, we begin by studying the change itself
according to the formula

change = future value − present value.

By collecting data over a period of time and plotting those data, we can discern patterns is the
data in the model and capture trends to model the change. If the behaviour is taking place over
discrete time periods, the preceding construct leads to a difference equation. Both are powerful
methodologies to study change. In the example below we will model the change using discrete
time periods.

13
2.1.1 Discrete Economic Models
A common and relevant problem in banking is the determination of the total value of an account
that earns interest over time. The following examples present some ideas that describe this
problem and identify the general patterns associated with compound interest calculations.

Example 2.1 (Compound Interest) Consider the value of savings account initially worth R1000,
that accumulates interest paid each month at 1% per month. Suppose we want to know the value
of the account after 6 months.
Below we consider the formal process of constructing a model for increase in value of the account.

Problem Identification: Given the value of a savings account that increases in value by some
fraction of its current value at each time step, we want to build a mathematical description of
the increase in the account value over n time steps.

Assumptions: The following sequence of numbers represents the month by month value of the
account,
A = (1000, 1010, 1020.10, 1030.30, . . .).

Suppose the value at month i is a i , then the first difference are as follows,

∆a 0 = a 1 − a 0 = 1010.00 − 1000.00 = 10.00

∆a 1 = a 2 − a 1 = 1020.10 − 1010.00 = 10.10
∆a 2 = a 3 − a 2 = 1030.30 − 1020.10 = 10.20

Note that the first differences represent the change in the sequence during one time period, or
the interest earned in the case.
The first difference is useful for modelling the change taking place in discrete intervals. In
this example, the change in the value of the account from one month to the next is merely
the interest paid during that month. If n is the number of months and a n the value of the
account after n months, then the change or interest growth in each month is represented by
the n t h difference
1
∆a n = a n+1 − a n =an .
100
This expression can be rewritten as the difference equation,
 ‹
1 1
a n+1 = a n + an = 1 + an .
100 100

Solution and Interpretation: We know the initial deposit of R1000 (initial value) that then gives
the dynamical system model,

a 0 = 1000 and a n +1 = 1.01a n for n = 0, 1, 2, 3, . . .

where a n represents the amount accrued after n months. Because n represents the non-
negative integers {0, 1, 2, 3, . . .}, (2.1) represents an infinite set of algebraic equations, called a

14
 dynamical system. Dynamical systems allow us to describe the change from one period to
the next. The difference equation formula computes the next term, given the immediately
previous term in the sequence, but it does not compute the value of a specific term directly
(e.g., the savings after n periods). We would iterate the sequence to a n to obtain that value.
Performing this explicitly, for an interest rate of α%, we find

a 1 = (1 + α)a 0
a 2 = (1 + α)a 1 = (1 + α)2 a 0
a 3 = (1 + α)a 2 = (1 + α)3 a 0
..
.
a n = (1 + α)a n−1 = (1 + α)n a 0

Clearly,
a n = (1 + α)n a 0 . (2.1)
The value of an account with initial value a 0 = R1000 that accrues compound interest at a
rate of 1% compounded monthly will have a value after 6 months of
 ‹  ‹
1 6 101 6
a6 = 1 + 1000 = 1000 = R1061.52.
100 100
This means that after 6 months of interest at 1% compounded monthly, an initial value of
R1000 increases in value to R1061.52.

Verification: As can be seen from a brute force computation,

a 0 = 1000.00
a 1 = 1010.00
a 2 = 1020.10
a 3 = 1030.30
a 4 = 1040.60
a 5 = 1051.01
a 6 = 1061.52

the value of the account is correctly described by the model.

Implementation: The functional output of this model is sufficiently simple to understand via the
direct, brute force check. As such, the computation can be achieved directly via brute force
when the number of steps is small. However, for a larger number of steps, it is clear that the
computational work load is vastly reduced by the simple exponential increase described by
(2.1).

Maintenance: Again, this model is sufficiently simple that there can be no updates to the
formulation without changing the problem statement. As such, this model requires no
maintenance.

15
 It is often change that we observe, so we can construct a difference equation by representing
or approximating the change from one period to the next. To modify, our example, if we were to
withdraw R50 from the account each month, the change during a period would be the interest
earned during that period minus the monthly withdrawal, or

1
∆a n = a n+1 − a n = a n − 50.
100

In most examples, the procedure of describing the change in a system mathematically is not as
simple as is illustrated here. It is often necessary to plot the change, observe a pattern and then
describe the change in mathematical terms. That is, we would try to find

change = ∆a n = some function f .

The change may be a function of previous terms in the sequence (as was the case with no monthly
withdrawals), or it may also involve some external terms (such as the amount of money withdrawn
in the current example or an expression involving the period n ). Thus, when constructing models
representing change for discrete intervals, it is a case of approximating the function f in the
following,

change = ∆a n = a n+1 − a n = f (terms in the sequence, external terms).

This can be generalised as is in the next example.

Example 2.2 (An Investment Annuity) Annuities are often planned for retirement purposes. They
are basically savings accounts that pay interest on the amount present and allows the investor
to withdraw a fixed amount each month until the account is depleted. An interesting issue is
to determine the amount one must save monthly to build an annuity allowing for withdrawals,
beginning at a certain age with a specified amount for a desired number of years, before the account’s
depletion. This is known as a retirement annuity. For now we consider only the value of the annuity
after retirement with an initial value of R90000 and after 6 monthly withdrawals of R1000 and
with interest of 1% compounded monthly.
Below we consider the formal process of constructing a model for the remaining value in the
annuity after n steps.

Problem Identification: The problem here is to determine the value of the annuity after a given
number of payouts subject to compound interest and step wise withdrawals from the account.

Assumption Specification: Consider 1% monthly interest accruing rate and a monthly

withdrawal of R1000. The problem is to determine the value of the annuity after it has paid
out 6 planned monthly instalments of R1000 with starting value R90000 and an interest rate
of 1% compounded monthly. At each step, the value contained within the annuity changes
while all other parametres remain fixed.

16
Variable Classification: Supposing that the above mentioned parameters are the only important
variables in this system, the problem is defined by the dynamical system,

a n+1 = 1.01a n − 1000.

where the value of the annuity at the next withdrawal is proportional to the current value
minus value of the current withdrawal.

Solution and Interpretation: To solve this problem, consider a rewriting of the equation that
describes the system,
a n = (1 + α)a n−1 − c .

where α is the interest rate and c is the value of the withdrawal. We can write out the values
of the annuity as a sequence following the first withdrawal,

a 1 = (1 + α)a 0 − c
a 2 = (1 + α)a 1 − c
= (1 + α)2 a 0 − c (1 + α) − c
a 3 = (1 + α)a 2 − c
= (1 + α)3 a 0 − c (1 + α)2 − c (1 + α) − c
a 4 = (1 + α)a 3 − c
= (1 + α)4 a 0 − c (1 + α)3 − c (1 + α)2 − c (1 + α) − c
..
.
a n = (1 + α)a n −1 − c
X
n−1
n
= (1 + α) a 0 − c (1 + i )m
m =0

This brute force evaluation can be somewhat simplified by noting the following identity


x n+1 − 1 = (x − 1) x n + x n −1 + x n−2 + · · · + x + 1 .

Rewriting this we find,

x n+1 − 1

= x n + x n−1 + x n−2 + · · · + x + 1 .
x −1

Then, we find,
(1 + α)n − 1
a n = (1 + α)n a 0 − c . (2.2)
α
From this we find that the annuity value is simply the expected value of the initial investment
value compounded monthly after n months minus the value of the withdrawals corrected for
interest.

17
Verification: Below is a brute force computation of a 6 ,

a 0 = 90000.00
a 1 = 89900.00
a 2 = 89799.00
a 3 = 89696.99
a 4 = 89593.96
a 5 = 89489.90
a 6 = 89384.80,

while using (2.2) we obtain

a 6 = 89384.80.
Clearly, the brute force computation of a 6 corresponds with the value obtained using (2.2).

Implementation: As in the case with Example 2.1 the brute force implementation works, but
requires much more computational effort to evaluate compared to evaluating (2.2).

Maintenance: Suppose we made the following initial investments

A: a 0 = 90000.
B: a 0 = 100000.
C: a 0 = 110000.

The solutions to these are tabulated below.

an
n A B C

0 90000 100000 110000

1 89900 100000 110100
2 89799 100000 110201
3 89697 100000 110303
4 89594 100000 110406
5 89490 100000 110510
.. .. .. ..
. . . .

Notice that the value R100000 is an equilibrium value and once it is reached, the system
remains fixed at this value. If we start above that equilibrium, there is growth without bound.
On the other hand if we start below R100000, the savings are used up at an increasing rate.
Note how drastically the long-term behaviours are even if we had to change the starting value
by only R0.01 in either direction. We would state that this is an unstable equilibrium.
Additionally, this model can be extended to include the value of inflation on the value of the
annuity at each payout.

18
 Example 2.1 and Example 2.2 present some of the basic elements of economic models that
are commonly used in every banking system. In each case, there are some commonalities that
are worth formally generalising so that we will recognise them when they appear later. We shall
consider this generalisation in the next section.

2.2 Formal Discrete Models

In Section 2.1, we constructed models using sequences of values, parameterised by the number
of time steps between a given initial value and the value that we are interested computing. This is
a general pattern of problem formulation in the process of mathematical modelling. It is worth
formalising some of these ideas. Let’s make some crisp and clear definitions to describe the
mathematical objects that have appeared thus far.

Definition 4 (Ordered Set) Consider a set of objects S = {a , b , c , d , . . .}. Suppose S is a set and for
any two distinct elements a , b ∈ S, either a < b or b < a . Suppose also that S has a least element 0
such that 0 < c for every other element c ∈ S. We may place each element in S in a sequence from
the least element to the greatest element. Now, S is an ordered set.

Definition 4 is useful to describe things that happen in a predefined order. For example, if
when the savings account in Example 2.1 was opened, the value of the account was R1000. After
the first round of interest, after one time period, the value of the account changed to R1010. The
sequence of account values is paired with the ordered time intervals of compounding interest,
starting at the initial time. Since we know how many time intervals have passed between the
opening of the account and some later date, we can associate to each increasing time an increasing
value of the account. So, the value of the account is a time-ordered set.

Remark 3 An ordered set S defines an ordering on a collection of objects. When objects in S are put
in one-to-one correspondence with objects in another set T, then T becomes an ordered set. Now,
the set S is an ordered index on T. Indexing elements of sets is a common task and the set of natural
numbers N = {1, 2, 3, . . .} is a convenient set to serve as an index. It is also common to use the set of
non-negative integers Z+ = {0, 1, 2, . . .}.

In Example 2.1 and Example 2.2, the value of each investment is evaluated at the end of some
time interval t . Since the time intervals are of fixed size, the total time that has passed between
the initial investment and the n-th interval is n δt , where δt is the length of the interval. In each
case, the value of each account depends on n, but not on δt , where n takes only discrete values, 0,
1, 2, 3 and so on. Therefore, the mathematical machine that determines the value of each account
should take as input the initial value of each account and the discrete labels for each time interval
and produce as output the value of the account at each interval. This leads us to the following
definition.

19
Definition 5 (Discrete Function) Suppose S is an ordered, discrete set of elements X and let f be
a function evaluated on elements k ∈ S. Define

fk = f (k ),

where k is the independent variable and fk is the dependent variable. The collection (k , fk ) is
now a discrete, ordered set of points and f is a discrete function.

Remark 4 It is common to write only the indexed dependent variable fk when referring to a discrete
function f indexed by k ∈ S.

It is possible to represent a discrete function graphically as a collection of points (k , xk ) in the

plane, see Figure 2.1. Relating the different step-wise defined values in a discrete equation leads
to the following class of functional equations.

4 xk = 12

k

xk = 5049 cos π
6k
1
+ 100
k
3 xk = (−1)

k
xk = 32 − 79
2
xk

−1

−2
0 5 10 15 20 25 30 35 40 45 50
k
Figure 2.1: A graph of discrete functions.

Definition 6 (Recurrence Relation) Let f be a discrete function. An equation relating dependent


variables fk in the (k , fk ) is called a recurrence relation (RR).

Remark 5 It is common to refer to recurrence relations as difference equations.

We aim to solve recurrence relations with the goal of solving associated differential equations
with series solutions.

Definition 7 (Order of a Recurrence Relation) The order of a recurrence relation is the largest
difference in index values appearing in the equation.

20
 Since recurrence relations are usually written in terms of xk , Definition 7 implies that the
order of a recurrence relation is the difference between the largest and smallest subscripts of the
dependent variable appearing in the equation.

Example 2.3 (Order of a Recurrence Relation) The recurrence relation

xk +1 − 2xk + 2xk −1 = 2k − 1,

has order (k + 1) − (k − 1) = 2.

The general n -th order linear recurrence relation with constant coefficients takes the form

a n xk +n + a n−1 xk +n−1 + a n −2 xk +n−2 + · · · + a 1 xn+1 + a 0 xn = f (k ), (2.3)

where {a i } are independent of xk , k = {0, 1, 2, . . . , n }, and f (k ) is any function of k . The theory of

solving linear recurrence relations such as (2.3) is similar to that for solving linear ODEs. Indeed,
there are equivalent definitions and methods in the theory of recurrence relations to those in the
theory of differential equations.

Definition 8 (Homogeneity of a Recurrence Relation) The recurrence relation is homogeneous

if and only if each term in the relation carries the same number of powers of the unknown function
that defines the relation. A recurrence relation that is not homogeneous is called
non-homogeneous.

Remark 6 (Homogeneity of a Recurrence Relation) We can test if a recurrence relation is

homogeneous by replacing each unknown function xi in the recurrence by a scalled multiple c xi ,
for any fixed constant c . If it is possible to simplify this rescalled recurrence relation to recover the
original recurrence relation then the recurrence relation is homogeneous, otherwise it is
non-homogeneous. A consequence of Definition 8 is that for linear recurrence relations, of the form
(2.3), necessarily f (k ) = 0.

The next section contains some general discussions of common discrete models.

2.3 General Discrete Models

It is a good point to discuss some generalisations of the patterns we have seen thus far. In this
section we consider a varied collection of related models and compare their differences and
similarities. To do this we consider first some general problems in Economics.

Remark 7 (Postulating Solutions) A common feature of the solutions to Example 2.1 and
Example 2.2 is the admission of a solution of the form xi = c λi where c and λ are fixed. Equipped
with this knowledge, we postulate a trial solution for systems of this kind, substitute this trial
solution into the corresponding equation we want to solve and determine c and λ. This is known
as the method of trial solution or method of undetermined coefficients. We shall use this pattern
of solution in the discussions to follow.

21
 In the discussions that follow, we omit the discussion on Verification, Implementation and
Maintenance, in favour of focusing on the assumptions, formulation and solution of each model
in question.

2.3.1 Discrete Population Growth

We now generalise the previous discrete models by considering two problems in Economics where
population growth is modelled using two different population growth patterns. The first pattern
is defined by a rate of population growth that is proportional to the current population size, and
a second is defined by a growth rate proportional to a change in population size. Although the
distinction seems small, these models will give markedly different predictions for the population
size at some time t . We discuss each below.

Example 2.4 (Discrete Population Growth Proportional to Current Population Size)

Consider a population where the change in population number over any time interval is
proportional to the size of the population at the start of that interval. What is the population size
after n intervals?

Problem Identification: We want to find the population function xn for the population after n
time intervals.

Assumptions: Suppose the initial population has size x0 . Let xi ≡ x (t i ) denote the size of the
population at time t i , the i -th time interval.

Variable Classification: The change in population size is given by the difference in population
between the start of one interval and the start of the next interval, so

xi +1 − xi = β xi ,

where β ∈ R+ is the positive constant of proportionality. This is a single interval relation

between the population sizes with constant rate of change.

Solution and Interpretation: It is convenient to rewrite the governing equation for this system as

xi +1 = (1 + β )xi = αxi ,

We solve this model using a brute force induction on i to yield,

x1 = αx0 ,
x2 = α2 x0 ,
x3 = α3 x0 ,
..
.
xn = αn x0 .

22
 So, after n successive intervals, starting at n = 0, the population is

xn = x0 αn .

The change in population size is proportional to the size of the original population. If 0 < α < 1
then the population size decreases for each iteration. If α = 1 then the population size is
constant. If 1 < α then the population grows with each iteration.

10
α=0
α = 0.5
8 α = 0.9
α = 1.0
α = 1.1
6
xn

0
0 1 2 3 4 5 6 7 8 9 10
n
Figure 2.2: Discrete population growth in a model where the rate of growth α is proportional to
the size of the population.

Next we consider a slight variation of Example 2.4.

Example 2.5 (Discrete Population Growth Proportional to Change in Population Size)

Consider a population where the change in population number over any time interval is
proportional to the change in population number over the previous interval. What is the
population size after n time intervals?

Problem Identification: We want to find the population function xn for the population after n
time intervals.

Assumptions: Suppose the initial population has size x0 . Let xi ≡ x (t i ) denote the size of the
population at time t i , the i -th time interval.

Variable Classification: Given this collection of assumptions, the change in population size across
successive time intervals is then

xi +1 − xi = α (xi − xi −1 ) .

23
 where α ∈ R is the constant of proportionality. Incrementing i → i + 1 allows to rewrite this
recurrence relation as
xi +2 − xi +1 = α (xi +1 − xi ) .

or,
xi +2 − (α + 1)xi +1 + αxi = 0.

This equation relates sequence values at three different index values, i , i + 1 and i + 2.

Solution and Interpretation: Suppose xi = c λi , then

xi +1 = λxi and xi +2 = λ2 xi

and substituting these expressions into the above recurrence relation yields


λ2 − (α + 1)λ + α xi = 0,
(λ − α) (λ − 1) xi = 0.

Since xi = 0 yields the trivial solution, the non-trivial solutions correspond to α = λ and λ = 1,
so
xn = c1 + c2 αn .

Note, that at n = 0, x (0) = c1 + c2 , so the initial population is x0 = c1 + c2 . Then the change in

population size at the n-th iteration is equal to the population at the initial iteration plus
some exponential correction. Adding as subtracting factors of c2 , we recover

xn = (c1 + c2 ) + c2 (αn − 1) .

Noting that
x1 − x0
x1 = c1 α + x0 − c1 or c1 = ,
α−1
and
x1 − x0 αx0 − x1
c0 = x0 − =
α α
we rewrite the solution as
x1 − x0 n αx0 − x1
xn = α + .
α−1 α−1
When α = 0 then the population size is c1 and remains so for all n . If 0 < α < 1, then the
population starts out at c1 + c2 and decreases with each iteration, approaching c1 . If 1 < α
then the population grows without bound.

As is demonstrated in Example 2.4 and Example 2.5, small changes in one recurrence relation
can bring about remarkable changes in the behaviour of a system. Next we consider a more
complicated problem.

24
 10
α=0
α = 0.5
8 α = 0.9
α = 1.0
α = 1.1
6
xn

0
0 1 2 3 4 5 6 7 8 9 10
n
Figure 2.3: Discrete population growth in a model where the rate of growth α is proportional to
the change in size of the population.

2.3.2 Discrete Supply and Demand

Example 2.6 (Discrete Supply and Demand) Consider an economy in which manufacturers
control production to meet demand, while retailers control demand to meet the supply. How do the
production and demand change over time given that the rate of production must adjust to meet
demand and that demand also depends on the level of production?

Problem Identification: We want to determine what the production and demand are as
mathematical functions

Assumptions: Suppose at time t the manufacturer production is P (t ) and retailer demand is D (t ).

Assume that at any time t the rate of production of some commodity is proportional to the
difference in production and demand, and the rate of consumption is proportional to the
difference in demand and production. Suppose that at time P (t = 0) = P0 and D (t = 0) = D0 .

Variable Classification: A suitable discrete model for this system is

Pn+1 − Pn = α (Dn − Pn ) and Dn+1 − Dn = β (Pn − Dn ) .

Note the symmetry of each model of this system under transposition of P and D , and α and
β . Next, we solve each model in turn and compare their relative behaviours.

Solution and Interpretation: Consider the discrete model,

Pn+1 − Pn = α(Dn − Pn ) and Dn+1 − Dn = β (Pn − Dn ),

25
which implies,

β Pn = Dn+1 + (β − 1)Dn and αDn = Pn+1 + (α − 1)Pn

and
β Pn +1 = Dn +2 + (β − 1)Dn+1 and αDn+1 = Pn+2 + (α − 1)Pn+1 .
From these equations we note that,

αβ Pn = αDn+1 + α(β − 1)Dn = Pn+2 + (α − 1)Pn+1 + (β − 1) (Pn+1 + (α − 1)Pn ) ,

β αDn = β Pn+1 + β (α − 1)Pn = Dn +2 + (β − 1)Dn+1 + (α − 1) Dn+1 + (β − 1)Dn .

or,

0 = Pn+2 + (α + β − 2)Pn+1 + (1 − α − β )Pn ,

0 = Dn+2 + (α + β − 2)Dn+1 + (1 − α − β )Dn .

Define, γ = α + β , then,

0 = Pn +2 + (γ − 2)Pn+1 + (1 − γ)Pn and 0 = Dn+2 + (γ − 2)Dn+1 + (1 − γ)Dn .

Assuming trial solutions Pn = a λn and Dn = a λn yield, in each case, λ = 1, or λ = 1 − γ, and

subsequent solutions,

P (n ) = a + b (1 − γ)n and D (n) = c + d (1 − γ)n .

Solving for a , b , c and d yields,

αD0 + β P0 αP0 − αD0 β D0 − β P0
a =c = , b= and d=
γ γ γ
to give,
αD0 + β P0 αP0 − αD0
P (n ) = + (1 − γ)n ,
γ γ
αD0 + β P0 β D0 − β P0
D (n ) = + (1 − γ)n .
γ γ
Note, again, the symmetry of each model of this system under transposition of P and D , and
α and β .
For the purpose of comparison, suppose α = β = k , then γ = 2k . The discrete model becomes
1 1
P (n ) = (D0 + P0 ) + (P0 − D0 ) (1 − 2k )n ,
2 2
1 1
D (n ) = (D0 + P0 ) + (D0 − P0 ) (1 − 2k )n .
2 2
The value of k in each exponential term determines the behaviour of each model. The discrete
models have two critical k values, namely 1 − 2k = 0, which implies k = 12 , and that splits the
behaviour of the discrete model according to k < 21 , k = 12 and k > 12 ; and 1 − 2k = 1, which
implies k = 0. These conditions on the k split the relative behaviour of the discrete model into
the following cases

26
 k = 0: The discrete model is constant.

0 < k < 12 : The discrete model decays exponentially to a constant value 21 (D0 + P0 ).

k = 12 : The discrete models are constant.

1
2 < k < 1: The discrete models oscillate with amplitude that decay exponentially about a
constant value.

k = 1: The discrete models oscillate indefinitely with fixed amplitude.

k > 1: The discrete models oscillate indefinitely with increasing amplitude. This model is
unstable.

The relative behaviours of this model at the critical values of k in Figure 2.4 and at the
non-critical values of k in Figure 2.5.

P0

D0 +P0
2

D0

0 1 2 3 4
n
Figure 2.4: The behaviour of P (square) and D (cross) at critical points k = 0 (red), k = 12 (blue)
and k = 1 (green).

2.3.3 Critical Values and Equilibria

Determining if equilibrium values exist, and classifying them as stable or unstable, assists us
immensely in analysing the long-term behaviour of a dynamical system. Again Example 2.1 and
Example 2.2 provide interesting study cases for the behaviour of systems undergoing change.
In particular, changing certain parameter values and initial conditions can result in decreasing,
unchanging or increasing model values. The collection of parametre values that result in constant
output values for each model are particularly interesting and provide a general classification
scheme to help study the kinds of change that can occur in different systems and the specific
conditions under which such changes, or lack thereof, occur.

27
 P0
D0 +P0
2
D0

0 1 2 3 4
n
1
Figure 2.5: The behaviour of P (square) and D (cross) between critical points 0 < k < 2 (red),
1
2 < k < 1 (blue) and 1 < k (green).

Consider a dynamical system of the form

a n +1 = r a n + b , (2.4)

and denote the equilibrium value by a , if one exists. Consider the given system and suppose that
the equilibrium value of this system is a . Then if the system starts at a , it must remain there for
all of n; that is, a n+1 = a n = a for all n. Substituting a for a n+1 and a n in (2.4) yields

a = r a + b. (2.5)

and solving for a in (2.5) we find

b
a= if r ̸= 1.
1−r
If r = 1 and b = 0, then every initial value results in a constant solution. In other words every
initial value would be an equilibrium value. If r = 1 and b ̸= 1, then no equilibrium value exists.
̸ 0, then we get the
If we are given a dynamical system of the form a n +1 = r a n + b where b =
following results,

Value of r Long-term behaviour observed

|r | < 1 Stable equilibrium

|r | > 1 Unstable equilibrium
r =1 No equilibrium

This classification is useful when the long term behaviour of the system it needed, without explicitly
knowing the general solution of the system.

28
2.4 Exercises
Exercise 2.1 Classify the following discrete equations

1. yk +2 − 2yk +1 + 3yk = 2k 2 .

2. 2yk +3 − k yk = 0.

3. yk yk +2 − 2yk3−2 = k 2 .

4. k 2 yk +5 + k !yk −4 = 5.

Exercise 2.2 Solve the following discrete equations

1. yk − 3yk −1 − 4yk −2 = 0.

2. xk +2 − 3xk +1 − 4xk = 0.

3. yk +2 − 2yk +1 − 2yk = 0.

In each case choose some collection of initial data and the determine the exact analytical solution
for each and then compare each against the corresponding numerically determined solution.

Exercise 2.3 Solve the following discrete equations numerically, first using y0 = 1 and then using
y0 = 2. Plot the corresponding solution and compare.

1. yk − yk −1 = 2k .

2. yk − yk −1 = k .

Exercise 2.4 Plot the solutions to the following discrete equations for k = 0, 1, 2, 3, . . . , 20

1. yk = (−1)k .

k
2. yk = − 43 .

k
3. yk = − 43 .
kπ


4. yk = sink 6 .

29
Chapter 3

Continuous Models

As was observed in Section 2.3 the behaviour of a discrete model changes when the interval size
of each step is changed. If we consider the limit where the step size becomes infinitely small, we
recover a continuum of intermediate states. The change in description from one state to the next
is given by an instantaneous rate of change,

∆f ∆f
= lim ,
∆t ∆t →0
∆t →0 ∆t

and the restriction ∆t → 0 corresponds to an infinitesimally small change in the value of the
parameter t for which there is a corresponding change ∆f . This chapter discusses some ideas
and examples of descriptions of infinitesimal change.

3.1 Modelling Continuous Change

Examples in Section 2.1 have an interesting behaviour when the number of steps is taken to be
infinitely large, while the interval that is stepped over in the model remains finite. The first such
interesting behaviour to note is the appearance of a maximum instantaneous rate of continuous
change. This is discussed next.

3.1.1 Euler’s Number

Euler’s number, e , is a constant that appears in many mathematical models and is used as the
base of the natural logarithm. The number can be defined as follows.
 ‹
1 n
e = lim 1 + .
n→∞ n
An interesting example of where this number appears is in compound interest. Consider the
scenario in which you wish to save your money in a bank which offers an interest rate of 100% per
annum. After one year of saving with a particular bank, the amount would double. That is, if we
had saved R1000 initially, we would end the year with R2000. Now, consider the case where half
way through the year we were paid 50% and were offered the chance to reinvest that amount with

31
the initial amount for the rest of the year at a rate of 50%. That is, if we initially saved R1000, half
way through the year we would receive R500, reinvesting it and earning a further 50%. We would
then earn R750, leaving us at a grand total of R2250 for the year. in the first scenario we increased
the amount by a multiplying factor of 2 or

(1 + 1)1 .

In the second scenario we increased the amount by a multiplying factor of 2.25 or

 ‹
1 2
1+ .
2

If we continue to divide the year and the interest rate into their relative proportions we notice
that the amount earned increases, but by a smaller amount each time. For 3 parts of the year we
get 2.3703 . . ., for 4 parts we get 2.4414 . . . and so on. If we were to break the year and the interest
rate up into infinitely many parts, we would find the multiplying factor would be equal to Euler’s
number.  ‹n
1
e = lim 1 + .
n→∞ n
This can also be described by the following equation:
1
e = lim (1 + n ) n . (3.1)
n→0

Both of these have an approximate solution of 2.7182818284 . . ..

In the next section we discuss Euler’s Number and its relationship with continuous models.

3.2 Formal Continuous Models

This section describes the technical aspect of describing continuous change and how it relates
to problems in mathematical modelling. General methods of solutions are deferred for a more
advanced course. Here we shall consider some basic definitions and ideas about differential
equations that will be useful later in this chapter.

3.2.1 Differential Equations

Before we formally construct any continuous models, it will be useful to discuss some definitions.

Definition 9 (Function) A function f taking a parameter t is a map that matches each parameter
value t to a new value f (t ). When f maps real numbers R into R then it is common to write f (t )
as an explicit mathematical expression in terms of t .

The explicit expression for f in terms of its parameter t allows us to analyse how f (t ) changes
with a given change of t , and gives rise to the following definition.

32
Definition 10 (Differentiation) The derivative of a function f (t ) with respect to the parameter t
is  ‹
df f (t + ε) − f (t )
= lim . (3.2)
dt ε→0 ε
The procedure of computing the derivative of f (t ) with respect to t is called differentiation.

The derivative of a function f with respect to t measures the rate-of-increase of f with respect
to an increase in t . We consider the difference in the function value f at points t + ε and t for
some positive value ε for increasing t . In this way it is more natural to state that differentiation
measures rate-of-increase rather than rate-of-change. The limit ε → 0 in (3.2) means that this
rate-of-increase is measured at the point t , rather than at the end points of some interval as is
commonly used when measuring a gradient. Therefore, the derivative of f with respect to t is the
instantaneous rate-of-increase of f at the point t . A combination of derivatives of functions in a
single expression leads to the following definition.

Definition 11 (Differential Equation) A differential equation is a mathematical equation that

relates a function with its derivatives.

There are many ways to relate a function with its derivatives. When this is done, we find the
following classification useful. If x is a function only of t , then a linear, n-th order, ordinary
differential equation (DE) of x and t with constant coefficients takes the form
 ‹  n−1   ‹
dn x d x dx
cn + cn−1 + · · · + c1 + c0 x = f (t ), (3.3)
dt n d t n−1 dt

where {ci } are independent of t , i = {0, 1, 2, . . . , n }, and cn ̸= 0.

Definition 12 (Order of a Differential Equation) The order of a differential equation is the

highest order of differentiation appearing in the equation.

dn
The order of (3.3) is n . This corresponds to the term with the highest derivative operator dt n .

Definition 13 (Linearity of a Differential Equation) The linearity of a differential equation is

determined by the presence or absence of non-linear terms. If any term in the DE is of the form


dm x k m
k

dt m ̸ 1, or x j dd t mx for any non-zero j or k then the DE is non-linear, otherwise

or x k where k =
it is linear.

Equation (3.3) is linear in x and the derivatives of x , therefore it is linear.

Remark 8 It is common to introduce a simplified notation

 ‹
d df
f ′ (t ) ≡ f (t ) =
dt dt

33
 df


for the derivative dt when f is a function of a single independent variable t . This is referred to as
prime notation. Higher order derivatives are signified by adding addtional primes
 2   3 
′′ d2 d f (t ) ′′′ d3 d f (t )
f (t ) ≡ f (t ) = and f (t ) ≡ f (t ) =
dt 2 dt 2 dt 3 dt 3

and so on. Clearly this notation becomes less convenient for more than four derivatives.

Differential equations appear naturally in mathematics and physics as the determining

equations governing the behaviour of physical systems. The property of a differential equation
that makes it different from an algebraic equation is that the solutions to algebraic equations are
numbers whereas the solutions to differential equations are themselves equations. This means
that the solution to a formal continuous model is not a number that might describe what a given
system is doing at any one point in time or place, but rather an equation that describes what the
system does. The particular behaviour of the system at a given point in time and space can be
determined by evaluating the solution to the differential equation at that point of interest.

Remark 9 We shall consider only linear, first order differential equations in a single parameter in
these notes.

In the next section we shall familiarise ourselves with some interesting aspects of the behaviour
of an important function and its derivatives that will appear in the models later in this chapter.

3.2.2 Euler’s Number and Differential Equations

Now let us consider the function y = e x . This function has many interesting characteristics. For
example, consider its derivative,

d x e x +h − e x
e = lim .
dx h →0 h
If we simply factor out the e x and note it is not affected by the limit, we get the following,

d x eh −1
e = e x lim .
dx h →0 h

With a simple substitution of n = e h − 1 and some algebraic manipulation we can show that this
is equivalent to
d x 1
e =ex € 1
Š.
dx lim ln (1 + n) n
n→0

At this point we note that

1
e = lim (1 + n ) n .
n→0

Substituting (3.1) into this we find a useful result.

d x
e = e x. (3.4)
dx

34
Using this knowledge we can consider evaluating the derivative of

y = ln (x ) . (3.5)

After a simple manipulation of (3.5) we get

e y = x.

Applying implicit differentiation to this we get

 ‹
y dy
e = 1.
dx
Rearranging we find  ‹
dy 1
= y.
dx e
By substituting (3.5) into this we get  ‹
dy 1
= ln(x ) ,
dx e
which simplifies to our next useful result.
d 1
ln (x ) = .
dx x
There is a special relationship between the instantaneous rate of change of a function and the
function itself. This is a general pattern that we observe when studying rates of change. In addition,
the number e shall play an important role in the solution to many models.

3.3 General Continuous Models

Following the development of the discrete and continuous model, we present an instructive
set of continuous models which explicate some generic patterns that appear in models of and
solutions to general continuous models. To do this, we shall consider first some general problems
in Economics and Physics.

3.3.1 Continuous Population Growth

The Malthusian Model is a simple population growth model. This model neglects some important
factors that are relevant to population growth, which limit its applicability, however it has some
interesting predictions and proposed consequences that are worth discussing.

Example 3.1 (Continuous Population Growth) Consider the problem of an initial population
which at some initial time is P0 . What might the population be at some later time t ?

Problem Identification: Suppose the population is P0 at the time t = t 0 , and we are interested in
predicting the population P at some future time t = t 1 . We want to find a population function
P (t ) for t 0 ≤ t ≤ t 1 satisfying P (t 0 ) = P0 .

35
Assumptions: Consider some factors that pertain to population growth. Two obvious ones are
the birth rate and the death rate. The birth rate and death rate are determined by different
factors. The birth rate is influenced by infant mortality rate, attitudes toward and availability
of contraceptives, attitudes toward abortion, health care during pregnancy, and so forth.
The death rate is affected by sanitation and public health, wars, pollution, medicines, diet,
psychological stress and anxiety, and so forth. Other factors that influence population growth
in a given region are immigration and emigration, living space restrictions, availability of food
and water, and epidemics. For our model, let’s neglect all these latter factors. We can revise
our model at a later stage if we are dissatisfied with the results. For now we’ll consider only
the birth rate and death rate. Because knowledge and technology have helped humankind
diminish the death rate below the birt hrate, human populations have tended to grow.
Let’s begin by assuming that during a small unit time period, a percentage b of the population
is newly born. Similarly, a percentage c of the population dies. In other words, the new
population P (t + ∆t ) is the old population P (t ) plus the number of births minus the number
of deaths during the time period ∆t . Symbolically,

P (t + ∆t ) = P (t ) + b P (t )∆t − c P (t )∆t ,

or
∆P
= b P − c P = (b − c )P = k P.
∆t
From our assumptions the average rate of change of the population over an interval of time
is proportional to the size of the population. Using the instantaneous rate of change to
approximate the average rate of change, we have the following differential equation model,
 ‹
dP
= k P, P (t 0 ) = P0 and t0 ≤ t ≤ t1, (3.6)
dt

where (for growth) k is a positive constant.

Solution and Interpretation: We can separate the variables and rewrite (3.6) by moving all terms
involving P and d P to one side of the equation and all terms in t and d t to the other. This
gives,
dP
= k dt .
P
Integration of both sides of this last equation yields

ln (P ) = k t + C , (3.7)

for some constant C . Applying the condition P (t 0 ) = P0 to (3.7) to find C results in

ln (P0 ) = k t 0 + C ,

or
C = ln (P0 ) − k t 0 .

36
 Then, substitution for C into (3.7) gives

ln (P ) = k t + ln (P0 ) − k t 0 .

or, simplifying algebraically,  ‹

P
ln = k (t − t 0 ).
P0
Finally, by exponentiating both sides of the preceding equation and multiplying the result by
P0 , we obtain the solution
P (t ) = P0 exp [k (t − t 0 )] (3.8)


Verifying the model: Note that ln PP0 = k (t − t 0 ), this model predicts a linear relationship between


ln PP0 and t − t 0 . This means that a graph showing the relationship between t − t 0 on the


horizontal axis and ln PP0 on the vertical axis would be given by a straight line, passing
through the origin with slope k . However, a plot the population data for an industrially
developed country over a time interval of several years is not well described by this model. For
example, the 1990 census for the population of the United States was 248 710 000.000 and in
1970 it was 203 211 926.000. Substituting these values into (3.8) and dividing the first result
by the second gives
248710000
= e k (1990−1970) .
203211926
Thus,  ‹
1 248710000
k= ln ≈ 0.01.
20 203211926
That is, during the 20-year period from 1970 to 1990, population in the United States increased
at an average rate of 1% per year. Using this information with (3.8) to predict the population
for the year 2000, where t 0 = 1990, P0 = 248710000 and k = 0.01 yields

P (2000) = 248710000e 0.01(2000−1990) = 303775080.

The 2000 census for the population of the United States was 281 400 000.000 (rounded to the
nearest thousand). This prediction is off by approximately 8%. An error of this magnitude
might not be worth correcting when only short time scales are considered, However a simple
calculation shows that this model predicts a population of 55 209.000 billion in the year
2300. This prediction is clearly wrong since it exceeds the sustainable carrying capacity of the
entire planet. Therefore, we conclude that this model is unreasonable over long time intervals.
Clearly, this model requires some refinement to make more sensible predictions. This can be
done by taking into account the factors that might reduce the population growth rate over
long time-intervals. These limiting factors might include limitations on food production and
other resources.

Refine the Model to Reflect Limited Growth: As a first attempt it is reasonable to incorporate some
limiting factors to the Malthusian population model to reflect the limitation on population
growth. Suppose that the proportionality factor k , which measures the rate of population

37
growth in (3.7), is not longer constant but a function of the population. As the population
increases and gets closer to the maximum population M , the rate k decreases. One simple
sub model for k is the linear one

k = r (M − P ) for r >0

where r is a constant. Substitution into (3.7) leads to

 ‹
dP
= r (M − P )P (3.9)
dt
or
dP
= r dt . (3.10)
P (M − P )
The model described in (3.9) is known as the logistic growth model. Again we will assume the
initial condition P (t 0 ) = P0 . It follows from elementary algebra that
 ‹
1 1 1 1
= + .
P (M − P ) M P M −P

Thus, (3.10) can be rewritten as

dP dP
+ = rM t +C
P M − P0

which integrates to
ln (P ) − ln (|M − P |) = r M t + C (3.11)

for some arbitrary constant C . Using the initial condition, we evaluate C in the case P < M ,
 ‹
P0
C = ln − r M t0.
M − P0

Substituting into (3.11) and simplifying gives

 ‹  ‹
P P0
ln − ln = r M (t − t 0 )
M −P M − P0
or  ‹
P (M − P0 )
ln = r M (t − t 0 ).
P0 (M − P )
Exponentiating both sides of this equation gives

P (M − P0 )
= e r M (t −t 0 )
P0 (M − P )
or
P0 (M − P )e r M (t −t 0 ) = P0 (M − P0 ).

Then,
P0 M e r M (t −t 0 ) = P0 (M − P0 ) + P0 P e r M (t −t 0 ) ,

38
 so solving for the population P gives

P0 M e r M (t −t 0 )
P (t ) = .
M − P0 + P0 e r M (t −t 0 )

To estimate P as t → ∞, we rewrite this last equation as

M P0
P (t ) = .
P0 + (M − P0 )e −r M (t −t 0 )

We can simplify this as follows. Define

 ‹
M
k = 1− ,
P0

then
M
P (t ) = . (3.12)
1 − k e −r M (t −t 0 )
This can be rewritten as either
P0
P (t ) = P0 P0

.
M + 1− M e −r M (t −t 0 )

or
M
P (t ) =
.
1− 1− M
P0 e
−r M (t −t 0 )

Notice from (3.12) that P (t ) approaches M as t tends to infinity. Moreover, from (3.9) we
calculate the second derivative

P ′′ = r M P ′ − 2r P P ′ = r P ′ (M − 2P )

so that P ′′ = 0 when P = M2 . This means that when the population P reaches half the limiting


population M , the growth ddPt is the most rapid and then starts to diminish toward zero.
M
One advantage of recognising that the maximum rate of growth occurs at P = 2 satisfied
that the growth involved is essentially logistic, if the point of maximum rate of growth has
M
been reached, then 2 can be estimated.

3.3.2 Continuous Supply and Demand

Example 3.2 (Continuous Supply and Demand) Consider an economy in which manufacturers
control production to meet demand, while retailers control demand to meet the supply. What might
be the supply and demand on retailers at some time t ? Note that in this example, we shall focus on
understanding the behaviour of the model and defer validation and refinement of the model until
later.

Problem Identification: Suppose at time t the manufacturer production is P and retailer demand
is D . We want to find the P (t ) and D (t ), the supply and demand functions, respectively.

39
 r = 0.10
r = 0.25
r = 0.50
r = 1.0
r =2
P0 + M

P0

P0 − M

0 1 2 3 4 5 6 7 8 9 10
t
Figure 3.1: The Malthusian limited continuous population growth in a model for various growth
factors r .

Assumptions: Since the retailers want to serve the demand by the customers, and their ability to
do so is limited by their ability to produce goods, we should expect that there is some relation
between the production P and the demand D . Clearly, it is not economical to produce more
than can be sold. Assume that at any time t the rate of production of some commodity is
proportional to the difference in production and demand, and the rate of consumption is
proportional to the difference in demand and production. Suppose that at time P (t = 0) = P0
and D (t = 0) = D0 . A suitable continuous model for this coupled system is
 ‹  ‹
dP dD
= α(D − P ) and = β (P − D ).
dt dt
Note the symmetry of this system under transposition of P and D , and α and β .

Solution and Interpretation: Consider first the continuous model, more conveniently rewritten as

P ′ = α(D − P ) and D ′ = β (P − D ),

which implies,
βP = D ′ + βD and αD = P ′ + αP,

and
β P ′ = D ′′ + β D ′ and αD ′ = P ′′ + αP ′ .

From these equations we note that,

αβ P = αD ′ + αβ D = (P ′′ + αP ′ ) + β (P ′ + αP ),
αβ D = β P ′ + αβ P = (D ′′ + β D ′ ) + α(D ′ + β D ).

40
 or,
P ′′ + (α + β )P ′ = 0 and D ′′ + (α + β )D ′ = 0.

Define, γ = α + β , then,

P ′′ + γP ′ = 0 and D ′′ + γD ′ = 0.

Since γ is a constant, we have

d
d

P ′′ + γP ′ = P ′ + γP = 0 and D ′′ + γD ′ = D ′ + γD = 0
dt dt
so P ′ + γP and D ′ + γD do not depend on t since their derivatives with respect to t are zero.
This implies that

P ′ + γP = cp and D ′ + γD = cd

for some constants cp and cd . It is now clear that P and D are determined by first order linear
differential equations, so our previous methods can be used again to solve these differential
equations. Since we do not know the values of cp and cd we will revert to solving the second
order equations using first order techniques. Assuming trial solutions P (t ) = a e λt and D (t ) =
a e λt yield, in each case, λ = 0, or λ = −γ, and subsequent solutions,

P (t ) = a + b e −γt and D (t ) = c + d e −γt ,

and so,

P ′ = −γb e −γt = α(c − a ) + α(d − b )e −γt ,

D ′ = −γd e −γt = β (a − c ) + β (b − d )e −γt .

Solving for a , b , c and d yields,

αD0 + β P0 αP0 − αD0 β D0 − β P0

a =c = , b= and d= ,
γ γ γ

to give,

αD0 + β P0 αP0 − αD0 −γt

P (t ) = + e ,
γ γ
αD0 + β P0 β D0 − β P0 −γt
D (t ) = + e .
γ γ

Note, again, the symmetry of each model of this system under transposition of P and D , and
α and β .

Discussion and Analysis: For the purpose of comparison, suppose α = β = k , then γ = 2k . Then
the continuous model becomes
1 1
P (t ) = (D0 + P0 ) + (P0 − D0 ) e −2k t ,
2 2

41
 1 1
D (t ) = (D0 + P0 ) + (D0 − P0 ) e −2k t .
2 2
The value of k in each exponential term determines the behaviour of each model. The
continuous models have a single critical k value, namely k = 0 that splits the behaviour of
the continuous model according to k < 0, k = 0 and k > 0. These conditions on the k split the
relative behaviour of the discrete model and continuous model into the following five cases.

k = 0: The continuous models give solutions that are identical to the discrete models.
0 < k < 12 : The continuous and discrete models are qualitatively the same, and decay
exponentially to a constant value 12 (D0 + P0 ), although the actual values given by each
model may differ.
k = 21 : The discrete models are constant, while the continuous model solutions tend to the
discrete model solutions after sufficiently large parameter values t .
1
2 < k < 1: The discrete models oscillate with amplitude that decay exponentially about a
constant value where as the continuous models are purely decaying and tend to the same
constant value.
k = 1: The discrete models oscillate indefinitely with fixed amplitude, while the continuous
models decay exponentially and tend to a constant value.
k > 1: The discrete models oscillate indefinitely with increasing amplitude, while the
continuous models decay exponentially and tend to a constant value. In this case the
discrete model is unstable.

The relative behaviours of this model at the critical values of k in Figure 3.2 and at the
non-critical values of k in Figure 3.3.

In the next section, we consider another important model of continuous change from Physics.

3.3.3 Newton’s Model of Cooling

Some physical systems have the interesting behaviour where by the rate of change of some part
of that system is proportional to the difference between some initial value and the current value.
An example of such a system is the familiar problem of the rate at which something cools by
losing heat to its environment. This feature of heat moving from a hot body to a cooler one can be
seen in the melting of ice, the radiance of the sun and, for a problem that is closer to home, the
cooling of tea at teatime. The following example discusses a solution to this problem that was
first developed by Isaac Newton.

Example 3.3 (Cooling a Cup of Tea) Consider the situation of a cup of hot tea on a table in a room
at room temperature. How would we model the cooling of the tea? The rate of change of the
temperature of an object is proportional to that of the temperature of the medium surrounding
that object. Two varying values are said to be proportional to one another if they are related by a

42
 P0

D0 +P0
2

D0

0 1 2 3 4
t
1
Figure 3.2: The behaviour of P (solid) and D (dotted) at critical points k = 0 (triangle), k = 2
(square) and k = 1 (diamond).

P0

D0 +P0
2

D0

0 1 2 3 4
t
Figure 3.3: The behaviour of P (solid) and D (dotted) between critical points 0 < k < 12 (triangle),
1
2 < k < 1 (square) and 1 < k (diamond).

multiplicative constant, this is known as the constant of proportionality. That is there is a constant
ratio between the two values. Keeping these concepts in mind, we can step through the modelling
process.

Problem Identification: Suppose the temperature of the tea was H at time t and the room
temperature was T . We want to find the temperature function at time t , H (t ).

43
Assumptions: We will assume that the temperature of the tea H is a differentiable function of t
and that the temperature of the room is constant, since the surrounding air is of such a large
volume that the effect of the temperature of the tea will be negligible. We will also assume
that the rate of change of the temperature of the tea is proportional to the difference in the
temperature of the tea and the temperature of the room. Since we want the temperature to go
down when the tea is warmer than its surrounding environment we will assign a negative
value to the coefficient. That is
 ‹
dH
= −k (H − T ). (3.13)
dt

Interpretation: When T < H , the gradient is negative, and when H < T the gradient is positive.
These both make sense since we would expect the tea to warm up if the room temperature is
higher than its temperature and to cool down when the room temperature is lower than its
temperature. When H = T we have reached an equilibrium point since the rate of change is 0.

Solution: We can solve this model by separating variables and integrating (3.13) to yield

H (t ) = C e −k t + T, (3.14)

where C is a constant of integration that is determined by some collection of conditions

specified by the initial state of the system. Several examples of (3.14) are plotted in Figure 3.4.

T +C
k = 0.1
k = 0.25
k = 0.5
k =1
k =2
H (t )

0
0 1 2 3 4 5 6 7 8 9 10
t
Figure 3.4: The cooling curves for a collection of systems for differing values of cooling coefficient
k.

44
3.4 Exercises
Exercise 3.1 Solve the following continuous equations analytically and the verify the solution by
comparing against the numerically solved solution on the interval x ∈ (1, 100]
dy


1. dx − 2y 2 = 0.
dy


2. dx − 2y = 0.
dy


3. dx + 5 = 10.
dy


4. dx + x5 y = 0.

In each case suppose y (1) = 1.

dx


Exercise 3.2 Suppose that the function x (t ) satisfies the equation dt = t + t x , and that x (0) = 2,
find the value of x (1) rounded to two decimal places. Determine the solution analytically and verify
the solution by comparing with the numerically determined solution.

Exercise 3.3 If the growth rate of bacteria in a culture medium is directly proportional to the
difference between the number of bacteria and a certain number N . Derive an expression for the
number of bacteria n as a function of the time t given that the initial number of bacteria placed in
the culture medium is N0 and the constant of proportionality is λ.

Exercise 3.4 Suppose that the growth rate of the South African population is proportional to the
number of people living at any given time. If in 1970 it was estimated that there were 15 million
people living in South Africa, and in 1980 it was estimated that this number had increased to 20
million, estimate the population of South Africa in 2019. Generate a plot of the population versus
time, starting in the year 1970 until 2025.

45
Appendices

47
Appendix A

Units and Dimensions

Table A.1: A collection of physical units and their associated dimension measured in units of Mass
(M ), Length (L ) and Time (T ).

Unit Symbol Dimension

Mass M M
Length L L
Time T T
Linear Velocity v LT −1
Linear Momentum p M LT −1
Linear Acceleration a LT −2
Angular Velocity ω T −1
Angular Momentum L M L 2 T −1
Angular Acceleration α T −2
Density ρ M L −3
Force F M LT −2
Power P M L 2 T −3
Energy E M L 2 T −2
Work W M L 2 T −2
Pressure P M L −1 T −2
Entropy S M L 2 T −2
Frequency f T −1
Heat Q M L 2 T −2
Rotational Inertia I M L2
Specific Weight γ M L −2 T −2
Surface Tension γ M T −2
Torque τ M L 2 T −2
Dynamic Viscosity µ M L −1 T −1
Kinematic Viscosity ν M L 2 T −1

49
Appendix B

Software

MATLAB is a powerful software package that is widely used in Science and Engineering. Here we
shall briefly discuss some instructive examples on its usage. The following note and discussion
on MATLAB are necessarily incomplete and narrowly focussed only on those aspects that are used
in this course. Additional information on MATLAB and the material presented in this appendix
can be found online at [Link] and in the MATLAB
documentation.

Remark 10 The purpose of this appendix is to present a minimal introduction to some useful
aspects of the MATLAB software package. It is necessarily an incomplete introduction to MATLAB. It is
neither a numerical analysis course nor general software development instructions set. The purpose
of this appendix is to present a small collection of useful pieces of MATLAB code as can be used in
computations related to mathematical modelling.

B.1 Getting Started

Below is a minimal collection of MATLAB code examples to get started using MATLAB to numerically
solve mathematical modelling problems and plot their solutions. Example B.1 lists examples of the
basic arithmetic operations provided in MATLAB. These operations can be applied to any primitive
datatype in MATLAB as part of a regular script file (.m), live script (.mlx) or on the command line
interface (CLI).

Remark 11 Each line of code is ended with a semi-colon to prevent the printing of each line of
code to the output. Each line of code is numbered with a corresponding description of the meaning
below. In each case, the reported outputs can be verified by executing the corresponding code.

Remark 12 Comments are preceded by any number of % characters and are included as notes to the
reader to help understand the code. The MATLAB interpreter ignores all comments during evaluation,
however the live script environment in uses %% to distinguish contiguous blocks of code. In the
example below, the numerical value of some assignments are included in the comments and can be
checked by evaluating the code explicitly.

51
 Example B.1 (Basic Operations) Basic arithmetic operations and variable assignment are
presented in the code listing below.

1 a = 3;
2 b = 2;
3
4 s = a + b; % 5
5 d = a - b; % 1
6 p = a * b; % 6
7 q = a / b; % 1.5

The variable a is assigned the value 3 in line 1, while the variable b is assigned the value 2 in
line 2. Then, s is assigned the value of the sum a + b in line 4, d is assigned the difference a - b in
line 5, while p is assigned the value the product of a * b in line 6 and q is assigned the value of the
quotient a / b in line 7. The values of a, b, s d, p and q can be used in any subsequent computation.
Evaluating the code in this listing will introduce the variables a, b, s d, p and q in to computation
session. The value of each of these variables can be retrieved and printed to the standard output by
evaluating the corresponding variable on the CLI or by omitting the semi-colon in the code.

A common task in any computing problem involves repeating one or more sequences of
mathematical operations. This introduces a flow control problem wherein some number of tasks
are completed with respect to some counting parameter that both tracks the sequence of steps
and may also enter as a numerically evaluated parameter in one or more of the tasks. An example
of such a task is presented in Example B.2 using the for loop control structure.

Example B.2 (For Loop) The following code computes the sum of the first ten natural numbers,
namely 1 + 2 + 3 + · · · + 9 + 10 by iteratively adding the next natural number to a total called s. To
perform this summation, the value of s before anything is added is set to zero, where after a loop
counter iterates over the first ten natural numbers, each time adding the corresponding natural
number to the running total. The numerical check and increment of the loop counter variable are
implicit.

1 s = 0; % initialize the sum counter

2
3 s % 0
4
5 for n = 1:10 % n is the loop counter
6 s = s + n;
7 end
8
9 s % 55

The variable s will hold the value of the sum and is initialized to a value of zero in line 1 and
the initial value of s is printed to the standard output on line 3. The loop control structure that

52
 establishes the iteration scheme begins on line 5 with the keyword for followed by the name of
the iteration variable n that is initialized to 1 on the first iteration of the loop and is incremented
by one on each subsequent iteration until it reaches a value of 10 and the loop terminates. On
each iteration of the loop, all the code contained between the keyword for and the keyword end is
executed; in this case, the code between line 5 and line 7 is executed each time the loop iterates. On
each iteration of this loop, the value of s is reassigned to a value equal to the present value of s + n
that is evaluated before assigning a new value to s. After the loop terminates, the new value of s is
printed on line 9.

Example B.3 gives an example implementation of the code in Example B.2 using the while
loop flow control structure.

Example B.3 (While Loop) The following code computes the same quantity as that computed
in Example B.2, this time using an explicit numerical check of the loop iteration counter and
incrementing the loop counter variable.

1 s = 0; % initialize the sum counter

2 n = 1; % initialize the loop counter
3
4 s % 0
5
6 while ( n < 10) % test counter
7 s = s + n;
8 n = n + 1; % increment counter
9 end
10
11 s % 55

The variable s will hold the value of the sum is initialized to zero in line 1 and the loop counter
n is initialized to 1 on line 2. The value of s is printed on line 4. The loop control structure that
establishes the iteration scheme begins on line 6 with the keywordwhile followed by the test statement
that returns true when the loop counter is less that 10 and false otherwisr. When the test statement
evaluates to true, then all the code between the keyword while and the keyword end is executed;
in this case, the code between line 6 and line 9 is executed each time the loop iterates. On each
iteration of this loop, the value of s is reassigned to a value equal to the present value of s + n
that is evaluated before assigning a new value to s on line 7 and then the value of the loop counter
variable is updateed on line 8. After the loop terminates, the new value of s is printed on line 11.

MATLAB includes a collection of standard functions that will be useful during the numerical
modelling process. Example B.4 gives some examples of these functions.

Example B.4 (Built-In Functions) Below is a short collection of built-in functions evaluated at
particular values of their input parameters and the corresponding numerical values.

53
 1 % % Standard Constants
2 pi % 3.1416
3
4 % % Powers and exponentials
5 e = exp (1) % 2.7183
6 y = log ( e ) % 1
7 z = e ^( pi ) % 23.1407
8 w = sqrt ( y ^2 + z ^2) % 23.3205
9
10 %% Standard Trigonometric functions
11 x1 = sin ( pi / 3) % 0.8660
12 x2 = cos ( pi / 3) % 0.5000
13 x3 = tan ( pi / 3) % 1.7321
14
15 %% Standard Inverse Trigonometric functions
16 y1 = asin ( x1 ) % 1.0472
17 y2 = acos ( x2 ) % 1.0472
18 y3 = atan ( x3 ) % 1.0472

The numerical value of π is evaluated and printed to standard output on line 1. Then tne
numerical value of Euler’s number e is assigned to the variable e by first evaluating the exponential
function exp at 1 on line 5, then the value of the logarithm of e is assigned to y on line 6. On line
7, the numerical value of the exponential of π is assigned to the variable z. On line 8, the built-in
function sqrt is used to compute and assign the value of the square root of the sum of the square of
y and z to w.
The numerical value of sin (π/3), cos (π/3) and tan (π/3) are evaluated and assigned to variables
x1, x2 and x3 on lines 11, 12 and 13, respectively. Finally, the trigonometric inverse functions are
evaluated on each of x1, x2 and x3 and assigned to y1, y2 and y3 on lines 16, 17 and 18, respectively.
The consistency of these assignments can be checked by explict evaluation of this code.

MATLAB is a software package that was originally designed to perform matrix computations.
One of the goals in designing MATLAB was to make setting-up and performing computations
involving matrices (of arbitrary dimensions) simple for the user. As such, MATLAB has several
useful functions that simplify the task of constructing matrices and accessing individual matrix
elements. Example B.5 contains example code to define and initialize matrices.

Example B.5 (Arrays) The construction of arrays is acheived using the following code.
X is defined as an empty array on line 1. X can be programmatically updated to modify the shape
and structure as needed by indexing into the array. MATLAB will reallocate memory to X as needed,
however, repeated reallocations may result in degraded performance. It is often more economical
to pre-allocate arrays as sequences of known size and shape and initialize them to some sensible
initial values. This is done on line 2 where the an array comprising ten rows and a single column of
zeros is assigned to the variable l̀stinlineZ. The values of both Z and X are updated in lines 5 and 6,
respectively. Notice the use of square brackets to denote and array in line 1 while individual entries

54
1 X = []; % define an uninitialized array of unspecified dimensions
2 Z = zeros (10 , 1) ; % define a column matrix initialized to all
zeros
3
4 for n = 1:10
5 Z(n) = n;
6 X(n) = 2 * Z(n);
7 end

in each array are accessed using parentheses, in the same way that parameters are passed to uild-in
functions. Arrays and matrices can be thought of as ’built-in’ functions in this way.

From time to time is will be necessary to clear the variable or symbol during a MATLAB session.
Example B.6 gives some example commands to clear the environment during a session.

Example B.6 (Clear Environment) The following code can be evaluated to clear part or all of the
computational environment during a running session. Comments are include after each line of
code for quick and easy reference.

1 clear all % remove all definitions from current session

2 close all % close unnecessary open windows
3 clc % clear the console window

All user defined values and symbols are cleared from the current session in line 1. Graphical
output windows are closed in line 2. The command line interface is cleared in line 3. Any
combination of these commands can be executed during a running MATLAB session.

B.2 Functions
Real world problems usually have mathematical models comprising numerous intricate
mathematical expressions. Moreover, the evaluation of these mathematical models often involve
a collection of carefully organized mathematical steps and evaluation of lengthy mathematical
expressions. It is often simpler to organize the numerical evaluation of such models into a
collection of well defined and understood functions that can be evaluated during the
computation process, storing previously computed quantities for later use and accumulating the
results of a collection of intermediate steps to generate the final result. The process of converting
the mathematical model in to a sequence of computational steps that can be communicated to
others and executed inside computers is an importand part of the modeling process. Since it is
common to repeat, reuse and combine mathematical expressions, it is instructive to consider a
few such cases explicitly. Example B.7 contains some of the basics on writing custom functions to
be used in the numerical evaluation of mathemaical models. Each custom function definition
starts with the keyword function followed by the return value list, written as a comma separated

55
 list of variable names inside square brackets, and set equal to the name of the function and the
comma separated parameter list in parenthesis. This function definition is followed by the
function body that contains the implementation detail of the function. The end of the function
implementastion is marked by the end keyword. Once the custom function is specified, it can be
used like any of the built-in functions.
Now consider the following functions


y1 (x ) = 2 sin πx 2 − 3 + x


y2 (x ) = −3 sin πx 2 + x − 1

where x is the name of the independent parameter in each case. Example B.7 and Example B.8
discuss two possible implementaitons of the functions y1 (x ) and y2 (x ), evaluated at a finite number
of evenly spaced points x on a given interval.

Example B.7 (Functions) Consider the following functions

y1 (x ) = 2 sin πx 2 − 3 + x


y2 (x ) = −3 sin πx 2 + x − 1

where x is the name of the independent parameter in each case. The code listing below gives example
implementations of these two functions.

1 function [ y ] = y1 ( x )
2 N = length ( x ) ;
3 y = zeros (N , 1) ;
4
5 for n = 1: N
6 y ( n ) = 2 * sin ( pi * x ( n ) * x ( n ) - 3) + x ( n ) ;
7 end
8 end
9
10 function [ y ] = y2 ( x )
11 N = length ( x ) ;
12 y = zeros (N , 1) ;
13
14 for n = 1: N
15 y ( n ) = -3 * sin ( x ( n ) * x ( n ) * pi ) + x ( n ) - 1;
16 end
17 end

The function y1 defined on line 1 contains the MATLAB implementation of the function y1 defined
above. The function y1, defined on line 1, returns an array of numerical values y1 (x ) for input value
x provided in the array x of length N. The first step is to determind the number of input parameters
in the array x in line 2 and pre-allocate an array with as many values for y1 as there are values x in
line 3. Then by iterating over each element in the input array x we assign the corresponding value
of the return value y, using a for loop to access individual entries in x and y. Similarly, the function
y2 defined on line 10 contains the implementation of the y2 function above.

56
 Example B.8 presents and alternate implementation of those functions in Example B.8 using
vectorization.

Example B.8 (Functions Alternate Implementation) When the points where a function is to be
evaluated is given in the form of a list or array, the evaluation process can be simplified by instructing
MATLAB to vectorize the evaluation. This is done by replacing the direct indexing of elements in the
input and output arrays with an implicit evaluation. Addition an subtraction are automatically
vectorized for arrays of the correct size, however multiplication is a more complicated operation
with requires the use of a special multiplication operation .* in place of the ususal * operator. In
the listing below presents this simplified, vectorized operation, which is both simpler than that
presented in Example B.7 and is also more like the mathematical definitions presented above.

1 function [ r ] = y1 ( x )
2 r = 2 * sin ( pi * x .* x - 3) + x ;
3 end
4
5 function [ r ] = y2 ( x )
6 r = -3 * sin ( x .* x * pi ) + x - 1;
7 end

In each of the alternate implementations for y1 and y2 presented above, the return value for
each expression is more closely related to that in the original function definitions.

Example B.9 gives and example of how multiple functions can be nested inside other functions,
where precomputed values are reused and return values are concatenated to generate new return
types with conveniently packaged data.

Example B.9 (Functions Combining Multiple Return Values) It is possible to concatenate the
return values of multiple functions to generate new functions with compound return values. This is
especially useful when constructing outputs with reused variables or when it is convenient to keep
related values together.

1 function [x , r1 , r2 ] = comparison ( xi , xf , N )
2 x = linspace ( xi , xf , N ) ;
3 r1 = y1 ( x ) ;
4 r2 = y2 ( x ) ;
5 end

The function comparison defined on line 1 takes as inputs the initial value xi, the final value
xf and the number of points N on the interval where the functions y1 and y2 should be evaluated
and returns an array of points x and the corresponding function valuies y1 and y2 packaged in an
array. The array of evenly spaced points x is computed from the inputs to the comparison function
on line 2 and constructed using the built-in function linspace. Then the corresponding values of
y1 and y2 are computed using x in line 3 and 4, respectively.

57
 B.3 Plotting
Plotting predictions from numerical models is an important tool in inderstanding the behavious
of a model. Example B.10 demonstrates how a plot of numerically generated values from a custon
function nay be generated in MATLAB.

Example B.10 (Plotting) Below is a comparison plot of the function ouput values of of the response
variables y1 and y2 versus the independent parameter x corresponding tothe functions y1 (x ) and
y2 (x ) defined in Example B.7 and Example B.8.

1 % % Evaluate Functions
2 [x , y1 , y2 ] = comparison (0 , 5 , 1000) ;
3
4 % % Generate Plot
5 figure (1) % set graphics window label for plot to 1
6 xlabel ( ’x ’) % set x - axis label
7 ylabel ( ’ response ’) % set y - label label
8 title (" Simple Function Plot ") % set plot tile
9 hold on % start accumulating multiple plots
10 plot (t , y1 , ’ -r ’) % plot y1 versus x
11 plot (t , y2 , ’ -b ’) % plot y2 versus x
12 hold off % end accumulation of multiple plots

The plots of the functions y1 (x ) and y2 (x ) are generated by constructing tables of values for each
on a given interval of x values and then ploting pairs of points. The array of point for plotting
is generated by the call to the function comparison in line 2, that was discussed previously in
Example B.7. This return a list of x -values and the corresponding y1 and y2 values as an array with
columns x, y1 and y2, respectively. A graphical plot of these functions is generated by executing a
call to the plot function an passing the array of x values x and response values y1 or y2 to build
an plot of x versus y1 (x ) or y2 (x ). Calls to the plot function are found in lines 10 and 11. To
generate a single graphical output containing both plots on a single set of axes, we first set the
output function window in line 5, specify the x -axis label in line 6 and the y -axis label in line 7, set
the plot title in line 8, and then collate the list of plots using hold on directive in line 9 to instruct
MATLAB to keep each plot, without overwriting the graphics window output. The end of the list of
accumulated graphics initiated using the hold command on line 9 is terminated using hold off
on line 12. Additional plot options ’-r’ and ’-b’ will generate line graphs in red (r) and blue (b),
respectively. Alternative choices for plot markers include the discrete point parker options ’o’ and
i’*’. Addtional plot options can be found in the docoumentation.

B.4 Scripts
As for the case of writing functions, it is often useful to collect all function definitions and numerical
evaluation code in to a collection of one or more script files. Standard script (.m) files are simple

58
human readable text files that contain all function defintions and numerical value assignemts. In
addtion to the standard script files, MATLAB also has interactive script (.mlx) files that allow users
to interact more easily with code.

59
Bibliography

[1] F.R. Giordano, W.P. Fox, and S.B. Horton. A First Course in Mathematical Modeling. Cengage
Learning, 2013. ISBN 9781285531762. URL [Link]
PYUWAAAAQBAJ.

[2] B. Jacobs. Introduction to Mathematical Modelling, Course Notes. University of the

Witwatersrand, Johannesburg, 2018.

[3] W.A. Carlson. Methods of Applied Mathematics to Engineers, Course Notes. University of the
Witwatersrand, Johannesburg, 2018.

61