Overview of ODE’s for

Computational Science

Stephen Roberts, Daniel MacKinley, and Peter Humburg

Contents

1 Introduction

1

2 A Population Problem

2

2.1 A Very Simple ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2.2 Testing the solution against the data . . . . . . . . . . . . . . . . . . . . . . .

3

2.3 Numerical solutions in Scilab . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.4 Numerical versus Analytic methods . . . . . . . . . . . . . . . . . . . . . . . .

6

2.5 Compartmental Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

3 A more complex ODE — Limited Population growth

7

3.1 Formulating the Differential Equation . . . . . . . . . . . . . . . . . . . . . .

7

3.2 Equilibrium Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

3.3 Logistic Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

4 Systems of ODEs - Modelling Predators

9

4.1 Getting a handle on coupled ODEs . . . . . . . . . . . . . . . . . . . . . . . .

10

4.2 Numerical solutions to the fox and rabbit problem . . . . . . . . . . . . . . .

10

4.3 Phase Plane Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

4.4 Equilibrium Populations and Null Clines . . . . . . . . . . . . . . . . . . . . .

13

4.5 Explaining the oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

4.6 Adding Logistic Growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

5 Extended Predator-Prey Model

16

5.1 Functional Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

5.2 Numerical Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

5.3 Implementing the Extended Model . . . . . . . . . . . . . . . . . . . . . . . .

17

6 Case Study: Competition, Predation and Diversity

18

1

Introduction

This is a brief introduction to simple Ordinary Differential Equations (ODEs), using

them to model physical process and to the considerations of approximating their solution
computationally.

Ordinary Differential Equations are equations that relate the rate of change of some

quantity (such as position, velocity, concentration) with respect to a quantity such as time,
with the value of the quantity. In particular the derivative of a quantity such as position
with respect to time may be proportional to the position. Unlike simple algebraic equations,
the solution to a differential equation involves not just finding a particular number which
satisfies the equation, but a particular function which satisfies the equation. ODEs arise
naturally in a huge number of areas, including physics, biology, chemistry and ecosystems.
They crop up in any area that the rate of change of something is related to how much of it
there already is.

In the language of calculus, differential equations relate the first or higher derivatives of

an unknown function to the function itself.

We examine what ODEs mean in one very particular class of problems: modelling the

populations of various species of living things in an ecosystem. As we all know, the more
animals there are making babies the faster the population of animals grows. We can thus
expect differential equations to be useful at describing this.

1

ODE population models are hopefully easy to understand yet behave in a fairly inter-

esting way. More importantly, the techniques that we develop are generally applicable to
very large classes of differential equations.

In this module we will first see how differential equations can be used to model a range

of population models. We will introduce the idea of a compartmental model. A number
of useful theoretical tools will be introduced, such as the phase plane diagrams. We will
use Scilab, an open source numerical computing environment to approximate the solu-
tion of our differential model, and also see how to use Scilab’s more sophisticated ODE
solvers. Scilab’s syntax is very similar to Matlab. If you know either of these software
packages you should have no problems to follow the examples in this module. The latest
Scilab version is available for free download from http://www.scilab.org. If you want
to learn more about Scilab, try the Scilab primer and the numerical methods tutorials at
http://comptlsci.anu.edu.au/tutorials.html.

Throughout this module you will find inserts like this. They contain
questions and and small exercises about the material covered in the section
above them. Try to answer them. They should help you to get some
practical experience with ODEs.

2

A Population Problem

Let’s look at some (fake) population figures. Say we wanted to predict what would happen
with an accidental release of rabbits into a previously rabbit-free National Park. Fortunately
for us, there are records available of what happened last time that rabbits were released into
a rabbit free property, and we can use these records to try and find the rule that they obey.

From a nearby national park, we get the following population figures:

Month

Rabbit population

Monthly increase

Percentage increase

January

200

—

—

February

220

20

10%

March

243

23

10.5%

April

268

25

10.3%

May

296

28

10.4%

June

326

30

10.1%

July

358

32

9.8%

August

395

37

10.3%

September

434

39

9.9%

October

478

44

10.1%

November

527

49

10.3%

December

577

50

9.5%

January

634

57

9.9%

What is the pattern here? How can we simulate it computationally? How do we begin

to tackle this?

One way might be to look at how the data changes each month. We look at the percentage

rate of change of population: It seems that the rabbit population increases by approximately
10% each month.

2.1

A Very Simple ODE

We write this using calculus: It seems that perhaps we can model this rabbit population
using some function which increases by 10% of its current value each unit of time. Let us
represent the population of rabbits at a given time with a function X(t). This would be
represented by

