Lab 5: Nonlinear Systems

Goals
In this lab you will use the pplane8 program to study two nonlinear systems by direct nu-
merical simulation. The first model, from population biology, displays interesting nonlinear
oscillations (so-called limit cycles). The second is a system whose solutions depend on a
parameter. Neither of these systems is described by exactly solvable systems of differen-
tial equations. Although much may be learned from strictly theoretical analyses, we must
ultimately rely on computational methods to extract quantitative predictions from these
systems.

Application 1: predator-prey species interactions

In class we considered a model of predator-prey species interactions known as the Lotka-
Volterra model (referred to in Section 6.3 as the predator-prey system). If x describes the
size of a population of rabbits and y describes a population of foxes (which feed on the
population of rabbits) then the Lotka-Volterra model of their interactions says that there
are positive constants a, b, c, d so that
dx
= x(a − by) , (1)
dt
dy
= y(−c + dx) . (2)
dt
That is, the exponential growth rate of rabbits is decreased by the presence of foxes and
the exponential death rate of foxes is decreased by the presence of rabbits. This model
predicts some unlikely behavior. In the absence of foxes (y = 0), equation (1) becomes
dx/dt = ax. In other words, without any foxes the rabbits will always grow exponentially
without bound. And even if the predator population is small, they will always eat the prey
at a rate proportional to their product. In other words, 10 foxes surrounded by 100,000
rabbits would each have to eat ten times more than 10 foxes surrounded by 10,000 rabbits.
If the rabbit population could be held at a fixed level x0 > c/d, equation (2) becomes
dy/dt = Cy where C = −c + dx0 > 0. In other words, if the rabbit population is maintained
at a given level, above some threshold, the fox population will always grow exponentially
without bound. None of these predictions are ecologically reasonable. The following model
addresses these problems.
For positive values of r, a more reasonable model of the two populations is the system
dx 5xy
= x(1 − x) − , (3)
dt 5x + 1

1
 dy y
 
= ry 1 − . (4)
dt x
In the absence of predators, the prey satisfies the logistic equation with equilibrium popu-
lation x = 1. In the presence of predators, prey is consumed at a rate 5xy/(5x + 1). That
is, if x is large compared to 1/5, then 5xy/(5x + 1) ≈ y, and the predators consume prey at
a rate proportional to the predator population. On the other hand, if x is small compared
to 1/5 then 5xy/(5x + 1) ≈ 5xy and only then do the predators consume prey at a rate
proportional to xy as in the Lotka-Volterra model. Furthermore, if the prey population x is
held fixed somehow, we ignore the differential equation governing x(t) and replace x(t) by a
constant, and we then see that the predator population satisfies a logistic growth equation
with equilibrium population x (the predators and prey are equal in number in equilibrium).
For this logistic growth model, the parameter r is the inverse relaxation time for the predator
population, that is, 1/r is proportional to the time it takes the predator population to equi-
librate. (Note: We have already scaled the variables and chosen some parameter values in
equations (3) and (4). The general version of the model would have many more parameters.)

Application 2: bifurcation
When one tries to understand the behavior of a nonlinear system one of the first things one
looks at is the set of equilibrium solutions. The number and type of equilibrium solutions
may well depend on some parameter(s) of the system: the mass of a component, the stiffness
of a spring, the length of a lever, the resistance of an electronic component, and so on. In
this section of the lab you will observe in a very simple case how the number of equilibrium
points of a system of differential equations changes as a parameter varies. A change in the
number of equilibrium points, or in the stability of the equilibria, usually leads to a dramatic
qualitative change in the behavior of the non-equilibrium solutions as well. Such a qualitative
change is called a bifurcation and the associated parameter value at which the number or
stability of equilibrium points changes is called a bifurcation point. A model nonlinear system
of differential equations that exhibits bifurcations as a parameter varies is
dx
= ax − y , (5)
dt
dy
= x + ay + x2 . (6)
dt
In this system a is the parameter.

Prelab assignment
Before arriving in the lab, answer the following questions. You will need your answers in
lab to work the problems, and your recitation instructor may check that you have brought
them. These problems are to be handed in as part of your lab report.

2
 1. Find all of the fixed points (critical points, or equilibria) of the modified predator-prey
system, equations (3) and (4). Calculate the numerical value of the coexistence point
corresponding to positive values of x and y.

2. By setting dx/dt = 0, find the curve in the phase plane where the trajectories of (3)
and (4) are vertical, and by setting dy/dt = 0, find the curve in the phase plane where
the trajectories are horizontal. Note that these two curves, called respectively the x-
nullcline and the y-nullcline, are not trajectories themselves, but since it is known that
any trajectory that crosses the x-nullcline does so vertically, and that any trajectory
that crosses the y-nullcline does so horizontally, sketching these non-solution curves
can help in sketching the phase portrait. Use the information from problem 1 along
with these curves to sketch possible phase portraits for both small and larger values of
r (to be concrete, think about r = 0.01 as being small, and r = 0.1 as being larger).
Your sketches should be qualitatively different for different values of r.

