11.4 Autonomous Systems as Mathematical Models
INTRODUCTION
Many applications from physics give rise to nonlinear autonomous second-order differential equations—that is, DEs of the form x″ = g(x, x′). For example, in the analysis of free, damped motion of a spring/mass system in Section 3.8 we assumed that the damping force was proportional to the velocity x′ and the resulting model mx″ = –βx′ – kx is a linear differential equation. But if the magnitude of the damping force is proportional to the square of the velocity, the new differential equation mx″ = –βx′ – kx is nonlinear. The corresponding plane autonomous system is nonlinear:
In this section we will also analyze the nonlinear pendulum, motion of a bead on a curve, the Lotka–Volterra predator–prey models, and the Lotka–Volterra competition model. Additional models are presented in these exercises.
Nonlinear Pendulum
In (6) of Section 3.11 we showed that the displacement angle θ for a simple pendulum satisfies the nonlinear second-order differential equation
When we let x = θ and y = θ′, this second-order differential equation may be rewritten as the dynamical system
The critical points are (±kπ, 0), and the Jacobian matrix is easily shown to be
If k = 2n + 1, Δ < 0, and so all critical points (±(2n + 1)π, 0) are saddle points. In particular, the critical point at (π, 0) is unstable as expected. See FIGURE 11.4.1. When k = 2n, the eigenvalues are pure imaginary, and so the nature of these critical points remains in doubt. Since we have assumed that there are no damping forces acting on the pendulum, we expect that all the critical points (±2nπ, 0) are centers. This can be verified using the phase-plane method. From
it follows that
If X(0) = (x0, 0), then
Note that y = 0 when x = –x0, and that (2g/l)(cos x – cos x0) > 0 for < < π. Thus each such x has two corresponding values of y, and so the solution X = X(t) that satisfies X(0) = (x0, 0) is periodic. We may conclude that (0, 0) is a center. Observe that x = θ increases for solutions, such as the one drawn in red in FIGURE 11.4.2, that correspond to large initial velocities. In this case the pendulum spins in complete circles about its pivot.
EXAMPLE 1 Periodic Solutions of the Pendulum DE
A pendulum in an equilibrium position with θ = 0 is given an initial angular velocity of ω0 rad/s. Determine under what conditions the resulting motion is periodic.
SOLUTION
We are asked to examine the solution of the plane autonomous system that satisfies X(0) = (0, ω0). From y2 = (2g/l) cos x + c it follows that
To establish that the solution X(t) is periodic it is sufficient to show that there are two x-intercepts x = ±x0 between –π and π and that the right-hand side is positive for |x| < |x0|. Each such x then has two corresponding values of y.
If y = 0, cos x = 1 – (l/2g), and this equation has two solutions x = ±x0 between –π and π, provided 1 – (l/2g) > –1. Note that (2g/l)(cos x – cos x0) is then positive for |x| < |x0|. This restriction on the initial angular velocity can be written as
≡
Nonlinear Oscillations: The Sliding Bead
Suppose, as shown in FIGURE 11.4.3, a bead with mass m slides along a thin wire whose shape is described by the function z = f(x). A wide variety of nonlinear oscillations can be obtained by changing the shape of the wire and by making different assumptions about the forces acting on the bead.
The tangential force F due to the weight W = mg has magnitude mg sin θ, and therefore the x-component of F is Fx = –mg sin θ cos θ. Since tan θ = f′(x), we can use the identities 1 + tan2 θ = sec2 θ and sin2 θ = 1 – cos2 θ to conclude that
Fx = –mg sin θ cos θ = –mg.
We assume (as in Section 3.8) that a damping force D, acting in the direction opposite to the motion, is a constant multiple of the velocity of the bead. The x-component of D is therefore
Dx = –β .
If we ignore the frictional force between the wire and the bead and assume that no other external forces are impressed on the system, it follows from Newton’s second law that
and the corresponding plane autonomous system is
If X1 = (x1, y1) is a critical point of the system, y1 = 0 and therefore f′(x1) = 0. The bead must therefore be at rest at a point on the wire where the tangent line is horizontal. When f is twice differentiable, the Jacobian matrix at X1 is
and so τ = – β/m, Δ = g f″(x1), and τ2 – 4Δ = β2/m2 – 4g f″(x1). Using the results of Section 11.3, we can make the following conclusions:
- f″(x1) < 0:
A relative maximum therefore occurs at x = x1, and since Δ < 0, an unstable saddle point occurs at X1 = (x1, 0). - f″(x1) > 0 and β > 0:
A relative minimum therefore occurs at x = x1, and since τ < 0 and Δ > 0, X1 = (x1, 0) is a stable critical point. If β2 > 4gm2f″(x1), the system is overdamped and the critical point is a stable node. If β2 < 4gm2f″(x1), the system is underdamped and the critical point is a stable spiral point. The exact nature of the stable critical point is still in doubt if β2 = 4gm2f″(x1). - f″(x1) > 0 and the system is undamped (β = 0):
In this case the eigenvalues are pure imaginary, but the phase-plane method can be used to show that the critical point is a center. Therefore solutions with X(0) = (x(0), x′(0)) near X1 = (x1, 0) are periodic.
EXAMPLE 2 Bead Sliding Along a Sine Wave
A 10-gram bead slides along the graph of z = sin x. According to conclusion (ii), the relative minima at x1 = –π/2 and 3π/2 give rise to stable critical points. See FIGURE 11.4.4. Since f″(–π/2) = f″(3π/2) = 1, the system will be underdamped provided β2 < 4gm2. If we use SI units, m = 0.01 kg and g = 9.8 m/s2, and so the condition for an underdamped system becomes β2 < 3.92 × 10–3.
If β = 0.01 is the damping constant, both of these critical points are stable spiral points. The two solutions corresponding to initial conditions X(0) = (x(0), x′(0)) = (–2π, 10) and X(0) = (–2π, 15), respectively, were obtained from a numerical solver and are shown in FIGURE 11.4.5. When x′(0) = 10, the bead has enough momentum to make it over the hill at x = –3π/2 but not over the hill at x = π/2. The bead then approaches the relative minimum based at x = –π/2. If x′(0) = 15, the bead has the momentum to make it over both hills, but then it rocks back and forth in the valley based at x = 3π/2 and approaches the point (3π/2, –1) on the wire. Experiment with other initial conditions using your numerical solver.
FIGURE 11.4.6 shows a collection of solution curves obtained from a numerical solver for the undamped case. Since β = 0, the critical points corresponding to x1 = –π/2 and 3π/2 are now centers. When X(0) = (–2π, 10), the bead has sufficient momentum to move over all hills. The figure also indicates that when the bead is released from rest at a position on the wire between x = –3π/2 and x = π/2, the resulting motion is periodic.≡
Lotka–Volterra Predator–Prey Model
A predator–prey interaction between two species occurs when one species (the predator) feeds on a second species (the prey). For example, the snowy owl feeds almost exclusively on a common arctic rodent called a lemming, while a lemming uses arctic tundra plants as its food supply. Interest in using mathematics to help explain predator–prey interactions has been stimulated by the observation of population cycles in many arctic mammals. In the MacKenzie River district of Canada, for example, the principal prey of the lynx is the snowshoe hare, and both populations cycle with a period of about 10 years.
There are many predator–prey models that lead to plane autonomous systems with at least one periodic solution. The first such model was constructed independently by the American mathematician and biophysicist Alfred James Lotka (1880–1949) and the Italian mathematician and physicist Vito Volterra (1860–1940). If x denotes the number of predators and y denotes the number of prey, then the Lotka–Volterra model takes the form
x′ = –ax + bxy = x(–a + by)
y′ = –cxy + dy = y(–cx + d),
where a, b, c, and d are positive constants.
Note that in the absence of predators (x = 0), y′ = dy, and so the number of prey grows exponentially. In the absence of prey, x′ = –ax, and so the predator population becomes extinct. The term –cxy represents the death rate due to predation. The model therefore assumes that this death rate is directly proportional to the number of possible encounters xy between predator and prey at a particular time t, and the term bxy represents the resulting positive contribution to the predator population.
The critical points of this plane autonomous system are (0, 0) and (d/c, a/b), and the corresponding Jacobian matrices are
The critical point at (0, 0) is a saddle point, and FIGURE 11.4.7 shows a typical profile of solutions that are in the first quadrant and near (0, 0).
Since the matrix A2 has pure imaginary eigenvalues λ = ±i, the critical point (d/c, a/b) may be a center. This possibility can be investigated using the phase-plane method. Since
we may separate variables and obtain
so that
or
The following argument establishes that all solution curves that originate in the first quadrant are periodic.
Typical graphs of the nonnegative functions F(x) = xde–cx and G(y) = yae–by are shown in FIGURE 11.4.8. It is not hard to show that F(x) has an absolute maximum at x = d/c, whereas G(y) has an absolute maximum at y = a/b. Note that, with the exception of 0 and the absolute maximum, F and G each take on all values in their range precisely twice.
These graphs can be used to establish the following properties of a solution curve that originates at a noncritical point (x0, y0) in the first quadrant.
- If y = a/b, the equation F(x)G(y) = c0 has exactly two solutions xm and xM that satisfy xm < d/c < xM.
- If xm < x1 < xM and x = x1, then F(x)G(y) = c0 has exactly two solutions y1 and y2 that satisfy y1 < a/b < y2.
- If x is outside the interval [xm, xM], then F(x)G(y) = c0 has no solutions.
We will give the demonstration of 1. and outline parts 2. and 3. in the exercises. Since (x0, y0) ≠ (d/c, a/b), F(x0)G(y0) < F(d/c)G(a/b). If y = a/b, then
Therefore F(x) = c0/G(a/b) has precisely two solutions xm and xM that satisfy xm < d/c < xM. The graph of a typical periodic solution is shown in FIGURE 11.4.9.
EXAMPLE 3 Predator–Prey Population Cycles
If we let a = 0.1, b = 0.002, c = 0.0025, and d = 0.2 in the Lotka–Volterra predator–prey model, the critical point in the first quadrant is (d/c, a/b) = (80, 50), and we know that this critical point is a center. See FIGURE 11.4.10 in which we have used a numerical solver to generate these cycles. The closer the initial condition X0 is to (80, 50), the more the periodic solutions resemble the elliptical solutions to the corresponding linear system. The eigenvalues of g′((80, 50)) are λ = ±i = ±(/10)i, and so the solutions near the critical point have period p ≈ 10π, or about 44.4.≡
Lotka–Volterra Competition Model
A competitive interaction occurs when two or more species compete for the food, water, light, and space resources of an ecosystem. The use of one of these resources by one population therefore inhibits the ability of another population to survive and grow. Under what conditions can two competing species coexist? A number of mathematical models have been constructed that offer insights into conditions that permit coexistence. If x denotes the number in species I and y denotes the number in species II, then the Lotka–Volterra model takes the form
(1)
Note that in the absence of species II (y = 0), x′ = (r1/K1)x(K1 – x), and so the first population grows logistically and approaches the steady-state population K1 (see Section 2.8 and Example 4 in Section 11.3). A similar statement holds for species II growing in the absence of species I. The term –α21xy in the second equation stems from the competitive effect of species I on species II. The model therefore assumes that this rate of inhibition is directly proportional to the number of possible competitive pairs xy at a particular time t.
This plane autonomous system has critical points at (0, 0), (K1, 0), and (0, K2). When α12α21 ≠ 0, the lines K1 – x – α12y = 0 and K2 – y – α21x = 0 intersect to produce a fourth critical point = (, ). FIGURE 11.4.11 shows the two conditions under which (, ) is in the first quadrant.
The trace and determinant of the Jacobian matrix at (, ) are, respectively,
In case (a), K1/α12 > K2 and K2/α21 > K1. It follows that α12α21 < 1, τ < 0, and Δ > 0. Since
τ2 – 4Δ > 0, and so (, ) is a stable node. Therefore if X(0) = X0 is sufficiently close to = (, ), limt→∞ X(t) = , and we may conclude that coexistence is possible. The demonstration that case (b) leads to a saddle point and the investigation of the nature of critical points at (0, 0), (K1, 0), and (0, K2) are left to the exercises.
When the competitive interactions between two species are weak, both of the coefficients α12 and α21 will be small, and so the conditions K1/α12 > K2 and K2/α21 > K1 may be satisfied. This might occur when there is a small overlap in the ranges of two predator species that hunt for a common prey.
EXAMPLE 4 A Lotka–Volterra Competition Model
A competitive interaction is described by the Lotka–Volterra competition model
x′ = 0.004x(50 – x – 0.75y)
y′ = 0.001y(100 – y – 3.0x).
Find and classify all critical points of the system.
SOLUTION
Critical points occur at (0, 0), (50, 0), (0, 100), and at the solution (20, 40) of the system
Since α12α21 = 2.25 > 1, we have case (b) in Figure 11.4.11, and so the critical point at (20, 40) is a saddle point. The Jacobian matrix is
and we obtain
Therefore (0, 0) is an unstable node, whereas both (50, 0) and (0, 100) are stable nodes. (Check this!) Since det A3 < 0, we have a second demonstration that (20, 40) is a saddle point.≡
Coexistence can also occur in the Lotka–Volterra competition model if there is at least one periodic solution lying entirely in the first quadrant. It is possible to show, however, that this model has no periodic solutions.
11.4 Exercises Answers to selected odd-numbered problems begin on page ANS-31.
The Nonlinear Pendulum
- A pendulum is released at θ = π/3 and is given an initial angular velocity of ω0 rad/s. Determine under what conditions the resulting motion is periodic.
-
- If a pendulum is released from rest at θ = θ0, show that the angular velocity is again 0 when θ = –θ0.
- The period T of the pendulum is the amount of time needed for θ to change from θ0 to –θ0 and back to θ0. Show that
The Sliding Bead
- A bead with mass m slides along a thin wire whose shape is described by the function z = f(x). If X1 = (x1, y1) is a critical point of the plane autonomous system associated with the sliding bead, verify that the Jacobian matrix at X1 is
- A bead with mass m slides along a thin wire whose shape is described by the function z = f(x). When f′(x1) = 0, f″(x1) > 0, and the system is undamped, the critical point X1 = (x1, 0) is a center. Estimate the period of the bead when x(0) is near x1 and x′(0) = 0.
- A bead is released from the position x(0) = x0 on the curve z = with initial velocity x′(0) = v0 cm/s.
- Use the phase-plane method to show that the resulting solution is periodic when the system is undamped.
- Show that the maximum height zmax to which the bead rises is given by
- Rework Problem 5 with z = cosh x.
Interaction Models
- Refer to Figure 11.4.9. If xm < x1 < xM and x = x1, show that F(x)G(y) = c0 has exactly two solutions y1 and y2 that satisfy y1 < a/b < y2. [Hint: First show that G(y) = c0/F(x1) < G(a/b).]
- From 1. and 3. on page 664, conclude that the maximum number of predators occurs when y = a/b.
- In many fishery science models the rate at which a species is caught is assumed to be directly proportional to its abundance. If both predator and prey are being exploited in this manner, the Lotka–Volterra differential equations take the form
x′ = – ax + bxy – ϵ1x
y′ = – cxy + dy – ϵ2 y,
where ϵ1 and ϵ2 are positive constants.
- When ϵ2 < d, show that there is a new critical point in the first quadrant that is a center.
- Volterra’s principle states that a moderate amount of exploitation increases the average number of prey and decreases the average number of predators. Is this fisheries model consistent with Volterra’s principle?
- A predator–prey interaction is described by the Lotka–Volterra model
x′ = –0.1x + 0.02xy
y′ = 0.2y – 0.025xy.
- Find the critical point in the first quadrant, and use a numerical solver to sketch some population cycles.
- Estimate the period of the periodic solutions that are close to the critical point in part (a).
- A competitive interaction is described by the Lotka–Volterra competition model
x′ = 0.08x(20 – 0.4x – 0.3y)
y′ = 0.06y(10 – 0.1y – 0.3x).
Find and classify all critical points of the system.
- In (1) show that (0, 0) is always an unstable node.
- In (1) show that (K1, 0) is a stable node when K1 > K2/α21 and a saddle point when K1 < K2/α21.
- Use Problems 12 and 13 to establish that (0, 0), (K1, 0), and (0, K2) are unstable when = (, ) is a stable node.
- In (1) show that = (, ) is a saddle point when
Miscellaneous Nonlinear Models
- If we assume that a damping force acts in a direction opposite to the motion of a pendulum and with a magnitude directly proportional to the angular velocity dθ/dt, the displacement angle θ for the pendulum satisfies the nonlinear second-order differential equation
- Write the second-order differential equation as a plane autonomous system, and find all critical points.
- Find a condition on m, l, and β that will make (0, 0) a stable spiral point.
- In the analysis of free, damped motion in Section 3.8 we assumed that the damping force was proportional to the velocity x′. Frequently the magnitude of this damping force is proportional to the square of the velocity, and the new differential equation becomes
- Write the second-order differential equation as a plane autonomous system, and find all critical points.
- The system is called overdamped when (0, 0) is a stable node and is called underdamped when (0, 0) is a stable spiral point. Physical considerations suggest that (0, 0) must be an asymptotically stable critical point. Show that the system is necessarily underdamped. [Hint: .]
- A bead with mass m slides along a thin wire whose shape may be described by the function z = f(x). Small stretches of the wire act like an inclined plane, and in mechanics it is assumed that the magnitude of the frictional force between the bead and wire is directly proportional to mg cos θ. See Figure 11.4.3.
- Explain why the new differential equation for the x-coordinate of the bead is
for some positive constant µ.
- Investigate the critical points of the corresponding plane autonomous system. Under what conditions is a critical point a saddle point? A stable spiral point?
- Explain why the new differential equation for the x-coordinate of the bead is
- An undamped oscillation satisfies a nonlinear second-order differential equation of the form x″ + f(x) = 0, where f(0) = 0 and xf(x) > 0 for x ≠ 0 and –d < x < d. Use the phase-plane method to investigate whether it is possible for the critical point (0, 0) to be a stable spiral point. [Hint: Let F(x) = f(u)du, and show that y2 + 2F(x) = c.]
- The Lotka–Volterra predator–prey model assumes that, in the absence of predators, the number of prey grows exponentially. If we make the alternative assumption that the prey population grows logistically, the new system is
where a, b, c, r, and K are positive and K > a/b.
- Show that the system has critical points at (0, 0), (0, K), and (, ), where = a/b and c = (K – ).
- Show that the critical points at (0, 0) and (0, K) are saddle points, whereas the critical point at (,) is either a stable node or a stable spiral point.
- Show that (, ) is a stable spiral point if < . Explain why this case will occur when the carrying capacity K of the prey is large.
- The nonlinear system
arises in a model for the growth of microorganisms in a chemostat, a simple laboratory device in which a nutrient from a supply source flows into a growth chamber. In the system, x denotes the concentration of the microorganisms in the growth chamber, y denotes the concentration of nutrients, and α > 1 and β > 0 are constants that can be adjusted by the experimenter. Find conditions on α and β that ensure that the system has a single critical point (, ) in the first quadrant, and investigate the stability of this critical point.
- Use the methods of this chapter together with a numerical solver to investigate stability in the nonlinear spring/mass system modeled by
x″ + 8x – 6x3 + x5 = 0.
See Problem 8 in Exercises 3.11.