2

dX

dt

= 0.1X,

X(0) = 200.

Writing our problem like this, relating a single unknown function of one variable, X(t)

and its derivatives, gives us what we call an Ordinary Differential Equation — an ODE.

The solution to this ODE is:

X(t) = 200e

0.1t

How do we get that?
It’s funny you should ask that: We notice that our equation can be written:

1

X

dX

dt

= 0.1

Integrate

Z

1

X

dX

dt

dt =

Z

0.1dt so

Z

1

X

dX = 0.1

Z

dt

So ln X = 0.1t + C for some arbitrary constant C.
Taking exponentials X(t) = Ae

0.1t

, where A = e

C

.

Using the fact that X(0) = 200 we get 200 = Ae

−kt

0

and A = 200e

−(0)(0.1)

Thus

X(t) = 200e

0.1(t)

Now we have solved the equations we set ourself. Does it model the data?

2.2

Testing the solution against the data

To see if our solution really fits the data we use Scilab to produce a plot of the data as

well as of the function X(t) = 200e

0.1(t)

. Take a look at this piece of Scilab code:

function x=f(t)

x=200*exp(0.1*t)

endfunction

tData = 0:12;
tFunction=0:0.1:12;

yData = [200;...

220;...
243;...
268;...
296;...
326;...
358;...
395;...
434;...
478;...
527;...
577;...
634];

yFunction = feval(tFunction,f);

xbasc();

plot2d(tData,yData,style=[-1,2]);

3

plot2d(tFunction,yFunction,style=[3,3]);
legend("Measured data","Analytical solution",2);

This script produces a plot that shows the measured data as black crosses and the

function X(t) as a green line. The result is shown in Figure 1.

0

2

4

6

8

10

12

200

250

300

350

400

450

500

550

600

650

700

Measured data

Analytical solution

Figure 1: Comparison of the data with the analytical solution of our equation X(t) = 200e

0.1(t)

.

Success! Or at least, close enough to it. We appear to have produced a model for a

rabbit population, and solved a differential equation along the way. Examining Figure 1 we
notice that our analytical solution is diverging slightly from the measured data. But the
result is still accurate enough for our purpose.

Note that our ‘solution’ in this case still has a variable in it — t. We haven’t found a

single specific number, but a function that describes the behaviour of the populations of
rabbits for all times t > 0.

Can you explain the difference between the analytical solution and the
data? Try to think of any assumptions we made on the way.

2.3

Numerical solutions in Scilab

There is another way to get a solution to the rabbit problem. Instead of solving the ODE

explicitly we can use Scilab to solve the differential equation for us using the following code:

function xDot=g(t,x)

xDot=0.1*x;

endfunction

tData = 0:12;
tFunction = 0:0.1:12;

xODE = ode(200,0,tFunction,g);

xbasc();

4

plot2d(tData,yData,style=[-1,2]);
plot2d(tFunction,xODE,style=[3,3]);
legend("Measured data","Numerical solution",2);

This produces the graph in Figure 2. It looks identical to Figure 1, and, in fact, it is to

0

2

4

6

8

10

12

200

250

300

350

400

450

500

550

600

650

700

Measured data

Numerical solution

Figure 2: Comparison of the data with the numerical solution of the differential equation

dX

dt

=

0.1X.

5 decimal places!

The trick here is: Scilab can estimate numerical solutions to the ordinary differential

equation without us having to go to the trouble of solving them.

How? Can Scilab tell us explicitly that it has found that the solution to our ODE is

X(t) = 200e

0.1(t)

? Well, no.

What Scilab can do is tell us how many rabbits there are at any time t using its internal

numerical ODE solver, without ever knowing what the explicit formula is.

Examining our program a little closer, this little snippet

function xDot=g(t,x)

xDot=0.1*x

endfunction

is Scilab-ese for

dX

dt

= 0.1X

and this command

xODE=ode(200,0,tFunction,g)

is enough to find the solution numerically for a given set of time values.

The numerical and analytical solutions to the rabbit problem are very
similar but not quite identical. Can you produce a plot in Scilab that
reveals the difference?

5

2.4

Numerical versus Analytic methods

We have just found a way of telling Scilab: “Give me an estimate of what will happen

to the population of rabbits that is constantly increasing at a rate of 10% per month, if
they start with a population of 200.” Without ever mentioning the explicit solution to our
problem, Scilab has done a very passable job.

What is the difference between the explicit formula solution and this computed numerical

solution then?

The thing is that the Scilab version is just as easy for more complicated differential

equations, and even for ones that are impossible to find numerical solutions for.

