Models, Mathematical Analysis and Simulation, Part I, Fall 2005
Homework 5, Phase Portraits.
(Score: 7.0)
Three problems are to be solved in this homework assignment. The rst consists in scaling
of a homogeneous and a nonhonogeneous differential equation. The second is a study of a
dynamical system with a simple bifurcation, and the third problem deals with predator-prey
models. Hints for making phaseplanes with Matlab are annexed.
Problem 1 (1.0).
The differential equation for damped free oscillations is a second order homogeneous linear
equation. It can be written in the following form:
m
d2u
du
+ c + ku = 0.
2
dt
dt
Here the state variable u is the deviation from a stationary position and t is the time. The equation
contains three parameters, i.e. the mass m, the damping constant c, and the spring constant k.
We rst introduce a dimensionless time by scaling of the original time variable t . After this
is done, we derive a differential equation for u() from the equation above for u(t ). One of the
results is that the equation for nondimensionalized variables will often contain fewer parameters
than the original equation. This fact will simplify the formal mathematical treatment of the
problem, and it will also lead to improved insight about both the problem and its solution.
We reparametrize by introducing two parameters as follows: The undamped (angular) fre
quency 0 is dened by 0 = k/m and the critical damping c0 is dened by c0 = 2 km.
Note that these new parameters are not dimensionless. We use them to dene the dimensionless
damping constant = c/c0 and the dimensionless time = 0 t .
Your rst task is to derive the differential equation for u() and to show that it equals
u + 2u + u = 0.
Observe that the scaling of the time has reduced the number of parameters from three to one.
With forced oscillations we study the differential equation
m
d2u
du
+ c + ku = F cos t .
dt 2
dt
Here there are two additional parameters, the amplitude F and the (angular) frequency of the
external force.
It is easy to nondimensionalize the driving frequency since we already have introduced the
reference frequency 0 . Thus we introduce a dimensionless driving frequency . Furthermore
1
we introduce a dimensionless state variable y by proper choice of scale factor (this should contain
F , of course).
Your second task is to derive a differential equation for y(t ) In this case the scaling of both
the dependent and the independent variables has reduced the number of parameters from ve to
two: and .
Problem 2 (3.0).
Consider the following linear system of differential equations:
u =
11
P1
u,
where P is a realvalued parameter.
(0.75) a) Determine the type of phase portrait as a function of P. (There are ve different types.)
Determine the slopes of the named manifolds (in those cases where they exist) as a function
of P.
(0.75) b) Determine the limit of u(t ) and of y(t )/x(t ) as t in all ve cases. Note that the
result may depend on the initial point u(0) and on the value of P, Determine the bifurcation
value P0 of P and describe the type of phase portrait on either side of it.
(0.75) c) Use Matlab to plot one phase portrait of each of the ve types, with your choice of
Pvalues. Include the named manifolds when they exist.
(0.75) d) There is a transition from one name to another of each named manifold when P passes
from a value below P0 to a value above P0 . Describe this transition for each manifold.
Describe in particular the boundary case of each of the manifolds when P = P0 .
Problem 3 (3.0).
(0.5) (a) Consider the two-species models of the form,
x = a(x, y)x,
y = b(x, y)y,
(1)
where x and y denote the populations of two species, and a and b denote the corresponding
growth rates. Assume that a and b are smooth functions. We are now interested in a
predator-prey model, where y denotes the predator and x denotes the prey population,
which have the following (meaningful) assumptions:
(i) If there is not enough prey, the predator population declines.
(ii) An increase in the prey population increases the growth rate of the predator.
(iii) If no predators are present, a small prey population will increase.
2
(iv) If the prey population goes beyond a certain size, it must decrease.
(v) If the predator population increases, the growth rate of the prey population declines.
Express the assumptions (i)(v) in terms of a, b, x and y. Verify that the simple model with
a(x, y) = y x and b(x, y) = x y (, , , , , > 0) satisfy (i)(v).
(2.2) (b) Consider the predator-prey model of a),
dx
dt
dy
dt
= ( y x)x
,
= (x y)y
where , , , , and are positive constants. Show that x(t ), y(t ) > 0 for t > 0 if
x(0), y(0) > 0. Determine the critical points of this system and determine their stability.
Hint: to reduce the number of parameters, introduce the scaling, := ct , u := dx and
v := ey.
Choose the scalings to make the transformed equations have three parameters which are
combinations of the original ones, and such that the Strang Predator-Prey model comes as
a special case:
u = (1 v Au)u
,
v = (u C Bv)v
What is the condition on A, B, and C that there be a strictly positive critical point? Choose
two sets of values, one to make the point a spiral, and one to make it a node. Hand in
plots of the phase portrait of the linearized equations and original equations around the
critical points which shows the difference between the solution of the linearized and nonlinearized equations.
(0.3) (c) In his proof that the predator-prey populations are cyclic, Strang says that the curve
a log u2 bu2 = cu1 d log u1 + constant,
(2)
is closed. Let u1 be xed. Show that (2) has either 0, 1 or 2 solutions for u2 (depending on
u1 ).
3
How to plot phase portraits with Matlab.
There are some orbits in phase portraits that are more important, and carry more information,
than others. Among the most important are the orbits on the so-called named manifolds. The
term named manifolds is used to refer to the slow and fast manifolds of nodes and the stable and
unstable manifolds of saddle points. We require of any plot of phase portraits that it includes
orbits on all named manifolds.
We show by treating a concrete example how you can use Matlab to plot the phase portrait
of a linear system in the plane.
The system of equations is written u = AA*u, where AA is a given 2 2-matrix and u is a
column vector. The matrix is denoted in Matlab with two letters instead of one since it will be
a so-called global variable, and it is customary in Matlab to denote global variables with long
names with capital letters.
The rst step is to write an m-le that computes the left-hand side of the system of equations
when the righthand side is known. Let the name of the le be e.g. lin.m. The le may look
like follows:
% lin.m Linear diff eq system in two variables.
function uprim = lin(t,u)
global AA
uprim = AA*u;
The variable AA is global. It is declared global both in this mle and in the main program. Its
value is set in the main program.
After you have saved this le, move to the main Matlab window and give the following
commands:
>>
>>
>>
>>
global AA
AA = [-1 1;1 -3];
tspan = [0 5];
Va = [-1 1 -1 1]; axis(Va), axis equal, hold on, grid
The second line sets the value AA. If you wish to with work another system you just feed in
a new matrix AA. The third line above sets the time span from initial value t0 = 0 to nal
value tf = 5 of the independent variable t. If you later wish to run backwards in time, then
all you have to do is to set tspan = [0,-5], that means t < 0. The command axis(Va)
in the fourth line denes a window in the x-y-plane, whose corners are determined by Va =
[xmin xmax ymin ymax]. axis equal sets equal units in both x- and y-directions. The
command hold on allows you to plot several curves in the same gure without removing
earlier plots. The command grid, nally, plots a grid in the plane.
The main tool for plotting phase portraits is Matlabs command ode23 or ode45. We
start plotting one orbit of the phase portrait by giving the following commands:
>> u0 = [1;1];
>> [t,u] = ode23(@lin,tspan,u0); plot(u(:,1),u(:,2))
4
The initial point is denoted u0. The semicolon makes it a column vector. The second line
contains the main command line. It will be used repeatedly in the sequel. Its rst command asks
Matlab to solve the system of equations in lin.m with the initial value u(t0) = u0. The
solution is given in the form of a column vector of t-values in the interval from t0 to tf, and a
matrix of corresponding values of u.
The number of t-values, and their spacing, is determined by the program. It is used here
with default values for tolerances etc. Read the help-le on ode23 (and ode45) to nd out
what they are and how to change them! You can check the number of steps the program has
taken by giving the command size(t). The matrix of u-values has two columns. The rst one
contains a vector of values of x, and the second one contains a column of values of y. Check with
size(u).
The second command on the main command line plots the orbit. Here, u(:,1) and u(:,2)
denote the rst and second columns of the matrix u.
When you wish to plot one more orbit, then you give a new initial point u0, after which you
repeat the main command line. You do not have to retype any commands; the old ones can be
recalled with the arrow-up key. Choose seven additional initial points uniformly distributed on
the circumference of the window determined by Va. You might want to experiment with a loop
like this1
n = 7;
st = exp(i*2*pi*(1:n)/n);
for s = st
u0 = [real(s);imag(s)];
[t,u] = ode23(@lin,tspan,u0); plot(u(:,1),u(:,2));
end
The phase portrait that you have produced so far shows eight orbits that all approach the origin
with the same slope. This behavior is typical of the stable node. To complete the phase portrait
you need to plot the orbits on the slow and fast manifolds.
You nd the orbits on the slow manifold by choosing initial points sufciently far away from
the origin. I suggest that you try u0 = [10;10] and then u0 = -u0.
The orbits on the fast manifold are found in a different way. One way would be to experiment
with a large number of initial points on the circumference of the window until you nd one that
gives you one of the orbits. An alternative is to start near the origin and run backwards in
time. The initial point should then be chosen in a direction from the origin determined by the
slope of the fast manifold. To determine the initial point, we need to determine eigenvalues and
eigenvectors of the matrix AA:
>> [V,D] = eig(AA)
The matrix V contains two eigenvectors as columns, and the matrix D has the corresponding
eigenvalues on its diagonal. The slopes of the two eigenvectors, sls and slf, are determined
with the following commands. The last letter in the slope designation is s for the slow manifold
and f for the fast one.
1A
more advanced possibility consists of using Matlab script les.
5
>> v1 = V(:,1);
>> v2 = V(:,2);
sls = v1(2)/v1(1)
slf = v2(2)/v2(1)
Change the numbering if the slow manifold eigenvector is the second column of V! To plot an
orbit on the fast manifold we choose an initial point whose coordinates are determined by the
fast manifold. In addition we run time backwards:
>> u0 = 0.02*v2; tspan = [0,-5];
After this we execute the main command line.
You get the other orbit on the same manifold by changing the sign of u0 and then executing
the main command line.
When the phase portrait is nished and you are satised with its esthetic appearance, you
may wish to equip it with a title. You can for example write as follows, replacing my name with
yours:
>> title(Phase portrait by my name)
>> xlabel(x), ylabel(y)
Finish off by giving the command print:
>> print -deps figure.eps
Even simpler, you may use the menu of the gure window (Export). Get your copy at the
printer, and equip the orbits with arrows that show the directions they are traversed as time
increases.
Let us plot one more phase portrait. The window is erased with the command clf (clear
gure). To prepare for a new phase portrait and reset the direction of the time axis to normal you
give the following command:
>> clf, axis(Va), axis equal, hold on, grid, tspan = [0,5];
Feed in a new matrix as follows:
>> AA = [-1 4;4 -3];
Now plot the orbits for a number of initial points on the boundary of the window determined
by Va. The orbits indicate that you are dealing with a saddle point. Use the command eig to
determine eigenvalues and eigenvectors of the matrix AA. Determine the slopes sls and slu
of the stable and unstable manifolds from the stable and unstable eigenvectors vs viz. vu. Use
initial points u0 = 0.02*v and u0 = -u0, where v is replaced by vs and vu. Let time run
forward to get the unstable manifold and backward to get the stable one.
It should be noted that the rules that we have followed for plotting orbits on the named
manifolds could be somewhat simplied if we are only concerned with linear systems. Our rules
have the virtue that they can also be used when we are dealing with nonlinear systems.
6
We summarize the rules for plotting orbits on the named manifolds of stable nodes and of
saddle points that have been used above:
The orbits on the slow manifold of a stable node are plotted with initial points far away from
the critical point, and with time running forward.
All other orbits are plotted with initial points close to the critical point. The initial point is
chosen on a straight line through the critical point whose slope is given by the corresponding
eigenvector.
The orbits on the unstable manifold of the saddle point are plotted with time running forward.
The orbits on both the stable manifold of the saddle point and on the fast manifold of the
stable node are plotted with time running backward.
An unstable node is transformed into a stable one by letting time run backward.
A note on using ode45 The initial value solver ode45 has a much higher order than
ode23. On one hand, this is nice: the computation time is much lower than with ode23. On the
other hand, the number of steps taken (e.g., size(t)) is much lower such that the phase plot
may become a polygon with visible corners. The latter can be avoided by using an alternative
main command line:
>>
>>
>>
>>
sol = ode45(@lin,tspan,u0);
t = linspace(tspan(1),tspan(end),101);
u = deval(sol,t);
plot(u(1,:),u(2,:))
% Note the transposition of u!
For details please consult Matlabs documentation.
G OOD
LUC K !
7