--- PDF PAGE 245 ---
230 Chapter 10 Closed Orbits and Limit Sets 2. Consider the three-dimensional system r ′ = r ( 1 − r ) θ ′ = 1 z ′ = − z . Compute the Poincar´ e map along the closed orbit lying on the unit circle given by r = 1 and show that this closed orbit is asymptotically stable. 3. Consider the three-dimensional system r ′ = r ( 1 − r ) θ ′ = 1 z ′ = z . Again compute the Poincar´ e map for this system. What can you now say about the behavior of solutions near the closed orbit. Sketch the phase portrait for this system. 4. Consider the system x ′ = sin x ( − 0.1 cos x − cos y ) y ′ = sin y ( cos x − 0.1 cos y ) . Show that all solutions emanating from the source at (π/ 2, π/ 2 ) have ω -limit sets equal to the square bounded by x = 0, π and y = 0, π . 5. The system r ′ = ar + r 3 − r 5 θ ′ = 1 depends on a parameter a . Determine the phase plane for representative a values and describe all bifurcations for the system. 6. Consider the system x ′ = − y − ( x 4 4 − x 2 2 + y 2 2 ) ( x 3 − x ) y ′ = x 3 − x − ( x 4 4 − x 2 2 + y 2 2 ) y . (a) Find all equilibrium points. (b) Determine the types of these equilibria.

--- PDF PAGE 246 ---
Exercises 231 A Figure 10.12 The region A is positively invariant. (c) Prove that all nonequilibrium solutions have ω -limit sets consisting of either one or two homoclinic solutions plus a saddle point. 7. Let A be an annular region in R 2 . Let F be a planar vector field that points inward along the two boundary curves of A . Suppose also that every radial segment of A is local section. See Figure 10.12. Prove there is a periodic solution in A . 8. Let F be a planar vector field and again consider an annular region A as in the previous problem. Suppose that F has no equilibria and that F points inward along the boundary of the annulus, as before. (a) Prove there is a closed orbit in A . (Notice that the hypothesis is weaker than in the previous problem.) (b) If there are exactly seven closed orbits in A , show that one of them has orbits spiraling toward it from both sides. 9. Let F be a planar vector field on a neighborhood of the annular region A above. Suppose that for every boundary point X of A , F ( X ) is a nonzero vector tangent to the boundary. (a) Sketch the possible phase portraits in A under the further assump- tion that there are no equilibria and no closed orbits besides the boundary circles. Include the case where the solutions on the boundary travel in opposite directions. (b) Suppose the boundary solutions are oppositely oriented and that the flow preserves area. Show that A contains an equilibrium. 10. Show that a closed orbit of a planar system meets a local section in at most one point. 11. Show that a closed and bounded limit set is connected (that is, not the union of two disjoint nonempty closed sets).

--- PDF PAGE 247 ---
232 Chapter 10 Closed Orbits and Limit Sets 12. Let X ′ = F ( X ) be a planar system with no equilibrium points. Suppose the flow φ t generated by F preserves area (that is, if U is any open set, the area of φ t ( U ) is independent of t ). Show that every solution is a closed set. 13. Let γ be a closed orbit of a planar system. Let λ be the period of γ . Let { γ n } be a sequence of closed orbits. Suppose the period of γ n is λ n . If there are points X n ∈ γ n such that X n → X ∈ γ , prove that λ n → λ . (This result can be false for higher dimensional systems. It is true, however, that if λ n → μ , then μ is an integer multiple of λ .) 14. Consider a system in R 2 having only a finite number of equilibria. (a) Show that every limit set is either a closed orbit or the union of equi- librium points and solutions φ t ( X ) such that lim t →∞ φ t ( X ) and lim t →−∞ φ t ( X ) are these equilibria. (b) Show by example (draw a picture) that the number of distinct solutions in ω( X ) may be infinite. 15. Let X be a recurrent point of a planar system; that is, there is a sequence t n → ±∞ such that φ t n ( X ) → X . (a) Prove that either X is an equilibrium or X lies on a closed orbit. (b) Show by example that there can be a recurrent point for a nonplanar system that is not an equilibrium and does not lie on a closed orbit. 16. Let X ′ = F ( X ) and X ′ = G ( X ) be planar systems. Suppose that F ( X ) · G ( X ) = 0 for all X ∈ R 2 . If F has a closed orbit, prove that G has an equilibrium point. 17. Let γ be a closed orbit for a planar system, and suppose that γ forms the boundary of an open set U . Show that γ is not simultaneously the ω - and α -limit set of points of U . Use this fact and the Poincar´ e–Bendixson Theorem to prove that U contains an equilibrium that is not a saddle. ( Hint : Consider the limit sets of points on the stable and unstable curves of saddles.)

--- PDF PAGE 248 ---
11 Applications in Biology In this chapter we make use of the techniques developed in the previous few chapters to examine some nonlinear systems that have been used as mathe- matical models for a variety of biological systems. In Section 11.1 we utilize preceding results involving nullclines and linearization to describe several bio- logical models involving the spread of communicable diseases. In Section 11.2 we investigate the simplest types of equations that model a predator–prey ecology. A more sophisticated approach is used in Section 11.3 to study the populations of a pair of competing species. Instead of developing explicit formulas for these differential equations, we instead make only qualitative assumptions about the form of the equations. We then derive geometric information about the behavior of solutions of such systems based on these assumptions. 11.1 Infectious Diseases The spread of infectious diseases such as measles or malaria may be modeled as a nonlinear system of differential equations. The simplest model of this type is the SIR model. Here we divide a given population into three disjoint groups. The population of susceptible individuals is denoted by S , the infected population by I , and the recovered population by R . As usual, each of these Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00011-7 c © 2013 Elsevier Inc. All rights reserved. 233

--- PDF PAGE 249 ---
234 Chapter 11 Applications in Biology are functions of time. We assume for simplicity that the total population is constant, so that ( S + I + R ) ′ = 0. In the most basic case we make the assumption that, once an individual has been infected and subsequently has recovered, that individual cannot be reinfected. This is the situation that occurs for such diseases as measles, mumps, and smallpox, among many others. We also assume that the rate of transmission of the disease is proportional to the number of encounters between susceptible and infected individuals. The easiest way to character- ize this assumption mathematically is to put S ′ = − β SI for some constant β > 0. We finally assume that the rate at which infected individuals recover is proportional to the number of infected. The SIR model is then S ′ = − β SI I ′ = β SI − ν I R ′ = ν I , where β and ν are positive parameters. As stipulated, we have ( S + I + R ) ′ = 0, so that S + I + R is a constant. This simplifies the system, for if we determine S ( t ) and I ( t ) , we then derive R ( t ) for free. Thus it suffices to consider the two-dimensional system S ′ = − β SI I ′ = β SI − ν I . The equilibria for this system are given by the S -axis ( I = 0). Linearization at ( S , 0 ) yields the matrix ( 0 − β S 0 β S − ν ) , so the eigenvalues are 0 and β S − ν . This second eigenvalue is negative if 0 < S < ν/β and positive if S > ν/β . The S -nullclines are given by the S - and I -axes. On the I -axis, we have I ′ = − ν I , so solutions simply tend to the origin along this line. The I -nullclines are I = 0 and the vertical line S = ν/β . Thus we have the nullcline diagram as shown in Figure 11.1. From this it appears that, given any initial population ( S 0 , I 0 ) with S 0 > ν/β and I 0 > 0, the susceptible population decreases mono- tonically, while the infected population at first rises, but eventually reaches a maximum and then declines to 0. We can actually prove this analytically, for we can explicitly compute a func- tion that is constant along solution curves. Note that the slope of the vector

--- PDF PAGE 250 ---
11.1 Infectious Diseases 235 S I S = ν / β Figure 11.1 Nullclines and direction field for the SIR model. field is a function of S alone: I ′ S ′ = β SI − ν I − β SI = − 1 + ν β S . Thus we have dI dS = dI / dt dS / dt = − 1 + ν β S , which we may immediately integrate to find I = I ( S ) = − S + ν β log S + constant. Therefore, the function I + S − (ν/β) log S is constant along solution curves. It then follows that there is a unique solution curve connecting each equi- librium point in the interval ν/β < S < ∞ to an equilibrium point in the interval 0 < S < ν/β as shown in Figure 11.2. A slightly more complicated model for infectious diseases arises when we assume that recovered individuals may lose their immunity and become rein- fected with the disease. Examples of this type of disease include malaria and tuberculosis. We assume that reinfection occurs at a rate proportional to the population of recovered individuals. This leads to the SIRS model (the extra S indicating that recovered individuals may reenter the susceptible group). The system becomes S ′ = − β SI + μ R I ′ = β SI − ν I R ′ = ν I − μ R .

--- PDF PAGE 251 ---
236 Chapter 11 Applications in Biology I S S = ν / β Figure 11.2 Phase portrait for the SIR system. Again we see that the total population S + I + R is a constant, which we denote by τ . We may eliminate R from this system by setting R = τ − S − I : S ′ = − β SI + μ(τ − S − I ) I ′ = β SI − ν I . Here β , μ , ν , and τ are all positive parameters. Unlike the SIR model, we now have at most two equilibria, one at (τ , 0 ) and the other at ( S ∗ , I ∗ ) = ( ν β , μ(τ − ν β ) ν + μ ) . The first equilibrium point corresponds to no disease whatsoever in the pop- ulation. The second equilibrium point only exists when τ ≥ ν/β . When τ = ν/β , we have a bifurcation as the two equilibria coalesce at (τ , 0 ) . The quantity ν/β is called the threshold level for the disease. The linearized system is given by Y ′ =   − β I − μ − β S − μ β I β S − ν   Y . At the equilibrium point (τ , 0 ) , the eigenvalues are − μ and βτ − ν , so this equilibrium point is a saddle provided that the total population exceeds the threshold level. At the second equilibrium point, a straightforward computa- tion shows that the trace of the matrix is negative, while the determinant is positive. It then follows from the results in Chapter 4 that both eigenvalues have negative real parts, and so this equilibrium point is asymptotically stable.

--- PDF PAGE 252 ---
11.2 Predator–Prey Systems 237 Biologically, this means that the disease may become established in the com- munity only when the total population exceeds the threshold level. We will only consider this case in what follows. Note that the SIRS system is only of interest in the region given by I , S ≥ 0 and S + I ≤ τ . Denote this triangular region by 1 (of course!). Note that the I -axis is no longer invariant, while on the S -axis, solutions increase up to the equilibrium at (τ , 0 ) . Proposition. The region 1 is positively invariant. Proof: We check the direction of the vector field along the boundary of 1 . The field is tangent to the boundary along the lower edge I = 0 as well as at ( 0, τ ) . Along S = 0 we have S ′ = μ(τ − I ) > 0, so the vector field points inward for 0 < I < τ . Along the hypotenuse, if 0 < S ≤ ν/β , we have S ′ = − β SI < 0 and I ′ = I (β S − ν) ≤ 0, so the vector field points inward. When ν/β < S < τ , we have − 1 < I ′ S ′ = − 1 + ν β S ≤ 0, so again the vector field points inward. This completes the proof.  The I -nullclines are given as in the SIR model by I = 0 and S = ν/β . The S -nullcline is given by the graph of the function I = I ( S ) = μ(τ − S ) β S + μ . A calculus student will compute that I ′ ( S ) < 0 and I ′′ ( S ) > 0 when 0 ≤ S < τ . So this nullcline is the graph of a decreasing and concave up function that passes through both (τ , 0 ) and ( 0, τ ) , as shown in Figure 11.3. Note that in this phase portrait, all solutions appear to tend to the equilibrium point ( S ∗ , I ∗ ) ; the proportion of infected to susceptible individuals tends to a “steady state.” To prove this, however, one would need to eliminate the possibility of closed orbits encircling the equilibrium point for a given set of parameters β , μ , ν , and τ . 11.2 Predator–Prey Systems We next consider a pair of species, one of which consists of predators with a population that is denoted by y and the other its prey with population x . We assume that the prey population is the total food supply for the predators. We also assume that, in the absence of predators, the prey population grows

--- PDF PAGE 253 ---
238 Chapter 11 Applications in Biology I S + I = τ S Figure 11.3 Nullclines and phase portrait in 1 for the SIRS system. Here β = ν = μ = 1 and τ = 2. at a rate proportional to the current population. That is, as in Chapter 1, when y = 0 we have x ′ = ax where a > 0. So in this case x ( t ) = x 0 exp ( at ) . When predators are present, we assume that the prey population decreases at a rate proportional to the number of predator–prey encounters. As in the pre- vious section, one simple model for this is bxy where b > 0. So the differential equation for the prey population is x ′ = ax − bxy . For the predator population, we make more or less the opposite assump- tions. In the absence of prey, the predator population declines at a rate proportional to the current population. So when x = 0 we have y ′ = − cy with c > 0, and thus y ( t ) = y 0 exp ( − ct ) . The predator species becomes extinct in this case. When there are prey in the environment, we assume that the predator population increases at a rate proportional to the predator–prey meetings, or dxy . We do not at this stage assume anything about overcrow- ing. Thus our simplified predator–prey system (also called the Volterra–Lotka system) is, x ′ = ax − bxy = x ( a − by ) y ′ = − cy + dxy = y ( − c + dx ) , where the parameters a , b , c , and d are all assumed to be positive. Since we are dealing with populations, we only consider x , y ≥ 0. As usual, our first job is to locate the equilibrium points. These occur at the origin and at ( x , y ) = ( c / d , a / b ) . The linearized system is X ′ = ( a − by − bx dy − c + dx ) X ,

--- PDF PAGE 254 ---
11.2 Predator–Prey Systems 239 x y y = a/b x = c/d Figure 11.4 Nullclines and direction field for the predator–prey system. so when x = y = 0 we have a saddle with eigenvalues a and − c . We know the stable and unstable curves: They are the y - and x -axes, respectively. At the other equilibrium point ( c / d , a / b ) , the eigenvalues are pure imaginary ± i √ ac , and so we cannot conclude anything at this stage about the stability of this equilibrium point. We next sketch the nullclines for this system. The x -nullclines are given by the straight lines x = 0 and y = a / b while the y -nullclines are y = 0 and x = c / d . The nonzero nullcline lines separate the region x , y > 0 into four basic regions in which the vector field points are as indicated in Figure 11.4. Thus the solutions wind in the counterclockwise direction about the equilibrium point. From this we cannot determine the precise behavior of solutions: They could possibly spiral in toward the equilibrium point, spiral toward a limit cycle, spiral out toward “infinity” and the coordinate axes, or else lie on closed orbits. To make this determination, we search for a Liapunov function L . Employing the trick of separation of variables , we look for a function of the form L ( x , y ) = F ( x ) + G ( y ) . Recall that ˙ L denotes the time derivative of L along solutions. We com- pute ˙ L ( x , y ) = d dt L ( x ( t ) , y ( t )) = dF dx x ′ + dG dy y ′ . Thus ˙ L ( x , y ) = x dF dx ( a − by ) + y dG dy ( − c + dx ) .

--- PDF PAGE 255 ---
240 Chapter 11 Applications in Biology We obtain ˙ L ≡ 0 provided x dF / dx dx − c ≡ y dG / dy by − a . Since x and y are independent variables, this is possible if and only if x dF / dx dx − c = y dG / dy by − a = constant. Setting the constant equal to 1, we obtain dF dx = d − c x , dG dy = b − a y . Integrating, we find F ( x ) = dx − c log x , G ( y ) = by − a log y . Thus the function L ( x , y ) = dx − c log x + by − a log y is constant on solution curves of the system when x , y > 0. By considering the signs of ∂ L /∂ x and ∂ L /∂ y , it is easy to see that the equilibrium point Z = ( c / d , a / b ) is an absolute minimum for L . It follows that L (or, more precisely, L − L ( Z )) is a Liapunov function for the system. Therefore, Z is a stable equilibrium. We note next that there are no limit cycles; this follows from corollary 5 in Chapter 10, Section 10.6, because L is not constant on any open set. We now prove the following. Theorem. Every solution of the predator–prey system is a closed orbit (except the equilibrium point Z and the coordinate axes). Proof: Consider the solution through W 6 = Z , where W does not lie on the x - or y -axis. This solution spirals around Z , crossing each nullcline infinitely often. Thus there is a doubly infinite sequence . . . < t − 1 < t 0 < t 1 < . . . such that φ t n ( W ) is on the line x = c / d , and t n → ±∞ as n → ±∞ . If W is not on a closed orbit, the points φ t n ( W ) are monotone along the line x = c / d , as discussed in the previous chapter. Since there are no limit cycles, either

--- PDF PAGE 256 ---
11.2 Predator–Prey Systems 241 x = c/d Figure 11.5 Nullclines and phase portrait for the Volterra–Lotka system. φ t n ( W ) → Z as n → ∞ or φ t n ( W ) → Z as n → −∞ . Since L is constant along the solution through W , this implies that L ( W ) = L ( Z ) . But this contradicts minimality of L ( Z ) . This completes the proof.  The phase portrait for this predator–prey system is shown in Figure 11.5. We conclude that, for any given initial populations ( x ( 0 ) , y ( 0 )) with x ( 0 ) 6 = 0, and y ( 0 ) 6 = 0, other than Z , the populations of predator and prey oscillate cyclically. No matter what the populations of prey and predator are, neither species will die out, nor will its population grow indefinitely. Now let us introduce overcrowding into the prey equation. As in the logistic model in Chapter 1, the equations for prey, in the absence of predators, may be written in the form x ′ = ax − λ x 2 . We also assume that the predator population obeys a similar equation, y ′ = − cy − μ y 2 , when x = 0. Incorporating the preceding assumptions yields the predator–prey equations for species with limited growth: x ′ = x ( a − by − λ x ) y ′ = y ( − c + dx − μ y ) . As before, the parameters a , b , c , d as well as λ and μ are all positive. When y = 0, we have the logistic equation x ′ = x ( a − λ x ) , which yields equilibria at

--- PDF PAGE 257 ---
242 Chapter 11 Applications in Biology the origin and at ( a /λ , 0 ) . As we saw in Chapter 1, all nonzero solutions on the x -axis tend to a /λ . When x = 0, the equation for y is y ′ = − cy − μ y 2 . Since y ′ < 0 when y > 0, it follows that all solutions on this axis tend to the origin. Thus we confine attention to the upper-right quadrant Q where x , y > 0. The nullclines are given by the x - and y -axes, together with the lines L : a − by − λ x = 0 M : − c + dx − μ y = 0. Along the lines L and M , we have x ′ = 0 and y ′ = 0, respectively. There are two possibilities, according to whether these lines intersect in Q or not. We first consider the case where the two lines do not meet in Q . In this case we have the nullcline configuration shown in Figure 11.6. All solutions to the right of M head upward and to the left until they meet M ; between the lines L and M solutions now head downward and to the left. Thus they either meet L or tend directly to the equilibrium point at ( a /λ , 0 ) . If solutions cross L , they then head right and downward, but they cannot cross L again. Thus they too tend to ( a /λ , 0 ) . All solutions in Q therefore tend to this equilibrium point. We conclude that, in this case, the predator population becomes extinct and the prey population approaches its limiting value a /λ . We may interpret the behavior of solutions near the nullclines as follows. Since both x ′ and y ′ are never both positive, it is impossible for both prey and predators to increase at the same time. If the prey population is above its limiting value, it must decrease. After a while the lack of prey causes the M L Figure 11.6 Nullclines and phase portrait for a predator–prey system with limited growth when the nullclines do not meet in Q .