However, the numerical solution is not exactly the same as the analytic solution. As we

pointed out before, in this case they agree to 5 decimal places- this is enough to be identical
to the human eye, but very far from identical! However, this sort of accuracy might be
enough — in the case of our rabbit population model, where the population figures are a
little imprecise anyway, we would probably never notice the difference. How do we know if
we have an extra 0.01 of one rabbit?

In any case, working out how accurate this ODE solver is, and choosing how accurate

to try to make it is a big subject — and one of the major goals of this course.

2.5

Compartmental Models

How can we relate these to the differential equation that we have constructed? Let’s

pick apart our assumptions. In particular, let’s look at what we assume is happening when
we blithely model ’growth’ in a rabbit population.

For a start, rabbits are not immortal. Obviously, rabbits are dying and being born all

the time. The population growth that we have described will only occur if the birth rate
outweighs the death rate. If we had a very high death rate for some reason, say disease or
starvation, we might expect to find a decreasing population of rabbits.

Our rabbit population model is one of the class of models using ordinary differential

equations known as compartmental models.







net rate of

change of

rabbit population







=

½

rate of

births

¾

−

½

death

rate

¾

½

rate of

births

¾

= βX(t).

½

death

rate

¾

= αX(t).

Adding these together, we get the following differential equation

dX

dt

= (β − α)X

If we set r = β − α then we may write this as

dX

dt

= rX

Our initial problem was precisely an equation describing a model with r = 0.1.
This technique of considering the inflow and outflows to the rabbit population will be

used later on when we consider more complex sets of equations.

inflow

−−−−−→

world

outflow

−−−−−→

Staying on our hypothetical rabbit population, inflow becomes rabbit births, and outflow

becomes rabbit deaths. Considering the model in this fashion is useful if we want to extend
it by having multiple compartments with complex interaction - for example, by looking

6

at the interaction of predators and prey, or even more complicated systems - say rabbits,
pasture, cows, foxes, native marsupials...

(

rate of

change of

population

)

=

(

rate

of

births

)

−

(

rate

of

deaths

)

3

A more complex ODE — Limited Population growth

It doesn’t take much thought to realise that with our current model the rabbit population

will increase forever, until eventually there are more rabbits than there can possibly fit in
one little national park.

If our model were correct then we should expect there to be exponentially expanding

masses of rabbits to be forcing their way into the room right now. This is clearly (hopefully?)
not the case.

Let’s say that, in the wild, the rabbit population exhibits a different sort of behaviour:

one where they start out growing exponentially, but as they approach the capacity of the
land to bear them, they begin to slow down their rate of population change - say, as they
start to consume all the vegetation, rabbits have a greater risk of staving to death, litters
are less likely to survive, rabbits are more likely to fight... If we extended our original
rabbit population data, that sort of pattern of growth might gradually approach some finite
number — the ‘carrying capacity’ — which denotes the largest number of rabbits that can
live on a piece of land indefinitely.

1

In short, we’re looking for something like the S-curve in Figure 3.
This graph reflects the following assumptions:

•

When rabbit populations are small, they will grow approximately exponentially.

•

As an (initially small) rabbit population grows, individuals eventually will compete
for limited resources.

•

A given environment can support only a limited number of individuals.

This limit number is called the carrying capacity (denoted K).
In our logistic equation, similarly, we are witnessing not a sudden freeze in the population

of rabbits, where a bunch of immortal bunnies hang about producing no offspring. Rather,
we are assuming that at the carrying capacity the birth rate has equalled the death rate,
and so the growth of the population stops.

3.1

Formulating the Differential Equation

With that in mind, we might try to construct an improved model, using differential

equations again:

As before we assume

n

rate of

births

o

= βX(t).

Now assume that the per capita death rate depends on population

(

per capita

death

rate

)

= α + γX(t).

A larger population results in a higher per capita death rate, for all the reasons we just

listed.

So

n

rate of

deaths

o

= (α + γX(t))X(t).

Now, for convenience we choose r = β − α. Then our final differential equation looks

like this:

dX

dt

= rX − γX

2

¡

= βX − αX − γX

2

¢

1

But actually, this is an ecological fallacy - There is in general no set carrying capacity for a given

ecosystem, and assuming there is one usually has disastrous consequences! But we’re here for the maths,
not the wildlife management, right?

7

0

10

20

30

40

50

60

70

80

0

2e3

4e3

6e3

8e3

10e3

12e3

Rabbit Population

0

10

20

30

40

50

60

70

80

0

2e3

4e3

6e3

8e3

10e3

12e3

Limit Population

Figure 3: Development of rabbit population with carrying capacity. The population is gradually
approaching the maximum number of individuals.