3. The system exhibits very different behavior depending on whether r > rc or r < rc ,
where rc = .053576 . . .. In one regime, the coexistence point is stable and all solutions
are attracted to it. In the other region, the coexistence point is unstable and population
levels starting near the point spiral outward. Which do you think happens for which
values of r? That is, do you think that a large or a small value of r ought to correspond
to the stable coexistence or to the oscillatory behavior? (Either provide a coherent
logical argument based on your solutions to the previous problems, or do a stability
analysis of the coexistence fixed point to justify your prediction.)

4. Find the equilibrium points for the second system; that is, for equations (5) and (6).

In the lab
We will be using Matlab with the program pplane8 to study phase portraits for nonlinear
systems. The program pplane8 is an analogue of dfield8 that we used in Lab 1 that is
specially adapted to systems of two differential equations. Using pplane8, we can let Matlab
handle the numerical approximation of solutions and instead focus on the meaning of the
solutions and the way that solutions change as parameters are varied.

Studying phase portraits using pplane8

To study the evolution of the fox and rabbit populations over time, you will want to generate
a phase portrait plotting x against y. The following describes how you are to use pplane8
to generate these phase portraits.
After you log onto the Macintosh system in the Instructional Computing Lab in the
basement of East Hall, or onto your own system, launch Matlab and type pplane8 in the
Command Window and hit return. As happened with the dfield8 program you used in Lab 1
a window with the title PPLANE8 Setup will open with lots of little boxes all filled in; ignore
them for now and click on the Proceed button. You will see a graph with a direction field

3
corresponding to a system of differential equations. Put the cursor on any point and click.
You have just chosen initial conditions for the system. Now you know what the solution to
the system of differential equations with your choice of initial conditions looks like. Plot a
few solutions in this way. If your graph is getting too cluttered go to the Edit menu on the
PPLANE8 Display window and select Erase all solutions.
After you have plotted a few solutions in the PPLANE8 Display window, go to the Graph
menu on that window and select y vs t. Your cursor will become cross hairs; center the cross
hairs on a solution curve and click. (Be sure to center the cross hairs on a solution curve
you have already plotted, otherwise no graph will appear and you will see an error message
at the bottom of the PPLANE8 Display window telling you that This object is not a solution
curve.) You will see the corresponding plot of y versus t in a new window. You can change
your mind and click on x vs t, or both, or 3D, and so on, using the Graph menu to the right
of the plot in the new window.
Finally, go to the Solutions menu on the PPLANE8 Display window and choose Find an
equilibrium point; your cursor will again become cross hairs, and if you position the cross
hairs near an equilibrium point and click you will get a red dot at the equilibrium point and
some additional useful information concerning the nature of the equilibrium point in a small
box. You can repeat the command and find another equilibrium point for this system.
When you use pplane8 to do this lab you will, of course, need to change the system in
the PPLANE8 Setup window to the system you want to study. In the first case, the modified
predator-prey system, there is a parameter r in the system. You can enter the equations
with the symbol r in them and then, below the equations box, enter r = 0.5 (or whatever
value you wish) in the parameters or expressions box. (A parameter is a constant of the
problem that may however change its value from one problem to the next.)
Some important points: when you set up your equation you also enter the minimum and
maximum values for x and y as you did with dfield8 but it is often more convenient to
zoom in or zoom back. You will find those commands under the Edit menu. In addition you
should have the solver evolve the solution forward in time. This can be done by changing
the solution direction in the Options menu for the PPLANE8 Display window. (Looking at
the solution only in the forward direction tells you whether solutions are moving towards or
away from an equilibrium point.) Now you have the tools to do the lab.

Lab problems
1. Use pplane8 to solve the system of equations (3) and (4) and print phase portaits.
Start from various initial conditions, and use r = 0.07, r = 0.05, and r = 0.03. You
should zoom in on important features you see in the phase portraits.

2. Check your prediction from Prelab problem 3. What, if anything, surprised you about
the behavior of the system?

3. What is different about the oscillatory state here compared to that of the (simpler,
but less realistic) Lotka-Volterra system of differential equations (1) and (2)? Discuss.

4
 (Hint: Consider the dependence of the steady-state oscillation amplitude on the initial
condition. How many different closed orbits do you see for each value of r?)

4. Classify the equilibrium points of the second system (equations (5) and (6)) when
a = −0.5, a = 0, and a = 0.5.

5. Print out graphs of the phase portrait showing the behavior of solutions near each
equilibrium point (by using zoom in) in each of the three cases of the parameter a
described in the previous problem. What is the bifurcation point for the parameter a?

Note: your lab report for Lab 5 should consist of your solutions to the prelab
problems, your solutions to lab problems 1–5 (including any graphs printed out),
and a brief, original, conclusions paragraph summarizing what you have learned.