--- PDF PAGE 258 ---
11.2 Predator–Prey Systems 243 predator population to begin to decrease (when the solution crosses M ). After that point the prey population can never increase past a /λ , and so the predator population continues to decrease. If the solution crosses L , the prey popula- tion increases again (but not past a /λ ), while the predators continue to die off. In the limit the predators disappear and the prey population stabilizes at a /λ . Suppose now that L and M cross at a point Z = ( x 0 , y 0 ) in the quadrant Q ; of course, Z is an equilibrium. The linearization of the vector field at Z is X ′ = ( − λ x 0 − bx 0 dy 0 − μ y 0 ) X . The characteristic polynomial has trace given by − λ x 0 − μ y 0 < 0 and deter- minant ( bd + λμ) x 0 y 0 > 0. From the trace–determinant plane of Chapter 4, we see that Z has eigenvalues that are either both negative or both complex with negative real parts. Thus Z is asymptotically stable. Note that, in addition to the equilibria at Z and ( 0, 0 ) , there is still an equi- librium at ( a /λ , 0 ) . Linearization shows that this equilibrium is a saddle; its stable curve lies on the x -axis. See Figure 11.7. It is not easy to determine the basin of Z ; nor do we know whether there are any limit cycles. Nevertheless, we can obtain some information. The line L meets the x -axis at ( a /λ , 0 ) and the y -axis at ( 0, a / b ) . Let 0 be a rectangle with corners that are ( 0, 0 ) , ( p , 0 ) , ( 0, q ) , and ( p , q ) with p > a /λ , q > a / b , and the point ( p , q ) lying in M . Every solution at a boundary point of 0 either enters 0 or is part of the boundary. Therefore, 0 is positively invariant. Every point in Q is contained in such a rectangle. q L Γ M p Figure 11.7 Nullclines and phase portrait for a predator–prey system with limited growth when the nullclines do meet in Q .

--- PDF PAGE 259 ---
244 Chapter 11 Applications in Biology By the Poincar´ e–Bendixson Theorem the ω -limit set of any point ( x , y ) in 0 , with x , y > 0, must be a limit cycle or contain one of the three equilibria ( 0, 0 ) Z , or ( a /λ , 0 ) . We rule out ( 0, 0 ) and ( a /λ , 0 ) by noting that these equilibria are saddles with stable curves that lie on the x - or y -axes. Therefore, ω( x , y ) is either Z or a limit cycle in 0 . By Corollary 4 of the Poincar´ e–Bendixson Theorem any limit cycle must surround Z . We observe further that any such rectangle 0 contains all limit cycles, for a limit cycle (like any solution) must enter 0 , and 0 is positively invari- ant. Fixing ( p , q ) as before, it follows that for any initial values ( x ( 0 ) , y ( 0 )) , there exists t 0 > 0 such that x ( t ) < p , y ( t ) < q if t ≥ t 0 . We conclude that in the long run, a solution either approaches Z or else spirals down to a limit cycle. From a practical standpoint, a solution that tends toward Z is indistinguish- able from Z after a certain time. Likewise, a solution that approaches a limit cycle γ can be identified with γ after it is sufficiently close. We conclude that any population of predators and prey that obeys these equations eventually settles down to either a constant or periodic population. Furthermore, there are absolute upper bounds that no population can exceed in the long run, no matter what the initial populations are. 11.3 Competitive Species We consider now two species that compete for a common food supply. Instead of analyzing specific equations, we follow a different procedure: We consider a large class of equations about which we assume only a few qualitative features. In this way considerable generality is gained, and little is lost, because specific equations can be very difficult to analyze. Let x and y denote the populations of the two species. The equations of growth of the two populations may be written in the form x ′ = M ( x , y ) x y ′ = N ( x , y ) y , where the growth rates M and N are functions of both variables. As usual, we assume that x and y are nonnegative. Thus the x -nullclines are given by x = 0 and M ( x , y ) = 0 and the y -nullclines are y = 0 and N ( x , y ) = 0. We make the following assumptions on M and N : 1. Because the species compete for the same resources, if the population of either species increases, then the growth rate of the other goes down.

--- PDF PAGE 260 ---
11.3 Competitive Species 245 Thus ∂ M ∂ y < 0 and ∂ N ∂ x < 0. 2. If either population is very large, both populations decrease. Thus there exists K > 0 such that M ( x , y ) < 0 and N ( x , y ) < 0 if x ≥ K or y ≥ K . 3. In the absence of either species, the other has a positive growth rate up to a certain population and a negative growth rate beyond it. Therefore, there are constants a , b > 0 such that M ( x , 0 ) > 0 for x < a and M ( x , 0 ) < 0 for x > a , N ( 0, y ) > 0 for y < b and N ( 0, y ) < 0 for y > b . By conditions (1) and (3) each vertical line { x } × R meets the set μ = M − 1 ( 0 ) exactly once if 0 ≤ x ≤ a and not at all if x > a . By (1) and the Implicit Function Theorem, μ is the graph of a nonnegative function f : [0, a ] → R such that f − 1 ( 0 ) = a . Below the curve μ , M is positive and above it, M is negative. In the same way, the set ν = N − 1 ( 0 ) is a smooth curve of the form { ( x , y ) | x = g ( y ) } , where g : [0, b ] → R is a nonnegative function with g − 1 ( 0 ) = b . The function N is positive to the left of ν and negative to the right. Suppose first that μ and ν do not intersect and that μ is below ν . Then the phase portrait can be determined immediately from the nullclines. The equilibria are ( 0, 0 ) , ( a , 0 ) , and ( 0, b ) . The origin is a source, while ( a , 0 ) is a saddle (assuming that (∂ M /∂ x )( a , 0 ) < 0). The equilibrium at ( 0, b ) is a sink (again assuming that (∂ N /∂ y )( 0, b ) < 0). All solutions with y 0 > 0 tend to the asymptotically stable equilibrium ( 0, b ) with the exception of solutions on the x -axis. See Figure 11.8. In case μ lies above ν , the situation is reversed and all solutions with x 0 > 0 tend to the sink that now appears at ( a , 0 ) . Suppose now that μ and ν intersect. We make the assumption that μ ∩ ν is a finite set and, at each intersection point, μ and ν cross transversely ; that is, they have distinct tangent lines at the intersection points. This assumption may be eliminated; we make it only to simplify the process of determining the flow.

--- PDF PAGE 261 ---
246 Chapter 11 Applications in Biology a b ν μ Figure 11.8 Phase portrait when μ and ν do not meet. The nullclines μ and ν and the coordinate axes bound a finite number of connected open sets in the upper right quadrant: These are the basic regions where x ′ 6 = 0 and y ′ 6 = 0. They are of four types: A : x ′ > 0, y ′ > 0, B : x ′ < 0, y ′ > 0, C : x ′ < 0, y ′ < 0, D : x ′ > 0, y ′ < 0. Equivalently, these are the regions where the vector field points northeast, northwest, southwest, or southeast, respectively. Some of these regions are indicated in Figure 11.9. The boundary ∂ R of a basic region R is made up of points of the following types: points of μ ∩ ν , called vertices ; points on μ or ν but not on both or on the coordinate axes, called ordinary boundary points; and points on the axes. A vertex is an equilibrium; the other equilibria lie on the axes at ( 0, 0 ) , ( a , 0 ) , and ( 0, b ) . At an ordinary boundary point Z ∈ ∂ R , the vector field is either vertical (if Z ∈ μ ) or horizontal (if Z ∈ ν ). This vector points either into or out of R since μ has no vertical tangents and ν has no horizontal tangents. We call Z an inward or outward point of ∂ R , accordingly. Note that, in Figure 11.9, the vector field either points inward at all ordinary points on the boundary of a basic region, or else points outward at all such points. This is no accident, for we have this proposition. Proposition. Let R be a basic region for the competitive species model. Then the ordinary boundary points of R are either all inward or all outward. Proof: There are only two ways that the curves μ and ν may intersect at a vertex P . As y increases along ν , the curve ν may either pass from below μ to above μ or from above to below μ . These two scenarios are illustrated

--- PDF PAGE 262 ---
11.3 Competitive Species 247 μ ν C D S R B Q b P a A Figure 11.9 Basic regions when the nullclines μ and ν intersect. μ ν ν (a) (b) μ Figure 11.10 In (a), ν passes from below μ to above μ as y increases. The situation is reversed in (b). in Figures 11.10(a) and (b). There are no other possibilities since we have assumed that these curves cross transversely. Since x ′ > 0 below μ and x ′ < 0 above μ , and since y ′ > 0 to the left of ν and y ′ < 0 to the right, we have the following configurations for the vector field in these two cases. See Figure 11.11. In each case we see that the vector field points inward in two opposite basic regions abutting P and outward in the other two basic regions. If we now move along μ or ν to the next vertex along this curve, we see that adjacent basic regions must maintain their inward or outward configuration. Therefore, at all ordinary boundary points on each basic region, the vector field either points outward or points inward, as required.  As a consequence of the proposition, it follows that each basic region and its closure is either positively or negatively invariant. What are the possible

--- PDF PAGE 263 ---
248 Chapter 11 Applications in Biology μ ν μ ν Figure 11.11 Configurations of the vector field near vertices. Γ μ ν Figure 11.12 All solutions must enter and then remain in 0 . ω -limit points of this system? There are no closed orbits because a closed orbit must be contained in a basic region, but this is impossible since x ( t ) and y ( t ) are monotone along any solution curve in a basic region. Therefore all ω -limit points are equilibria. We note also that each solution is defined for all t ≥ 0 because any point lies in a large rectangle 0 with corners at ( 0, 0 ) , ( x 0 , 0 ) , ( 0, y 0 ) , and ( x 0 , y 0 ) with x 0 > a and y 0 > b ; such a rectangle is positively invariant. See Figure 11.12. Thus we have shown the following. Theorem. The flow φ t of the competitive species system has the following property: for all points ( x , y ) , with x ≥ 0, y ≥ 0 , the limit lim t →∞ φ t ( x , y ) exists and is one of a finite number of equilibria.  We conclude that the populations of two competing species always tend to one of a finite number of limiting populations.

--- PDF PAGE 264 ---
11.3 Competitive Species 249 P μ ν Figure 11.13 This configuration of μ and ν leads to an asymptotically stable equilibrium point. Examining the equilibria for stability, one finds the following results. A ver- tex where μ and ν each have negative slope but μ is steeper is asymptotically stable. See Figure 11.13. One sees this by drawing a small rectangle with sides parallel to the axes around the equilibrium, putting one corner in each of the four adjacent basic regions. Such a rectangle is positively invariant; since it can be arbitrarily small, the equilibrium is asymptotically stable. This may also be seen as follows. We have slope of μ = − M x M y < slope of ν = − N x N y < 0, where M x = ∂ M /∂ x , M y = ∂ M /∂ y , and so on, at the equilibrium. Now recall that M y < 0 and N x < 0. Therefore, at the equilibrium point, we also have M x < 0 and N y < 0. Linearization at the equilibrium point yields the matrix ( xM x xM y yN x yN y ) . The trace of this matrix is xM x + yN y < 0 while the determinant is xy ( M x N y − M y N x ) > 0. Thus the eigenvalues have negative real parts, and so we have a sink. A case-by-case study of the different ways μ and ν can cross shows that the only other asymptotically stable equilibrium in this model is ( 0, b ) when ( 0, b ) is above μ , or ( a , 0 ) when ( a , 0 ) is to the right of ν . All other equilibria are unstable. There must be at least one asymptotically stable equilibrium. If ( 0, b ) is not one, then it lies under μ ; and if ( a , 0 ) is not one, it lies over μ . In that case μ and ν cross, and the first crossing to the left of ( a , 0 ) is asymptotically stable. For example, this analysis tells us that, in Figure 11.14, only P and ( 0, b ) are asymptotically stable; all other equilibria are unstable. In particular, assuming

--- PDF PAGE 265 ---
250 Chapter 11 Applications in Biology Z R P S Q b B ∞ B ∞ μ ν Figure 11.14 Note that solutions on either side of the point Z in the stable curve of Q have very different fates. that the equilibrium Q in Figure 11.14 is hyperbolic, then it must be a saddle because certain nearby solutions tend toward it, while others tend away. The point Z lies on one branch of the stable curve through Q . All points in the region denoted B ∞ to the left of Z tend to the equilibrium at ( 0, b ) , while points to the right go to P . Thus, as we move across the branch of the stable curve containing Z , the limiting behavior of solutions changes radically. Since solutions just to the right of Z tend to the equilibrium point P , it follows that the populations in this case tend to stabilize. On the other hand, just to the left of Z , solutions tend to an equilibrium point where x = 0. Thus, in this case, one of the species becomes extinct. A small change in initial conditions has led to a dramatic change in the fate of populations. Ecologically, this small change could have been caused by the introduction of a new pesticide, the importation of additional members of one of the species, a forest fire, or the like. Mathematically, this event is a jump from the basin of P to that of ( 0, b ) . 11.4 Exploration: Competition and Harvesting In this exploration we investigate the competitive species model where we allow either harvesting (emigration) or immigration of one of the species. We

--- PDF PAGE 266 ---
11.5 Exploration: Adding Zombies to the SIR Model 251 consider the system x ′ = x ( 1 − ax − y ) y ′ = y ( b − x − y ) + h . Here a , b , and h are parameters. We assume that a , b > 0. If h < 0, then we are harvesting species y at a constant rate, whereas if h > 0, we add to the popula- tion y at a constant rate. The goal is to understand this system completely for all possible values of these parameters. As usual, we only consider the regime where x , y ≥ 0. If y ( t ) < 0 for any t > 0, then we consider this species to have become extinct. 1. First assume that h = 0. Give a complete synopsis of the behavior of this system by plotting the different behaviors you find in the ab -parameter plane. 2. Identify the points or curves in the ab -plane where bifurcations occur when h = 0 and describe them. 3. Now let h < 0. Describe the ab -parameter plane for various (fixed) h - values. 4. Repeat the previous exploration for h > 0. 5. Describe the full three-dimensional parameter space using pictures, flip- books, 3D models, movies, or whatever you find most appropriate. 11.5 Exploration: Adding Zombies to the SIR Model The recent uptick in the number of zombie movies suggests that we might revise the SIR model so that the infected population now consists of zombies. In this situation the zombie population (again denoted by I ) is much more active in infecting the susceptible population. One way to model this is as follows. Assume that S , I , and R denote the fraction of the total population that is susceptible, infected, and recovered. So S + I + R = 1. Then √ I is much larger than I when I is small, so the zombies are much more likely to cause a problem. The new equations are S ′ = − β S √ I I ′ = β S √ I − ν I . Another possible scenario is that zombies continue to infect susceptibles until they are destroyed by a susceptible. This would mean that the infected

--- PDF PAGE 267 ---
252 Chapter 11 Applications in Biology population would now decrease at a rate that is proportional to the number of susceptibles present, or S ′ = − β SI I ′ = β SI − γ S . 1. In each case, describe the behavior of solutions of the corresponding system. 2. What happens if we revise the SIRS model so that either of the preceding assumptions hold? E X E R C I S E S 1. For the SIRS model, prove that all solutions in the triangular region 1 tend to the equilibrium point (τ , 0 ) when the total population does not exceed the threshold level for the disease. 2. Sketch the phase plane for the following variant of the predator–prey system x ′ = x ( 1 − x ) − xy y ′ = y ( 1 − y x ) . 3. A modification of the predator–prey equations is given by x ′ = x ( 1 − x ) − axy x + 1 y ′ = y ( 1 − y ) , where a > 0 is a parameter. (a) Find all equilibrium points and classify them. (b) Sketch the nullclines and the phase portraits for different values of a . (c) Describe any bifurcations that occur as a varies. 4. Another modification of the predator–prey equations is given by x ′ = x ( 1 − x ) − xy x + b y ′ = y ( 1 − y ) , where b > 0 is a parameter.

--- PDF PAGE 268 ---
Exercises 253 (a) Find all equilibrium points and classify them. (b) Sketch the nullclines and the phase portraits for different values of b . (c) Describe any bifurcations that occur as b varies. 5. The equations x ′ = x ( 2 − x − y ) , y ′ = y ( 3 − 2 x − y ) satisfy conditions (1) through (3) in Section 11.3 for competing species. Determine the phase portrait for this system. Explain why these equa- tions make it mathematically possible, but extremely unlikely, for both species to survive. 6. Consider the competing species model x ′ = x ( a − x − ay ) y ′ = y ( b − bx − y ) , where the parameters a and b are positive. (a) Find all equilibrium points for this system and determine their stability type. These types will, of course, depend on a and b . (b) Use the nullclines to determine the various phase portraits that arise for different choices of a and b . (c) Determine the values of a and b for which there is a bifurcation in this system and describe the bifurcation that occurs. (d) Record your findings by drawing a picture of the ab -plane and indi- cating in each open region of this plane the qualitative structure of the corresponding phase portraits. 7. Two species x , y are in symbiosis if an increase of either population leads to an increase in the growth rate of the other. Thus we assume x ′ = M ( x , y ) x y ′ = N ( x , y ) y with ∂ M ∂ y > 0 and ∂ N ∂ x > 0 and x , y ≥ 0. We also suppose that the total food supply is limited; thus

--- PDF PAGE 269 ---
254 Chapter 11 Applications in Biology for some A > 0, B > 0 we have M ( x , y ) < 0 if x > A , N ( x , y ) < 0 if y > B . If both populations are very small, they both increase; thus M ( 0, 0 ) > 0 and N ( 0, 0 ) > 0. Assuming that the intersections of the curves M − 1 ( 0 ) , N − 1 ( 0 ) are finite, and that all are transverse, show that (a) Every solution tends to an equilibrium in the region 0 < x < A , 0 < y < B . (b) There are no sources. (c) There is at least one sink. (d) If ∂ M /∂ x < 0 and ∂ N /∂ y < 0, there is a unique sink Z , and Z is the ω -limit set for all ( x , y ) with x > 0, y > 0. 8. Give a system of differential equations for a pair of mutually destruc- tive species. Then prove that, under plausible hypotheses, two mutually destructive species cannot coexist in the long run. 9. Let y and x denote predator and prey populations. Let x ′ = M ( x , y ) x , y ′ = N ( x , y ) y , where M and N satisfy the following conditions. (a) If there are not enough prey, the predators decrease. Thus, for some b > 0, N ( x , y ) < 0 if x < b . (b) An increase in the prey improves the predator growth rate; thus ∂ N /∂ x > 0. (c) In the absence of predators a small prey population will increase; thus M ( 0, 0 ) > 0. (d) Beyond a certain size, the prey population must decrease; thus there exists A > 0 with M ( x , y ) < 0 if x > A . (e) Any increase in predators decreases the rate of growth of prey; thus ∂ M /∂ y < 0. (f) The two curves M − 1 ( 0 ) , N − 1 ( 0 ) intersect transversely, and at only a finite number of points.