Does this give us a limited population size? Has it produced a model of rabbit populations
on land with a limited carrying capacity?

We answer that by looking at the equilibrium solutions.

3.2

Equilibrium Solutions

If we have formulated the equations we set out to, we would like to know that at some

number of rabbits, the population will not be growing, i.e.

dX

dt

= 0.

Finding static population levels like this is looking for equilibrium solutions to the

differential equation.

We look for solutions to the equation

dX

dt

= 0

i.e. rX − γX

2

= 0

Solving

rX − γX

2

= 0

⇒ X(r − γX) = 0
⇒ X = 0 or X = r/γ

There are two equilibrium solutions.
The first, X(t) = 0 correspond to no rabbits at all - and since it takes rabbits to make

rabbits, it makes sense that one possible way of having a steady population of rabbits is to
have no rabbits at all.

8

The second equilibrium solution, X = r/γ is more interesting. This is precisely the

carrying capacity that we sought earlier. That is K = r/γ.

Hence our differential equation becomes:

dX

dt

= rX

³

1 −

γ

r

X

´

dX

dt

= rX(1 − X/K)

These fixed points are clearly depending on the initial populations of rabbits. However,

the solution as a whole does not. That is, we can have the same differential equation
describing our rabbits’ population for a different starting population. This is important!

One of the common characteristics of ODEs is that solutions depend on both the dif-

ferential equations, but also on the constraints that we set on the problem - the initial
conditions and the boundary values. When we find a particular function that satisfies both
the initial conditions of our problem and the differential equation that defines it, we say
that we have solved the initial value problem (the IVP) for that ODE.

In general certain characteristics of ODE solutions are given by the ODE itself, and

certain ones are given by the initial conditions. For example, in this case, the equilibrium
points of our solutions are given by the ODE that defines it. But whether the rabbit
population ever approaches this equilibrium value requires us to know the initial value.

To see this, ignore for a moment the fact that rabbits are supposed to
have non-negative populations and look at what happens to the function
if we set X(0) = −10.

3.3

Logistic Equation

If

dX

dt

= rX(1 − X/K)

Then we can solve the equation by separation of variables, much as before! The solution to
this ODE is given by

X(t) =

K

1 + m e

−rt

where m =

K

x

0

− 1

4

Systems of ODEs - Modelling Predators

We can take this technique further, solving several related functions simultaneously. The

quintessential example of this are the Predator-Prey ODEs. We will consider a particular
example of these, the Lotka-Volterra equations.

Let’s examine the consequences of introducing a new species whose growth is determined

by how many rabbits there are to eat...

n

Rate of change

of rabbit population

o

= { Rabbit birth rate } − { Rate of fox rabbit kills }

n

Rate of change

of fox population

o

= { Fox birth rate } − { Rate of fox deaths }

A little thought reveals that we will need not one but two differential equations to reflect

this. Since there are two populations, there are two related functions, one for each.

9

Let X(t) be the number of prey per unit area and Y (t) be the number of predator per

unit area.

Let β

1

> 0 be the prey per-capita birth rate and α

2

be the predator per-capita death

rate.

Let the prey deaths per-capita be c

1

Y (t) (i.e. dependent on the number of predators).

Let the rate of predator births per capita be c

2

X(t).

So, rate of prey births = β

1

X(t), rate of predator deaths = α

2

Y (t). Rate of prey killed

by predator = c

1

X(t)Y (t) and rate of predator births = c

2

X(t)Y (t). The term X(t)Y (t)

can be thought of s representing the chance of an encounter between prey and predators
at a given moment. Each ’encounter’ entails a fox eating a rabbit, and a decrease in the
rabbit population, and a smaller increase in the fox population (because a fox surely eats
more rabbits than it has babies in its life...)

dX

dt

= β

1

X − c

1

XY.

(1)

dY

dt

= c

2

XY − α

2

Y.

(2)

This is the Lotka-Volterra predator-prey system, and was proposed independently in the

1920s by two ecologists.

This equation has, in fact, no analytic solution - as much as we might play with calculus,

we will not find an equation which encapsulates this without an ordinary differential term.

4.1

Getting a handle on coupled ODEs

What will the solution look like?
If X = 0 (meaning, there is no prey),

dY

dt

= −α

2

Y , which implies exponential decay for

the predators. That is, the predators die out.

And if Y = 0 (no predators) then

dX

dt

= β

1

X, which implies that the prey grow expo-

nentially - exactly as in our very first rabbit population model. In fact, in the case that
Y (t) = 0 and β

1

= 0.1 we have found exactly the same model.

In the general case, though, with arbitrary values for the constants, if we wish to observe

the behaviour of the equations over time, we require Scilab to come to our aid.

4.2

Numerical solutions to the fox and rabbit problem

By now it should be fairly routine to produce Scilab code to calculate and display a

numerical solution for some arbitrarily chosen constants and initial conditions.

//define our constants

beta1 = 1;
alpha2 = 0.5;

c1 = 0.01;
c2 = 0.005;

//give the initial values
y0 = [200;80];

//define the predator-prey system of equations
function PopsDot=predprey(t,y)

PopsDot(1) = beta1*y(1)-c1*y(1)*y(2);
PopsDot(2) = c2*y(1)*y(2)-alpha2*y(2);

endfunction

//evaluate at the given time steps

10

t = [1.d-5:0.01:60];

y = ode(y0,0,t,predprey);

//plot the results (prey = blue curve)
xbasc();

plot2d(t,y(1,:),style = 4);
plot2d(t,y(2,:),style = 6);
legend("Rabbit Population","Fox Population");

0

10

20

30

40

50

60

20

40

60

80

100

120

140

160

180

200

220

Rabbit Population

Fox Population

Figure 4: Lotka-Volterra predator-pray system. The coupled rabbit and fox populations show
oscillating behaviour.

This produces the graph shown in Figure 4. Notice the oscillations in both populations.

This makes some sense intuitively and reflects the behaviour observed in some natural
systems.

This model has four parameters. Instead of the values we used above try
α

2

= 1.5. What happened? Can you explain this behaviour? Try other

sets of parameters. How do the different parameters influence the result?

Can we quantify these oscillations analytically?

4.3

Phase Plane Analysis

We get an idea of how X and Y behave with respect to one another by plotting the

‘phase plane’, which shows the direction of motion of the system at a given moment. This
is perhaps best demonstrated by example. Figure 5 shows the phase plane for the Lotka-
Volterra system discussed above.

The ‘arrows’ in the graph are produced by the following command in Scilab:

xbasc();
plot2d(y(1,:),y(2,:),style=4);
fchamp(predprey,0,[0:20:240],[0:20:200]);

11

0

50

100

150

200

250

0

20

40

60

80

100

120

140

160

180

200

Rabbits

Foxes

Figure 5: Phase plane for the Lotka-Volterra system.

a = get("current_axes");
a.x_label.text = "Rabbits";
a.y_label.text = "Foxes";

The first command, plot2d we have seen before. In this case it plots X versus Y rather

than X versus t.

fchamp is a new one. It plots, for a selection of points, the direction that the system

of populations of creatures is evolving. For each particular pair of possible values of the
populations, X and Y , the rate of change of the system along each axis is given by

dX

dt

and

dY

dt

respectively.

Each of the arrowed vectors in our plot is hence

µ

dX

dt

dY

dt

¶

. If both

dX

dt

and

dY

dt

are ‘smooth’

enough, this plot will show us how these values evolve in time.

Indeed, the blue path depicts a system of populations whose curve follows those vectors.

Note that the periodic behaviour we observed in Figure 4 corresponds to a closed curve in
the phase plane diagram. Such a curve is called a limit cycle.

For each different set of initial conditions, we will get a (possibly different) path through

the phase plane diagram - for example, Figure 6 is a similar plot, but with starting fox
populations of 1, 4, 16, and 64.

Can you work out which starting population was used for each curve in
Figure 6?

12

0

100

200

300

400

500

600

700

800

0

100

200

300

400

500

600

0

100

200

300

400

500

600

700

800

0

100

200

300

400

500

600

0

100

200

300

400

500

600

700

800

0

100

200

300

400

500

600

0

100

200

300

400

500

600

700

800

0

100

200

300

400

500

600

0

100

200

300

400

500

600

700

800

0

100

200

300

400

500

600

Figure 6: Phase diagram for starting fox populations of 1, 4, 16, and 64.

4.4

Equilibrium Populations and Null Clines

Now lets take a closer look at the phase plane diagram. You can see that some of the

arrows are parallel to one axes. This means that either

dX

dt

= 0 or

dY

dt

= 0. The points

where at least one of our differential equations is 0 all occur along a few lines. These are
called null clines. In Figure 7 the null clines were added to a schematic representation of
the phase plane diagram for our predator-prey model. They separate the plane into four
regions. In each region the two differential equations show a different behaviour. This is
indicated by the arrows in Figure 7.

Region I:

X decreases and Y increases.

Region II:

X decreases and Y decreases.

Region III:

X increases and Y decreases.

Region IV:

X increases and Y increases.