--- PDF PAGE 270 ---
Exercises 255 Show that if there is some ( u , v ) with M ( u , v ) > 0 and N ( u , v ) > 0, then there is either an asymptotically stable equilibrium or an ω -limit cycle. Moreover, show that, if the number of limit cycles is finite and positive, one of them must have orbits spiraling toward it from both sides. 10. Consider the following modification of the predator–prey equations: x ′ = x ( 1 − x ) − axy x + c y ′ = by ( 1 − y x ) , where a , b , and c are positive constants. Determine the region in the parameter space for which this system has a stable equilibrium with both x , y 6 = 0. Prove that, if the equilibrium point is unstable, this system has a stable limit cycle.

--- PDF PAGE 271 ---
This page intentionally left blank

--- PDF PAGE 272 ---
12 Applications in Circuit Theory In this chapter we first present a simple but very basic example of an electri- cal circuit and then derive the differential equations governing this circuit. Certain special cases of these equations are analyzed using the techniques developed in Chapters 8, 9, and 10 in the next two sections; these are the classical equations of Li´ enard and van der Pol. In particular, the van der Pol equation could perhaps be regarded as one of the fundamental examples of a nonlinear ordinary differential equation. It possesses an oscillation or peri- odic solution that is a periodic attractor. Every nontrivial solution tends to this periodic solution; no linear system has this property. Whereas asymptoti- cally stable equilibria sometimes imply death in a system, attracting oscillators imply life. We give an example in Section 12.4 of a continuous transition from one such situation to the other. 12.1 An RLC Circuit In this section, we present our first example of an electrical circuit. This cir- cuit is the simple but fundamental series RLC circuit shown in Figure 12.1. We begin by explaining what this diagram means in mathematical terms. The circuit has three branches : one resistor marked by R , one inductor marked by Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00012-9 c © 2013 Elsevier Inc. All rights reserved. 257

--- PDF PAGE 273 ---
258 Chapter 12 Applications in Circuit Theory R L C α γ β Figure 12.1 RLC circuit. L , and one capacitor marked by C . We think of a branch of this circuit as a certain electrical device with two terminals. For example, in this circuit, the branch R has terminals α and β and all of the terminals are wired together to form the points or nodes α , β , and γ . In the circuit there is a current flowing through each branch that is mea- sured by a real number. More precisely, the currents in the circuit are given by the three real numbers i R , i L , and i C , where i R measures the current through the resistor, and so on. Current in a branch is analogous to water flowing in a pipe; the corresponding measure for water would be the amount flowing in unit time or, better, the rate at which water passes by a fixed point in the pipe. The arrows in the diagram that orient the branches tell us how to read which way the current (read water!) is flowing; for example, if i R is positive, then according to the arrow, current flows through the resistor from β to α (the choice of the arrows is made once and for all at the start). The state of the currents at a given time in the circuit is thus represented a point i = ( i R , i L , i C ) ∈ R 3 . But Kirchhoff ’s current law (KCL) says that in reality there is a strong restriction on which i can occur. KCL asserts that the total current flowing into a node is equal to the total current flowing out of that node. (Think of the water analogy to make this plausible.) For our circuit this is equivalent to KCL : i R = i L = − i C . This defines the one-dimensional subspace K 1 of R 3 of physical current states . Our choice of orientation of the capacitor branch may seem unnatural. In fact, these orientations are arbitrary; in this example they were chosen so that the equations eventually obtained relate most directly to the history of the subject. The state of the circuit is characterized by the current i = ( i R , i L , i C ) together with the voltage (or, more precisely, the voltage drop) across each branch. These voltages are denoted by v R , v L , and v C for the resistor branch, inductor

--- PDF PAGE 274 ---
12.1 An RLC Circuit 259 branch, and capacitor branch, respectively. In the water analogy one thinks of the voltage drop as the difference in pressures at the two ends of a pipe. To measure voltage one places a voltmeter (imagine a water pressure meter) at each of the nodes α , β , and γ that reads V (α) at α , and so on. Then v R is the difference in the reading at α and β : V (β) − V (α) = v R . The direction of the arrow tells us that v R = V (β) − V (α) rather than V (α) − V (β) . An unrestricted voltage state of the circuit is then a point v = ( v R , v L , v C ) ∈ R 3 . This time the Kirchhoff voltage law (KVL) puts a physical restriction on v : KVL : v R + v L − v C = 0. This defines a two-dimensional subspace K 2 of R 3 . KVL follows immediately from our definition of the v R , v L , and v C in terms of voltmeters; that is, v R + v L − v C = ( V (β) − V (α)) + ( V (α) − V (γ )) − ( V (β) − V (γ )) = 0. The product space R 3 × R 3 is called the state space for the circuit. Those states ( i , v ) ∈ R 3 × R 3 satisfying Kirchhoff ’s laws form a three-dimensional subspace of the state space. Now we give a mathematical definition of the three kinds of electrical devices in the circuit. First consider the resistor element. A resistor in the R branch imposes a “functional relationship” on i R and v R . In our exam- ple we take this relationship to be defined by a function f : R → R , so that v R = f ( i R ) . If R is a conventional linear resistor, then f is linear and so f ( i R ) = ki R . This relation is known as Ohm’s law . Nonlinear functions yield a generalized Ohm’s law. The graph of f is called the characteristic of the resistor. A couple of examples of characteristics are given in Figure 12.2. (A characteristic like that in Figure 12.2(b) occurs in a tunnel diode.) A physical state ( i , v ) ∈ R 3 × R 3 is a point that satisfies KCL, KVL, and also f ( i R ) = v R . These conditions define the set of physical states 6 ⊂ R 3 × R 3 . Thus 6 is the set of points ( i R , i L , i C , v R , v L , v C ) in R 3 × R 3 that satisfy 1. i R = i L = − i C (KCL) 2. v R + v L − v C = 0 (KVL) 3. f ( i R ) = v R (generalized Ohm’s law) Now we turn to the differential equations governing the circuit. The induc- tor (which we think of as a coil; it is hard to find a water analogy) specifies

--- PDF PAGE 275 ---
260 Chapter 12 Applications in Circuit Theory (a) (b) ν R iR ν R iR Figure 12.2 Several possible characteristics for a resistor. that L di L ( t ) dt = v L ( t ) ( Faraday’s law ) , where L is a positive constant called the inductance . On the other hand, the capacitor (which may be thought of as two metal plates separated by some insulator; in the water model it is a tank) imposes the condition C dv C ( t ) dt = i C ( t ) , where C is a positive constant called the capacitance . Let’s summarize the development so far: A state of the circuit is given by the six numbers ( i R , i L , i C , v R , v L , v C ) , that is, a point in the state space R 3 × R 3 . These numbers are subject to three restrictions: Kirchhoff ’s current law, Kirchhoff ’s voltage law, and the resistor characteristic or generalized Ohm’s law. Therefore the set of physical states is a certain subset 6 ⊂ R 3 × R 3 . The way a state changes in time is determined by the preceding two differential equations. Next, we simplify the set of physical states 6 by observing that i L and v C determine the other four coordinates. This follows since i R = i L and i C = − i L by KCL, v R = f ( i R ) = f ( i L ) by the generalized Ohm’s law, and v L = v C − v R = v C − f ( i L ) by KVL. Therefore, we can use R 2 as the state space, with coordi- nates given by ( i L , v C ) . Formally, we define a map π : R 3 × R 3 → R 2 , which sends ( i , v ) ∈ R 3 × R 3 to ( i L , v C ) . Then we set π 0 = π | 6 , the restriction of π to 6 . The map π 0 : 6 → R 2 is one to one and onto; its inverse is given by the map ψ : R 2 → 6 , where ψ( i L , v C ) = ( i L , i L , − i L , f ( i L ) , v C − f ( i L ) , v C ) .

--- PDF PAGE 276 ---
12.2 The Li ´ enard Equation 261 It is easy to check that ψ( i L , v C ) satisfies KCL, KVL, and the generalized Ohm’s law, so ψ does map R 2 into 6 . It is also easy to see that π 0 and ψ are inverse to each other. We therefore adopt R 2 as our state space. The differential equations gov- erning the change of state must be rewritten in terms of our new coordinates ( i L , v C ) : L di L dt = v L = v C − f ( i L ) C dv C dt = i C = − i L . For simplicity, and since this is only an example, we set L = C = 1. If we write x = i L and y = v C , we then have a system of differential equations in the plane of the form dx dt = y − f ( x ) dy dt = − x . This is one form of the equation known as the Li´ enard equation . We analyze this system in the following section. 12.2 The Li ´ enard Equation In this section we begin the study of the phase portrait of the Li´ enard system from the circuit of the previous section: dx dt = y − f ( x ) dy dt = − x . In the special case where f ( x ) = x 3 − x , this system is called the van der Pol equation . First consider the simplest case where f is linear. Suppose f ( x ) = kx , where k > 0. Then the Li´ enard system takes the form Y ′ = AY where A = ( − k 1 − 1 0 ) .

--- PDF PAGE 277 ---
262 Chapter 12 Applications in Circuit Theory The eigenvalues of A are given by λ ± = ( − k ± ( k 2 − 4 ) 1 / 2 )/ 2. Since λ ± is either negative or else has a negative real part, the equilibrium point at the origin is a sink. It is a spiral sink if k < 2. For any k > 0, all solutions of the system tend to the origin; physically, this is the dissipative effect of the resistor. Note that we have y ′′ = − x ′ = − y + kx = − y − ky ′ , so that the system is equivalent to the second-order equation y ′′ + ky ′ + y = 0, which is often encountered in elementary differential equations courses. Next we consider the case of a general characteristic f . There is a unique equilibrium point for the Li´ enard system that is given by ( 0, f ( 0 )) . Lineariza- tion yields the matrix ( − f ′ ( 0 ) 1 − 1 0 ) , with eigenvalues given by λ ± = 1 2 ( − f ′ ( 0 ) ± √ ( f ′ ( 0 )) 2 − 4 ) . We conclude that this equilibrium point is a sink if f ′ ( 0 ) > 0 and a source if f ′ ( 0 ) < 0. In particular, for the van der Pol equation where f ( x ) = x 3 − x , the unique equilibrium point is a source. To analyze the system further, we define the function W : R 2 → R 2 by W ( x , y ) = 1 2 ( x 2 + y 2 ) . Then we have ˙ W = x ( y − f ( x )) + y ( − x ) = − xf ( x ) . In particular, if f satisfies f ( x ) > 0 if x > 0, f ( x ) < 0 if x < 0, and f ( 0 ) = 0, then W is a strict Liapunov function on all of R 2 . It follows that, in this case, all solutions tend to the unique equilibrium point lying at the origin. In circuit theory, a resistor is called passive if its characteristic is contained in the set consisting of ( 0, 0 ) and the interior of the first and third quadrant. Therefore, in the case of a passive resistor, − xf ( x ) is negative except when x = 0, and so all solutions tend to the origin. Thus the word passive correctly describes the dynamics of such a circuit.

--- PDF PAGE 278 ---
12.3 The van der Pol Equation 263 12.3 The van der Pol Equation In this section we continue the study of the Li´ enard equation in the special case where f ( x ) = x 3 − x . This is the van der Pol equation : dx dt = y − x 3 + x dy dt = − x . Let φ t denote the flow of this system. In this case we can give a fairly complete phase portrait analysis. Theorem. There is one nontrivial periodic solution of the van der Pol equa- tion, and every other solution (except the equilibrium point at the origin) tends to this periodic solution. “The system oscillates.”  We know from the previous section that this system has a unique equilib- rium point at the origin, and that this equilibrium is a source, since f ′ ( 0 ) < 0. The next step is to show that every nonequilibrium solution “rotates” in a cer- tain sense around the equilibrium in a clockwise direction. To see this, note that the x -nullcline is given by y = x 3 − x and the y -nullcline is the y -axis. We subdivide each of these nullclines into two pieces given by v + = { ( x , y ) | y > 0, x = 0 } v − = { ( x , y ) | y < 0, x = 0 } g + = { ( x , y ) | x > 0, y = x 3 − x } g − = { ( x , y ) | x < 0, y = x 3 − x } . These curves are disjoint; together with the origin they form the boundaries of the four basic regions A , B , C , D shown in Figure 12.3. From the configuration of the vector field in the basic regions, it appears that all nonequilibrium solutions wind around the origin in the clockwise direction. This is indeed the case. Proposition. Solution curves starting on v + cross successively through g + , v − , and g − before returning to v + . Proof: Any solution starting on v + immediately enters the region A since x ′ ( 0 ) > 0. In A we have y ′ < 0, so this solution must decrease in the y -direction. Since the solution cannot tend to the source, it follows that this

--- PDF PAGE 279 ---
264 Chapter 12 Applications in Circuit Theory D C B A ν + ν − g + g − Figure 12.3 Basic regions and nullclines for the van der Pol system. solution must eventually meet g + . On g + we have x ′ = 0 and y ′ < 0. Conse- quently, the solution crosses g + and then enters the region B . Once inside B , the solution heads southwest. Note that the solution cannot reenter A since the vector field points straight downward on g + . There are thus two possibilities: Either the solution crosses v − , or the solution tends to −∞ in the y -direction and never crosses v − . We claim that the latter cannot happen. Suppose that it does. Let ( x 0 , y 0 ) be a point on this solution in region B and consider φ t ( x 0 , y 0 ) = ( x ( t ) , y ( t )) . Since x ( t ) is never 0, it follows that this solution curve lies for all time in the strip S given by 0 < x ≤ x 0 , y ≤ y 0 and we have y ( t ) → −∞ as t → t 0 for some t 0 . We first observe that, in fact, t 0 = ∞ . To see this, note that y ( t ) − y 0 = ∫ t 0 y ′ ( s ) ds = ∫ t 0 − x ( s ) ds . But 0 < x ( s ) ≤ x 0 , so we may only have y ( t ) → −∞ if t → ∞ . Now consider x ( t ) for 0 ≤ t < ∞ . We have x ′ = y − x 3 + x . Since the quan- tity − x 3 + x is bounded in the strip S and y ( t ) → −∞ as t → ∞ , it follows that x ( t ) − x 0 = ∫ t 0 x ′ ( s ) ds → −∞ as t → ∞ as well. But this contradicts our assumption that x ( t ) > 0. Thus this solution must cross v − . Now the vector field is skew-symmetric about the origin. That is, if G ( x , y ) is the van der Pol vector field, then G ( − x , − y ) = − G ( x , y ) . Exploiting this symmetry, it follows that solutions must then pass through the regions C and D in similar fashion.  As a consequence of this result, we may define a Poincar´ e map P on the half-line v + . Given ( 0, y 0 ) ∈ v + , we define P ( y 0 ) to be the y -coordinate of the

--- PDF PAGE 280 ---
12.3 The van der Pol Equation 265 y g + g − P (y) ν + ν − Figure 12.4 The Poincar ´ e map on ν + . first return of φ t ( 0, y 0 ) to v + with t > 0. See Figure 12.4. As in Chapter 10, Section 10.3, P is a C ∞ function that is one to one. The Poincar´ e map is also onto. To see this, simply follow solutions starting on v + backward in time until they reintersect v + , as they must by the proposition. Let P n = P ◦ P n − 1 denote the n -fold composition of P with itself. Our goal now is to prove this theorem. Theorem. The Poincar´ e map has a unique fixed point in v + . Further- more, the sequence P n ( y 0 ) tends to this fixed point as n → ∞ for any nonzero y 0 ∈ v + .  Clearly, any fixed point of P lies on a periodic solution. On the other hand, if P ( y 0 ) 6 = y 0 , then the solution through ( 0, y 0 ) can never be periodic. Indeed, if P ( y 0 ) > y 0 , then the successive intersections of φ t ( 0, y 0 ) with v + is a mono- tone sequence as in Chapter 10, Section 10.4. Thus the solution crosses v + in an increasing sequence of points and so the solution can never meet itself. The case P ( y 0 ) < y 0 is analogous. We may define a “semi-Poincar´ e map” α : v + → v − by letting α( y ) be the y -coordinate of the first point of intersection of φ t ( 0, y ) with v − where t > 0. Also define δ( y ) = 1 2 ( α( y ) 2 − y 2 ) . Note for later use that there is a unique point ( 0, y ∗ ) ∈ v + and time t ∗ such that 1. φ t ( 0, y ∗ ) ∈ A for 0 < t < t ∗ 2. φ t ∗ ( 0, y ∗ ) = ( 1, 0 ) ∈ g + See Figure 12.5.

--- PDF PAGE 281 ---
266 Chapter 12 Applications in Circuit Theory g + (1, 0) y* y α ( y ) ν + ν − Figure 12.5 Semi-Poincar ´ e map. The theorem will now follow directly from the following rather delicate result. Proposition. The function δ( y ) satisfies the following: 1. δ( y ) > 0 if 0 < y < y ∗ 2. δ( y ) decreases monotonically to −∞ as y → ∞ for y > y ∗  We will prove the proposition shortly; first we use it to complete the proof of the theorem. We exploit the fact that the vector field is skew-symmetric about the origin. This implies that if ( x ( t ) , y ( t )) is a solution curve, then so is ( − x ( t ) , − y ( t )) . Proof: Part of the graph of δ( y ) is shown schematically in Figure 12.6. The intermediate value theorem and the proposition imply that there is a unique y 0 ∈ v + with δ( y 0 ) = 0. Consequently, α( y 0 ) = − y 0 and it follows from the skew-symmetry that the solution through ( 0, y 0 ) is periodic. Since δ( y ) 6 = 0 except at y 0 , we have α( y ) 6 = y for all other y -values, and so it follows that φ t ( 0, y 0 ) is the unique periodic solution. We next show that all other solutions (except the equilibrium pont) tend to this periodic solution. Toward that end we define a map β : v − → v + , sending each point of v − to the first intersection of the solution (for t > 0) with v + . By symmetry we have β( y ) = − α( − y ) . Note also that P ( y ) = β ◦ α( y ) . We identify the y -axis with the real numbers in the y -coordinate. Thus if y 1 , y 2 ∈ v + ∪ v − , we write y 1 > y 2 if y 1 is above y 2 . Note that α and β reverse this ordering while P preserves it.