If both differential equations are 0 there is no change in the system. These points are
equilibrium solutions, just like the one we encountered in section 3.2. These equilibrium
solutions lie at the intersection of two null clines.

Before you continue to the next section take another look at Figure 7.
There are four intersections between horizontal and vertical lines. Only
two of these intersections are equilibrium solutions. Can you tell which
ones? Why is there no equilibrium at the other two points?

4.5

Explaining the oscillations

If initially the populations are x

0

and y

0

and we travel along the curve we eventually

return to x

0

, y

0

.

13

III

IV

I

II

X = 0

Y = β

1

/c

1

Y = 0

H

1

X = α

2

/c

2

H

2

V

1

V

2

@

@

@

I



@

@

@

R

6

?

6

?

-



-



6

?

?

6



-



-

Figure 7: The phase plane diagram is partitioned into four regions, each representing a different
general behaviour of the two differential equations.

To find the equilibrium points set

dX

dt

=

dY

dt

= 0. Then

X(β

1

− c

1

Y ) = 0

(3)

Y (c

2

X − α

2

) = 0

(4)

If X = 0 , then from equation 4, Y = 0. If β

1

− c

1

Y = 0, Y = β

1

/c

1

and X = α

2

/c

2

. The

equilibrium points are (0, 0) and (α

2

/c

2

, β

1

/c

1

). Figure 8 shows a phase plane diagram with

added null clines and marked equilibrium points.

Lets look at the equilibrium point (α

2

/c

2

, β

1

/c

1

). At this point the amount of rabbits

that is eaten by foxes during a time interval is equal to the amount of rabbits born dur-
ing that time. And the number of foxes born is equal the the number of foxes that die.
Both populations are in an equilibrium. But what happens when we move away from this
equilibrium point?

Consider region II We know that X <

α

2

c

2

and Y >

β

1

c

1

. We can write this as X =

α

2

c

2

−δ

1

,

Y =

β

1

c

1

+ δ

2

for δ

1/2

> 0. Using equations 1 and 2 we get

dX

dt

= X(β

1

− c

1

(

β

1

c

1

+ δ

2

)) = −Xc

1

δ

2

< 0

dY

dt

= Y (c

2

(

α

2

c

2

− δ

1

) − α

2

) = −Y c

2

δ

1

< 0

Thus, both the predator and the prey population decrease in region II. Using similar
calculations it is easy to see what happens in the other regions:

X > α

2

/c

2

, Y > β

1

/c

1

:

From Equation 2 Y increases, From Equation 1 X decreases.

X < α

2

/c

2

, Y > β

1

/c

1

:

From Equation 2 Y decreases, From Equation 1 X decreases.

X < α

2

/c

2

, Y < β

1

/c

1

:

From Equation 2 Y decreases, From Equation 1 X increases.

X > α

2

/c

2

, Y < β

1

/c

1

:

From Equation 2 Y increases, From Equation 1 X increases.

14

Examine the phase plane diagram and describe the behaviour of the
predator-prey system moving along trajectories in anti-clockwise direc-
tion. What happens when the system enters a new region? Try to identify
the four regions from the phase plane diagram in Figure 4.

4.6

Adding Logistic Growth

Recall

dX

dt

= β

1

X − c

1

XY. = X(β

1

− c

1

Y )

dY

dt

= c

2

XY − α

2

Y = Y (c

2

X − α

2

).

Previously we considered the case where β

1

= 1, α

2

= 0.5, c

1

= 0.01, c

2

= 0.005. These

equations produced oscillations in time for various conditions and parameter combinations.

The model assumes the prey grows exponentially in the absence of the predator. We

could replace the growth term for the prey population with a density dependent growth
with carrying capacity K. Then

dX

dt

= β

1

X

µ

1 −

X
K

¶

− c

1

XY.

dY

dt

= c

2

XY − α

2

Y.

We represent this in Scilab as

K = 800;

function PopsDot=predprey(t,y)

PopsDot(1)=beta1*y(1)*(1-y(1)/K)-c1*y(1)*y(2);
PopsDot(2)=c2*y(1)*y(2)-alpha2*y(2);

endfunction

Which gives us a slightly different set of results. Figure 9 shows the development of the

predator and prey populations over time. We see a damped oscillation where the changes
in the rabbit and fox populations get less extreme as they approach an equilibrium. This
kind of behaviour is easy to observe on our phase plane diagram (Figure 10). We can see
how the rabbit-fox system ‘spirals’ towards an equilibrium point. It’s a little more difficult
to show that the solutions approach an equilibrium point, rather than circling it. We need a
little more experience with partial derivatives, and a neat tool called the Jacobian matrix.
An introduction to that is a little beyond the scope of this document!