--- PDF PAGE 282 ---
12.3 The van der Pol Equation 267 y* y 0 y δ ( y ) Figure 12.6 Graph of δ ( y ). Now choose y ∈ v + with y > y 0 . Since α( y 0 ) = − y 0 , we have α( y ) < − y 0 and P ( y ) > y 0 . On the other hand, δ( y ) < 0 which implies that α( y ) > − y . Therefore P ( y ) = β(α( y )) < y . We have shown that y > y 0 implies y > P ( y ) > y 0 . Similarly P ( y ) > P ( P ( y )) > y 0 and by induction P n ( y ) > P n + 1 ( y ) > y 0 for all n > 0. The decreasing sequence P n ( y ) has a limit y 1 ≥ y 0 in v + . Note that y 1 is a fixed point of P , for, by continuity of P , we have P ( y 1 ) − y 1 = lim n →∞ P ( P n ( y )) − y 1 = y 1 − y 1 = 0. Since P has only one fixed point, we have y 1 = y 0 . This shows that the solu- tion through y spirals toward the periodic solution as t → ∞ . The same is true if y < y 0 ; the details are left to the reader. Since every solution except the equilibrium meets v + , the proof of the theorem is complete. (See Figure 12.7.)  Figure 12.7 Phase portrait of the van der Pol equation.

--- PDF PAGE 283 ---
268 Chapter 12 Applications in Circuit Theory Finally, we turn to the proof of the proposition. We adopt the following notation. Let γ : [ a , b ] → R 2 be a smooth curve in the plane and let F : R 2 → R . We write γ ( t ) = ( x ( t ) , y ( t )) and define ∫ γ F ( x , y ) = b ∫ a F ( x ( t ) , y ( t )) dt . If it happens that x ′ ( t ) 6 = 0 for a ≤ t ≤ b , then along γ , y is a function of x , so we may write y = y ( x ) . In this case we can change variables: b ∫ a F ( x ( t ) , y ( t )) dt = x ( b ) ∫ x ( a ) F ( x , y ( x )) dt dx dx . Therefore, ∫ γ F ( x , y ) = x ( b ) ∫ x ( a ) F ( x , y ( x )) dx / dt dx . We have a similar expression if y ′ ( t ) 6 = 0. Now recall the function W ( x , y ) = 1 2 ( x 2 + y 2 ) introduced in the previous section. Let p ∈ v + . Suppose α( p ) = φ τ ( p ) . Let γ ( t ) = ( x ( t ) , y ( t )) for 0 ≤ t ≤ τ = τ ( p ) be the solution curve joining p ∈ v + to α( p ) ∈ v − . By definition δ( p ) = 1 2 ( y (τ ) 2 − y ( 0 ) 2 ) = W ( x (τ ) , y (τ )) − W ( x ( 0 ) , y ( 0 )) . Thus δ( p ) = τ ∫ 0 d dt W ( x ( t ) , y ( t )) dt . Recall from Section 12.2 that we have ˙ W = − xf ( x ) = − x ( x 3 − x ) .

--- PDF PAGE 284 ---
12.3 The van der Pol Equation 269 Thus we have δ( p ) = τ ∫ 0 − x ( t )( x ( t ) 3 − x ( t )) dt = τ ∫ 0 x ( t ) 2 ( 1 − x ( t ) 2 ) dt . This immediately proves part (1) of the proposition because the integrand is positive for 0 < x ( t ) < 1. We may rewrite the last equality as δ( p ) = ∫ γ x 2 ( 1 − x 2 ) . We restrict attention to points p ∈ v + with p > y ∗ . We divide the correspond- ing solution curve γ into three curves γ 1 , γ 2 , γ 3 as shown in Figure 12.8. The curves γ 1 and γ 3 are defined for 0 ≤ x ≤ 1, while the curve γ 2 is defined for y 1 ≤ y ≤ y 2 . Then δ( p ) = δ 1 ( p ) + δ 2 ( p ) + δ 3 ( p ) , where δ i ( p ) = ∫ γ i x 2 ( 1 − x 2 ) , i = 1, 2, 3. g + y 1 y 2 y 0 γ 3 γ 2 γ 1 Figure 12.8 Curves γ 1 , γ 2 , and γ 3 shown on the closed orbit through y 0 .

--- PDF PAGE 285 ---
270 Chapter 12 Applications in Circuit Theory Notice that, along γ 1 , y ( t ) may be regarded as a function of x . Thus we have δ 1 ( p ) = 1 ∫ 0 x 2 ( 1 − x 2 ) dx / dt dx = 1 ∫ 0 x 2 ( 1 − x 2 ) y − f ( x ) dx , where f ( x ) = x 3 − x . As p moves up the y -axis, y − f ( x ) increases (for ( x , y ) on γ 1 ). Thus δ 1 ( p ) decreases as p increases. Similarly δ 3 ( p ) decreases as p increases. On γ 2 , x ( t ) may be regarded as a function of y that is defined for y ∈ [ y 1 , y 2 ] and x ≥ 1. Therefore, since dy / dt = − x , we have δ 2 ( p ) = y 1 ∫ y 2 − x ( y ) ( 1 − x ( y ) 2 ) dy = y 2 ∫ y 1 x ( y ) ( 1 − x ( y ) 2 ) dy , so that δ 2 ( p ) is negative. As p increases, the domain [ y 1 , y 2 ] of integration becomes steadily larger. The function y → x ( y ) depends on p , so we write it as x p ( y ) . As p increases, the curves γ 2 move to the right; thus x p ( y ) increases and so x p ( y )( 1 − x p ( y ) 2 ) decreases. It follows that δ 2 ( p ) decreases as p increases, and evidently lim p →∞ δ 2 ( p ) = −∞ . Consequently, δ( p ) also decreases and tends to −∞ as p → ∞ . This completes the proof of the proposition. 12.4 A Hopf Bifurcation We now describe a more general class of circuit equations where the resistor characteristic depends on a parameter μ and is denoted by f μ . (Perhaps μ is the temperature of the resistor.) The physical behavior of the circuit (see Figure12.9) is then described by the system of differential equations on R 2 : dx dt = y − f μ ( x ) , dy dt = − x .

--- PDF PAGE 286 ---
12.4 A Hopf Bifurcation 271 μ = − 0.5 μ = − 0.1 μ = 0.05 μ = 0.2 Figure 12.9 Hopf bifurcation in the system x ′ = y − x 3 + μ x , y ′ = − x . Consider as an example the special case where f μ is described by f μ ( x ) = x 3 − μ x and the parameter μ lies in the interval [ − 1, 1]. When μ = 1 we have the van der Pol system from the previous section. As before, the only equilibrium point lies at the origin. The linearized system is Y ′ = ( μ 1 − 1 0 ) Y , and the eigenvalues are λ ± = 1 2 ( μ ± √ μ 2 − 4 ) . Thus the origin is a spiral sink for − 1 ≤ μ < 0 and a spiral source for 0 < μ ≤ 1. Indeed, when − 1 ≤ μ ≤ 0, the resistor is passive as the graph of f μ lies in the first and third quadrants. Therefore, all solutions tend to the origin in this case. This holds even in the case where μ = 0 and the linearization yields

--- PDF PAGE 287 ---
272 Chapter 12 Applications in Circuit Theory a center. The circuit is physically dead in that, after a period of transition, all the currents and voltages stay at 0 (or as close to 0 as we want). However, as μ becomes positive, the circuit becomes alive. It begins to oscil- late. This follows from the fact that the analysis of Section 12.3 applies to this system for all μ in the interval ( 0, 1]. We therefore see the birth of a (unique) periodic solution γ μ as μ increases through 0 (see Exercise 4 at the end of this chapter). As just shown, this solution attracts all other nonzero solutions. As in Chapter 8, Section 8.5, this is an example of a Hopf bifurcation . Further elaboration of the ideas in Section 12.3 can be used to show that γ μ → 0 as μ → 0 with μ > 0. Review Figure 12.9 for some phase portraits associated with this bifurcation. 12.5 Exploration: Neurodynamics One of the most important developments in the study of the firing of nerve cells or neurons was the development of a model for this phenomenon in giant squid in the 1950s by Hodgkin and Huxley [23]. They developed a four-dimensional system of differential equations that described the electro- chemical transmission of neuronal signals along the cell membrane, a work for which they later received the Nobel Prize. Roughly speaking, this system is similar to systems that arise in electrical circuits. The neuron consists of a cell body, or soma , that receives electrical stimuli. These stimuli are then conducted along the axon , which can be thought of as an electrical cable that connects to other neurons via a collection of synapses. Of course, the motion is not really electrical, as the current is not really made up of electrons but rather ions (predominantly sodium and potassium). See Edelstein-Keshet [15] or Murray [34] for a primer on the neurobiology behind these systems. The four-dimensional Hodgkin–Huxley system is difficult to deal with primarily because of the highly nonlinear nature of the equations. An impor- tant breakthrough from a mathematical point of view was achieved by Fitzhugh [18] and Nagumo et. al. [35], who produced a simpler model of the Hodgkin–Huxley model. Although this system is not as biologically accu- rate as the original system, it nevertheless does capture the essential behavior of nerve impulses, including the phenomenon of excitability alluded to in the following. The Fitzhugh–Nagumo system of equations is given by x ′ = y + x − x 3 3 + I y ′ = − x + a − by ,

--- PDF PAGE 288 ---
Exercises 273 where a and b are constants satisfying 0 < 3 2 ( 1 − a ) < b < 1, and I is a parameter. In these equations x is similar to the voltage and rep- resents the excitability of the system; the variable y represents a combination of other forces that tend to return the system to rest. The parameter I is a stimulus parameter that leads to excitation of the system; I is like an applied current. Note the similarity of these equations with the van der Pol equation of Section 12.3. 1. First assume that I = 0. Prove that this system has a unique equilibrium point ( x 0 , y 0 ) . Hint: Use the geometry of the nullclines for this rather than explicitly solving the equations. Also remember the restrictions placed on a and b . 2. Prove that this equilibrium point is always a sink. 3. Now suppose that I 6 = 0. Prove that there is still a unique equilibrium point ( x I , y I ) and that x I varies monotonically with I . 4. Determine values of x I for which the equilibrium point is a source and show that there must be a stable limit cycle in this case. 5. When I 6 = 0, the point ( x 0 , y 0 ) is no longer an equilibrium point. Nonetheless, we can still consider the solution through this point. Describe the qualitative nature of this solution as I moves away from 0. Explain in mathematical terms why biologists consider this phenomenon the “excitement” of the neuron. 6. Consider the special case where a = I = 0. Describe the phase plane for each b > 0 (no longer restricted to b < 1) as completely as possible. Describe any bifurcations that occur. 7. Now let I vary as well and again describe any bifurcations that occur. Describe in as much detail as possible the phase portraits that occur in the Ib -plane, with b > 0. 8. Extend the analysis of the previous problem to the case b ≤ 0. 9. Now fix b = 0 and let a and I vary. Sketch the bifurcation plane (the Ia - plane) in this case. E X E R C I S E S 1. Find the phase portrait for the differential equation x ′ = y − f ( x ) , f ( x ) = x 2 , y ′ = − x . Hint: Exploit the symmetry about the y -axis.

--- PDF PAGE 289 ---
274 Chapter 12 Applications in Circuit Theory 2. Let f ( x ) =    2 x − 3 if x > 1 − x if − 1 ≤ x ≤ 1 2 x + 3 if x < − 1. Consider the system x ′ = y − f ( x ) y ′ = − x . (a) Sketch the phase plane for this system. (b) Prove that this system has a unique closed orbit. 3. Let f a ( x ) =    2 x + a − 2 if x > 1 ax if − 1 ≤ x ≤ 1 2 x − a + 2 if x < − 1. Consider the system x ′ = y − f a ( x ) y ′ = − x . (a) Sketch the phase plane for this system for various values of a . (b) Describe the bifurcation that occurs when a = 0. 4. Consider the system described in Section 12.4: x ′ = y − ( x 3 − μ x ) y ′ = − x , where the parameter μ satisfies 0 ≤ μ < 1. Fill in the details of the proof that a Hopf bifurcation occurs at μ = 0. 5. Consider the system x ′ = μ( y − ( x 3 − x )) , μ > 0 y ′ = − x . Prove that this system has a unique nontrivial periodic solution γ μ . Show that as μ → ∞ , γ μ tends to the closed curve consisting of two horizontal line segments and two arcs on y = x 3 − x as shown in Figure 12.10. This

--- PDF PAGE 290 ---
Exercises 275 Figure 12.10 type of solution is called a relaxation oscillation . When μ is large, there are two quite different time scales along the periodic solution. When moving horizontally, we have x ′ very large, and so the solution makes this transit very quickly. On the other hand, near the cubic nullcline, x ′ = 0 while y ′ is bounded. Thus this transit is comparatively much slower. 6. Find the differential equations for the network shown in Figure 12.11, where the resistor is voltage controlled; that is, the resistor characteristic is the graph of a function g : R → R , i R = g ( v R ) . 7. Show that the LC circuit consisting of one inductor and one capacitor wired in a closed loop oscillates. 8. Determine the phase portrait of the following differential equation and in particular show there is a unique nontrivial periodic solution: x ′ = y − f ( x ) , y ′ = − g ( x ) , where all of the following are assumed: (a) g ( − x ) = − g ( x ) and xg ( x ) > 0 for all x 6 = 0 (b) f ( − x ) = − f ( x ) and f ( x ) < 0 for 0 < x < a R C L Figure 12.11

--- PDF PAGE 291 ---
276 Chapter 12 Applications in Circuit Theory (c) for x > a , f ( x ) is positive and increasing (d) f ( x ) → ∞ as x → ∞ 9. Consider the system x ′ = y y ′ = a ( 1 − x 4 ) y − x (a) Find all equilibrium points and classify them. (b) Sketch the phase plane. (c) Describe the bifurcation that occurs when a becomes positive. (d) Prove that there exists a unique closed orbit for this system when a > 0. (e) Show that all nonzero solutions of the system tend to this closed orbit when a > 0.

--- PDF PAGE 292 ---
13 Applications in Mechanics We turn our attention in this chapter to the earliest important examples of differential equations that, in fact, are connected with the origins of calculus. These equations were used by Newton to derive and unify the three laws of Kepler. In this chapter we give a brief derivation of two of Kepler’s laws and then discuss more general problems in mechanics. The equations of Newton, our starting point, have retained importance throughout the history of modern physics and lie at the root of that part of physics called classical mechanics. The examples here provide us with con- crete examples of historical and scientific importance. Furthermore, the case we consider most thoroughly here, that of a particle moving in a central force gravitational field, is simple enough so that the differential equations can be solved explicitly using exact, classical methods (just calculus!). However, with an eye toward the more complicated mechanical systems that cannot be solved in this way, we also describe a more geometric approach to this problem. 13.1 Newton’s Second Law We will be working with a particle moving in a force field F . Mathematically, F is just a vector field on the (configuration) space of the particle, which in our case will be R n . From the physical point of view, F ( X ) is the force exerted on a particle located at position X . Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00013-0 c © 2013 Elsevier Inc. All rights reserved. 277

--- PDF PAGE 293 ---
278 Chapter 13 Applications in Mechanics The example of a force field we will be most concerned with is the gravita- tional field of the sun: F ( X ) is the force on a particle located at X which attracts the particle to the sun. We go into details of this system in Section 13.3. The connection between the physical concept of force field and the mathe- matical concept of differential equation is Newton’s second law : F = ma . This law asserts that a particle in a force field moves in such a way that the force vector at the location X of the particle, at any instant, equals the accelera- tion vector of the particle times the mass m . That is, Newton’s law gives the second-order differential equation mX ′′ = F ( X ) . As a system, this equation becomes X ′ = V V ′ = 1 m F ( X ) , where V = V ( t ) is the velocity of the particle. This is a system of equations on R n × R n . This type of system is often called a mechanical system with n degrees of freedom . A solution X ( t ) ⊂ R n of the second-order equation is said to lie in config- uration space . The solution of the system ( X ( t ) , V ( t )) ⊂ R n × R n lies in the phase space or state space of the system. Example. Recall the simple undamped harmonic oscillator from Chapter 2. In this case the mass moves in one dimension and its position at time t is given by a function x ( t ) , where x : R → R . As we saw, the differential equation governing this motion is mx ′′ = − kx for some constant k > 0. That is, the force field at the point x ∈ R is given by − kx .  Example. The two-dimensional version of the harmonic oscillator allows the mass to move in the plane, so the position is now given by the vector X ( t ) = ( x 1 ( t ) , x 2 ( t )) ∈ R 2 . As in the one-dimensional case, the force field is F ( X ) = − kX so the equations of motion are the same mX ′′ = − kX ,

--- PDF PAGE 294 ---
13.1 Newton’s Second Law 279 with solutions in configuration space given by x 1 ( t ) = c 1 cos ( √ k / m t ) + c 2 sin ( √ k / m t ) , x 2 ( t ) = c 3 cos ( √ k / m t ) + c 4 sin ( √ k / m t ) for some choices of the c j , as is easily checked using the methods in Chapter 6.  Before dealing with more complicated cases of Newton’s Law, we need to recall a few concepts from multivariable calculus. Recall that the dot product (or inner product ) of two vectors X , Y ∈ R n is denoted by X · Y and defined by X · Y = n ∑ i = 1 x i y i , where X = ( x 1 , . . . , x n ) . Thus X · X = | X | 2 . If X , Y : I → R n are smooth func- tions, then a version of the product rule yields ( X · Y ) ′ = X ′ · Y + X · Y ′ , as can be easily checked using the coordinate functions x i and y i . Recall also that if g : R n → R , the gradient of g , denoted grad g , is defined by grad g ( X ) = ( ∂ g ∂ x 1 ( X ) , . . . , ∂ g ∂ x n ( X ) ) . As we saw in Chapter 9, grad g is a vector field on R n . Next, consider the composition of two smooth functions g ◦ F , where F : R → R n and g : R n → R . The chain rule applied to g ◦ F yields d dt g ( F ( t )) = grad g ( F ( t )) · F ′ ( t ) = n ∑ i = 1 ∂ g ∂ x i ( F ( t )) dF i dt ( t ) . We will also use the cross product (or vector product) U × V of vectors U , V ∈ R 3 . By definition, U × V = ( u 2 v 3 − u 3 v 2 , u 3 v 1 − u 1 v 3 , u 1 v 2 − u 2 v 1 ) ∈ R 3 .