Find the null clines and equilibrium points for this new model. You can
check your results by comparing them to Figure 12

We can plot multiple starting populations with comparative ease within Scilab.

K=400

xbasc();
y0 = [2;200];
y = ode(y0,0,t,predprey2);
plot2d(y(1,:),y(2,:),style=1);

y0 = [10;200];
y = ode(y0,0,t,predprey2);

15

plot2d(y(1,:),y(2,:),style=4);

y0 = [50;200];
y = ode(y0,0,t,predprey2);
plot2d(y(1,:),y(2,:),style=3);

fchamp(predprey2,0,[0:20:500],[0:20:240]);

You can examine the output of this script in Figure 11. Note that we changed the carrying
capacity for the prey population to 400 to get a nicer plot.

Scilab can’t easily compute null-clines for us, but we have already done enough algebra

to persuade it to at least plot them for us. The lines marked in red are null-clines, and the
crosses mark equilibrium points, on the intersection of an X and a Y null cline. The result
is quite revealing. We see that the behaviour of solutions is, in some sense, conditioned by
the equilibrium solutions and null-clines.

5

Extended Predator-Prey Model

In section 4 we modelled a predator-prey system with a simple Lotka-Volterra model

using the two differential equations

dX

dt

= r

1

X − c

1

XY

dY

dt

= c

2

XY − α

2

Y

This model resulted in oscillatory behaviour of the system. We then refined our model by
adding logistic growth for the prey population:

dX

dt

= r

1

µ

1 −

X
K

¶

− c

1

XY

dY

dt

= c

2

XY − α

2

Y

This led to a loss of the oscillatory behaviour. Instead we observed damped oscillations while
the system approached an equilibrium. But in nature we observe oscillatory behaviour! Now
we want to extend our model even further. Maybe we can get the oscillatory behaviour back.

In practice, field ecologists measure

•

functional response

F (X)

•

numerical response

N (X, Y )

Functional response is the rate at which a single predator kills prey, as a function of prey

population. Numerical response is the per capita growth rate of the predators which may
depend on prey and predator populations.

dX

dt

= r

1

X

µ

1 −

X
K

¶

− F (X)Y

dY

dt

= N (X, Y )Y

5.1

Functional Response

Assume a single predator takes time t

h

to handle each prey. Consider a time interval [0, T ]

with

T = T

h

+ T

s

16

where T

h

is the time taken to handle all prey encountered before time T and T

s

is the time

left to search for new prey.

(

number of prey

eaten in time T

by one predator

)

= c

1

T

s

X(t) = N

e

(

rate of prey

eaten by one

predator

)

=

c

1

T

s

X(t)

T

Now T = T

h

+ T

s

= N

e

t

h

+ T

s

= c

1

T

s

X(t)t

h

+ T

S

.

(

rate of prey

eaten by one

predator

)

=

c

1

X

1 + t

h

c

1

X

Ã

F (x) → 0

as x → 0

F (x) →

1

t

h

as x → ∞

!

This is known as Holling type II response function.

5.2

Numerical Response

Now we assume logistic growth for predator population

r

2

µ

1 −

Y

K

2

¶

.

But assume the carrying capacity K

2

is proportional to prey density, ie K

2

= c

3

Y . Hence

N (X, Y ) = r

2

µ

1 −

Y

c

3

X

¶

Putting everything together we get an improved model:

dX

dt

= r

1

X

µ

1 −

X
K

¶

−

c

1

XY

1 + t

h

c

1

X

dY

dt

= r

2

Y

µ

1 −

Y

c

3

X

¶

This is known as Holling-Tanner predator-prey model. This extended model now has six
parameters:

r

1

> 0:

Prey birth rate.

K > 0:

Prey carrying capacity.

c

1

> 0:

Prey per capita death rate.

0 < t

h

< T :

Time a single predator takes to handle each prey.

r

2

> 0:

Predator birth rate.

K

2

= c

3

X > 0:

Predator carrying capacity.

5.3

Implementing the Extended Model

Again we can implement this in Scilab.

//define our constants
r1 = 0.16;
r2 = 0.05;

c1 = 0.006;
c3 = 0.6;

17

th = 2;

K = 800;

//give the initial values
y0 = [200;92];

//define the predator-prey system of equations
function PopsDot=predprey3(t,y)

PopsDot(1) = r1*y(1)*(1-y(1)/K)-(c1*y(1)*y(2)/(1+th*c1*y(1)));
PopsDot(2) = r2*y(2)*(1-y(2)/(c3*y(1)));