--- PDF PAGE 295 ---
280 Chapter 13 Applications in Mechanics Recall from multivariable calculus that we have U × V = − V × U = | U || V | N sin θ , where N is a unit vector perpendicular to U and V , with the orientations of the vectors U , V , and N given by the “right-hand rule.” Here θ is the angle between U and V . Note that U × V = 0 if and only if one vector is a scalar multiple of the other. Also, if U × V 6 = 0, then U × V is perpendicular to the plane containing U and V . If U and V are functions of t in R , then another version of the product rule asserts that d dt ( U × V ) = U ′ × V + U × V ′ , as one can again check by using coordinates. 13.2 Conservative Systems Many force fields appearing in physics arise in the following way. There is a smooth function U : R n → R such that F ( X ) = − ( ∂ U ∂ x 1 ( X ) , ∂ U ∂ x 2 ( X ) , . . . , ∂ U ∂ x n ( X ) ) = − grad U ( X ) . (The negative sign is traditional.) Such a force field is called conservative . The associated system of differential equations X ′ = V V ′ = − 1 m grad U ( X ) is called a conservative system. The function U is called the potential energy of the system. (More properly, U should be called a potential energy since adding a constant to it does not change the force field − grad U ( X ) .) Example. The preceding planar harmonic oscillator corresponds to the force field F ( X ) = − kX . This field is conservative, with potential energy U ( X ) = 1 2 k | X | 2 .

--- PDF PAGE 296 ---
13.2 Conservative Systems 281 For any moving particle X ( t ) of mass m , the kinetic energy is defined to be K = 1 2 m | V ( t ) | 2 . Note that the kinetic energy depends on velocity, while the potential energy is a function of position. The total energy (or sometimes simply energy ) is defined on phase space by E = K + U . The total energy function is important in mechanics because it is constant along any solution curve of the system. That is, in the language of Chapter 9, Section 9.4, E is a constant of the motion or a first integral for the flow.  Theorem. (Conservation of Energy) Let ( X ( t ) , V ( t )) be the solution curve of a conservative system. Then the total energy E is constant along this solution curve. Proof: To show that E ( X ( t )) is constant in t , we compute ˙ E = d dt ( 1 2 m | V ( t ) | 2 + U ( X ( t )) ) = mV · V ′ + ( grad U ) · X ′ = V · ( − grad U ) + ( grad U ) · V = 0.  We remark that we may also write this type of system in Hamiltonian form. Recall from Chapter 9 that a Hamiltonian system on R n × R n is a system of the form x ′ i = ∂ H ∂ y i y ′ i = − ∂ H ∂ x i , where H : R n × R n → R is the Hamiltonian function . As we have seen, the function H is constant along solutions of such a system. To write the conser- vative system in Hamiltonian form, we make a simple change of variables. We introduce the momentum vector Y = mV and then set H = K + U = 1 2 m ∑ y 2 i + U ( x 1 , . . . , x n ) . This puts the conservative system in Hamiltonian form, as is easily checked.

--- PDF PAGE 297 ---
282 Chapter 13 Applications in Mechanics 13.3 Central Force Fields A force field F is called central if F ( X ) points directly toward or away from the origin for every X . In other words, the vector F ( X ) is always a scalar multiple of X : F ( X ) = λ( X ) X , where the coefficient λ( X ) depends on X . We often tacitly exclude from con- sideration a particle at the origin; many central force fields are not defined (or are “infinite”) at the origin. We deal with these types of singularities in Sec- tion 13.7. Theoretically, the function λ( X ) could vary for different values of X on a sphere given by | X | = constant. However, if the force field is conservative, this is not the case. Proposition. Let F be a conservative force field. Then the following state- ments are equivalent: 1. F is central 2. F ( X ) = f ( | X | ) X 3. F ( X ) = − grad U ( X ) and U ( X ) = g ( | X | ) Proof: Suppose (3) is true. To prove (2) we find, from the Chain Rule: ∂ U ∂ x j = g ′ ( | X | ) ∂ ∂ x j ( x 1 2 + x 2 2 + x 3 2 ) 1 / 2 = g ′ ( | X | ) | X | x j . This proves (2) with f ( | X | ) = − g ′ ( | X | )/ | X | . It is clear that (2) implies (1). To show that (1) implies (3) we must prove that U is constant on each sphere: S α = { X ∈ R n | | X | = α > 0 } . Since any two points in S α can be connected by a curve in S α , it suffices to show that U is constant on any curve in S α . Thus if J ⊂ R is an interval and γ : J → S α is a smooth curve, we must show that the derivative of the composition U ◦ γ is identically 0. This derivative is d dt U (γ ( t )) = grad U (γ ( t )) · γ ′ ( t ) ,

--- PDF PAGE 298 ---
13.3 Central Force Fields 283 as in Section 13.1. Now grad U ( X ) = − F ( X ) = − λ( X ) X since F is central. Thus we have d dt U (γ ( t )) = − λ(γ ( t ))γ ( t ) · γ ′ ( t ) = − λ(γ ( t )) 2 d dt | γ ( t ) | 2 = 0 because | γ ( t ) | ≡ α .  Consider now a central force field, not necessarily conservative, defined on R 3 . Suppose that, at some time t 0 , P ⊂ R 3 denotes the plane containing the position vector X ( t 0 ) , the velocity vector V ( t 0 ) , and the origin (assuming, for the moment, that the position and velocity vectors are not collinear). Suppose that the force vector F ( X ( t 0 )) also lies in P . This makes it plausible that the particle stays in the plane P for all time. In fact, this is true: Proposition. A particle moving in a central force field in R 3 always moves in a fixed plane. Proof: Suppose X ( t ) is the path of a particle moving under the influence of a central force field. We have d dt ( X × V ) = V × V + X × V ′ = X × X ′′ = 0 because X ′′ is a scalar multiple of X . Therefore, Y = X ( t ) × V ( t ) is a constant vector. If Y 6 = 0, this means that X and V always lie in the plane orthogonal to Y , as asserted. If Y = 0, then X ′ ( t ) = g ( t ) X ( t ) for some real function g ( t ) . This means that the velocity vector of the moving particle is always directed along the line through the origin and the particle, as is the force on the particle. This implies that the particle always moves along the same line through the origin. To prove this, let ( x 1 ( t ) , x 2 ( t ) , x 3 ( t )) be the coordinates of X ( t ) . Then we have three separable differential equations dx k dt = g ( t ) x k ( t ) , for k = 1, 2, 3.

--- PDF PAGE 299 ---
284 Chapter 13 Applications in Mechanics Integrating, we find x k ( t ) = e h ( t ) x k ( 0 ) , where h ( t ) = t ∫ 0 g ( s ) ds . Therefore, X ( t ) is always a scalar multiple of X ( 0 ) and so X ( t ) moves in a fixed line and thus in a fixed plane.  The vector m ( X × V ) is called the angular momentum of the system, where m is the mass of the particle. By the proof of the preceding proposition, this vector is also conserved by the system. Corollary. (Conservation of Angular Momentum) Angular momentum is constant along any solution curve in a central force field.  We now restrict attention to a conservative central force field. Because of the previous proposition, the particle remains for all time in a plane, which we may take to be x 3 = 0. In this case angular momentum is given by the vector ( 0, 0, m ( x 1 v 2 − x 2 v 1 )) . Let ` = m ( x 1 v 2 − x 2 v 1 ) . Thus the function ` is also constant along solutions. In the planar case we also call ` the angular momentum. Introducing polar coordinates x 1 = r cos θ and x 2 = r sin θ , we find v 1 = x ′ 1 = r ′ cos θ − r sin θ θ ′ v 2 = x ′ 2 = r ′ sin θ + r cos θ θ ′ . Then x 1 v 2 − x 2 v 1 = r cos θ( r ′ sin θ + r cos θ θ ′ ) − r sin θ( r ′ cos θ − r sin θ θ ′ ) = r 2 ( cos 2 θ + sin 2 θ)θ ′ = r 2 θ ′ . Thus, in polar coordinates, ` = mr 2 θ ′ . We can now prove one of Kepler’s laws. Let A ( t ) denote the area swept out by the vector X ( t ) in the time from t 0 to t . In polar coordinates we have dA = 1 2 r 2 d θ . We define the areal velocity to be A ′ ( t ) = 1 2 r 2 ( t )θ ′ ( t ) ,

--- PDF PAGE 300 ---
13.4 The Newtonian Central Force System 285 the rate at which the position vector sweeps out the area. Kepler observed that the line segment joining a planet to the sun sweeps out equal areas in equal times, which we interpret to mean A ′ = constant. We have therefore proved more generally that this is true for any particle moving in a conservative central force field. We now have found two constants of the motion or first integrals for a con- servative system generated by a central force field: total energy and angular momentum. In the nineteenth century, the idea of solving a differential equa- tion was tied to the construction of a sufficient number of such constants of the motion. In the twentieth century, it became apparent that first integrals do not exist for differential equations very generally; the culprit here is chaos, which we will discuss in the next two chapters. Basically, chaotic behavior of solutions of a differential equation in an open set precludes the existence of first integrals in that set. 13.4 The Newtonian Central Force System We now direct the discussion to the Newtonian central force system. This sys- tem deals with the motion of a single planet orbiting around the sun. We assume that the sun is fixed at the origin in R 3 and that the relatively small planet exerts no force on the sun. The sun exerts a force on a planet given by Newton’s law of gravitation , which is also called the inverse square law . This law states that the sun exerts a force on a planet located at X ∈ R 3 with a magni- tude that is gm s m p / r 2 , where m s is the mass of the sun, m p is the mass of the planet, and g is the gravitational constant. The direction of the force is toward the sun. Therefore, Newton’s law yields the differential equation m p X ′′ = − gm s m p X | X | 3 . For clarity, we change units so that the constants are normalized to one and so the equation becomes more simply X ′′ = F ( X ) = − X | X | 3 , where F is now the force field. As a system of differential equations, we have X ′ = V V ′ = − X | X | 3 .

--- PDF PAGE 301 ---
286 Chapter 13 Applications in Mechanics This system is called the Newtonian central force system . Our goal in this sec- tion is to describe the geometry of this system; in the next section we derive a complete analytic solution of this system. Clearly, this is a central force field. Moreover, it is conservative, since X | X | 3 = grad U ( X ) , where the potential energy U is given by U ( X ) = − 1 | X | . Observe that F ( X ) is not defined at 0; indeed, the force field becomes infinite as the moving mass approaches collision with the stationary mass at the origin. As in the previous section we may restrict attention to particles moving in the plane R 2 . Thus we look at solutions in the configuration space C = R 2 − { 0 } . We denote the phase space by P = ( R 2 − { 0 } ) × R 2 . We visualize phase space as the collection of all tangent vectors at each point X ∈ C . Let T X = { ( X , V ) | V ∈ R 2 } . T X is the tangent plane to the configuration space at X . Then P = ⋃ X ∈ C T X is the tangent space to the configuration space, which we may naturally identify with a subset of R 4 . The dimension of phase space is four. However, we can cut this dimension in half by making use of the two known first integrals, total energy and angular momentum. Recall that energy is constant along solutions and is given by E ( X , V ) = K ( V ) + U ( X ) = 1 2 | V | 2 − 1 | X | . Let 6 h denote the subset of P consisting of all points ( X , V ) with E ( X , V ) = h . 6 h is called an energy surface with total energy h . If h ≥ 0, then 6 h meets each T X in a circle of tangent vectors satisfying | V | 2 = 2 ( h + 1 | X | ) . The radius of these circles in the tangent planes at X tends to ∞ as X → 0 and decreases to 2 h as | X | tends to ∞ .

--- PDF PAGE 302 ---
13.4 The Newtonian Central Force System 287 When h < 0, the structure of the energy surface 6 h is different. If | X | > − 1 / h , then there are no vectors in T X ∩ 6 h . When | X | = − 1 / h , only the zero vector in T X lies in 6 h . The circle r = − 1 / h in configuration space is therefore known as the zero velocity curve . If X lies inside the zero velocity curve, then T X meets the energy surface in a circle of tangent vectors as before. Figure 13.1 gives a caricature of 6 h in the case h < 0. We now introduce polar coordinates in configuration space and new variables ( v r , v θ ) in the tangent planes via V = v r ( cos θ sin θ ) + v θ ( − sin θ cos θ ) . We have V = X ′ = r ′ ( cos θ sin θ ) + r θ ′ ( − sin θ cos θ ) so that r ′ = v r and θ ′ = v θ / r . Differentiating once more, we find − 1 r 2 ( cos θ sin θ ) = − X | X | 3 = V ′ = ( v ′ r − v θ θ ′ ) ( cos θ sin θ ) + ( v r v θ r + v ′ θ ) ( − sin θ cos θ ) . Zero velocity curve Figure 13.1 Over each nonzero point inside the zero velocity curve, T X meets the energy surface 6 h in a circle of tangent vectors.

--- PDF PAGE 303 ---
288 Chapter 13 Applications in Mechanics Therefore, in the new coordinates ( r , θ , v r , v θ ) , the system becomes r ′ = v r θ ′ = v θ / r v ′ r = − 1 r 2 + v 2 θ r v ′ θ = − v r v θ r . In these coordinates, total energy is given by 1 2 ( v 2 r + v 2 θ ) − 1 r = h , and angular momentum is given by ` = rv θ . Let 6 h , ` consist of all points in phase space with total energy h and angular momentum ` . For simplicity, we will restrict attention to the case where h < 0. If ` = 0, we must have v θ = 0. So if X lies inside the zero velocity curve, the tangent space at X meets 6 h ,0 in precisely two vectors of the form ± v r ( cos θ sin θ ) , both of which lie on the line connecting 0 and X , one pointing toward 0, the other pointing away. On the zero velocity curve, only the zero vector lies in 6 h ,0 . Thus we see immediately that each solution in 6 h ,0 lies on a straight line through the origin. The solution leaves the origin and travels along a straight line until reaching the zero velocity curve, after which time it recedes back to the origin. In fact, since the vectors in 6 h ,0 have magnitude tending to ∞ as X → 0, these solutions reach the singularity in finite time in both directions. Solutions of this type are called collision-ejection orbits . When ` 6 = 0, a different picture emerges. Given X inside the zero velocity curve, we have v θ = `/ r , so that, from the total energy formula, r 2 v 2 r = 2 hr 2 + 2 r − ` 2 . ( ∗ ) The quadratic polynomial in r on the right in Equation ( ∗ ) must therefore be nonnegative, so this puts restrictions on which r -values can occur for X ∈ 6 h , ` . The graph of this quadratic polynomial is concave down since h < 0. It has no real roots if ` 2 > − 1 / 2 h . Therefore, the space 6 h , ` is empty in this case. If ` 2 = − 1 / 2 h , we have a single root that occurs at r = − 1 / 2 h . Thus this is the only allowable r -value in 6 h , ` in this case. In the tangent plane at ( r , θ) ,

--- PDF PAGE 304 ---
13.4 The Newtonian Central Force System 289 we have v r = 0, v θ = − 2 h ` , so this represents a circular closed orbit (traversed clockwise if ` < 0, counterclockwise if ` > 0). If ` 2 < − 1 / 2 h , then this polynomial has a pair of distinct roots at α , β with α < − 1 / 2 h < β . Note that α > 0. Let A α , β be the annular region α ≤ r ≤ β in configuration space. We therefore have that motion in configuration space is confined to A α , β . Proposition. Suppose h < 0 and ` 2 < − 1 / 2 h. Then 6 h , ` ⊂ P is a two- dimensional torus. Proof: We compute the set of tangent vectors lying in T X ∩ 6 h , ` for each X ∈ A α , β . If X lies on the boundary of the annulus, the quadratic term on the right of Equation ( ∗ ) vanishes, and so v r = 0 while v θ = `/ r . Thus there is a unique tangent vector in T X ∩ 6 h , ` when X lies on the boundary of the annulus. When X is in the interior of A α , β , we have v ± r = ± 1 r √ 2 hr 2 + 2 r − ` 2 , v θ = `/ r , so that we have a pair of vectors in T X ∩ 6 h , ` in this case. Note that these vectors all point either clockwise or counterclockwise in A α , β , since v θ has the same sign for all X . See Figure 13.2. Thus we can think of 6 h , ` as being given by a pair of graphs over A α , β : a positive graph given by v + r and a negative graph given by v − r that are joined together along the boundary circles r = α and r = β . (Of course, the “real” picture is a subset of R 4 .) This yields the torus.  It is tempting to think that the two curves in the torus given by r = α and r = β are closed orbits for the system, but this is not the case. This follows A α , β Figure 13.2 A selection of vectors in 6 h , ` .

--- PDF PAGE 305 ---
290 Chapter 13 Applications in Mechanics A α , β Figure 13.3 Solutions in 6 h , ` that meet r = α or r = β . since, when r = α , we have v ′ r = − 1 α 2 + v 2 θ α = 1 α 3 ( − α + ` 2 ) . However, since the right side of Equation ( ∗ ) vanishes at α , we have 2 h α 2 + 2 α − ` 2 = 0, so that − α + ` 2 = ( 2 h α + 1 )α . Since α < − 1 / 2 h , it follows that r ′′ = v ′ r > 0 when r = α , so the r -coordinate of solutions in 6 h , ` reaches a minimum when the curve meets r = α . Similarly, along r = β , the r -coordinate reaches a maximum. Thus solutions in A α , β must behave as shown in Figure 13.3. One can eas- ily show, incidentally, that these curves are preserved by rotations about the origin, so all of these solutions behave symmetrically. More, however, can be said. Each of these solutions actually lies on a closed orbit that traces out an ellipse in configuration space. To see this, we need to turn to analysis. 13.5 Kepler’s First Law For most nonlinear mechanical systems, the geometric analysis of the previous section is just about all we can hope for. In the Newtonian central force system, however, we get lucky: As has been known for centuries, we can write down explicit solutions for this system.

--- PDF PAGE 306 ---
13.5 Kepler’s First Law 291 Consider a particular solution curve of the differential equation. We have two constants of the motion for this system, namely the angular momentum ` and total energy E . The case ` = 0 yields collision–ejection solutions, as we saw before. Thus we assume ` 6 = 0. We will show that in polar coordinates in configuration space, a solution with nonzero angular momentum lies on a curve given by r ( 1 +  cos θ) = κ , where  and κ are constants. This equation defines a conic section, as can be seen by rewriting this equation in Cartesian coordinates. This fact is known as Kepler’s first law . To prove this, recall that r 2 θ ′ is constant and nonzero. Thus the sign of θ ′ remains constant along each solution curve, and so θ is always increasing or always decreasing in time. Therefore, we may also regard r as a function of θ along the curve. Let W ( t ) = 1 / r ( t ) ; then W is also a function of θ . Note that W = − U . The following proposition gives a convenient formula for kinetic energy. Proposition. The kinetic energy is given by K = ` 2 2 ( ( dW d θ ) 2 + W 2 ) . Proof: In polar coordinates, we have K = 1 2 ( ( r ′ ) 2 + ( r θ ′ ) 2 ) . Since r = 1 / W , we also have r ′ = − 1 W 2 dW d θ θ ′ = − ` dW d θ . Finally, r θ ′ = ` r = ` W . Substitution into the formula for K then completes the proof.  Now we find a differential equation relating W and θ along the solution curve. Observe that K = E − U = E + W . From the proposition we get ( dW d θ ) 2 + W 2 = 2 ` 2 ( E + W ) . ( ∗∗ )

--- PDF PAGE 307 ---
292 Chapter 13 Applications in Mechanics Differentiating both sides with respect to θ , dividing by 2 dW / d θ , and using dE / d θ = 0 (conservation of energy), we obtain d 2 W d θ 2 + W = 1 ` 2 , where 1 /` 2 is a constant. Note that this equation is just the equation for a harmonic oscillator with constant forcing 1 /` 2 . From elementary calculus, solutions of this second- order equation may be written in the form W (θ) = 1 ` 2 + A cos θ + B sin θ or, equivalently, W (θ) = 1 ` 2 + C cos (θ + θ 0 ) , where the constants C and θ 0 are related to A and B . If we substitute this expression into Equation ( ∗∗ ) and solve for C (at, say, θ + θ 0 = π/ 2 ) , we find C = ± 1 ` 2 √ 1 + 2 ` 2 E . Inserting this into the preceding solution, we find W (θ) = 1 ` 2 ( 1 ± √ 1 + 2 E ` 2 cos (θ + θ 0 ) ) . There is no need to consider both signs in front of the radical since cos (θ + θ 0 + π) = − cos (θ + θ 0 ) . Moreover, by changing the variable θ to θ − θ 0 we can put any particular solution into the form 1 ` 2 ( 1 + √ 1 + 2 E ` 2 cos θ ) .