endfunction

//evaluate at the given time steps
t = [0d-5:0.1:600];

y = ode(y0,0,t,predprey3);

//plot the results (prey = blue curve)
xbasc();

plot2d(t,y(1,:),style = 4);
plot2d(t,y(2,:),style = 6);
legend("Rabbit Population","Fox Population");

The result of this is shown in Figure 13. Examining the Scilab code above you can

see that we used r

1

= 0.16, r

2

= 0.05, c

1

= 0.006, c

3

= 0.6, t

h

= 2, K = 800. This set

of parameters gives us a nice periodic behaviour. Unfortunately these parameter values are
not making much sense biologically. It seems unlikely that a fox needs two month to hunt
down and eat one rabbit. Also note that the oscillations are now much slower than before.
One period takes approximately 100 months!

Try different sets of parameters. How does the system react to different
parameter choices? How does this compare to the different parameter
sets for Lotka-Volterra model?
Can you find a more realistic parameter set that results in periodic be-
haviour?

Figure 14 shows the phase plane diagram for this system. Again you can see the limit

cycle. Note that we have adjusted the starting population of foxes from the previous example
to get a starting condition on the limit cycle.

Try different starting populations for this model. Examine the resulting
phase plane diagrams. What do you notice?

6

Case Study: Competition, Predation and Diversity

By now you should be fairly comfortable with handling ODEs. We have considered different
population models, implemented them in Scilab, and analysed the results. Now it is time
to try your skills on a more complex example.

In the last few sections we have considered systems with one prey and one predator

species. Now we will add a second prey species into the mix.

Consider two prey species N

1

, N

2

and one predator N

3

. In our case, we will be looking

at foxes and rabbits again. The rabbits are competing for the same resources and we should

18

consider the total number of rabbits for the per capita growth rate of both prey species.
The foxes are eating both rabbit species but they may have a preference for one. We allow
for this possibility by using different parameters for the interaction between foxes and the
two different rabbit species.

The general form for a system like this is

dN

1

dt

= N

1

(²

1

− α

11

N

1

− α

12

N

2

− α

13

N

3

)

dN

2

dt

= N

2

(²

2

− α

21

N

1

− α

22

N

2

− α

23

N

3

)

dN

3

dt

= N

3

(−²

3

+ α

31

N

1

+ α

32

N

2

− α

33

N

3

)

This has a simple competition model between N

1

and N

2

, and a predator-prey model be-

tween N

1

and N

3

, and N

2

and N

3

.

Using the above equations formulate a predator-prey model using a com-
bined carrying capacity K

1

for both rabbit species, i.e. K

1

≥ N

1

+ N

2

,

and a carrying capacity K

2

for the foxes that depends on the amount of

available prey.
Implement this in Scilab. Choose a suitable set of parameters and plot
the result. Produce a time as well as a phase plane diagram. Interpret
the results.
Try different parameter sets to see how they influence the result.

Two species competition models often lead to extinction of one species. Nature seems

to be more stable. Stability is enhanced with more interaction of species.

19

−50

0

50

100

150

200

250

−50

0

50

100

150

200

250

Rabbits

Foxes

Figure 8: Phase plane diagram for the Lotka-Volterra predator-prey model with added null clines.
Dark blue lines indicate points where

dY

dt

= 0 and red lines indicate points where

dX

dt

= 0. Equilib-

rium points are marked with crosses.

0

10

20

30

40

50

60

50

70

90

110

130

150

170

190

210

0

10

20

30

40

50

60

50

70

90

110

130

150

170

190

210

Figure 9: Introducing a carrying capacity for the prey population results in a damped oscillation.

20

0

40

80

120

160

200

240

0

20

40

60

80

100

120

140

160

0

40

80

120

160

200

240

0

20

40

60

80

100

120

140

160

Figure 10: Phase plane diagram for the Lotka-Volterra model with carrying capacity for the prey
population. The system approaches an equilibrium.

0

50

100

150

200

250

300

350

400

450

500

0

50

100

150

200

250

Figure 11: Phase plane diagram for different starting populations.

21

−100

0

100

200

300

400

500

−50

0

50

100

150

200

250

Figure 12: Phase plane diagram for different starting population. Null clines are shown by red
lines, the equilibrium points are marked with crosses.

0

100

200

300

400

500

600

0

50

100

150

200

250

Rabbit Population

Fox Population

Figure 13: Development of a Holling-Tanner predator-prey system over time.

22

0

50

100

150

200

250

0

20

40

60

80

100

120

140

160

Figure 14: Phase plane diagram for a Holling-Tanner predator-prey model.

23