--- PDF PAGE 308 ---
13.6 The Two-Body Problem 293 This looks pretty complicated. However, recall from analytic geometry (or from Exercise 2 at the end of this chapter) that the equation of a conic in polar coordinates is 1 r = 1 κ ( 1 +  cos θ) . Here κ is the latus rectum and  ≥ 0 is the eccentricity of the conic. The origin is a focus, and the three cases  > 1,  = 1, and  < 1 correspond respectively to a hyperbola, parabola, and ellipse. The case  = 0 is a circle. In our case we have  = √ 1 + 2 E ` 2 , so the three different cases occur when E > 0, E = 0, or E < 0. We have proved the following: Theorem. (Kepler’s First Law) The path of a particle moving under the influence of Newton’s law of gravitation is a conic of eccentricity √ 1 + 2 E ` 2 . This path lies along a hyperbola, parabola, or ellipse according to whether E > 0 , E = 0 , or E < 0 .  13.6 The Two-Body Problem We now turn our attention briefly to what at first appears to be a more difficult problem, the two-body problem . In this system we assume that we have two masses that move in space acording to their mutual graviational attraction. Let X 1 , X 2 denote the positions of particles of mass m 1 , m 2 in R 3 . So X 1 = ( x 1 1 , x 1 2 , x 1 3 ) and X 2 = ( x 2 1 , x 2 2 , x 2 3 ) . From Newton’s law of gravitation, we find the equations of motion m 1 X ′′ 1 = gm 1 m 2 X 2 − X 1 | X 2 − X 1 | 3 m 2 X ′′ 2 = gm 1 m 2 X 1 − X 2 | X 1 − X 2 | 3 .

--- PDF PAGE 309 ---
294 Chapter 13 Applications in Mechanics Let’s examine these equations from the perspective of a viewer living on the first mass. Let X = X 2 − X 1 . We then have X ′′ = X ′′ 2 − X ′′ 1 = gm 1 X 1 − X 2 | X 1 − X 2 | 3 − gm 2 X 2 − X 1 | X 1 − X 2 | 3 = − g ( m 1 + m 2 ) X | X | 3 . But this is just the Newtonian central force problem, with a different choice of constants. So, to solve the two-body problem, we first determine the solution of X ( t ) of this central force problem. This then determines the right side of the differ- ential equations for both X 1 and X 2 as functions of t , and so we may simply integrate twice to find X 1 ( t ) and X 2 ( t ) . Another way to reduce the two-body problem to the Newtonian central force is as follows. The center of mass of the two-body system is the vector X c = m 1 X 1 + m 2 X 2 m 1 + m 2 . A computation shows that X ′′ c = 0. Therefore, we must have X c = At + B where A and B are fixed vectors in R 3 . This says that the center of mass of the system moves along a straight line with constant velocity. We now change coordinates so that the origin of the system is located at X c . That is, we set Y j = X j − X c for j = 1, 2. Therefore, m 1 Y 1 ( t ) + m 2 Y 2 ( t ) = 0 for all t . Rewriting the differential equations in terms of the Y j , we find Y ′′ 1 = − gm 3 2 ( m 1 + m 2 ) 3 Y 1 | Y 1 | 3 Y ′′ 2 = − gm 3 1 ( m 1 + m 2 ) 3 Y 2 | Y 2 | 3 which yields a pair of central force problems. However, since we know that m 1 Y 1 ( t ) + m 2 Y 2 ( t ) = 0, we need only solve one of them. 13.7 Blowing up the Singularity The singularity at the origin in the Newtonian central force problem is the first time we have encountered such a situation. Usually our vector fields have been

--- PDF PAGE 310 ---
13.7 Blowing up the Singularity 295 well defined on all of R n . In mechanics, such singularities can sometimes be removed by a combination of judicious changes of variables and time scalings. In the Newtonian central force system, this may be achieved using a change of variables introduced by McGehee [32]. We first introduce scaled variables u r = r 1 / 2 v r u θ = r 1 / 2 v θ . In these variables the system becomes r ′ = r − 1 / 2 u r θ ′ = r − 3 / 2 u θ u ′ r = r − 3 / 2 ( 1 2 u 2 r + u 2 θ − 1 ) u ′ θ = r − 3 / 2 ( − 1 2 u r u θ ) . We still have a singularity at the origin, but note that the last three equations are all multiplied by r − 3 / 2 . We can remove these terms by simply multiplying the vector field by r 3 / 2 . In doing so, solution curves of the system remain the same but are parametrized differently. More precisely, we introduce a new time variable τ via the rule dt d τ = r 3 / 2 . By the Chain Rule we have dr d τ = dr dt dt d τ and similarly for the other variables. In this new time scale the system becomes ˙ r = ru r ˙ θ = u θ ˙ u r = 1 2 u 2 r + u 2 θ − 1 ˙ u θ = − 1 2 u r u θ ,

--- PDF PAGE 311 ---
296 Chapter 13 Applications in Mechanics where the dot now indicates differentiation with respect to τ . Note that, when r is small, dt / d τ is close to zero, so “time” τ moves much more slowly than time t near the origin. This system no longer has a singularity at the origin. We have “blown up” the singularity and replaced it with a new set given by r = 0 with θ , u r , u θ arbi- trary. On this set the system is now perfectly well defined. Indeed, the set r = 0 is an invariant set for the flow since ˙ r = 0 when r = 0. We have thus introduced a fictitious flow on r = 0. Although solutions on r = 0 mean nothing in terms of the real system, by continuity of solutions, they can tell us a lot about how solutions behave near the singularity. We need not concern ourselves with all of r = 0 since the total energy relation in the new variables becomes hr = 1 2 ( u 2 r + u 2 θ ) − 1. On the set r = 0, only the subset 3 defined by u 2 r + u 2 θ = 2, θ arbitrary matters. 3 is called the collision surface for the system; how solutions behave on 3 dictate how solutions move near the singularity since any solution that approaches r = 0 necessarily comes close to 3 in our new coordinates. Note that 3 is a two-dimensional torus: It is formed by a circle in the θ -direction and a circle in the u r u θ -plane. On 3 the system reduces to ˙ θ = u θ ˙ u r = 1 2 u 2 θ ˙ u θ = − 1 2 u r u θ , where we have used the energy relation to simplify ˙ u r . This system is easy to analyze. We have ˙ u r > 0 provided u θ 6 = 0. Thus the u r -coordinate must increase along any solution in 3 with u θ 6 = 0. On the other hand, when u θ = 0, the system has equilibrium points. There are two circles of equilibria, one given by u θ = 0, u r = √ 2, and θ arbitrary, the other by u θ = 0, u r = − √ 2, and θ arbitrary. Let C ± denote these two circles with u r = ± √ 2 on C ± . All other solutions must travel from C − to C + since v θ increases along solutions.

--- PDF PAGE 312 ---
13.7 Blowing up the Singularity 297 To fully understand the flow on 3 , we introduce the angular variable ψ in each u r u θ -plane via u r = √ 2 sin ψ u θ = √ 2 cos ψ . The torus is now parametrized by θ and ψ . In θψ –coordinates, the system becomes ˙ θ = √ 2 cos ψ ˙ ψ = 1 √ 2 cos ψ . The circles C ± are now given by ψ = ± π/ 2. Eliminating time from this equation, we find d ψ d θ = 1 2 . Thus all nonequilibrium solutions have constant slope 1 / 2 when viewed in θψ -coordinates. See Figure 13.4. Now recall the collision–ejection solutions described in Section 13.4. Each of these solutions leaves the origin and then returns along a ray θ = θ ∗ in configuration space. The solution departs with v r > 0 (and so u r > 0) and returns with v r < 0 ( u r < 0). In our new four-dimensional coordinate system, C + C − ψ θ Figure 13.4 Solutions on 3 in θψ -coordinates. Recall that θ and ψ are both defined mod 2 π , so opposite sides of this square are identified to form a torus.

--- PDF PAGE 313 ---
298 Chapter 13 Applications in Mechanics Λ θ u θ ur θ ∗ Figure 13.5 A collision–ejection solution in the region r > 0 leaving and returning to 3 and a connecting orbit on the collision surface. it follows that this solution forms an unstable curve associated with the equi- librium point ( 0, θ ∗ , √ 2, 0 ) and a stable curve associated with ( 0, θ ∗ , − √ 2, 0 ) . See Figure 13.5. What happens to nearby noncollision solutions? Well, they come close to the “lower” equilibrium point with θ = θ ∗ , u r = − √ 2, then follow one of two branches of the unstable curve through this point up to the “upper” equilibrium point θ = θ ∗ , u r = + √ 2, and then depart near the unstable curve leaving this equilibrium point. Interpreting this motion in con- figuration space, we see that each near-collision solution approaches the origin and then retreats after θ either increases or decreases by 2 π units. Of course, we know this already, since these solutions whip around the origin in tight ellipses. 13.8 Exploration: Other Central Force Problems In this exploration, we consider the (non-Newtonian) central force problem where the potential energy is given by U ( X ) = − 1 | X | ν , where ν > 1. The primary goal is to understand near-collision solutions.

--- PDF PAGE 314 ---
13.9 Exploration: Classical Limits of QMSs 299 1. Write this system in polar coordinates ( r , θ , v r , v θ ) and state explicitly the formulas for total energy and angular momentum. 2. Using a computer, investigate the behavior of solutions of this system when h < 0 and ` 6 = 0. 3. Blow up the singularity at the origin via the change of variables u r = r ν/ 2 v r , u θ = r ν/ 2 v θ and an appropriate change of time scale; write down the new system. 4. Compute the vector field on the collision surface 3 determined in r = 0 by the total energy relation. 5. Describe the bifurcation that occurs on 3 when ν = 2. 6. Describe the structure of 6 h , ` for all ν > 1. 7. Describe the change of structure of 6 h , ` that occurs when ν passes through the value 2. 8. Describe the behavior of solutions in 6 h , ` when ν > 2. 9. Suppose 1 < ν < 2. Describe the behavior of solutions as they pass close to the singularity at the origin. 10. Using the fact that solutions of this system are preserved by rotations about the origin, describe the behavior of solutions in 6 h , ` when h < 0 and ` 6 = 0. 13.9 Exploration: Classical Limits of Quantum Mechanical Systems In this exploration we investigate the anisotropic Kepler problem. This is a classical mechanical system with two degrees of freedom that depends on a parameter μ . When μ = 1 the system reduces to the Newtonian central force system discussed in Section 13.4. When μ > 1 some anisotropy is introduced into the system, so that we no longer have a central force field. We still have some collision–ejection orbits, as in the central force system, but the behavior of nearby orbits is quite different from those when μ = 1. The anisotropic Kepler problem was first introduced by Gutzwiller as a clas- sical mechanical approximation to certain quantum mechanical systems. In particular, this system arises naturally when one looks for bound states of an electron near a donor impurity of a semiconductor. Here the potential is due to an ordinary Coulomb field, while the kinetic energy becomes anisotropic because of the electronic band structure in the solid. Equivalently, we can view this system as having an anisotropic potential energy function. Gutzwiller suggests that this situation is akin to an electron with mass in one direction

--- PDF PAGE 315 ---
300 Chapter 13 Applications in Mechanics that is larger than in other directions. For more background on the quantum mechanical applications of this work, we refer to Gutzwiller [20]. The anisotropic Kepler system is given by x ′′ = − μ x (μ x 2 + y 2 ) 3 / 2 y ′′ = − y (μ x 2 + y 2 ) 3 / 2 , where μ is a parameter that we assume is greater than 1. 1. Show that this system is conservative with potential energy given by U ( x , y ) = − 1 √ μ x 2 + y 2 ; write down an explicit formula for total energy. 2. Describe the geometry of the energy surface 6 h for energy h < 0. 3. Restricting to the case of negative energy, show that the only solutions that meet the zero velocity curve and are straight-line collision–ejection orbits for the system lie on the x - and y -axes in configuration space. 4. Show that angular momentum is no longer an integral for this system. 5. Rewrite this system in polar coordinates. 6. Using a change of variables and time rescaling as for the Newtonian central force problem (Section 13.7), blow up the singularity and write down a new system without any singularities at r = 0. 7. Describe the structure of the collision surface 3 (the intersection of 6 h with r = 0 in the scaled coordinates). In particular, why would someone call this surface a “bumpy torus?” 8. Find all equilibrium points on 3 and determine the eigenvalues of the linearized system at these points. Determine which equilibria on 3 are sinks, sources, and saddles. 9. Explain the bifurcation that occurs on 3 when μ = 9 / 8. 10. Find a function that is nondecreasing along all nonequilibrium solutions in 3 . 11. Determine the fate of the stable and unstable curves of the saddle points in the collision surface. Hint: Rewrite the equation on this surface to eliminate time and estimate the slope of solutions as they climb up 3 . 12. When μ > 9 / 8, describe in qualitative terms what happens to solutions that approach collision close to the collision–ejection orbits on the x -axis. In particular, how do they retreat from the origin in the con- figuration space? How do solutions approach collision when traveling near the y -axis?

--- PDF PAGE 316 ---
Exercises 301 13.10 Exploration: Motion of a Glider In this exploration we investigate the motion of a glider moving in an xy -plane, where x measures the horizontal direction and y the vertical direc- tion. Let v > 0 be the velocity and θ the angle of the nose of the plane, with θ = 0 indicating the horizontal direction. Besides gravity, there are two other forces that determine the motion of the glider: the drag (which is parallel to the velocity vector but in the opposite direction) and the lift (which is perpendicular to the velocity vector). 1. Assuming that both the drag and the lift are proportional to v 2 , use Newton’s second law to show that the system of equations governing the motion of the glider may be written as θ ′ = v 2 − cos θ v v ′ = − sin θ − Dv 2 , where D ≥ 0 is a constant. 2. Find all equilibrium solutions for this system and use linearization to determine their types. 3. When D = 0 show that v 3 − 3 v cos θ is constant along solution curves. Sketch the phase plane in this case and describe the corresponding motion of the glider in the xy -plane. 4. Describe what happens when D becomes positive. E X E R C I S E S 1. Which of the following force fields on R 2 are conservative? (a) F ( x , y ) = ( − x 2 , − 2 y 2 ) (b) F ( x , y ) = ( x 2 − y 2 , 2 xy ) (c) F ( x , y ) = ( x , 0 ) 2. Prove that the equation 1 r = 1 h ( 1 +  cos θ) determines a hyperbola, parabola, and ellipse when  > 1,  = 1, and  < 1, respectively.

--- PDF PAGE 317 ---
302 Chapter 13 Applications in Mechanics 3. Consider the case of a particle moving directly away from the origin at time t = 0 in the Newtonian central force system. Find a specific for- mula for this solution and discuss the corresponding motion of the particle. For which initial conditions does the particle eventually reverse direction? 4. In the Newtonian central force system, describe the geometry of 6 h , ` when h > 0 and h = 0. 5. Let F ( X ) be a force field on R 3 . Let X 0 , X 1 be points in R 3 and let Y ( s ) be a path in R 3 with s 0 ≤ s ≤ s 1 , parametrized by arc length s , from X 0 to X 1 . The work done in moving a particle along this path is defined to be the integral ∫ s 1 s 0 F ( y ( s )) · y ′ ( s ) ds , where Y ′ ( s ) is the (unit) tangent vector to the path. Prove that the force field is conservative if and only if the work is independent of the path. In fact, if F = − grad V , then the work done is V ( X 1 ) − V ( X 0 ) . 6. Describe solutions to the non-Newtonian central force system given by X ′′ = − X | X | 4 . 7. Discuss solutions of the equation X ′′ = X | X | 3 . This equation corresponds to a repulsive rather than attractive force at the origin. 8. The following three problems deal with the two-body problem. Let the potential energy be U = gm 1 m 2 | X 2 − X 1 | and grad j ( U ) = ( ∂ U ∂ x j 1 , ∂ U ∂ x j 2 , ∂ U ∂ x j 3 ) . Show that the equations for the two-body problem may be written m j X ′′ j = − grad j ( U ) .

--- PDF PAGE 318 ---
Exercises 303 9. Show that the total energy K + U of the system is a constant of the motion, where K = 1 2 ( m 1 | V 1 | 2 + m 2 | V 2 | 2 ) . 10. Define the angular momentum of the system by ` = m 1 ( X 1 × V 1 ) + m 2 ( X 2 × V 2 ) , and show that ` is also a first integral.

--- PDF PAGE 319 ---
This page intentionally left blank

--- PDF PAGE 320 ---
14 The Lorenz System So far, in all of the differential equations we have studied, we have not encoun- tered any “chaos.” The reason is simple: The linear systems of the first few chapters always have straightforward, predictable behavior. (OK, we may see solutions wrap densely around a torus as in the oscillators of Chapter 6, but this is not chaos.) Also, for the nonlinear planar systems of the last few chap- ters, the Poincar´ e–Bendixson Theorem completely eliminates any possibility of chaotic behavior. So, to find chaotic behavior, we need to look at nonlinear, higher-dimensional systems. In this chapter we investigate the system that is, without doubt, the most famous of all chaotic differential equations, the Lorenz system from meteorol- ogy. First formulated in 1963 by E. N. Lorenz as a vastly oversimplified model of atmospheric convection, this system possesses what has come to be known as a strange attractor. Before the Lorenz model started making headlines, the only types of stable attractors known in differential equations were equilibria and closed orbits. The Lorenz system truly opened up new horizons in all areas of science and engineering, as many of the phenomena present in the Lorenz system have later been found in all of the areas we have previously investigated (biology, circuit theory, mechanics, and elsewhere). In the ensuing nearly 50 years, much progress has been made in the study of chaotic systems. Be forewarned, however, that the analysis of the chaotic behavior of particular systems, such as the Lorenz system, is usually extremely difficult. Most of the chaotic behavior that is readily understandable arises from geometric models for particular differential equations, rather than from Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00014-2 c © 2013 Elsevier Inc. All rights reserved. 305

--- PDF PAGE 321 ---
306 Chapter 14 The Lorenz System the actual equations themselves. Indeed, this is the avenue we pursue here. We shall present a geometric model for the Lorenz system which can be com- pletely analyzed using tools from discrete dynamics. Although this model has been known for some 30 years, it is interesting to note the fact that this model was only shown to be equivalent to the Lorenz system in the year 1999. 14.1 Introduction In 1963, E. N. Lorenz [29] attempted to set up a system of differential equa- tions that would explain some of the unpredictable behavior of the weather. Most viable models for weather involve partial differential equations; Lorenz sought a much simpler and easier-to-analyze system. The Lorenz model may be somewhat inaccurately thought of as follows. Imagine a planet with an “atmosphere” that consists of a single fluid particle. As on earth, this particle is heated from below (and thus rises) and cooled from above (so then falls back down). Can a weather expert predict the “weather” on this planet? Sadly, the answer is no, which raises a lot of questions about the possibility of accurate weather prediction down here on earth, where we have quite a few more particles in our atmosphere. A little more precisely, Lorenz looked at a two-dimensional fluid cell that was heated from below and cooled from above. The fluid motion can be described by a system of differential equations involving infinitely many vari- ables. Lorenz made the tremendous simplifying assumption that all but three of these variables remained constant. The remaining independent variables then measured, roughly speaking, the rate of convective “overturning” ( x ), and the horizontal and vertical temperature variation ( y and z , respectively). The resulting motion led to a three-dimensional system of differential equa- tions which involved three parameters: the Prandtl number σ , the Rayleigh number r , and another parameter b that is related to the physical size of the system. When all of these simplifications were made, the system of differential equations involved only two nonlinear terms and was given by x ′ = σ ( y − x ) y ′ = rx − y − xz z ′ = xy − bz . In this system all three parameters are assumed to be positive and, more- over, σ > b + 1. We denote this system by X ′ = L ( X ) . In Figure 14.1, we have displayed the solution curves through two different initial conditions

--- PDF PAGE 322 ---
14.1 Introduction 307 x y z P 1 P 2 Figure 14.1 The Lorenz attractor. Two solutions with initial conditions P 1 = (0, 2, 0) and P 2 = (0, − 2, 0). P 1 = ( 0, 2, 0 ) and P 2 = ( 0, − 2, 0 ) when the parameters are σ = 10, b = 8 / 3, and r = 28. These are the original parameters that led to Lorenz’ discovery. Note how both solutions start out very differently, but eventually have more or less the same fate: They both seem to wind around a pair of points, alter- nating at times which point they encircle. This is the first important fact about the Lorenz system: All nonequilibrium solutions tend eventually to the same complicated set, the so-called Lorenz attractor . There is another important ingredient lurking in the background here. In Figure 14.1, we started with two relatively far apart initial conditions. Had we started with two very close initial conditions, we would not have observed the “transient behavior” apparent in Figure 14.1. Rather, more or less the same picture would have resulted for each solution. This, however, is mislead- ing. When we plot the actual coordinates of the solutions, we see that these two solutions actually move quite far apart during their journey around the Lorenz attractor. This is illustrated in Figure 14.2, where we have graphed the x -coordinates of two solutions that start out nearby, one at ( 0, 2, 0 ) , the other (in gray) at ( 0, 2.01, 0 ) . These graphs are nearly identical for a certain time period, but then they dif- fer considerably as one solution travels around one of the lobes of the attractor while the other solution travels around the other. No matter how close two solutions start, they always move apart in this manner when they are close to the attractor. This is sensitive dependence on initial conditions , one of the main features of a chaotic system.

--- PDF PAGE 323 ---
308 Chapter 14 The Lorenz System x t Figure 14.2 The x ( t ) graphs for two nearby initial conditions P 1 = (0, 2, 0) and P 2 = (0, 2.01, 0). We will describe in detail the concept of an attractor and chaos in this chap- ter. But first we need to investigate some of the more familiar features of the system. 14.2 Elementary Properties of the Lorenz System As usual, to analyze this system, we begin by finding the equilibria. Some easy algebra yields three equilibrium points, the origin, and Q ± = ( ± √ b ( r − 1 ) , ± √ b ( r − 1 ) , r − 1 ) . The last two equilibria only exist when r > 1, so already we see that we have a bifurcation when r = 1. Linearizing, we find the system Y ′ =   − σ σ 0 r − z − 1 − x y x − b   Y . At the origin, the eigenvalues of this matrix are − b and λ ± = 1 2 ( − (σ + 1 ) ± √ (σ + 1 ) 2 − 4 σ ( 1 − r ) ) .

--- PDF PAGE 324 ---
14.2 Elementary Properties of the Lorenz System 309 Note that both λ ± are negative when 0 ≤ r < 1. Thus the origin is a sink in this case. The Lorenz vector field L ( X ) possesses a symmetry. If we let S ( x , y , z ) = ( − x , − y , z ) , then we have DS ( L ( X )) = L ( S ( X )) . That is, reflection through the z -axis preserves the vector field. In particular, if ( x ( t ) , y ( t ) , z ( t )) is a solution of the Lorenz equations, then so is ( − x ( t ) , − y ( t ) , z ( t )) . When x = y = 0, we have x ′ = y ′ = 0, so the z -axis is invariant. On this axis, we have simply z ′ = − bz , so all solutions tend to the origin on this axis. In fact, the solution through any point in R 3 tends to the origin when r < 1, for we have the following. Proposition. Suppose r < 1 . Then all solutions of the Lorenz system tend to the equilibrium point at the origin. Proof: We construct a strict Liapunov function on all of R 3 . Let L ( x , y , z ) = x 2 + σ y 2 + σ z 2 . Then we have ˙ L = − 2 σ ( x 2 + y 2 − ( 1 + r ) xy ) − 2 σ bz 2 . We therefore have ˙ L < 0 away from the origin provided that g ( x , y ) = x 2 + y 2 − ( 1 + r ) xy > 0 for ( x , y ) 6 = ( 0, 0 ) . This is clearly true along the y -axis. Along any other straight-line y = mx in the plane we have g ( x , mx ) = x 2 ( m 2 − ( 1 + r ) m + 1 ) . But the quadratic term m 2 − ( 1 + r ) m + 1 is positive for all m if r < 1, as is easily checked. Thus g ( x , y ) > 0 for ( x , y ) 6 = ( 0, 0 ) .  When r increases through 1, two things happen. First, the eigenvalue λ + at the origin becomes positive, so the origin is now a saddle with a two- dimensional stable surface and an unstable curve. Second, the two equilibria Q ± are born at the origin when r = 1 and move away as r increases. Proposition. The equilibrium points Q ± are sinks provided 1 < r < r ∗ = σ ( σ + b + 3 σ − b − 1 ) .

--- PDF PAGE 325 ---
310 Chapter 14 The Lorenz System Proof: From the linearization, we calculate that the eigenvalues at Q ± satisfy the cubic polynomial f r (λ) = λ 3 + ( 1 + b + σ )λ 2 + b (σ + r )λ + 2 b σ ( r − 1 ) = 0. When r = 1 the polynomial f 1 has distinct roots at 0, − b , and − σ − 1. These roots are distinct since σ > b + 1, so − σ − 1 < − σ + 1 < − b < 0. Thus, for r close to 1, f r has three real roots close to these values. Note that f r (λ) > 0 for λ ≥ 0 and r > 1. Looking at the graph of f r , it follows that, at least for r close to 1, the three roots of f r must be real and negative. We now let r increase and ask what is the lowest value of r for which f r has an eigenvalue with zero real part. Note that this eigenvalue must in fact be of the form ± i ω with ω 6 = 0, since f r is a real polynomial that has no roots equal to 0 when r > 1. Solving f r ( i ω) = 0 by equating both real and imaginary parts to zero then yields the result (recall that we have assumed σ > b + 1).  We remark that a Hopf bifurcation is known to occur at r ∗ , but proving this is beyond the scope of this book. When r > 1 it is no longer true that all solutions tend to the origin. How- ever, we can say that solutions that start far from the origin do at least move closer in. To be precise, let V ( x , y , z ) = rx 2 + σ y 2 + σ ( z − 2 r ) 2 . Note that V ( x , y , z ) = ν > 0 defines an ellipsoid in R 3 centered at ( 0, 0, 2 r ) . We will show the following. Proposition. There exists ν ∗ such that any solution that starts outside the ellipsoid V = ν ∗ eventually enters this ellipsoid and then remains trapped therein for all future time. Proof: We compute ˙ V = − 2 σ ( rx 2 + y 2 + b ( z 2 − 2 rz ) ) = − 2 σ ( rx 2 + y 2 + b ( z − r ) 2 − br 2 ) . The equation rx 2 + y 2 + b ( z − r ) 2 = μ

--- PDF PAGE 326 ---
14.2 Elementary Properties of the Lorenz System 311 also defines an ellipsoid when μ > 0. When μ > br 2 we have ˙ V < 0. Thus we may choose ν ∗ large enough so that the ellipsoid V = ν ∗ strictly contains the ellipsoid rx 2 + y 2 + b ( z − r ) 2 = br 2 in its interior. Then ˙ V < 0 for all ν ≥ ν ∗ .  As a consequence, all solutions starting far from the origin are attracted to a set that sits inside the ellipsoid V = ν ∗ . Let 3 denote the set of all points with solutions that remain for all time (forward and backward) in this ellip- soid. Then the ω -limit set of any solution of the Lorenz system must lie in 3 . Theoretically, 3 could be a large set, perhaps bounding an open region in R 3 . However, for the Lorenz system, this is not the case. To see this, recall from calculus that the divergence of a vector field F ( X ) on R 3 is given by div F = 3 ∑ i = 1 ∂ F i ∂ x i ( X ) . The divergence of F measures how fast volumes change under the flow φ t of F . Suppose D is a region in R 3 with a smooth boundary, and let D ( t ) = φ t ( D ) , the image of D under the time t map of the flow. Let V ( t ) be the volume of D ( t ) . Then Liouville’s Theorem asserts that dV dt = ∫ D ( t ) div F dx dy dz . For the Lorenz system, we compute immediately that the divergence is the constant − (σ + 1 + b ) so that volume decreases at a constant rate: dV dt = − (σ + 1 + b ) V . Solving this simple differential equation yields V ( t ) = e − (σ + 1 + b ) t V ( 0 ) , so that any volume must shrink exponentially fast to 0. In particular, we have this proposition. Proposition. The volume of 3 is zero.  The natural question is what more can we say about the structure of the “attractor” 3 . In dimension 2, such a set would consist of a collection of

--- PDF PAGE 327 ---
312 Chapter 14 The Lorenz System limit cycles, equilibrium points, and solutions connecting them. In higher dimensions, these attractors may be much “stranger,” as we show in the next section. 14.3 The Lorenz Attractor The behavior of the Lorenz system as the parameter r increases is the subject of much contemporary research; we are decades (if not centuries) away from rig- orously understanding all of the fascinating dynamical phenomena that occur as the parameters change. Sparrow has written an entire book devoted to this subject [44]. In this section we will deal with one specific set of parameters where the Lorenz system has an attractor. Roughly speaking, an attractor for the flow is an invariant set that “attracts” all nearby solutions. The following definition is more precise. Definition Let X ′ = F ( X ) be a system of differential equations in R n with flow φ t . A set 3 is called an attractor if 1. 3 is compact and invariant. 2. There is an open set U containing 3 such that for each X ∈ U , φ t ( X ) ∈ U and ∩ t ≥ 0 φ t ( U ) = 3 . 3. (Transitivity) Given any points Y 1 , Y 2 ∈ 3 and any open neighbor- hoods U j about Y j in U , there is a solution curve that begins in U 1 and later passes through U 2 . The transitivity condition in this definition may seem a little strange. Basi- cally, we include it to guarantee that we are looking at a single attractor rather than a collection of dynamically different attractors. For example, the transitivity condition rules out situations such as that given by the planar system x ′ = x − x 3 . y ′ = − y . The phase portrait of this system is shown in Figure 14.3. Note that any solution of this system enters the set marked U and then tends to one of the three equilibrium points: either to one of the sinks at ( ± 1, 0 ) or to the sad- dle ( 0, 0 ) . The forward intersection of the flow applied to U is the interval

--- PDF PAGE 328 ---
14.3 The Lorenz Attractor 313 U Figure 14.3 The interval on the x -axis between the two sinks is not an attractor for this system, despite the fact that all solutions enter U . − 1 ≤ x ≤ 1. This interval meets conditions (1) and (2) in the definition, but condition (3) is violated, as there are no solution curves passing close to points in both the left and right half of this interval. We choose not to consider this set an attractor since most solutions tend to one of the two sinks. We really have two distinct attractors in this case. As a remark, there is no universally accepted definition of an attractor in mathematics; some people choose to say that a set 3 that meets only condi- tions (1) and (2) is an attractor, while if 3 also meets condition (3), it is called a transitive attractor. For planar systems, condition (3) is usually easily veri- fied; in higher dimensions, however, this can be much more difficult, as we shall see. For the rest of this chapter, we restrict attention to the very special case of the Lorenz system where the parameters are given by σ = 10, b = 8 / 3, and r = 28. Historically, these are the values Lorenz used when he first encoun- tered chaotic phenomena in this system. Thus, the specific Lorenz system we consider is X ′ = L ( X ) =   10 ( y − x ) 28 x − y − xz xy − ( 8 / 3 ) z   . As in the previous section, we have three equilibria: the origin and Q ± = ( ± 6 √ 2, ± 6 √ 2, 27 ) . At the origin we find eigenvalues λ 1 = − 8 / 3 and λ ± = − 11 2 ± √ 1201 2 .

--- PDF PAGE 329 ---
314 Chapter 14 The Lorenz System x z y Figure 14.4 Linearization at the origin for the Lorenz system. For later use, note that these eigenvalues satisfy λ − < − λ + < λ 1 < 0 < λ + . The linearized system at the origin is then Y ′ =   λ − 0 0 0 λ + 0 0 0 λ 1 Y   . The phase portrait of the linearized system is shown in Figure 14.4. Note that all solutions in the stable plane of this system tend to the origin tangentially to the z -axis. At Q ± a computation shows that there is a single negative real eigenvalue and a pair of complex conjugate eigenvalues with positive real parts. In Figure 14.5, we have displayed a numerical computation of a portion of the left and right branches of the unstable curve at the origin. Note that the right portion of this curve comes close to Q − and then spirals away. The left portion behaves symmetrically under reflection through the z -axis. In Figure 14.6, we have displayed a significantly larger portion of these unstable curves. Note that they appear to circulate around the two equilibria, some- times spiraling around Q + , sometimes around Q − . In particular, these curves continually reintersect the portion of the plane z = 27 containing Q ± in which the vector field points downward. This suggests that we may construct a Poincar´ e map on a portion of this plane.

--- PDF PAGE 330 ---
14.3 The Lorenz Attractor 315 Q − Q + (0, 0, 0) Figure 14.5 Unstable curve at the origin. (0,0,0) Figure 14.6 More of the unstable curve at the origin. As we have seen before, computing a Poincar´ e map is often impossible, and this case is no different. So we will content ourselves with building a simplified model that exhibits much of the behavior we find in the Lorenz system. As we shall see in the following section, this model provides a computable means to assess the chaotic behavior of the system.

--- PDF PAGE 331 ---
316 Chapter 14 The Lorenz System 14.4 A Model for the Lorenz Attractor In this section we describe a geometric model for the Lorenz attractor origi- nally proposed by Guckenheimer and Williams [21]. Tucker [46] showed that this model does indeed correspond to the Lorenz system for certain parame- ters. Rather than specify the vector field exactly, we give instead a qualitative description of its flow, much as we did in Chapter 11. The specific numbers we use are not that important; only their relative sizes matter. We will assume that our model is symmetric under the reflection ( x , y , z ) → ( − x , − y , z ) , as is the Lorenz system. We first place an equilibrium point at the origin in R 3 and assume that, in the cube S given by | x | , | y | , | z | ≤ 5, the system is linear. Rather than use the eigenvalues λ 1 and λ ± from the actual Lorenz system, we simplify the computations a bit by assuming that the eigenvalues are − 1, 2, and − 3, and that the system is given in the cube by x ′ = − 3 x y ′ = 2 y z ′ = − z . Note that the phase portrait of this system agrees with that in Figure 14.4 and that the relative magnitudes of the eigenvalues are the same as in the Lorenz case. We need to know how solutions make the transit near ( 0, 0, 0 ) . Consider a rectangle R 1 in the plane z = 1 given by | x | ≤ 1, 0 < y ≤  < 1. As time moves forward, all solutions that start in R 1 eventually reach the rectangle R 2 in the plane y = 1 defined by | x | ≤ 1, 0 < z ≤ 1. Thus we have a function h : R 1 → R 2 defined by following solution curves as they pass from R 1 to R 2 . We leave it as an exercise to check that this function assumes the form h ( x y ) = ( x 1 z 1 ) = ( xy 3 / 2 y 1 / 2 ) . It follows that h takes lines y = c in R 1 to lines z = c 1 / 2 in R 2 . Also, since x 1 = xz 3 1 , we have that h maps lines x = c to curves of the form x 1 = cz 3 1 . Each of these image curves meet the xy -plane perpendicularly, as shown in Figure 14.7. Mimicking the Lorenz system, we place two additional equilibria in the plane z = 27, one at Q − = ( − 10, − 20, 27 ) and the other at Q + = ( 10, 20, 27 ) . We assume that the lines given by y = ± 20, z = 27 form portions of the stable lines at Q ± and that the other two eigenvalues at these points are complex with positive real parts.

--- PDF PAGE 332 ---
14.4 A Model for the Lorenz Attractor 317 h ( x , y ) ( x , y ) x z y 1 2 Figure 14.7 Solutions making the transit near (0, 0, 0). Let 6 denote the square | x | , | y | ≤ 20, z = 27. We assume that the vector field points downward in the interior of this square. Thus solutions spiral away from Q ± in the same manner as in the Lorenz system. We also assume that the stable surface of ( 0, 0, 0 ) first meets 6 in the line of intersection of the xz -plane and 6 . Let ζ ± denote the two branches of the unstable curve at the origin. We assume that these curves make a passage around 6 and then enter this square as shown in Figure 14.8. We denote the first point of intersection of ζ ± with 6 by ρ ± = ( ± x ∗ , ∓ y ∗ ) . Now consider a straight line y = ν in 6 . If ν = 0, all solutions beginning at points on this line tend to the origin as time moves forward. Thus these solutions never return to 6 . We assume that all other solutions originating in 6 do return to 6 as time moves forward. How these solutions return leads to our major assumptions about this model. These assumptions are 1. Return condition: Let 6 + = 6 ∩ { y > 0 } and 6 − = 6 ∩ { y < 0 } . We assume that the solutions through any point in 6 ± return to 6 in forward time. Thus we have a Poincar´ e map 8 : 6 + ∪ 6 − → 6 . We assume that the images 8(6 ± ) are as shown in Figure 14.8. By symmetry, we have 8( x , y ) = − 8( − x , − y ) . 2. Contracting direction: For each ν 6 = 0 we assume that 8 maps the line y = ν in 6 into the line y = g (ν) for some function g . Moreover, we assume that 8 contracts this line in the x -direction. 3. Expanding direction: We assume that 8 stretches 6 + and 6 − in the y - direction by a factor greater than √ 2, so that g ′ ( y ) > √ 2.

--- PDF PAGE 333 ---
318 Chapter 14 The Lorenz System Q − Q + F ( Σ + ) F ( Σ − ) Σ + Σ− ρ + ζ + ζ − ρ - z y x Figure 14.8 Solutions ζ ± and their intersection with 6 in the model for the Lorenz attractor. 4. Hyperbolicity condition: Besides the expansion and contraction, we assume that D 8 maps vectors tangent to 6 ± with slopes that are ± 1 to vectors with slopes of a magnitude larger than μ > 1. Analytically, these assumptions imply that the map 8 assumes the form 8( x , y ) = ( f ( x , y ) , g ( y )) , where g ′ ( y ) > √ 2 and 0 < ∂ f /∂ x < c < 1. The hyperbolicity condition implies that g ′ ( y ) > μ ∣ ∣ ∣ ∣ ∂ f ∂ x ± ∂ f ∂ y ∣ ∣ ∣ ∣ . Geometrically, this condition implies that the sectors in the tangent planes given by | y | ≥ | x | are mapped by D 8 strictly inside a sector with steeper slopes. Note that this condition holds if | ∂ f /∂ y | and c are sufficiently small throughout 6 ± . Technically, 8( x , 0 ) is not defined, but we do have lim y → 0 ± , 8( x , y ) = ρ ± where we recall that ρ ± is the first point of intersection of ζ ± with 6 . We call ρ ± the tip of 8(6 ± ) . In fact, our assumptions on the eigenvalues guarantee that g ′ ( y ) → ∞ as y → 0 (see Exercise 3 at the end of this chapter).

--- PDF PAGE 334 ---
14.4 A Model for the Lorenz Attractor 319 − 20 20 y* − y* R + R − Φ ( R + ) Φ ( R − ) Φ ( y = − y* ) Φ ( y = y* ) Figure 14.9 Poincar ´ e map 8 on R . To find the attractor, we may restrict attention to the rectangle R ⊂ 6 given by | y | ≤ y ∗ , where we recall that ∓ y ∗ is the y -coordinate of the tips ρ ± . Let R ± = R ∩ 6 ± . It is easy to check that any solution starting in the interior of 6 ± but outside R must eventually meet R , so it suffices to consider the behavior of 8 on R . A planar picture of the action of 8 on R is displayed in Figure 14.9. Note that 8( R ) ⊂ R . Let 8 n denote the n th iterate of 8 , and let A = ∞ ⋂ n = 0 8 n ( R ) . Here U denotes the closure of the set U . The set A will be the intersection of the attractor for the flow with R . That is, let A = ( ⋃ t ∈ R φ t ( A ) ) ∪ { ( 0, 0, 0 ) } . We add the origin here so that A will be a closed set. We will prove the following theorem. Theorem. A is an attractor for the model Lorenz system. Proof: The proof that A is an attractor for the flow follows immediately from the fact that A is an attractor for the mapping 8 (where an attractor for a map- ping is defined completely analogously to that for a flow). Clearly, A is closed. Technically, A itself is not invariant under 8 since 8 is not defined along y = 0. However, for the flow, the solutions through all such points do lie in A and so A is invariant. If we let O be the open set given by | x | < 20, | y | < 20 − 

--- PDF PAGE 335 ---
320 Chapter 14 The Lorenz System for an  with y ∗ < 20 −  , then for any ( x , y ) ∈ O , there is an n such that 8 n ( x , y ) ∈ R . Thus, lim n →∞ 8 n ( x , y ) ⊂ A for all ( x , y ) ∈ O . By definition, A = ∩ n ≥ 0 8 n ( R ) , and so A = ∩ n ≥ 0 8 n ( O ) as well. Therefore, conditions (1) and (2) in the definition of an attractor hold for 8 . It remains to show the transitivity property. We need to show that if P 1 and P 2 are points in A , and W j are open neighborhoods of P j in O , then there exists an n ≥ 0 such that 8 n ( W 1 ) ∩ W 2 6 = ∅ . Given a set U ⊂ R , let 5 y ( U ) denote the projection of U onto the y -axis. Also let ` y ( U ) denote the length of 5 y ( U ) , which we call the y -length of U . In the following, U will be a finite collection of connected sets, so ` y ( U ) is well defined. We need a lemma.  Lemma. For any open set W ⊂ R, there exists n > 0 such that 5 y (8 n ( W )) is the open interval ( − y ∗ , y ∗ ) . Equivalently, 8 n ( W ) meets each line y = c in the interior of R. Proof: First suppose that W contains a connected piece W ′ that extends from one of the tips to y = 0. Thus ` y ( W ′ ) = y ∗ . Then 8( W ′ ) is connected and we have ` y (8( W ′ )) > √ 2 y ∗ . Moreover, 8( W ′ ) also extends from one of the tips, but now crosses y = 0 since its y -length exceeds y ∗ . Now apply 8 again. 8 2 ( W ′ ) contains two pieces, one of which extends to ρ + , the other to ρ − . Moreover, ` y (8 2 ( W ′ )) > 2 y ∗ . Thus it follows that 5 y (8 2 ( W ′ )) = [ − y ∗ , y ∗ ] and so we are done in this special case. For the general case, suppose first that W is connected and does not cross y = 0. Then we have ` y (8( W )) > √ 2 ` y ( W ) as before, so the y -length of 8( W ) grows by a factor of more than √ 2. If W does cross y = 0, then we may find a pair of connected sets W ± with W ± ⊂ { R ± ∩ W } and ` y ( W + ∪ W − ) = ` y ( W ) . The images 8( W ± ) extend to the tips ρ ± . If either of these sets also meets y = 0, then we are done according to the preceding. If neither 8( W + ) nor 8( W − ) crosses y = 0, then we may apply 8 again. Both 8( W + ) and 8( W − ) are connected sets, and we have ` y (8 2 ( W ± )) > 2 ` y ( W ± ) . Thus, for one of W + or W − , we have ` y (8 2 ( W ± )) > ` y ( W ) and again the y -length of W grows under iteration. Thus, if we continue to iterate 8 or 8 2 and choose the appropriate largest subset of 8 j ( W ) at each stage as above, then we see that the y -lengths of these images grow without bound. This completes the proof of the lemma.  We now complete the proof of the theorem.

--- PDF PAGE 336 ---
14.5 The Chaotic Attractor 321 Proof: We must find a point in W 1 with an image under an iterate of 8 that lies in W 2 . Toward that end, note that ∣ ∣ ∣ 8 k ( x 1 , y ) − 8 k ( x 2 , y ) ∣ ∣ ∣ ≤ c k | x 1 − x 2 | since 8 j ( x 1 , y ) and 8 j ( x 2 , y ) lie on the same straight line parallel to the x -axis for each j and since 8 contracts distances in the x -direction by a factor of c < 1. We may assume that W 2 is a disk of diameter  . Recalling that the width of R in the x -direction is 40, we choose m such that 40 c m <  . Consider 8 − m ( P 2 ) . Note that 8 − m ( P 2 ) is defined since P 2 ∈ ∩ n ≥ 0 8 n ( R ) . Say 8 − m ( P 2 ) = (ξ , η) . From the lemma, we know that there exists n such that 5 y (8 n ( W 1 )) = [ − y ∗ , y ∗ ]. Thus we may choose a point (ξ 1 , η) ∈ 8 n ( W 1 ) . Say (ξ 1 , η) = 8 n ( ˜ x , ˜ y ) , where ( ˜ x , ˜ y ) ∈ W 1 , so that 8 n ( ˜ x , ˜ y ) and 8 − m ( P 2 ) have the same y -coordinate. Then we have | 8 m + n ( ˜ x , ˜ y ) − P 2 | = | 8 m (ξ 1 , η) − P 2 | = | 8 m (ξ 1 , η) − 8 m (ξ , η) | ≤ 40 c m <  . We have found a point ( ˜ x , ˜ y ) ∈ W 1 with a solution that passes through W 2 . This concludes the proof.  Note that, in the preceding proof, the solution curve that starts near P 1 and comes close to P 2 need not lie in the attractor. However, it is possible to find such a solution that does lie in A (see Exercise 4 at the end of this chapter). 14.5 The Chaotic Attractor In the previous section we reduced the study of the behavior of solutions of the Lorenz system to the analysis of the dynamics of the Poincar´ e map 8 . In the process, we dropped from a three-dimensional system of differential equations to a two-dimensional mapping. But we can do better. According to our assumptions, two points that share the same y -coordinate in 6 are mapped to two new points with y -coordinates given by g ( y ) and thus are again the same.

--- PDF PAGE 337 ---
322 Chapter 14 The Lorenz System Moreover, the distance between these points is contracted. It follows that, under iteration of 8 , we need not worry about all points on a line y = con- stant; we need only keep track of how the y -coordinate changes under iteration of g . Then, as we shall see, the Poincar´ e map 8 is completely determined by the dynamics of the one-dimensional function g defined on the interval [ − y ∗ , y ∗ ]. Indeed, iteration of this function completely determines the behav- ior of all solutions in the attractor. In this section, we begin to analyze the dynamics of this one-dimensional discrete dynamical system . In Chapter 15, we plunge more deeply into this topic. Let I be the interval [ − y ∗ , y ∗ ]. Recall that g is defined on I , except at y = 0, and satisfies g ( − y ) = − g ( y ) . From the results of the previous section, we have g ′ ( y ) > √ 2, 0 < g ( y ∗ ) < y ∗ , and − y ∗ < g ( − y ∗ ) < 0. Also, lim y → 0 ± = ∓ y ∗ . Thus the graph of g resembles that shown in Figure 14.10. Note that all points in the interval [ g ( − y ∗ ) , g ( y ∗ ) ] have two preimages, while points in the inter- vals ( − y ∗ , g ( − y ∗ )) and ( g ( y ∗ ) , y ∗ ) have only one. The endpoints of I , namely ± y ∗ , have no preimages in I since g ( 0 ) is undefined. Let y 0 ∈ I . We will investigate the structure of the set A ∩ { y = y 0 } . We define the (forward) orbit of y 0 to be the set ( y 0 , y 1 , y 2 , . . .) , where y n = g ( y n − 1 ) = g n ( y 0 ) . For each y 0 , the forward orbit of y 0 is uniquely determined, though it terminates if g n ( y 0 ) = 0. A backward orbit of y 0 is a sequence of the form ( y 0 , y − 1 , y − 2 . . .) , where g ( y − k ) = y − k + 1 . Unlike forward orbits of g , there are infinitely many distinct backward orbits for a given y 0 except in the case where y 0 = ± y ∗ (since these two points have no preimages in I ). To see this, suppose first that y 0 does not lie on the forward orbit of ± y ∗ . Then each point y − k must have either one or two distinct preimages since y − k 6 = ± y ∗ . If y − k has only one preimage y − k − 1 , − y* y* g ( − y* ) g ( y* ) Figure 14.10 Graph of the one-dimensional function g on I = [ − y ∗ , y ∗ ].

--- PDF PAGE 338 ---
14.5 The Chaotic Attractor 323 then y − k lies in either ( − y ∗ , g ( − y ∗ ) or ( g ( y ∗ ) , y ∗ ) . But then the graph of g shows that y − k − 1 must have two preimages, so no two consecutive points in a given backward orbit can have only one preimage. This shows that y 0 has infinitely many distinct backward orbits. If we happen to have y − k = ± y ∗ for some k > 0, then this backward orbit stops since ± y ∗ has no preimage in I . However, y − k + 1 must have two preim- ages, one of which is the endpoint and the other is a point in I that does not equal the other endpoint. Thus we can continue taking preimages of this second backward orbit as before, thereby generating infinitely many distinct backward orbits as in the preceding. We claim that each of these infinite backward orbits of y 0 corresponds to a unique point in A ∩ { y = y 0 } . To see this, consider the line J − k given by y = y − k in R . Then 8 k ( J − k ) is a closed subinterval of J 0 for each k . Note that 8( J − k − 1 ) ⊂ J − k , since 8( y = y − k − 1 ) is a proper subinterval of y = y − k . Thus the nested intersection of the sets 8 k ( J − k ) is nonempty, and any point in this intersection has backward orbit ( y 0 , y − 1 , y − 2 , . . .) by construction. Further- more, the intersection point is unique, since each application of 8 contracts the intervals y = y − k by a factor of c < 1. In terms of our model, we therefore see that the attractor A is a complicated set. We have proved this proposition. Proposition. The attractor A for the model Lorenz system meets each of the lines y = y 0 6 = y ∗ in R at infinitely many distinct points. In forward time all of the solution curves through each point on this line either 1. Meet the line y = 0 , in which case the solution curves all tend to the equilibrium point at ( 0, 0, 0 ) , or 2. Continually reintersect R, and the distances between these intersection points on the line y = y k tends to 0 as time increases.  Now we turn to the dynamics of 8 in R . We first discuss the behavior of the one-dimensional function g , and then use this information to understand what happens for 8 . Given any point y 0 ∈ I , note that nearby forward orbits of g move away from the orbit of y 0 since g ′ > √ 2. More precisely, we have the following proposition. Proposition. Let 0 < ν < y ∗ . Let y 0 ∈ I = [ − y ∗ , y ∗ ] . Given any  > 0 , we may find u 0 , v 0 ∈ I with | u 0 − y 0 | <  and | v 0 − y 0 | <  and n > 0 such that | g n ( u 0 ) − g n ( v 0 ) | ≥ 2 ν . Proof: Let J be the interval of length 2  centered at y 0 . Each iteration of g expands the length of J by a factor of at least √ 2, so there is an iteration for which g n ( J ) contains 0 in its interior. Then g n + 1 ( J ) contains points arbitrarily

--- PDF PAGE 339 ---
324 Chapter 14 The Lorenz System close to both ± y ∗ , and thus there are points in g n + 1 ( J ) where the distance from each other is at least 2 ν . This completes the proof.  Let’s interpret the meaning of this proposition in terms of the attractor A . Given any point in the attractor, we may always find points arbitrarily nearby with forward orbits that move apart just about as far as they possibly can. This is the hallmark of a chaotic system: We call this behavior sensitive dependence on initial conditions. A tiny change in the initial position of the orbit may result in drastic changes in the eventual behavior of the orbit. Note that we must have a similar sensitivity for the flow in A ; certain nearby solution curves in A must also move far apart. This is the behavior we witnessed earlier in Figure 14.2. This should be contrasted with the behavior of points in A that lie on the same line y = constant with − y ∗ < y < y ∗ . As we saw before, there are infinitely many such points in A . Under iteration of 8 , the successive images of all of these points move closer together rather than separating. Recall now that a subset of I is dense if its closure is all of I . Equivalently, a subset of I is dense if there are points in the subset arbitrarily close to any point whatsoever in I . Also, a periodic point for g is a point y 0 for which g n ( y 0 ) = y 0 for some n > 0. Periodic points for g correspond to periodic solutions of the flow. Proposition. The periodic points of g are dense in I. Proof: As in the proof in the last section that A is an attractor, given any subinterval J of I − { 0 } , we may find n so that g n maps some subinterval J ′ ⊂ J in one-to-one fashion over either ( − y ∗ , 0 ) or ( 0, y ∗ ) . Thus either g n ( J ′ ) con- tains J ′ or the next iteration, g n + 1 ( J ′ ) , contains J ′ . In either case, the graphs of g n or g n + 1 cross the diagonal line y = x over J ′ . This yields a periodic point for g in J .  Now let us interpret this result in terms of the attractor A . We claim that periodic points for 8 are also dense in A . To see this, let P ∈ A and U be an open neighborhood of P . We assume that U does not cross the line y = 0 (otherwise just choose a smaller neighborhood nearby that is disjoint from y = 0). For small enough  > 0, we construct a rectangle W ⊂ U centered at P and having width 2  (in the x -direction) and height  (in the y -direction). Let W 1 ⊂ W be a smaller square centered at P with sidelength / 2. By the transitivity result of the previous section, we may find a point Q 1 ∈ W 1 such that 8 n ( Q 1 ) = Q 2 ∈ W 1 . By choosing a subset of W 1 if necessary, we may assume that n > 4 and furthermore that n is so large that c n < / 8. It follows that the image of 8 n ( W ) (not 8 n ( W 1 ) ) crosses through the interior of W nearly vertically and extends beyond its top and bottom boundaries, as shown in Figure 14.11. This fact uses the hyperbolicity condition.

--- PDF PAGE 340 ---
14.5 The Chaotic Attractor 325 W W 1 Q 1 Q 2 P y = c 0 Φ ( y = c 0 ) Φ n ( W ) Figure 14.11 8 maps W across itself. Now consider the lines y = c in W . These lines are mapped to other such lines in R by 8 n . Since the vertical direction is expanded, some of the lines must be mapped above W and some below. It follows that one such line y = c 0 must be mapped inside itself by 8 n , and therefore there must be a fixed point for 8 n on this line. Since this line is contained in W , we have produced a periodic point for 8 in W . This proves density of periodic points in A . In terms of the flow, a solution beginning at a periodic point of 8 is a closed orbit. Thus the set of points lying on closed orbits is a dense subset of A . The structure of these closed orbits is quite interesting from a topological point of view, as many of these closed curves are actually “knotted.” See Birman and Williams [10] and exercise 10 at the end of this chapter. Finally, we say that a function g is transitive on I if, for any pair of points y 1 and y 2 in I and neighborhoods U i of y i , we can find ˜ y ∈ U 1 and n such that g n ( ˜ y ) ∈ U 2 . Just as in the proof of density of periodic points, we may use the fact that g is expanding in the y -direction to prove the following. Proposition. The function g is transitive on I. We leave the details to the reader. In terms of 8 , we almost proved the corre- sponding result when we showed that A was an attractor. The only detail we did not provide was the fact that we could find a point in A with an orbit that made the transit arbitrarily close to any given pair of points in A. For this detail, we refer to Exercise 4 at the end of this chapter. Thus we can summarize the dynamics of 8 on the attractor A of the Lorenz model as follows.  Theorem. (Dynamics of the Lorenz Model) The Poincar´ e map 8 restricted to the attractor A for the Lorenz model has the following properties: 1. 8 has sensitive dependence on initial conditions 2. Periodic points of 8 are dense in A 3. 8 is transitive on A