odechap4

Chapter 4 Nonlinear Systems of Di®erential Equations Introduction This ¯nal chapter brings to bear all the material presented before and pushes on to the heart of the subject, nonlinear systems of di®erential equations. Section 1 begins with a demonstration of existence and uniqueness (for t close to t0) of solutions to (0.1)dx dt=F(t; x); x(t0) =x0: Here x(t) is a path in ­ ½RnandFis bounded and continuous on I£­ (with t02I), and satis¯es a Lipschitz condition in x. (See (1.2) for a de¯nition.) We study the issue of global existence, including positive results when F(t; x) is linear in x. Section 2 studies the smoothness of the solution to (0.1) as a function of x0, given various additional hypotheses on F, and related issues. Section 3 reveals a geometric °avor to (0.1), described in the language of vector ¯elds and the °ows they generate. A vector ¯eld on ­½Rnis a map F: ­!Rn. This is a special case of (0.1), where Fis independent of t. The path x(t) in ­ satisfying (0.1) for such Fis called the orbit ofFthrough x0; denote it ©t(x0). This gives rise to the family of maps ©t, called the °ow generated by F. The phase portrait is introduced as a tool to understand the orbits and °ow, from a visual perspective. We pay particular attention to how phase portraits look near critical points of a vector ¯eld F(which are points where Fvanishes), including special types known as sources, sinks, saddles, and centers. 229 230 4. Nonlinear Systems of Di®erential Equations Section 4 discusses a particular class of vector ¯elds, gradient vector ¯elds, on a domain ­ ½Rn. In case n= 2, this relates to the topic of exact equations, discussed in many texts early on. We have broken with tradition and moved the discussion of exactness to here, to see it in a broader context. We move from generalities about nonlinear systems to settings in which they arise. Section 5 introduces a class of di®erential equations arising from Newton's law F=ma. This resumes the study introduced in x5 of Chapter 1. This time we are studying the interaction of several bodies, each moving inn-dimensional space. We concentrate on central force problems. We show how a two-body central force problem (for motion in Rn) gives rise to a second order n£nsystem, in \center of mass coordinates." We look at this two-body problem in more detail in x6, and derive Newton's epoch-making analysis of the planetary motion problem. Inx7 we introduce another (though ultimately related) class of problems that lead to di®erential equations, namely variational problems. The general setup is to consider (0.2) I(u) =Zb aL(u(t); u0(t))dt; for paths u: [a; b]!­½Rn, given smooth Lon ­£Rn, and ¯nd conditions under which Ihas a minimum, or maximum, or more generally a stationary point u. We produce a di®erential equation known as the Lagrange equa- tion for u. This method has many important rami¯cations. One of the most important is to produce di®erential equations for physical problems, providing an alternative to the method discussed in x5. We illustrate this inx7 by obtaining a derivation of the pendulum equation, alternative to that given in x6 of Chapter 1. We proceed to more sophisticated uses of the variational method. In x8 we discuss the \brachistochrone problem," tossed about by the early leading lights of calculus, one of the foundational varia- tional problems. In x9 we discuss the double pendulum, a physical problem that is confounding when one uses the F=maapproach, and which well illustrates the \Lagrangian" approach. An alternative to Lagrangian dif- ferential equations is the class of Hamiltonian di®erential equations. The passage from Lagrangian to Hamiltonian equations is previewed (in special cases) in xx7 and 9, and developed further in x10. The majority of the systems studied in this chapter are not amenable to solution in terms of explicit formulas. In x11 we introduce a tool that has revolutionized the study of these equations, namely numerical approxima- tion. Behind this revolution is the availability of personal computers. In x11 we present several techniques that allow for accurate approximation of so- lutions to (0.1), the most important being Runge-Kutta di®erence schemes. 1. Existence and uniqueness of solutions 231 Inx12 we return to the study of qualitative features of phase portraits, initiated in x3. We de¯ne limit sets of orbits, and establish a result known as the Poincar¶ e-Bendixson theorem, which provides a condition under which a limit set for an orbit of a planar vector ¯eld can be shown to be a closed curve, called a limit cycle. Sections 13{14 are devoted to some systems of di®erential equations arising to model the populations of interacting species. In x13 we study \predator-prey" equations. We study several models. In some, all the or- bits are periodic, except for one critical point. In others, there is a limit cycle, arising via the mechanism examined in x12. In x14 we look at other interacting species equations, namely equations modeling competing species. One phenomenon behind the Poincar¶ e-Bendixson theorem is that an orbit of a vector ¯eld Fin the plane locally divides the plane into two parts, one to the left of the orbit and one to the right. Since another orbit of F cannot cross it, this tends to separate the plane into pieces, in each of which the phase portrait has a fairly simple appearance. In dimension three and higher, this mechanism to enforce simplicity does not work, and far more complicated scenarios are possible. This leads to the occurrence of \chaos" forn£nsystems of di®erential equations when n¸3. We explore some aspects of this in the last section of this chapter, x15. This chapter has six appendices. In Appendix A we give basic informa- tion on the derivative of functions of several variables, reviewing material typically covered in third semester calculus and setting up notation that is used in the chapter. Appendix B discusses some basic results about conver- gence, including the notion of compactness . In Appendix C we show that if the linearization of a vector ¯eld Fat a critical point behaves like a saddle, so does F. Appendix D discusses an approximation procedure for comput- ing the periods of orbits, for a certain family of planar vector ¯elds, with reference to how Einstein's correction of Newton's equations for planetary motion yields a calculation of the precession of the planet's perihelion. In Appendix E we show that a spherically symmetric planet produces the same gravitational ¯eld as if all its mass were concentrated at its center. In Ap- pendix F we prove the Brouwer ¯xed-point theorem (in dimension 2), a use of which arises in x15. The proof we give makes use of material developed inx4. 1. Existence and uniqueness of solutions We investigate existence and uniqueness of solutions to a ¯rst order nonlinear n£nsystem of di®erential equations, (1.1)dx dt=F(t; x); x(t0) =x0: 232 4. Nonlinear Systems of Di®erential Equations We assume Fis bounded and continuous on I£­, where Iis an open interval about t0and ­ is an open subset of Rn, containing x0. We also assume Fsatis¯es a Lipschitz condition in x: (1.2) kF(t; x)¡F(t; y)k ·Lkx¡yk; for all t2I; x; y 2­, with L2(0;1). Such an estimate holds if ­ is convex and FisC1inxand satis¯es (1.3) kDxF(t; x)k ·L; for all t2I; x2­. At this point, the reader might want to review the concept of the derivative of a function of nvariables, by looking in Appendix A. The estimate (1.3) follows readily from (A.9). Our ¯rst goal is to prove the following. Proposition 1.1. Assume F:I£­!Rnis bounded and continuous and satis¯es the Lipschitz condition (1.2), and let x02­. Then there exists T0>0and a unique C1solution to (1.1) for jt¡t0j< T0. The ¯rst step in proving this is to rewrite (1.1) as an integal equation: (1.4) x(t) =x0+Zt t0F(s; x(s))ds: The equivalence of (1.1) and (1.4) follows from the Fundamental Theorem of Calculus. It su±ces to ¯nd a continuous solution xto (1.4) on [ t0¡T0; t0+ T0], since then the right side of (1.4) will be C1int. We will apply a technique known as Picard iteration to construct a solution to (1.4). We set x0(t)´x0and then de¯ne xn(t) inductively by (1.5) xn+1(t) =x0+Zt t0F(s; xn(s))ds: We show that this converges uniformly to a solution to (1.4), for jt¡t0j ·T0, ifT0is taken small enough. To get this, we quantify some hypotheses made above. We assume (1.6) BR(x0) =fx2Rn:kx¡x0k ·Rg ½­ and (1.7) kF(s; x)k ·M;8s2I; x2BR(x0): 1. Existence and uniqueness of solutions 233 Clearly x0(t)´x0takes values in BR(x0) for all t. Suppose that xn(t) has been constructed, taking values in BR(x0), and xn+1(t) is de¯ned by (1.5). We have (1.8) kxn+1(t)¡x0k ·Zt t0kF(s; xn(s))kds·Mjt¡t0j; soxn+1(t) also takes values in BR(x0) provided jt¡t0j ·T0and (1.9) T0·R M: As long as (1.9) holds and [ t0¡T0; t0+T0]½I, we get an in¯nite sequence xn(t) of functions, related by (1.5). We produce one more constraint on T0, which will guarantee conver- gence. Note that, for n¸1, (1.10)kxn+1(t)¡xn(t)k=°°°Zt t0£ F(s; xn(s))¡F(s; xn¡1(s))¤ ds°°° ·Zt t0kF(s; xn(s))¡F(x; xn¡1(s))kds ·LZt t0kxn(s)¡xn¡1(s)kds; the last inequality by (1.2). Hence (1.11) max jt¡t0j·T0kxn+1(t)¡xn(t)k ·LT0max js¡t0j·T0kxn(s)¡xn¡1(s)k: The additional constraint we impose on T0is (1.12) T0·1 2L: Noting that (1.13) max jt¡t0j·T0kx1(t)¡x0k ·R; we see that (1.14) max jt¡t0j·T0kxn+1(t)¡xn(t)k ·2¡nR: Consequently, the in¯nite series (1.15) x(t) =x0+1X n=0¡ xn+1(t)¡xn(t)¢ 234 4. Nonlinear Systems of Di®erential Equations is absolutely and uniformly convergent for jt¡t0j ·T0, with a continuous sum, satisfying (1.16) max jt¡t0j·T0kx(t)¡xn(t)k ·21¡nR: It readily follows that (1.17)Zt t0F(s; xn(s))ds¡!Zt t0F(s; x(s))ds; so (1.4) follows from (1.5) in the limit n! 1 . To ¯nish the proof of Proposition 1.1, we establish uniqueness. Suppose y(t) also satis¯es (1.4) for jt¡t0j ·T0. Then kx(t)¡y(t)k=°°°Zt t0£ F(s; x(s))¡F(s; y(s))¤ ds°°° ·Zt t0kF(s; x(s))¡F(s; y(s))kds ·LZt t0kx(s)¡y(s)kds; and hence (1.18) max jt¡t0j·T0kx(t)¡y(t)k ·T0Lmax js¡t0j·T0kx(s)¡y(s)k: As long as (1.12) holds, T0L·1=2, so (1.18) clearly implies max jt¡t0j·T0kx(t)¡ y(t)k= 0, which gives the asserted uniqueness. Note that the Lipschitz hypothesis (1.2) was needed only for x; y2 BR(x0). Thus we can extend Proposition 1.1 to the following setting: (1.19)For each closed, bounded K½­;there exists LK<1such that kF(t; x)¡F(t; y)k ·LKkx¡yk;8x; y2K; t2I: We can also replace the bound on Fby (1.20)For each Kas above, there exists MK<1such that kF(t; x)k ·MK;8x2K; t2I: Results of Appendix B imply that there exists RK>0 such that eK=[ x2KBRK(x) is a compact subset of ­ : 1. Existence and uniqueness of solutions 235 It follows that for each x02K, the solution to (1.1) exists on the interval (1.21) ft2I:jt¡t0j ·min(RK=MeK;1=2LeK)g: Now that we have local solutions to (1.1), it is of interest to investigate when global solutions exist. Here is an example of breakdown: (1.22)dx dt=x2; x(0) = 1 : Here I=R; n= 1;­ =R, and F(x) =x2is smooth, satisfying the local bounds (1.20){(1.21). The equation (1.22) has the unique solution (1.23) x(t) =1 1¡t; t2(¡1;1); which blows up as t%1. It is useful to know that \blowing up" is the only way a solution can fail to exist globally. We have the following result. Proposition 1.2. LetFbe as in Proposition 1.1, but with the Lipschitz and boundedness hypotheses relaxed to (1.19){(1.20). Assume [a; b]is contained in the open interval Iand assume x(t)solves (1.1) for t2(a; b). Assume there exists a closed, bounded set K½­such that x(t)2Kfor all t2(a; b). Then there exist a1< aandb1> bsuch that x(t)solves (1.1) for t2(a1; b1). Proof. We deduce from (1.21) that there exists ± >0 such that for each x12K; t 12[a; b], the solution to (1.24)dx dt=F(t; x); x(t1) =x1 exists on the interval [ t1¡±; t1+±]. Now, under the current hypotheses, take t12(b¡±=2; b); x1=x(t1), with x(t) solving (1.1). Then solving (1.24) continues x(t) past t=b. Similarly one can continue x(t) past t=a. Here is an example of a global existence result that can be deduced from Proposition 1.2. Consider the 2 £2 system for x= (y; v): (1.25)dy dt=v; dv dt=¡y3: Here we have ­ = R2; F(t; x) =F(t; y; v ) = ( v;¡y3). If (1.25) holds for t2(a; b), we have (1.26)d dt³v2 2+y4 4´ =vdv dt+y3dy dt= 0; 236 4. Nonlinear Systems of Di®erential Equations so each x(t) = (y(t); v(t)) solving (1.25) lies in a level curve y4=4+v2=2 =C, hence is con¯ned to a closed, bounded subset of R2, yielding global existence of solutions to (1.25). We can also apply Proposition 1.2 to establish global existence of solu- tions to linear systems, (1.27)dx dt=A(t)x; x (0) = x0; given A(t) continuous in t2I(an interval about 0), with values in M(n;C). It su±ces to establish the following. Proposition 1.3. IfkA(t)k · Kfort2I, then the solution to (1.27) satis¯es (1.28) kx(t)k ·eKjtjkx0k: Proof. It su±ces to prove (1.28) for t¸0. Then y(t) =e¡Ktx(t) satis¯es (1.29)dy dt=C(t)y; y (0) = x0; with C(t) =A(t)¡K. Hence C(t) satis¯es (1.30) Re ( C(t)u; u)·0;8u2Cn: Then (1.28) is a consequence of the following estimate, of interest in its own right. Lemma 1.4. Ify(t)solves (1.29) and (1.30) holds for C(t), then (1.31) ky(t)k · k y(0)kfort¸0: Proof. We have (1.32)d dtky(t)k2= (y0(t); y(t)) + ( y(t); y0(t)) = 2 Re ( C(t)y(t); y(t)) ·0: Thanks to Proposition 1.3, we have for s; t2I, the solution operator for (1.27), (1.33) S(t; s)2M(n;C); S(t; s)x(s) =x(t); 1. Existence and uniqueness of solutions 237 introduced in x8 of Chapter 3. As noted there, we have the Duhamel formula (1.34) x(t) =S(t; t0) +Zt t0S(t; s)f(s)ds; for the solution to (1.35)dx dt=A(t)x+f(t); x(t0) =x0: IfF(t; x) depends explicitly on t, we call (1.1) a non-autonomous system. IfFdoes not depend explicitly on t, we say (1.1) is autonomous. The following device converts a non-autonomous system to an autonomous one. Take the n£nsystem (1.1). Then the ( n+ 1)£(n+ 1) system (1.36)dx dt=F(y; x);dy dt= 1; x(t0) =x0; y(t0) =t0 has the autonomous form (1.37)dz dt=G(z); z(t0) = (x0; t0); forz= (x; y), with G(z) = (F(y; x);1), and the solution to (1.36) is ( x(t); t), where x(t) solves (1.1). Thus for many purposes it su±ces to consider au- tonomous sytems. To close this section, we note how a higher order n£nsystem, such as (1.38)dkx dtk=F(t; x; : : : ; x(k¡1)); x(t0) =x0; : : : ; x(k¡1)(t0) =xk¡1; can be converted to a ¯rst order nk£nksystem, for (1.39) y=0 @y0... yk¡11 A; y j(t)2Rn;0·j·k¡1: The system is (1.40)dy0 dt=y1; ... dyk¡2 dt=yk¡1; dyk¡1 dt=F(t; y0; : : : ; y k¡1); with initial data (1.41) yj(t0) =xj;0·j·k¡1: Ify(t) solves (1.40){(1.41), then x(t) =y0(t) solves (1.38), and we have (1.42) x(j)(t) =yj(t);0·j·k¡1: Note how this construction is parallel to that done in the linear case in Chapter 3, x3. 238 4. Nonlinear Systems of Di®erential Equations Exercises 1. Apply the Picard iteration method to dx dt=ax; x (0) = 1 ; given a2C. Taking x0(t)´1, show that xn(t) =nX k=0ak k!tk: 2. Discuss the matrix analogue of Exercise 1. 3. Consider the initial value problem dx dt=x2; x(0) = 1 : Take x0´1 and use the Picard iteration method (1.5) to write out xn(t); n = 1;2;3: Compare the results with the formula (1.23). 4. Given A0; A12M(n;C), consider the initial value problem dx dt= (A0+A1t)x; x (0) = x0: Take x0(t)´x0and use the Picard iteration (1.5) to write out xn(t); n = 1;2;3: Compare and contrast the results with calculations from x10 of Chapter 3. 5. Modify the system (1.25) to dy dt=v;dv dt=¡y3¡v: 1. Existence and uniqueness of solutions 239 Show that solutions satisfy d dt³v2 2+y4 4´ ·0; and use this to establish global existence for t¸0. 6. Consider the initial value problem dx dt=jxj1=2; x(0) = 0 : Note that x(t)´0 is a solution, and x(t) =1 4t2; t¸0; 0; t·0 is another solution, on t2(¡1;1). Why does this not contradict the uniqueness part of Proposition 1.1? Can you produce other solutions to this initial value problem? 7. Take ¯2(0;1) and consider the initial value problem dx dt=x¯; x(0) = 1 : Show that this has a solution for all t¸0 if and only if ¯·1. 8. Let F:Rn!RnbeC1and suppose x(t) solves (1.43)dx dt=F(x); x(t0) =x0; fort2I, an open interval containing t0. Show that, for t2I, (1.44)d dtkx(t)k2= 2x(t)¢F(x(t)): Show that, if ® >0 and x(t)6= 0, (1.45)d dtkx(t)k®=®kx(t)k®¡2x(t)¢F(x(t)): 9. In the setting of Exercise 8, suppose Fsatis¯es an estimate (1.46) kF(x)k ·C(1 +kxk)¯;8x2Rn; C < 1; ¯ < 1: Show that there exists ® >0 and K <1such that, if kx(t)k ¸1 for t2I, d dtkx(t)k®·K;8t2I: Use this to establish that the solution to (1.43) exists for all t2R. 240 4. Nonlinear Systems of Di®erential Equations Exercises 10{12 below will extend the conclusion of Exercise 9 to the case ¯= 1 in (1.46). One approach is via the following result, known as Gron- wall's inequality . Proposition 1.5. Assume (1.47) g2C1(R); g0¸0: Letuandvbe real valued, continuous functions on Isatisfying (1.48)u(t)·A+Zt t0g(u(s))ds; v(t)¸A+Zt t0g(v(s))ds: Then (1.49) u(t)·v(t);fort2I; t¸t0: Proof. Setw(t) =u(t)¡v(t). Then (1.50)w(t)·Zt t0£ g(u(s))¡g(v(s))¤ ds =Zt t0M(s)w(s)ds; where (1.51) M(s) =Z1 0g0(¿u(s) + (1 ¡¿)v(s))d¿: Hence we have (1.52) w(t)·Zt t0M(s)w(s)ds; M (s)¸0; M 2C(I); and we claim this implies (1.53) w(t)·0;8t2I; t¸t0: In other words, we claim that w(t)·0 on [ t0; b] whenever [ t0; b]½I. To see this, let t1be the largest number in [ t0; b] with the property that w·0 on [t0; t1]. We claim that t1=b. 1. Existence and uniqueness of solutions 241 Assume to the contrary that t1< b. Noting thatRt1 t0M(s)w(s)ds·0, we deduce from (1.52) that (1.54) w(t)·Zt t1M(s)w(s)ds;8t2[t1; b]: Hence, with (1.56) K= max [t1;b]M(s)<1; we have, for a2(t1; b), (1.57) max [t1;a]w(t)·(a¡t1)Kmax [t1;a]w(s): If we pick a2(t1; b) such that ( a¡t1)K < 1, this implies (1.58) w(t)·0;8t2[t1; a]; contradicting the maximality of t1. Hence actually t1=b, and we have the implication (1.52) )(1.53), completing the proof of Proposition 1.5. 10. Assume v¸0 is a C1function on I= (a; b), satisfying (1.59)dv dt·Cv; v (t0) =v0; where C2(0;1) and t02I. Using Proposition 1.5, show that (1.60) v(t)·eC(t¡t0)v0;8t2[t0; b): 11. In the setting of Exercise 10, avoid use of Proposition 1.5 as follows. Write (1.59) as (1.61)dv dt=Cv¡g(t); v(t0) =v0; g¸0; with solution (1.62) v(t) =eC(t¡t0)v0¡Zt t0eC(t¡s)g(s)ds: Deduce (1.60) from this. 12. Return to the setting of Exercise 8, and replace the hypothesis (1.46) by (1.63) kF(x)k ·C(1 +kxk);8x2Rn: Show that the solution to (1.43) exists for all t2R. Hint. Take v(t) = 1 + kx(t)k2and use (1.44). Show that Exercise 10 (or 11) applies. 242 4. Nonlinear Systems of Di®erential Equations 2. Dependence of solutions on initial data and other parameters We study how the solution to a system of di®erential equations (2.1)dx dt=F(x); x(0) = y depends on the initial condition y. As shown in x1, there is no loss of generality in considering the autonomous system (2.1). We will assume F: ­!Rnis smooth, ­ ½Rnopen and convex, and denote the solution to (2.1) by x=x(t; y). We want to examine smoothness in y. Let DF(x) denote the n£nmatrix valued function of partial derivatives of F. (See Appendix A for more on this derivative.) To start, we assume Fis of class C1, i.e., DFis continuous on ­, and we want to show x(t; y) is di®erentiable in y. Let us recall what this means. Take y2­ and pick R >0 such that BR(y), de¯ned as in (1.6), is contained in ­. We seek an n£nmatrix W(t; y) such that, for w02Rn;kw0k ·R, (2.2) x(t; y+w0) =x(t; y) +W(t; y)w0+r(t; y; w 0); where (2.3) r(t; y; w 0) =o(kw0k); which means (2.4) lim w0!0r(t; y; w 0) kw0k= 0: When this holds, x(t; y) is di®erentiable in y, and (2.5) Dyx(t; y) =W(t; y): In other words, (2.6) x(t; y+w0) =x(t; y) +Dyx(t; y)w0+o(kw0k): In the course of proving this di®erentiability, we also want to produce an equation for W(t; y) =Dyx(t; y). This can be done as follows. Suppose x(t; y) were di®erentiable in y. (We do not yet know that it is, but that is okay.) Then F(x(t; y)) is di®erentiable in y, so we can apply Dyto (2.1). Using the chain rule, we get the following equation, (2.7)dW dt=DF(x)W; W (0; y) =I; 2. Dependence of solutions on initial data and other parameters 243 called the linearization of (2.1). Here, Iis the n£nidentity matrix. Equiv- alently, given w02Rn, (2.8) w(t; y) =W(t; y)w0 is expected to solve (2.9)dw dt=DF(x)w; w (0) = w0: Now, we do not yet know that x(t; y) is di®erentiable, but we do know from results of x1 that (2.7) and (2.9) are uniquely solvable. It remains to show that, with such a choice of W(t; y), (2.2){(2.3) hold. To rephrase the task, set (2.10) x(t) =x(t; y); x 1(t) =x(t; y+w0); z(t) =x1(t)¡x(t); and let w(t) solve (2.9). The task of verifying (2.2){(2.3) is equivalent to the task of verifying (2.11) kz(t)¡w(t)k=o(kw0k): To show this, we will obtain for z(t) an equation similar to (2.9). To begin, (2.10) implies (2.12)dz dt=F(x1)¡F(x); z(0) = w0: Now the fundamental theorem of calculus gives (2.13) F(x1)¡F(x) =G(x1; x)(x1¡x); with (2.14) G(x1; x) =Z1 0DF¡ ¿x1+ (1¡¿)x¢ d¿: IfFisC1, then Gis continuous. Then (2.12){(2.13) yield (2.15)dz dt=G(x1; x)z; z (0) = w0: Given that (2.16) kDF(u)k ·L;8u2­; 244 4. Nonlinear Systems of Di®erential Equations which we have by continuity of DF, after possibly shrinking ­ slightly, we deduce from Proposition 1.3 that (2.17) kz(t)k ·ejtjLkw0k; that is, (2.18) kx(t; y)¡x(t; y+w0)k ·ejtjLkw0k: This establishes that x(t; y) isLipschitz iny. To proceed, since Gis continuous and G(x; x) =DF(x), we can rewrite (2.15) as (2.19)dz dt=G(x+z; x)z=DF(x)z+R(x; z); z(0) = w0; where (2.20) F2C1(­) =) kR(x; z)k=o(kzk) =o(kw0k): Now comparing (2.19) with (2.9), we have (2.21)d dt(z¡w) =DF(x)(z¡w) +R(x; z);(z¡w)(0) = 0 : Then Duhamel's formula gives (2.22) z(t)¡w(t) =Zt 0S(t; s)R(x(s); z(s))ds; where S(t; s) is the solution operator for d=dt¡B(t), with B(t) =G(x1(t); x(t)), which as in (2.17), satis¯es (2.23) kS(t; s)k ·ejt¡sjL: We hence have (2.11), i.e., (2.24) kz(t)¡w(t)k=o(kw0k): This is precisely what is required to show that x(t; y) is di®erentiable with respect to y, with derivative W=Dyx(t; y) satisfying (2.7). Hence we have: Proposition 2.1. IfF2C1(­)and if solutions to (2.1) exist for t2 (¡T0; T1), then, for each such t; x(t; y)isC1iny, with derivative Dyx(t; y) satisfying (2.7). 2. Dependence of solutions on initial data and other parameters 245 We have shown that x(t; y) is both Lipschitz and di®erentiable in y. The continuity of W(t; y) inyfollows easily by comparing the di®erential equations of the form (2.7) for W(t; y) and W(t; y+w0), in the spirit of the analysis of z(t) done above. IfFpossesses further smoothness, we can establish higher di®erentia- bility of x(t; y) inyby the following trick. Couple (2.1) and (2.7), to get a system of di®erential equations for ( x; W): (2.25)dx dt=F(x); dW dt=DF(x)W; with initial conditions (2.26) x(0) = y; W (0) = I: We can reiterate the preceding argument, getting results on Dy(x; W), hence onD2 yx(t; y), and continue, proving: Proposition 2.2. IfF2Ck(­), then x(t; y)isCkiny. Similarly, we can consider dependence of the solution to (2.27)dx dt=F(¿; x); x(0) = y on a parameter ¿, assuming Fsmooth jointly in ( ¿; x). This result can be deduced from the previous one by the following trick. Consider the system (2.28)dx dt=F(z; y);dz dt= 0; x(0) = y; z(0) = ¿: Then we get smoothness of x(t; ¿; y ) jointly in ( ¿; y). As a special case, let F(¿; x) =¿F(x). In this case x(t0; ¿; y) =x(¿t0; y), so we can improve the conclusion in Proposition 2.2 to the following: (2.29) F2Ck(­) =)x2Ckjointly in ( t; y): Exercises 1. Suppose ¿2Rin (2.27). Show that »=@x=@¿ satis¯es d» dt=DxF(¿; x)»+@ @¿F(¿; x); »(0) = 0 : 246 4. Nonlinear Systems of Di®erential Equations 2. Consider the family of di®erential equations for x¿(t), dx dt=x+¿x2; x(0) = 1 : Write down the di®erential equations satis¯ed by »=@x=@¿ and by ´=@2x=@¿2. 3. Let x=x¿(t); y=y¿(t) solve (2.30)dx dt=¡y+¿(x2+y2);dy dt=x; x (0) = 1 ; y(0) = 0 : Knowing smooth dependence on ¿, ¯nd di®erential equations for the coe±cients Xj(t); Yj(t) in power series expansions (2.31)x¿(t) =X0(t) +¿X1(t) +¿2X2(t) +¢¢¢; y¿(t) =Y0(t) +¿Y1(t) +¿2Y2(t) +¢¢¢: Note that X0(t) = cos t; Y 0(t) = sin t. 4. Using the substitution »(t) =¡x(¡t); ´(t) =y(¡t), show that, for ¿ su±ciently small, solutions to (2.30) are periodic in t. 5. Let p(¿) denote the period of the solution to (2.30). Using (2.31), show that p(¿) is smooth in ¿forj¿jsmall. Note that p(0) = 2 ¼. Compute p0(0). Compare results in Appendix C. 6. Suppose yin (2.1) is a critical point of F, i.e., F(y) = 0. Show that (2.7) becomes dW dt=LW; W (0) = I;where L=DF(y); hence F(y) = 0 = )Dyx(t; y) =etL: 3. Vector ¯elds, orbits, and °ows Let ­ ½Rnbe an open set. A vector ¯eld on ­ is simply a map (3.1) F: ­¡!Rn; 3. Vector ¯elds, orbits, and °ows 247 such as encountered in (2.1). We say Fis aCkvector ¯eld if Fis aCkmap. AC1vector ¯eld is said to be smooth. By convention, if we simply call F a vector ¯eld, we mean it is a smooth vector ¯eld. In this section we always assume Fis at least C1. One can also look at time-dependent vector ¯elds (cf. (1.1)), but in this section we restrict attention to the autonomous case. The solution to (2.1), i.e., to (3.2)dx dt=F(x); x(0) = y; will be denoted (3.3) x(t) = ©t F(y): Results of xx1{2 imply that for each closed bounded K½­ there exists an interval I= (¡T0; T1) about 0 such that, for each t2I, (3.4) ©t F:K¡!­; and this is a Ckmap if Fis aCkvector ¯eld. The family of maps ©t Ffrom Kto ­ is called the °owgenerated by F. We have (3.5) ©0 F(y)´y; i.e., ©0 Fis the identity map. We also have (3.6) ©s+t F(y) = ©t F±©s F(y); provided all these maps are well de¯ned. Given y2­, the path (3.7) t7!©t F(y) is called the orbit through y. Another way to state the de¯ning property of ©t Fis that (3.5) holds and (3.8)d dt©t F(x) =F(©t F(x)): We next obtain interesting information on the t-derivative of (3.9) vt(x) =v(©t F(x)); given v2C1 0(­), i.e., vis of class C1and vanishes outside some closed bounded K½­. The chain rule (cf. Appendix A, especially (A.8)) plus (3.8) yields (3.10)d dtvt(x) =F(©t F(x))¢ rv(©t F(x)): 248 4. Nonlinear Systems of Di®erential Equations In particular, (3.11)d dsv(©s F(x))¯¯¯ s=0=F(x)¢ rv(x): Herervis the gradient of v, given by rv= (@v=@x 1; : : : ; @v=@x n)t. A useful alternative formula to (3.10) is (3.12)d dtvt(x) =d dsvt(©s F(x))¯¯¯ s=0 =F(x)¢ rvt(x); the ¯rst equality following from (3.6) and the second from (3.11), with v replaced by vt. One signi¯cant consequence of (3.12), which will lead to the important result (3.17) below, is that, for v2C1 0(­), (3.13)d dtZ ­v(©t F(x))dx=Z ­F(x)¢ rvt(x)dx =¡Z ­divF(x)v(©t F(x))dx: Here div F(x) is the divergence of the vector ¯eld F(x) = (F1(x); : : : ; F n(x))t, de¯ned by (3.14) div F(x) =@F1 @x1(x) +¢¢¢+@Fn @xn(x): The last equality in (3.13) follows by integration by parts, Z ­Fk(x)@vt @xkdx=¡Z ­@Fk @xkvt(x)dx; followed by summation over k. We reiterate the content of (3.13): (3.15)d dtZ ­v(©t F(x))dx=¡Z ­divF(x)v(©t F(x))dx: So far, we have (3.15) for v2C1 0(­). We can extend this by noting that (3.15) implies (3.16)Z ­v(©t F(x))dx¡Z ­v(x)dx =¡Zt 0Z ­divF(x)v(©s F(x))dx ds: 3. Vector ¯elds, orbits, and °ows 249 Basic results on the integral allow one to pass from v2C1 0(­) in (3.16) to more general v, including v=ÂB(the characteristic function of B, de¯ned to be equal to 1 on Band 0 on ­ nB), for smoothly bounded closed B½­, amongst other functions. In more detail, if B½­ is a smoothly bounded, closed set, let B±=fx2 Rn: dist( x; B)·±g. There exists ±0>0 such that B±½­ for ±2(0; ±0]. For such ±, one can produce v±2C1 0(­) such that v±= 1 on B; 0·v±·1; v ±= 0 onRnnB±: Then ¯¯¯Z ÂB(x)dx¡Z v±(x)dx¯¯¯·vol(B±nB)!0;as±!0; so, as ±!0, Z ­v±(x)dx¡!Z ­ÂB(x)dx: Similar arguments give Z ­v±(©t F(x))dx¡!Z ­ÂB(©t F(x))dx; and Zt 0Z ­divF(x)v±(©s F(x))dx ds¡!Zt 0Z ­divF(x)ÂB(©s F(x))dx ds: These results allow one to take v=ÂBin (3.16). Now one can pass from (3.16) back to (3.15), via the fundamental the- orem of calculus. Note that Vol ©t F(B) =Z ÂB(©¡t F(x))dx: We can apply (3.15) with treplaced by ¡t, and vbyÂB, and deduce the following. Proposition 3.1. IfFis aC1vector ¯eld, generating the °ow ©t F, well de¯ned on ­fort2I, and B½­is smoothly bounded, then, for t2I, (3.17)d dtVol ©t F(B) =Z ©t F(B)divF(x)dx: 250 4. Nonlinear Systems of Di®erential Equations This result is behind the notation div F, i.e., the divergence ofF. Vector ¯elds Fwith positive divergence generate °ows ©t Fthat magnify volumes as tincreases, while vector ¯elds with negative divergence generate °ows that shrink volumes as tincreases. We say the °ow generated by a vector ¯eld Fiscomplete provided ©t F(y) is de¯ned for all t2R; y2­. We say it is forward complete if ©t F(y) is de¯ned for all t2[0;1); y2­. The °ow is backward complete if ©t F(y) is de¯ned for all t2(¡1;0]; y2­. Here is an occasionally useful criterion for forward completeness. Proposition 3.2. LetFbe a C1vector ¯eld on ­ =Rn. Assume there exists R <1and a function V2C1(Rn)such that (3.18) V(x)!+1askxk ! 1 and (3.19) kxk ¸R=) rV(x)¢F(x)·0: Then the °ow ©t Fis forward complete. Proof. Letx(t) = ©t F(x0) be an orbit, de¯ned for t2I, some interval about 0. Then (3.20) kx(t)k ¸R=)d dtV(x(t)) =rV(x(t))¢F(x(t))·0: Hence, for t2I; t¸0; x(t) is con¯ned to the closed bounded set (3.21)n x2Rn:V(x)·maxV(y); y2BR(0)[ fx0go : From here, Proposition 1.2 yields forward completeness. One way to display the behavior of the °ow generated by a vector ¯eld Fon a domain ­ is to draw a \phase portrait." This consists of graphs of selected integral curves of F, with arrows indicating the direction of F along each integral curve. Such portraits are particularly revealing when dim ­ = 2, and also of considerable use when dim ­ = 3. As an example, consider Fig. 3.1, the phase portrait of the °ow associated to the 2 £2 system (3.22)dµ dt=Ã; dà dt=¡g `sinµ; 3. Vector ¯elds, orbits, and °ows 251 Figure 3.1 which arises from the pendulum equation (cf. Chapter 1, (6.6)) (3.23)d2µ dt2+g `sinµ= 0; by adding the variable Ã=dµ=dt . Here g; ` > 0. The system (3.22) has the form (3.2) with x= (µ; Ã) and (3.24) F(µ; Ã) =µÃ ¡g `sinµ¶ : Note that Fig. 3.1 looks like Fig. 6.2 of Chapter 1, except that here we have added arrows, to indicate the direction of the °ow. As noted in Chapter 1, the orbits of this °ow are level curves of the function (3.25) E(µ; Ã) =Ã2 2¡g `cosµ; since if ( µ(t); Ã(t)) solves (3.22), (3.26)d dtE(µ; Ã) =ÃÃ0+g `(sinµ)µ0= 0: It is instructive to expand on this last calculation. In general, if ( µ0; Ã0) = F(µ; Ã), (3.27)d dtE(µ; Ã) =rE(µ; Ã)¢F(µ; Ã);where rE(µ; Ã) =µ@E=@µ @E=@ö : 252 4. Nonlinear Systems of Di®erential Equations Figure 3.2 Now the formula (3.24) gives (3.28) F(µ; Ã) =¡JrE(µ; Ã); where (3.29) J=µ0¡1 1 0¶ ; so the vanishing of dE(µ; Ã)=dtfollows from (3.27){(3.28) and the skew- symmetry of J, which implies (3.30) v¢Jv= 0;8v2R2: A vector ¯eld of the form (3.28) is a special case of a Hamiltonian vector ¯eld, a class of vector ¯elds that will be discussed further in xx5, 7, and 10. We mention some noteworthy features of the phase portrait in Fig. 3.1, features to look for in other such portraits. First, there are the critical points of F, i.e., the points where Fvanishes. In case (3.24), the set of critical points in f(k¼;0) :k2Zg: Fig. 3.1 indicates di®erent natures of the orbits near these critical points, depending on whether kis even or odd. For keven, the orbits near ( k¼;0) consist of closed curves. We say these critical points are centers ; cf. Fig. 3.2. Forkodd, the orbits near p= (k¼;0) consist of curves of the following 3. Vector ¯elds, orbits, and °ows 253 Figure 3.3 nature: (3.31)(a) two orbits that approach past!+1; (b) two orbits that approach past! ¡1 ; (c) orbits that miss p;looking like saddles. We say these critical points are saddles . Cf. Fig. 3.3. Sometimes one calls them hyperbolic critical points. Considerable insight is obtained from the study of the linearization ofF at each critical point. Generally, if Fis aC1vector ¯eld on ­ ½Rn; x02­, andF(x0) = 0, the linearization of Fatx0is given by (3.32) L=DF(x0)2 L(Rn): This construction extends the notion of linearization given in x8 of Chapter 1. We expect that ©t F(x0+y)¼x0+etLy; forkyksmall. Cf. Exercise 6 of x2 (but mind the change in notation). Going further, we expect some important qualitative features of the °ow ©t Fnear x0to be captured by the behavior of etL, and this is born out, with some exceptions. If DF(x0) has zero as an eigenvalue (we say x0is a degenerate critical point) this approximation is not typically useful. It has a better chance if det DF(x0)6= 0. We then say x0is a nondegenerate critical point forF. 254 4. Nonlinear Systems of Di®erential Equations In case Fis given by (3.24), with critical points at pk= (k¼;0), we have (3.33) L0=DF(0;0) =µ0 1 ¡g `0¶ : The eigenvalues of this matrix are §ip g=`, and the orbits of etL0are ellipses, with qualitative features like Fig. 3.2, a center. Meanwhile, (3.34) L1=DF(§¼;0) =µ0 1 g `0¶ : The eigenvalues of this matrix are §p g=`, with corresponding eigenvectors (1;§p g=`)t, and the orbit structure for etL1has qualitative features like Fig. 3.3, a saddle. In general, if Fis a planar vector ¯eld with a nondegenerate critical point at x0, and if all the eigenvalues of DF(x0) are purely imaginary, F itself might not have a center at x0, i.e., the orbits of Fnear x0might not be closed orbits surrounding x0. Here is an example. Take (3.35) F(x) =Jx¡ kxk2x; x 2R2; with Jas in (3.29). Then x0= 0 is a critical point, and DF(0) = J. Thus the linearization has a center. However, if x(t) is an orbit for this vector ¯eld, then (3.36)d dtkx(t)k2= 2x¢x0 = 2x¢(Jx¡ kxk2x) =¡2kxk4; i.e.,½(t) =kx(t)k2satis¯es (3.37)d½ dt=¡2½2: This is separable and we have (3.38) ½(0) = ½0=)½(t) =½0 1 + 2 t½0!0 as t%+1; so the orbits of this vector ¯eld spiral into the origin as t%+1, though much more slowly than they do in the case of spiral sinks, a type of critical point that we will encounter shortly. 3. Vector ¯elds, orbits, and °ows 255 Despite the existence of such examples as (3.35), the fact that (0 ;0) is a center for F, given by (3.24), is no accident, but rather a consequence of the fact that Fhas the form (3.28), (3.39) F(x) =¡JrE(x); so that, as derived in (3.27){(3.30), orbits of Flie on level curves of E. Generally, if Eis a smooth real-valued function on a planar domain ­ ½R2 and the vector ¯eld Fis given by (3.39), (nondegenerate) critical points of Fand (nondegenerate) critical points of Ecoincide. If x02­ is such a point (3.40) DF(x0) =¡JD2E(x0); where D2E(x0) is the matrix of second-order partial derivatives of Eatx0, i.e., D2E=µ@2E=@µ2@2E=@Ã@µ @2E=@µ@à @2E=@Ã2¶ : We recall the following result, established in basic multivariable calculus. Letx0be a nondegenerate critical point of E, soD2E(x0) is an invertible, real symmetric matrix. Then (3.41)D2E(x0) positive de¯nite , E has a local minimum at x0; D2E(x0) negative de¯nite , E has a local maximum at x0; D2E(x0) inde¯nite , E has a saddle at x0; We also note that, whenever A2M(2;R) is symmetric and invertible, (3.42)Apositive de¯nite ,detA >0 and Tr A >0; Anegative de¯nite ,detA >0 and Tr A <0; Ainde¯nite ,detA <0: Furthermore, if Ais such a matrix and (3.43) B=¡JA; then (3.44) det B= det A; and, for such B2M(2;R), (3.45)Bhas 2 real eigenvalues of opposite signs ,detB <0; Bhas 2 purely imaginary eigenvalues ,detB >0 and Tr B= 0; 256 4. Nonlinear Systems of Di®erential Equations Figure 3.4 Putting these observations together (cf. also Exercise 5 below), we have: Proposition 3.3. LetEbe a smooth function on ­½R2, with a nonde- generate critical point at x0. Let Fbe given by (3.39). Then (3.46)DF(x0)has 2 purely imaginary eigenvalues , E has a local max or local min at x0; and (3.47)DF(x0)has 2 real eigenvalues of opposite sign , E has a saddle at x0: We move on to the 2 £2 system (3.48)dµ dt=Ã; dà dt=¡® mág `sinµ; which arises from the damped pendulum equation (cf. Chapter 1, (7.6)), (3.49)d2µ dt2+® mdµ dt+g `sinµ= 0; by adding the variable Ã=dµ=dt . Here g; `; ®; m > 0. The system (3.48) has the form (3.2) with x= (µ; Ã) and (3.50) F(µ; Ã) =µÃ ¡® mág `sinµ¶ : The phase portrait for this system is illustrated in Fig. 3.4. We compare and contrast this portrait with that depicted in Fig. 3.1. 3. Vector ¯elds, orbits, and °ows 257 To start, the vector ¯eld (3.50) has the same critical points as the ¯eld given by (3.24), namely f(k¼;0) : k2Zg. The ¯rst striking di®erence is in the behavior near the critical points ( k¼;0) with keven. Fig. 3.4 depicts orbits spiraling into these critical points, as opposed to the picture in Fig. 3.1 of closed orbits circling these critical points. Let us consider the linearizations about these critical points. For Fas in (3.50), we have (3.51) L=DF(0;0) =µ0 1 ¡g `¡® m¶ ; with characteristic polynomial ¸(¸+®=m) +g=`, hence with eigenvalues (3.52) ¸§=¡® 2m§r ®2 m2¡4g `: There are three cases: Case I. ®2=m2<4g=`. Then ¸§are complex conjugates, each with real part¡®=2m. Case II. ®2=m2= 4g=`. Then ¸+=¸¡=¡®=2m. Case III. ®2=m2>4g=`. Then ¸+and¸¡are distinct real numbers, each negative. In all three cases, we have etLv!0 ast%+1, for each v2R2. In Case I, there is also spiraling, and the orbits look like those in Fig. 3.5(a). Fig. 3.4 depicts such behavior. In Case III, the orbits look like those in Fig. 3.5(c). In Case II, the orbits look like a cross between Fig. 3.5(b) and Fig. 3.5(c). These critical points are all called sinks. (Reverse the sign on F, and the associated orbits are called sources ; cf. Fig. 3.6.) The three cases described above correspond to damped oscillatory, critically damped, and overdamped motion, as discussed in x9 of Chapter 1. Further information on the nature of these orbits spiraling in toward these sinks can be obtained from a computation of the rate of change along the orbits of E(µ; Ã), given by (3.25), i.e., (3.53) E(µ; Ã) =Ã2 2¡g `cosµ: 258 4. Nonlinear Systems of Di®erential Equations Figure 3.5 Figure 3.6 This time, instead of (3.26), we have (3.54)d dtE(µ; Ã) =ÃÃ0+g `(sinµ)µ0 =¡Ã³® mÃ+g `sinµ´ +g `(sinµ)à =¡® mÃ2: While this calculation applies nicely to the problem at hand, it is useful to note the following general phenomenon. Proposition 3.4. LetFbe a smooth vector ¯eld on ­½Rn, with a critical point at x02­. Assume (3.55) all the eigenvalues of DF(x0)have negative real part. Then there exists ± >0such that (3.56) kx¡x0k ·±=)lim t!+1©t Fx=x0: 3. Vector ¯elds, orbits, and °ows 259 To prove this, we bring in the following linear algebra result. Lemma 3.5. LetL2M(n;R)and assume all the eigenvalues of Lhave real part<0. Then there exists a (symmetric, positive de¯nite) inner product h;ionRnand a positive constant Ksuch that (3.57) hLv; vi · ¡ Khv; vi;8v2Rn: We show how Lemma 3.5 allows us to prove Proposition 3.4. Apply the lemma to L=DF(x0). Note that there exist a; b2(0;1) such that (3.58) akvk2· hv; vi ·bkvk2;8v2Rn; where hv; viis as in (3.57) and, as usual, kvk2=v¢v. Since Fis smooth, (3.59) F(x0+y) =Ly+R(y); with Rsmooth on a ball about 0 and DR(0) = 0. Hence (3.60) kR(y)k ·Ckyk2·C0hy; yi: Fory(t) = ©t F(x0+y0)¡x0, we have (3.61)d dthy(t); y(t)i= 2hy0(t); y(t)i = 2hF(x0+y); yi = 2hLy; yi+ 2hR(y); yi: Now (3.57) applies to the ¯rst term in the last line of (3.61), while Cauchy's inequality plus (3.60) yields (3.62)jhR(y); yij · h R(y); R(y)i1=2hy; yi1=2 ·Chy; yi3=2: Hence (3.63)d dthy; yi · ¡ 2Khy; yi+Chy; yi3=2 · ¡Khy; yi; the last inequality holding provided hy; yi1=2·K=C . As long as ±in (3.56) is small enough that fx2Rn:kx¡x0k ·±gis contained in ­ and kvk · 260 4. Nonlinear Systems of Di®erential Equations ±) hv; vi1=2·K=C , ifx=x0+y0andky0k ·±, then (3.63) holds for y(t) = ©t F(x0+y0)¡x0for all t¸0, and yields (3.64) hy(t); y(t)i ·e¡Kthy0; y0i; which in turn gives (3.56). We now prove Lemma 3.5. As shown in x8 of Chapter 2, Cnhas a basis fv1; : : : ; v ngwith respect to which Lis upper triangular, i.e., (3.65) Lvj=¸jvj+X k<jajkvk: Alternatively, Appendix B of Chapter 2 shows that Cnhas an orthonormal basisfvjgfor which (3.65) holds. The eigenvalues of Lare¸j, so by hy- pothesis there exists K12(0;1) such that Re ¸j· ¡K1for all j. Now if we take " >0 and set wj="jvj, we get (3.66) Lwj=¸jwj+X k<j"j¡kajkwk: Then setting (3.67)DX ajwj;X bkwkE = ReX ajbj de¯nes a positive de¯nite inner product (depending on " >0) onCn, hence by restriction on Rn, and if " >0 is taken su±ciently small, the desired conclusion (3.57) follows, with K=K1=2, from (3.66). Having discussed the critical points of the vector ¯eld (3.50) at (0 ;0) and related issues, we now consider the critical points at ( §¼;0). We have (3.68) DF(§¼;0) =µ0 1 g `¡® m¶ : This matrix has eigenvalues (3.69) ¸§=¡® 2m§r ®2 m2+4g `; one positive and one negative. These critical points are saddles. The orbits near these critical points have a behavior such as described in (3.31). Unlike the case of Fgiven by (3.24), where the orbits are level curves of E, the proof of this is more subtle in the present situation. See Appendix C for a proof. 3. Vector ¯elds, orbits, and °ows 261 Having studied the various critical points depicted in Figs. 3.1 and 3.4, we point out some special orbits that appear in these phase portraits, namely orbits connecting two critical points. Generally, if Fis aC1vector ¯eld on ­½Rnwith critical points p1; p22­, an orbit x(t) of ©t Fsatisfying (3.70) lim t!¡1x(t) =p1;lim t!+1x(t) =p2 is called a heteroclinic orbit , from p1top2, ifp16=p2. Ifp1=p2, such an orbit is called a homoclinic orbit . In Fig. 3.1, we see heteroclinic orbits connecting p1= (¡¼;0) and p2= (¼;0), one from p1top2and one from p2 top1. These lie on level curves where E(µ; Ã) =g=`. Such a heteroclinic orbit describes the motion of a pendulum that is heading towards pointing vertically upward. As time goes on, the pendu- lum ascends more and more slowly, never quite reaching the vertical position. With a little less energy, the pendulum would stop a bit short of vertical and fall back, swinging back and forth. With a little more energy, the pendulum would swing past the vertical position. Recall that Fig. 3.1 portrays the mo- tion of an idealized pendulum, without friction. The motion of a pendulum with friction is portrayed in Fig. 3.4. In Fig. 3.4, we see a heteroclinic orbit from ( ¡¼;0) to (0 ;0), another from ( ¡¼;0) to ( ¡2¼;0), another from ( ¼;0) to (0 ;0), another from ( ¼;0) to (2 ¼;0), etc. Given that there is an orbit x(t) = ( µ(t); Ã(t)) here such that lim t!¡1 x(t) = (¡¼;0) and Ã(t)>0 for large negative t, the fact that limt!+1x(t) = (0 ;0) can be deduced from (3.54), i.e., (3.71)d dtE(µ; Ã) =¡® mÃ2: We end this section with a look at the phase portrait for one more vector ¯eld, namely (3.72) F(µ; Ã) =µg `sinµ ö : See Fig. 3.7. In this case, (3.73) F=rE=µ@E=@µ @E=@ö ; withEgiven by (3.25), i.e., (3.74) E(µ; Ã) =Ã2 2¡g `cosµ: 262 4. Nonlinear Systems of Di®erential Equations Figure 3.7 Such a vector ¯eld is called a gradient vector ¯eld , and its °ow ©t Fis called agradient °ow . Note that if F(x) =rE(x) and x(t) is an orbit of ©t F, then (3.75)d dtE(x(t)) =rE(x(t))¢ rE(x(t)) =krE(x(t))k2: The critical points of Fagain consist of f(k¼;0) :k2Zg, and again they behave di®erently for even kthan for odd k. This time (3.76) DF(0;0) =µg `0 0 1¶ ; which is positive de¯nite. The origin is a (non-spiraling) source ; cf. Fig. 3.6. In particular, if x(t) = (µ(t); Ã(t)) is an orbit and x(0) is close to (0 ;0), then (3.77) lim t!¡1x(t) = (0 ;0): This can be deduced from Proposition 3.4 by reversing time. It also follows directly from (3.75). For kodd, we have saddles: (3.78) DF(§¼;0) =µ¡g `0 0 1¶ : In this case, segments of the real axis provide heteroclinic orbits, from (0 ;0) to (¡¼;0), from (0 ;0) to ( ¼;0), etc. 3. Vector ¯elds, orbits, and °ows 263 Exercises 1. If Fgenerates the °ow ©t Fandvt(x) =v(©t F(x)), show that (3.79) D©t F(x)F(x) =F(©t F(x)) and (3.80) Dvt(x) =Dv(©t F(x))D©t F(x): Relate these identities to the simultaneous validity of (3.10) and (3.12). Hint. To get (3.79), use (3.81)d dt©t F(x) =d ds©t F±©s F(x)¯¯¯ s=0=D©t F(x)F(x); and compare (3.8). 2. Extend Proposition 3.2 as follows. Replace hypothesis (3.19) by rV(x)¢F(x)·K;8x2Rn; for some K <1. Show that the °ow ©t Fis forward complete. 3. Let ­ =Rnand assume Fis aC1vector ¯eld on ­. Show that if kF(x)k ·C(1 +kxk); then the °ow generated by Fis complete. (Hint. Recall Exercise 12 of x1.) Show that the °ow is forward complete if F(x)¢x·C(1 +kxk2): 4. Let ­ ½Rnbe open and Fbe a C1vector ¯eld on ­. Let U½­ be an open set whose closure Uis a compact subset of ­, and whose boundary @Uis smooth. Let n:@U!Rndenote the outward pointing unit normal to @U. Assume (3.82) n(x)¢F(x)<0;8x2@U: 264 4. Nonlinear Systems of Di®erential Equations Show that ©t F(x)2Uifx2Uandt¸0, and deduce that ©t Fis forward complete on U, and also on U. 5. In the setting of Exercise 4, relax the hypothesis (3.82) to (3.83) n(x)¢F(x)·0;8x2@U: Show that ©t F(x)2Uifx2Uandt¸0, and deduce that ©t Fis forward complete on U. Hint. Find a smooth family F¿ofC1vector ¯elds on ­ such that F0=F and, for 0 < ¿ < 1,F¿has the property given in (3.82). Then make use of Exercise 4 and of results of x2. 6. In the setting of Exercise 5, replace the hypothesis (3.83) by (3.84) n(x)¢F(x) = 0 ;8x2@U: Show that ©t Fis complete on U, and that ©t F(x)2@Uwhenever x2@U andt2R. 7. Show that if Fis given by (3.24), then divF= 0; while if Fis given by (3.50), then divF=¡® m; and if Fis given by (3.72), then divF=g `cosµ+ 1: 8. Let A2M(2;R), and take Jas in (3.29). Show that if Ais positive de¯nite then A=P2with Ppositive de¯nite. Show that ¡JAand¡PJP are similar ; and deduce that A2M(2;R) positive de¯nite = )B=¡JAhas 2 purely imaginary eigenvalues : Relate this to Proposition 3.3. 3. Vector ¯elds, orbits, and °ows 265 9. Give an alternative proof of Proposition 3.4, avoiding use of Lemma 3.5, starting with the representation of y(t) = ©t F(x0+y0)¡x0as y(t) =etLy0+Zt 0e(t¡s)LR(y(s))ds; and using the fact that (3.55) implies ketLk ·Ce¡Kt; for some C; K2(0;1). 10. Consider the system (3.85)dx dt=y;dy dt= 1¡x2: Take E(x; y) =y2=2 +x3=3¡x. Show that if ( x(t); y(t)) solves (3.85), then dE(x(t); y(t)) = 0. Show that the associated vector ¯eld has two critical points, one a center and the other a saddle. Sketch level curves ofEand put in arrows to show the phase space portrait of F. Show that there is a homoclinic orbit connecting the saddle to itself. 11. Returning to the context of Exercise 1, show that (2.2) gives (3.86)d dtD©t F(x) =DF(©t F(x))D©t F(x); D ©0 F(x) =I: Recall from (8.6){(8.10) of Chapter 3 that, for an n£nmatrix function M(t), d dtM(t) =A(t)M(t) =)d dtdetM(t) = (Tr A(t)) det M(t): Deduce that (3.87)d dtdetD©t F(x) = Tr DF(©t F(x)) det D©t F(x) = div F(©t F(x)) det D©t F(x)): Relate this to (3.13), using the change of variable formula (3.88)Z u(x)dx=Z u(©t F(x)) det D©t F(x)dx: 266 4. Nonlinear Systems of Di®erential Equations 12. Use (8.10) of Chapter 3 to conclude from (3.87) that (3.89) det D©t F(x) = expnZt 0divF(©s F(x))dso : 13. Let U½­½Rnbe a smoothly bounded domain. The divergence theorem says that if Fis aC1vector ¯eld on ­, (3.90)Z UdivF(x)dx=Z @Un(x)¢F(x)dS(x); where n(x) is the outward pointing unit normal to @UanddS(x) is (n¡1)-dimensional surface area on @U(arc length if n= 2). Given this identity, we see that, in the setting of Proposition 3.1, (3.17) is equivalent to (3.91)d dtVol ©t F(B) =Z @©t F(B)n(x)¢F(x)dS(x): Show that this holds if and only if for each smoothly bounded U½­, (3.92)d dtVol ©t F(U)¯¯ t=0=Z @Un(x)¢F(x)dS(x): Try to provide a direct demonstration of (3.92) (at least for n= 2). 4. Gradient vector ¯elds As mentioned in x3, a vector ¯eld Fon an open subset ­ ½Rnis a gradient vector ¯eld provided there exists u2C1(­) such that (4.1) F=ru; i.e.,F= (F1; : : : ; F n)twith Fk=@u=@x k. It is of interest to characterize which vector ¯elds are gradient ¯elds. Here is one necessary condition. Suppose u2C2(­) and (4.1) holds. Then (4.2)@Fk @xj=@ @xj@u @xk; and (4.3)@ @xj@u @xk=@ @xk@u @xj; 4. Gradient vector ¯elds 267 so if (4.1) holds then (4.4)@Fk @xj=@Fj @xk;8j; k2 f1; : : : ; n g: We will establish the following converse. Proposition 4.1. Assume ­½Rnis a connected open set satisfying the condition (4.13) given below. Let Fbe aC1vector ¯eld on ­. If (4.4) holds on­, then there exists u2C2(­)such that (4.1) holds. We will construct uas a line integral. Namely, ¯x p2­, and for each x2­ let °be a smooth path from ptox: (4.5) °: [0;1]¡!­; °(0) = p; °(1) = x: We propose that, under the hypotheses of Proposition 4.1, we can take (4.6) u(x) =Z °F(y)¢dy: Here the line integral is de¯ned by (4.7)Z °F(y)¢dy=Zt 0F(°(t))¢°0(t)dt: For this to work, we need to know that (4.6) is independent of the choice of such a path. A key step to getting this is to consider a smooth 1-parameter family of paths °sfrom ptox: (4.8)°s(t) =°(s; t); ° : [0;1]£[0;1]¡!­; °(s;0) = p; ° (s;1) = x: Lemma 4.2. IfFis aC1vector ¯eld satisfying (4.4) and °sis a smooth family satisfying (4.8), then (4.9)Z °sF(y)¢dyis independent of s2[0;1]: 268 4. Nonlinear Systems of Di®erential Equations Proof. We compute the s-derivative of this family of line integrals, i.e., of (4.10)Z1 0F(°(s; t))¢@° @t(s; t)dt =Z1 0X jFj(°(s; t))@°j @t(s; t)dt: Thes-derivative of the integrand is obtained via the product rule and the chain rule. We obtain (4.11)d dsZ °sF(y)¢dy=Z1 0X j;k@Fj @xk(°(s; t))@ @s°k(s; t)@ @t°j(s; t)dt +Z1 0X jFj(°(s; t))@ @s@ @t°j(s; t)dt: We can apply the identity @ @s@ @t°j(s; t) =@ @t@ @s°j(s; t) to the second integrand on the right side of (4.11) and then integrate by parts. This involves applying @=@t toFj(°(s; t)), and hence another appli- cation of the chain rule. When this is done, the second integral on the right side of (4.11) becomes (4.12) ¡Z1 0X j;k@Fj @xk(°(s; t))@ @t°k(s; t)@ @s°j(s; t)dt: Now if we interchange the roles of jandkin (4.12), we cancel the ¯rst integral on the right side of (4.11), provided (4.4) holds. This proves the lemma. Given ­ ½Rnopen and connected, we say ­ is simply connected pro- vided it has the following property: (4.13)Given p; x2­;if°0and°1are smooth paths from ptox; they are connected by a smooth family °sof paths from ptox: Here is a class of such domains. Lemma 4.3. If­½Rnis an open convex domain, then ­is simply con- nected. 4. Gradient vector ¯elds 269 Proof. If ­ is convex, two paths °0and°1from p2­ to x2­ are connected by (4.14) °s(t) = (1 ¡s)°0(t) +s°1(t);0·s·1: Of course there are many other simply connected domains, as the reader is invited to explore. Now that we have Lemma 4.2, under the hypotheses of Proposition 4.1 we simply write (4.15) u(x) =Zx pF(y)¢dy: Note that if qis another point in ­, we can take a smooth path from pto x, passing through q, and write (4.16) u(x) =Zq pF(y)¢dy+Zx qF(y)¢dy: Again using the path independence, we see we can independently choose paths from ptoqand from qtoxin (4.16); these paths need not match up smoothly at q. We are now in a position to complete the proof of Proposition 4.1. Take ± >0 so that fy2Rn:kx¡yk ·±g ½­. Take k2 f1; : : : ; n g, ¯xqksuch thatjqk¡xkj< ±, and write (4.17) u(x) =Z(x1;:::;qk;:::;x n) pF(y)¢dy+Zx (x1;:::;qk;:::;x n)F(y)¢dy: Here the intermediate point is obtained by replacing xkinx= (x1; : : : ; x n) byqk. The ¯rst term on the right side of (4.17) is independent of xk, so (4.18)@u @xk(x) =@ @xkZx (x1;:::;qk;:::;x n)F(y)¢dy =@ @xkZxk qkFk(x1; : : : ; x k¡1; s; x k+1; : : : ; x n)ds =Fk(x); the last identity by the fundamental theorem of calculus. This proves Propo- sition 4.1. 270 4. Nonlinear Systems of Di®erential Equations An example of a domain that is not simply connected is the punctured planeR2n0. Consider on this domain the vector ¯eld (4.19) F(x) =Jx kxk2; J =µ0¡1 1 0¶ ; with components (4.20) F1(x) =¡x2 x2 1+x2 2; F 2(x) =x1 x2 1+x2 2: We have (4.21)@F1 @x2=x2 2¡x2 1 kxk4=@F2 @x1;onR2n0: However, Fis not a gradient vector ¯eld on R2n0. Up to an additive constant, the only candidate for uin (4.1) is the angular coordinate µ: (4.22) F(x) =rµ(x); and this identity is true on any region ­ formed by removing from R2a ray starting from the origin. However, µcannot be de¯ned as a smooth, single valued function on R2n0. Let us linger on the case n= 2 and make contact with the concept of \exact equations." Consider a 2 £2 system (4.23)dx dt=f1(x; y);dy dt=f2(x; y): We take ( x; y)2­½R2and assume fj2C1(­). This system turns into a single di®erential equation for yas a function of x: (4.24)dy dx=f2(x; y) f1(x; y); which we rewrite as (4.25)g1(x; y)dx+g2(x; y)dy= 0; g1(x; y) =f2(x; y); g 2(x; y) =¡f1(x; y): The equation (4.25) is called exact if there exists u2C2(­) such that (4.26) g1=@u @x; g 2=@u @y: 4. Gradient vector ¯elds 271 If there is such a u, solutions to (4.24) or (4.25) are given by (4.27) u(x; y) =C: Now (4.26) is the condition that G= (g1; g2)tbe a gradient vector ¯eld on ­. Note that the relation between F= (f1; f2)tandG= (g1; g2)t, with components given by (4.25), is (4.28) G=¡JF; where (4.29) J=µ0¡1 1 0¶ : As we have seen, when ­ is simply connected, (4.26) holds for some uif and only if (4.30)@g1 @y=@g2 @x: Note that this is equivalent to (4.31) div F= 0: Remark. IfF= (F1; F2; F3)tis a vector ¯eld on ­ ½R3, its curl is de¯ned as (4.32)curlF=r £F = det0 @i j k @=@x @=@y @=@z F1 F2 F31 A =³@F3 @y¡@F2 @z´ i+³@F1 @y¡@F2 @x´ j+³@F2 @x¡@F1 @y´ k: We see that (4.33) (4.4) holds ()curlF= 0: We conclude with some remarks on how to construct u(x), satisfying (4.34)@u @xj(x) =Fj(x);1·j·n; 272 4. Nonlinear Systems of Di®erential Equations given the compatibility conditions (4.4), without evaluating line integrals. We start with (4.35) un(x) =Z Fn(x)dxn;so@un @xn=Fn(x): Then @(u¡un)=@xn= 0, so (4.36) u(x) =un(x) +v(x0); x0= (x1; : : : ; x n¡1): It remains to ¯nd v, a function of fewer variables. It must solve (4.37)@v @xj=Fj(x)¡@un @xj;1·j·n¡1: Note that the left side is independent of xn, which requires that the right side have this property. To check this, we calculate (4.38)@ @xn³ Fj(x)¡@un @xj´ =@Fj @xn¡@ @xn@un @xj =@Fn @xj¡@ @xj@un @xn = 0; the second identity by (4.4) (and (4.3)). Thus (4.37) takes the form (4.39)@v @xj=Gj(x0);1·j·n¡1; with Gj(x0) =Fj(x)¡@un=@xj. Note that, for 1 ·j; k·n¡1, (4.40)@Gj @xk=@Fj @xk¡@ @xk@un @xj =@Fk @xj¡@ @xj@un @xk =@Gk @xj; so the task of solving (4.39) is just like that in (4.34), but with one fewer variable. An iteration yields the solution to (4.34). Example. Take (4.41) F(x; y; z ) = (y; x+z2;2yz)t: 4. Gradient vector ¯elds 273 One readily ver¯es (4.4), or equivalently that curl F= 0. Here (4.35) gives (4.42) u3(x; y; z ) =Z 2yz dz =yz2; so (4.43) u(x; y; z ) =yz2+v(x; y): Next, requiring @u=@y =x+z2means (4.44)@v @y=x; so (4.45) v(x; y) =xy+w(x): Then, requiring @u=@x =ymeans @w=@x = 0, so we get (4.46) u(x; y; z ) =yz2+xy; as the unique function on R3such that ru=F, up to an additive constant. One can turn the method given by (4.35){(4.40) into an alternative proof of Proposition 4.1, at least if ­ is an n-dimensional box. The reader is invited to look into what happens when this method is applied to Fgiven onR2n0 by (4.19). Exercises For (1){(4), identify which vector ¯elds are gradient ¯elds. If the ¯eld is a gradient ¯eld ru, ¯nd u. (yz; xz; xy ); (1) (xy; yz; xz ); (2) (2x; z; y ); (3) (2x; y; z ): (4) 274 4. Nonlinear Systems of Di®erential Equations For (5){(8), identify which equations are exact. If the equation is exact, write down the solution, in implicit form (4.27). (2x+y)dx+x dy= 0; (5) x dx+ (2x+y)dy= 0; (6) dx+x dy= 0; (7) eydx+xeydy= 0: (8) Given f(x; y)dx+g(x; y)dy, a function u(x; y) is called an integrating factor ifuf dx +ug dy is exact. For example, eyis an integrating factor fordx+x dy. Find integrating factors for the left sides of (9){(12), and use them to ¯nd solutions, in implicit form. (x2+y2¡1)dx¡2xy dy = 0; (9) x2y3dx+x(1 +y2)dy= 0; (10) y dx+ (2x¡yey)dy= 0; (11) dx+ 2xy dy = 0: (12) 13. Establish the following variant of Lemma 4.2: Lemma 4.2A. IfFis aC1vector ¯eld on ­satisfying (4.4) and °sis a smooth family satisfying °s(t) =°(s; t); ° : [0;1]£[0;1]!­; °(s;0)´°(s;1); then Z °sF(y)¢dyis independent of s2[0;1]: 5. Newtonian equations In Chapter 1 we saw how Newton's law F=maleads to a second order di®erential equation for the motion on a line of a single particle, acted on by a force. Newton's laws also apply to a system of minteracting particles, moving in n-dimensional space, to give a second order system of the form (5.1) mkd2xk dt2=X fj:j6=kgFjk(xk¡xj);1·k·m: 5. Newtonian equations 275 Each xktakes values in Rn, sox= (x1; : : : ; x m) takes values in Rmn. Here xk is the location of a particle of mass mk. The law that each action produces an equal and opposite reaction translates to (5.2) Fjk(xk¡xj) =¡Fkj(xj¡xk): A particularly important class of forces Fjk(xk¡xj) are those parallel (or antiparallel) to the line from xjtoxk: (5.3) Fjk(xk¡xj) =fjk(kxk¡xjk)(xk¡xj): In such a case, (5.2) is equivalent to (5.4) fjk(r) =fkj(r): A force ¯eld of the form (5.3) is a gradient vector ¯eld: (5.5)fjk(kuk)u=¡rVjk(u); Vjk(u) =vjk(kuk); v0 jk(r) =¡rfjk(r): If (5.4) holds, (5.6) Vjk(u) =Vkj(u): The total energy of this system of interacting particles is (5.7) E=1 2X kmk°°°dxk dt°°°2 +1 2X j6=kVjk(xk¡xj): The ¯rst sum is the total kinetic energy and the second sum is the total potential energy. The following calculations yield conservation of energy. First, (5.8)dE dt=X kmkd2xk dt2¢dxk dt+1 2X j6=krVjk(xk¡xj)¢³dxk dt¡dxj dt´ : Next, (5.1) implies that the ¯rst sum on the right side of (5.8) is equal to (5.9)X j6=kFjk(xk¡xj)¢dxk dt; and (5.3){(5.5) imply that the second sum on the right side of (5.8) is equal to (5.10) ¡1 2X j6=kFjk(xk¡xj)¢³dxk dt¡dxj dt´ =¡X j6=kFjk(xk¡xj)¢dxk dt: 276 4. Nonlinear Systems of Di®erential Equations Comparing (5.9) and (5.10), we have energy conservation: (5.11)dE dt= 0: We can convert the second order system (5.1) for mnvariables into a ¯rst order system for 2 mnvariables. One way would be to introduce the velocities vk=x0 k, but we get a better mathematical structure by instead using the momenta : (5.12) pk=mkdxk dt;1·k·m: We can express the energy Ein (5.7) as a function of position x= (x1; : : : ; x m) and momentum p= (p1; : : : ; p m): (5.13) E(x; p) =X k1 2mkkpkk2+1 2X j6=kVjk(xk¡xj): Recall that xk= (xk1; : : : ; x kn)2Rnandpk= (pk1; : : : ; p kn)2Rn. We have (5.14)@E @pk`=1 mkpk`; and (5.15)@E @xk`=X fj:j6=kg@Vjk @u`(xk¡xj); invoking (5.6). Let us write (5.14){(5.15) in vector form, (5.16)@E @pk=1 mkpk;@E @xk=X fj:j6=kgrVjk(xk¡xj); where @E=@p k= (@E=@p k1; : : : ; @E=@p kn)t, etc. Now the system (5.1) yields the ¯rst order system (5.17)dxk dt=1 mkpk;dpk dt=X fj:j6=kgFjk(xk¡xj); which in turn, given (5.3){(5.5), gives (5.18)dxk dt=@E @pk;dpk dt=¡@E @xk: 5. Newtonian equations 277 The system (5.18) is said to be in Hamiltonian form . We can place the study of Hamiltonian equations in a more general framework, as follows. Let R2Khave points ( x; p); x= (x1; : : : ; x K); p= (p1; : : : ; p K). Let ­ ½R2Kbe open and E2C1(­). Consider the system (5.19)dxk dt=@E @pk; dpk dt=¡@E @xk; for 1·k·K. This is called a Hamiltonian system. It is of the form (5.20)d dtµx p¶ =XE(x; p); where XEis a vector ¯eld on ­, called the Hamiltonian vector ¯eld asso- ciated to E. In this general setting, Eis constant on each solution curve (x(t); p(t)) of (5.19). Indeed, in such a case, (5.21)d dtE(x(t); p(t)) =X k@E @xk¢dxk dt+X k@E @pk¢dpk dt =X k@E @xk¢@E @pk¡X k@E @pk¢@E @xk = 0: Returning to the setting (5.1){(5.2), we next discuss the conservation of the total momentum (5.22) P=X kpk=X kmkdxk dt: Indeed, (5.23)dP dt=X kmkd2xk dt2 =X j6=kFjk(xk¡xj) = 0; the last identity by (5.2). Thus, for each solution x(t) to (5.1), there exist a; b2Rnsuch that (5.24)1 MX kmkxk(t) =a+bt; M =X kmk: 278 4. Nonlinear Systems of Di®erential Equations The left side is the center of mass of the system of interacting particles. The vectors a; b2Rnare given by the initial data for (5.1): (5.25) a=1 MX kmkxk(0); b =1 MX kmkx0 k(0): Given this, we can obtain a system similar to (5.1) for the variables (5.26) yk(t) =xk(t)¡(a+bt): We have y00 k=x00 kandyk¡yj=xk¡xj, so (5.1) gives (5.27) mkd2yk dt2=X fj:j6=kgFjk(yk¡yj);1·k·m: In this case we have the identity (5.28)X kmkyk(t)´0; as a consequence of (5.24). We can use this to reduce the size of (5.27), from a system of mnequations to a system of ( m¡1)nequations, by substituting (5.29) ym=¡1 mmm¡1X `=1m`y` into (5.27), for 1 ·k·m¡1. One calls ( y1; : : : ; y m) center of mass coordinates. In case m= 2, this substitution works out quite nicely. We have (5.30) y2=¡m1 m2y1; and the system (5.30) reduces to (5.31) m1d2y1 dt2=F21³³ 1 +m1 m2´ y1´ ; the equation of motion of a single particle in an external force ¯eld. Form > 2, the resulting equations are not so neat. For example, for m= 3, we have (5.32) y3=¡m1 m3y1¡m2 m3y2; and the system (5.27) reduces to (5.33)m1y00 1=F21(y1¡y2) +F31³³ 1 +m1 m2´ y1+m2 m1y2´ ; m2y00 2=F12(y2¡y1) +F32³m1 m3y1+³ 1 +m2 m3´ y2´ : 6. Central force problems and two-body planetary motion 279 Exercises In (1){(5), take n= 1; m = 3, and m1=m2=m3= 1. Set up the equations of motion in center of mass coordinates and analyze the solution. Fjk(x) =x (1) Fjk(x) =¡x (2) F12(x) =F13(x) =x; F 23(x) =¡x (3) F12(x) =F13(x) =¡x; F 23(x) =x: (4) F12(x) =F23(x) =¡x; F 23(x) = 1 : (5) In all cases, (5.2) must be enforced. 6. Central force problems and two-body planetary motion As seen in x5, one can transform the m-body problem (5.1) to center of mass coordinates, under the hypothesis (5.2), and obtain a smaller system, which form= 2 is given by (5.31). Changing notation, we rewrite (5.31) as (6.1) md2x dt2=F(x): Here x2Rn. We assume F2C1(Rnn0) but allow blowup at x= 0. Under hypotheses (5.3){(5.4) for the two body problem, we have (6.2) F(x) =f(kxk)x: In such a case, (6.1) is called a central force problem . Parallel to (5.5), we have (6.3)F(x) =¡rV(x); V(x) =v(kxk); v0(r) =rf(r): The total energy is given by (6.4) E=1 2m°°°dx dt°°°2 +V(x); and if x(t) solves (6.1), then (6.5)dE dt=md2x dt¢dx dt+rV(x)¢dx dt= 0; 280 4. Nonlinear Systems of Di®erential Equations yielding conservation of energy. There are further conservation laws, starting with the following. Proposition 6.1. Assume x(0)6= 0and let W½Rnbe the linear span of x(0)andx0(0). Ifx(t)solves (6.1) for t2Iand (6.2) holds, we have (6.6) x(t)2W;8t2I: Proof. De¯ne A2 L(Rn) by (6.7)Av=v; 8v2W; Av=¡v;8v2W?: Note that Ais an orthogonal transformation. Let y(t) =Ax(t). The hy- pothesis on the initial data gives (6.8) y(0) = x(0); y0(0) = x0(0): Also, given F(x) of the form (6.2), we have AF(x) =F(y), so y(t) solves (6.1). The basic uniqueness result proven in x1 implies y´xonI, which in turn gives (6.6). Thus each path x(t) solving (6.1) lies in a plane, and we can take n= 2. For the next step, it is actually convenient to take n= 3. Thus x(t) solves (6.1) and x(t) is a path in R3. We de¯ne the angular momentum (6.9) ®(t) =mx(t)£x0(t): We then have, under hypothesis (6.2), (6.10)®0(t) =mx(t)£x00(t) =x(t)£F(x) =f(kxk)x(t)£x(t) = 0: This yields conservation of angular momentum: (6.11) x(t)£x0(t)´L; where L=x(0)£x0(0)2R3. In case x(t) = (x1(t); x2(t);0), we have (6.12) x(t)£x0(t) = (0 ;0; x1(t)x0 2(t)¡x0 1(t)x2(t)); 6. Central force problems and two-body planetary motion 281 so the conservation law (6.11) gives (6.13) x1(t)x0 2(t)¡x0 1(t)x2(t)´L3: Let's return to the planar setting, and also use complex notation: (6.14) x(t) =x1(t) +ix2(t) =r(t)eiµ(t): A computation gives (6.15)x0= (r0+irµ0)eiµ; x00= [r00¡r(µ0)2+i(2r0µ0+rµ00)]eiµ; so (6.1){(6.2) becomes (6.16) m£ r00¡r(µ0)2+i(2r0µ0+rµ00)¤ =f(r)r: Equating real and imaginary parts separately, we get (6.17)r00¡r(µ0)2=f(r)r m; 2r0µ0+rµ00= 0: Note that (6.18)d dt(r2µ0) =r(2r0µ0+rµ00); so the second equation in (6.17) says r2µ0is independent of t. This is actually equivalent to the conservation of angular momentum, (6.13). In fact, we have x1=rcosµ; x 2=rsinµ, hence (6.19) x0 1=r0cosµ¡rµ0sinµ; x0 2=r0sinµ+rµ0cosµ; and hence (6.20) x1x0 2¡x0 1x2=r2µ0: Thus we have in two ways derived the identity (6.21) r2µ0=L: (For notational simplicity, we drop the subscript 3 from (6.13).) There is the following geometrical interpretation of (6.21). The (signed) areaA(t) swept out by the ray from 0 to x(s), assruns from t0tot, is given by (6.22) A(t) =1 2Zµ(t) µ(t0)r2dµ=1 2Zt t0r(s)2µ0(s)ds; 282 4. Nonlinear Systems of Di®erential Equations so (6.23) A0(t) =1 2r2µ0=L 2: This says (6.24) Equal areas are swept out in equal times, which, as we will discuss below, is Kepler's second law. Next, we can plug µ0=L=r2into the ¯rst equation of (6.17), obtaining (6.25)d2r dt2=f(r)r m+L2 r3: This has the form (6.26)d2r dt2=g(r); treated in Chapter 1, x5. We recall that treatment. Take w(r) such that g(r) =¡w0(r), so (6.26) becomes (6.27)d2r dt2=¡w0(r): Then form the \energy" (6.28) E=1 2³dr dt´2 +w(r); and compute that if r(t) solves (6.27) then (6.29)dE dt=d2r dt2dr dt+w0(r)dr dt= 0; so for each solution to (6.27), there is a constant Esuch that (6.30)dr dt=§p 2E¡2w(r): Separation of variables gives (6.31)Zdrp 2E¡2w(r)=§t+C: This integral can be quite messy. 6. Central force problems and two-body planetary motion 283 Note that dividing (6.30) by (6.21) yields a di®erential equation for ras a function of µ: (6.32)dr dµ=§r2 Lp 2E¡2w(r); which separates to (6.33) LZdr r2p 2E¡2w(r)=§µ+C: Let us recall that (6.34) w0(r) =¡f(r)r m¡L2 r3: Typically the integral in (6.33) is as messy as the one in (6.31). These integrals do turn out to be tractable in one very important case, the Kepler problem, to which we now turn. This problem is named after the astronomer Johannes Kepler, who from observations formulated the following three laws for planetary motion. 1. The planets move on ellipses with the sun at one focus. 2. The line segment from the sun to a planet sweeps out equal areas in equal time intervals. 3. The period of revolution of a planet is proportional to a3=2, where ais the semi-major axis of its ellipse. The Kepler problem is to provide a theoretical framework in which to derive these three laws. This was solved by Isaac Newton, who formulated his universal law of gravitation, used it to derive a di®erential equation for the position of a planet, and solved the di®erential equation. Newton's law of gravitation speci¯es the force between two objects, of mass m1andm2, located at points x1andx2inR3. Let us say the center of the planet is at x1and the center of the sun is at x2. In the framework of (5.1), this means specifying the vector ¯eld F21onR3. The formula is (6.35) F21(x) =¡Gm1m2x kxk3: Here Gis the universal gravitational constant. If we go to center of mass coordinates, the motion of the planet is governed by (5.31), and changing notation, replacing y1byx, yields (6.1) with (6.36) F(x) =¡Kmx kxk3; K =G(m+m2): 284 4. Nonlinear Systems of Di®erential Equations Here m=m1is the mass of the planet and m2is the mass of the sun. Consequently we have (6.17) with (6.37)f(r)r m=¡K r2; and (6.25) becomes (6.38)d2r dt2=¡K r2+L2 r3: Thus w(r) in (6.27){(6.34) is given by (6.39) w(r) =¡K r+L2 2r2: Thus the integral in (6.31) is (6.40)Zr drp 2Er2+ 2Kr¡L2; and the integral in (6.33) is (6.41)Zdr rp 2Er2+ 2Kr¡L2: The integral (6.40) can be evaluated by completing the square for 2 Er2+ 2Kr¡L2. The integral (6.41) can also be evaluated, but rather than tackling this directly, we instead produce a di®erential equation for u, de¯ned by (6.42) u=1 r: By the chain rule, (6.43)dr dt=¡r2du dt=¡r2du dµdµ dt=¡Ldu dµ; the last identity by (6.21). Taking another t-derivative gives (6.44)d2r dt2=¡Ld dtdu dµ=¡Ld2u dµ2dµ dt=¡L2u2d2u dµ2; again using (6.21). Comparing this with (6.38), we get (6.45) ¡L2u2d2u dµ2=L2u3¡Ku2; 6. Central force problems and two-body planetary motion 285 or equivalently (6.46)d2u dµ2+u=K L2: Miraculously, we have obtained a linear equation! The general solution to (6.46) is (6.47) u(µ) =Acos(µ¡µ0) +K L2; which by (6.42) gives (6.48) rh Acos(µ¡µ0) +K L2i = 1: This is equivalent to (6.49) r£ 1 +ecos(µ¡µ0)¤ =p; p =L2 K; e=AL2 K: Ife= 0, this is the equation of a circle. If 0 < e < 1, it is the equation of an ellipse. If e= 1, it is the equation of a parabola, and if e >1, it is the equation of one branch of a hyperbola. Among these curves, those that are bounded are the ellipses, and the circle, which we regard as a special case of an ellipse. Since planets move in bounded orbits, this establishes Kepler's ¯rst law (with caveats, which we discuss below). Kepler's second law holds for general central force problems, as noted already in (6.24). To establish the third law, recall from (6.23) that L=2 is the rate at which such area is swept out, so the period Tof the orbit satis¯es (6.50)L 2T= area enclosed by the ellipse =¼ab; where ais the semi-major axis and bthe semi-minor axis. For an ellipse given by (6.49), we have (6.51) a=p 1¡e2; b =pp 1¡e2=p1=2a1=2; which yields (6.52) T=2¼ab L= 2¼pp La3=2=2¼p Ka3=2: 286 4. Nonlinear Systems of Di®erential Equations This establishes Kepler's third law. We now discuss some caveats. Our solar system has nine planets, plus numerous other satellites. In the calculations above, all but one planet was ignored. One can expect this approximation to work best for Jupiter. Jupiter has about 10¡3the sun's mass, and its distance from the sun is about 400 times the sun's radius. Hence the center of mass of Jupiter and the sun is located about 0.4 times the sun's radius from the center of the sun. The sun and Jupiter engage in a close to circular elliptical orbit with a focus at this center of mass. Clearly this motion is going to in°uence the orbits of the other planets. In fact, each planet in°uences all the others, including Jupiter, in ways not captured by the calculations of this section. Realization of this situation led to a vigorous development of the subject known as celestial mechanics, from Newton's time on. Material on this can be found in [ AM] and [ Gr], and references given there. Advances in celestial mechanics led to the discovery of the planet Nep- tune. By the early 1900s, this subject was su±ciently well developed that astronomers were certain that an observed anomaly in the motion of Mer- cury could not be explained by the Newtonian theory. This discrepancy was accounted for by Einstein's theory of general relativity, which provided a new foundation for the theory of gravity. This is discussed in [ ABS ] and also in Chapter 18 of [ T]. While a derivation is well outside the scope of this book, we mention that the relativistic treatment leads to the following variant of (6.46): (6.53)d2u dµ2+u=A+"u2; where A¼K=L2and"is a certain (small) positive constant, determined by the mass of the sun. This can be converted into the ¯rst order system (6.54)du dµ=v;dv dµ=¡u+A+"u2: In analogy with (6.26){(6.29), we can form (6.55) F(u; v) =1 2v2+1 2u2¡Au¡" 3u3; and check that if ( u(µ); v(µ)) solves (6.54), then (6.56)d dµF(u; v) = 0 ; so the orbits for (6.54) lie on level curves of F. As long as A"2(0;1=4), Fhas two critical points, a minimum and a saddle. Thus (6.54) has some 6. Central force problems and two-body planetary motion 287 solutions periodic in µ. However, the period is generally not equal to 2 ¼. (See Appendix D for results related to computing this period.) This fact leads to the precession of the perihelion of the planet orbiting the sun, where the perihelion is the place where uis maximal, so ris minimal. In the non- relativistic situation covered by (6.46), all the solutions in (6.47) are periodic inµof period 2 ¼. Exercises 1. Solve explicitly w00(t) =¡w(t); forwtaking values in R2=C. Show that jw(t)j2+jw0(t)j2= 2E is constant on each orbit. 2. For w(t) taking values in C, de¯ne a new curve by z(s) =w(t)2;ds dt=jw(t)j2: Show that if w00(t) =¡w(t), then z00(s) =¡4Ez(s) jz(s)j3; soz(s) solves the Kepler problem. 3. Take u= 1=ras in (6.42), and generalize the calculations (6.43){(6.46) to obtain a di®erential equation for uas a function of µ, for more general central forces. Consider particularly f(x) =¡rV(x) in the cases V(x) =¡Kkxk2; V (x) =¡Kkxk: 4. Take the following steps to show that if p >0 and 0 < e < 1, then (6.57) r(1 +ecosµ) =p 288 4. Nonlinear Systems of Di®erential Equations is the equation in polar coordinates of an ellipse. (a) Show that (6.57) describes a closed, bounded curve, since 1+ ecosµ > 0 for all µif 0< µ < 1, and cos µis periodic in µof period 2 ¼. Denote the curve by °(µ) = (x(µ); y(µ)), in Cartesian coordinates. (b) Show that this curve is symmetric about the x-axis and cuts the axis at two points, whose distance apart is 2a=r(0) + r(¼); so (6.58) a=p 1¡e2: (c) Show that the midpoint between °(0) and °(¼) is given by x0=¡ea; y 0= 0: (d) For °(µ) = (x(µ); y(µ)), as in part (a), show that (6.59)(x+ea)2 a2+y2 b2= 1; i.e., that (rcosµ+ea)2 a2+r2(1¡cos2µ) b2= 1; provided (6.57) holds, when ais given by (6.58) and (6.60) b=pp 1¡e2: 5. As an approximation, assume that the earth has a circular orbit about the sun with a radius (6.61) a= 1:496£1011m; and its period is one year, i.e., (6.62) T= 31:536£106sec: 6. Central force problems and two-body planetary motion 289 The gravitational constant Ghas been measured as (6.63) G= 6:674£10¡11m3=(kg sec2): With this information, use (6.36) and (6.52) to calculate the mass m2 of the sun. Assume the mass of the earth is negligible compared to m2. You should get (6.64) m2=®£1030kg; with ®between 1 and 10. Remark. Historically, Twas measured by the position of the \¯xed stars." Modern methods to measure ainvolve bouncing a radar signal o® Venus to measure its distance, given that we have an accurate mea- surement of the speed of light. Then trigonometry is used to determine a. See [ GM] for a discussion of how Ghas been measured; this is the most di±cult issue. 6. The force of gravity the earth exerts on a body of mass mat the earth's surface is (6.65) ¡Gmm er¡2; where Gis given in Exercise 5, (6.66) r= 6:38£106m is the radius of the earth, and meis the mass of the earth. It is observed that the earth's gravity accelerates objects at its surface downward at 9:8 m/sec2, so we have (6.67) 9 :8 m/sec2=Gmer¡2: Use this to compute me. You should get (6.68) me=¯£1024kg; with ¯between 1 and 10. Remark. See Appendix E for more on (6.65). 7. As an approximation, assume that the moon has a circular orbit about the earth, of radius a= 3:8£108m; 290 4. Nonlinear Systems of Di®erential Equations and its period is 27.3 days, i.e., T= 2:359£106sec: Assume the mass of the moon is negligible compared to the mass of the earth. Use the method of Exercise 5 to calculate the mass of the earth. Compare your result with that of Exercise 6. 8. Use the data presented in Exercises 5 and 7 to calculate the ratio of the masses of the earth and the sun, irrespective of the knowledge of G. 9. Jupiter has a moon, Ganymede, which orbits the planet at a distance 1:07£109m, with a period of 7.15 earth days. Using the method of Exercise 5 (or 8), compute the mass mJof Jupiter. You should get mJ¼318me: 7. Variational problems and the stationary action principle A rich source of second order systems of di®erential equations is provided by variational problems, which we will consider here. Let ­ ½Rnbe open, and let L2C2(­£Rn), say L=L(x; v). For a path u: [a; b]!­, consider (7.1) I(u) =Zb aL(u(t); u0(t))dt: We desire to ¯nd equations for a path that minimizes I(u), among all such paths for which the endpoints u(a) = pandu(b) = qare ¯xed. More generally, we desire to specify when uis a stationary path, meaning that (7.2)d dsI(us)¯¯¯ s=0= 0; for all smooth families of paths ussuch that u0=u; u s(a) =p, and us(b) = q. Let us write (7.3)@ @sus(t)¯¯¯ s=0=w(t); sow: [a; b]!Rnis an arbitrary smooth function such that w(a) =w(b) = 0. To compute ( d=ds)I(us), let us denote (7.4) Lxk=@L @xk; L vk=@L @vk: 7. Variational problems and the stationary action principle 291 Then (7.5)d dsI(us)¯¯¯ s=0=Zb aX kLxk(u(t); u0(t))wk(t)dt +Zb aX kLvk(u(t); u0(t))w0 k(t)dt: We can apply integration by parts to the last integral. The condition that wk(a) =wk(b) = 0 implies that there are no endpoint contributions, so (7.6)d dtI(us)¯¯¯ s=0=Zb aX kh Lxk(u(t); u0(t))¡d dtLvk(u(t); u0(t))i wk(t)dt: For this to vanish for all smooth wkthat vanish at t=aandb, it is necessary and su±cient that (7.7)d dtLvk(u(t); u0(t))¡Lxk(u(t); u0(t)) = 0 ;8k: This system is called the Lagrange equation for stationarity of (7.1). Ap- plying the chain rule to the ¯rst sum, we can expand this out as (7.8)X `Lvkv`(u(t); u0(t))u00 `(t) +X `Lvkx`(u(t); u0(t))u0 `(t) ¡Lxk(u(t); u0(t)) = 0 ;8k: This can be converted to a ¯rst order system for ( u(t); u0(t)), to which the results of x1 apply, provided the n£nmatrix (7.9)³ Lvkv`(x; v)´ of second order partial derivatives of L(x; v) with respect to vis invertible. The Newtonian equations of motion can be put into this Lagrangian framework, as follows. A particle of mass m, position x, and velocity v, moving in a force ¯eld F(x) =¡rV(x), has kinetic energy and potential energy (7.10) T=1 2mkvk2;and V=V(x); respectively. The Lagrangian L(x; v) is given by the di®erence : (7.11) L(x; v) =T¡V=1 2mkvk2¡V(x): 292 4. Nonlinear Systems of Di®erential Equations Figure 7.1 In such a case, (7.12) Lvk(x; v) =mvk; L xk(x; v) =¡@V @xk; and the Lagrange system (7.7) becomes the standard Newtonian system (7.13) md2u dt2=¡rV(u): In this setting, the integral (7.1) is called the action . The assertion that the laws of motion are given by the stationary condition for (7.1) where Lis the Lagrangian (7.11) is the stationary action principle. The Lagrangian approach can be particularly convenient in situations where coordinates other than Cartesian coordinates are used. As an exam- ple, we consider the simple pendulum problem, and give a treatment that can be compared and contrasted with that given in x6 of Chapter 1. As there, we have a rigid rod, of length `, suspended at one end. We assume the rod has negligible mass, except for an object of mass mat the other end. See Fig. 7.1. The rod makes an angle µwith the downward vertical. We seek a di®erential equation for µas a function of t. The end with the mass mtraces out a path in a plane, which, as in Chapter 1, we identify with the complex plane, with the origin at the point where the pendulum is suspended and the real axis pointing vertically down. We can write the path as (7.14) z(t) =`eiµ(t): 7. Variational problems and the stationary action principle 293 The velocity is (7.15) v(t) =z0(t) =i`µ0(t)eiµ(t); so the kinetic energy is (7.16) T=1 2mkv(t)k2=m`2 2µ0(t)2: Meanwhile the potential energy, due to the force of gravity, is (7.17) V=¡mg`cosµ: Taking Ã=µ0, we have the Lagrangian (7.18)L(µ; Ã) =m`2 2Ã2+mg`cosµ; LÃ(µ; Ã) =m`2Ã; L µ(µ; Ã) =¡mg`sinµ; and Lagrange's equation (7.19)d dtLÃ(µ(t); µ0(t))¡Lµ(µ(t); µ0(t)) = 0 yields the pendulum equation (7.20)d2µ dt2+g `sinµ= 0; in agreement with (6.6) of Chapter 1. The approach above avoided a computation of the force acting on the pendulum (cf. (6.4) of Chapter 1), and is arguably a bit simpler than the approach given in Chapter 1. The Lagrangian approach can be very much simpler in more complex situations, such as the double pendulum, which we will discuss in x9. An important variant of these variational problems is the class of con- strained variational problems , which we now discuss. For the sake of de¯- niteness, let Mbe either a smooth curve in ­ ½R2or a smooth surface in ­½R3, and let n(x) be a smooth unit normal to M, forx2M. Again, let L2C2(­£Rn); n= 2 or 3, and de¯ne I(u) by (7.1). We look for equations for (7.21) u: [a; b]¡!M; satisfying the stationary condition (7.2), not for all smooth families of paths ussuch that u0=uandus(0) = p; u s(b) =q, but rather for all such paths satisfying the constraint (7.22) us: [a; b]¡!M: 294 4. Nonlinear Systems of Di®erential Equations Again we take w(t) as in (7.3), and this time we obtain an arbitrary smooth function w: [a; b]!Rn, satisfying w(a) =w(b) = 0, and the additional constraint (7.23) w(t)¢n(u(t))´0: The calculations (7.4){(7.6) still apply, but from here we get a conclusion di®erent from (7.7). Since (7.6) holds for all w(t) described as just above, the conclusion is (7.24)d dtLv(u(t); u0(t))¡Lx(u(t); u0(t)) is parallel to n(u(t)); where Lv= (Lv1; : : : ; L vn)tandLx= (Lx1; : : : ; L xn)t. In case n= 3, an equivalent formulation of (7.24) is (7.25)hd dtLv(u(t); u0(t))¡Lx(u(t); u0(t))i £n(u(t)) = 0 : Let's specialize this constrained variational problem to the case (7.26) L(x; v) =1 2kvk2: The associated integral (7.27) E(u) =1 2Zb aku0(t)k2dt is called the energy ofu: [a; b]!M. In this case, Lv=vandLx= 0, so (7.24) becomes (7.28) u00(t) is parallel to n(u(t)): That is, u00(t) = a(t)n(u(t)). Taking the inner product with n(t) gives a(t) =n(u(t))¢u00(t), so (7.28) yields (7.29) u00(t) =n(u(t))¢u00(t)n(u(t)): An equation with a better form can be obtained by di®erentiating (7.30) u0(t)¢n(u(t))´0; to get (7.31) u00¢n(u(t)) =¡u0(t)¢d dtn(u(t)): 7. Variational problems and the stationary action principle 295 Plugging this into the right side of (7.29) gives the di®erential equation (7.32) u00(t) +u0(t)¢³d dtn(u(t))´ n(u(t)) = 0 : Note by (7.28) that u00is orthogonal to u0(t), so (7.33)d dtku0(t)k2= 2u0(t)¢u00(t)´0: Thus stationary paths u: [a; b]!Mfor the energy have constant speed. Such curves on Maregeodesics . These curves are also constant speed curves on Mthat are stationary curves for the arclength: (7.34) `(u) =Zb aku0(t)kdt: We will not go further into this here. The reader can consult texts on elementary di®erential geometry, such as [ DoC ], [Hen] ], or [ Op], or see [T], Chapter 1, x11. We next present another approach to ¯nding equations for stationary paths of (7.27). Suppose ­ = O £RandMis the graph of a function z='(x1; x2), for x= (x1; x2)2 O. Then a curve u: [a; b]!Mhas the form (7.35) u(t) =¡ x(t); '(x(t))¢ ; and (7.36) u0(t) = (x0(y);r'(x(t))¢x0(t)); so (7.37)ku0(t)k2=kx0(t)k2+ (r'(x(t))¢x0(t))2 =x0(t)¢G(x(t))x0(t); where (7.38) G(x) =µ1 +'1(x)2'1(x)'2(x) '1(x)'2(x) 1 + '2(x)2¶ ; ' j(x) =@' @xj: Thus the problem of ¯nding a constrained stationary path u(t) for the energy (7.27) is equivalent to the problem of ¯nding an unconstrained stationary path x(t) for (7.39) E(x) =1 2Zb ax0(t)¢G(x(t))x(t)dt: 296 4. Nonlinear Systems of Di®erential Equations In this case, (7.40)L(x; v) =1 2v¢G(x)v; Lv(x; v) =G(x)v;and; Lx(x; v) =1 2v¢ rG(x)v; where the last identity means (7.41) Lxk(x; v) =1 2v¢@G @xkv: In this setting, the Lagrange equation (7.7) becomes (7.42)d dth G(x(t))x0(t)i ¡1 2x0(t)¢ rG(x(t))x0(t) = 0 ; i.e., (7.43)d dtX jGkj(x(t))x0 j(t)¡1 2X i;jx0 i(t)@Gij @xkx0 j(t) = 0 ;8k: Exercises 1. Given a Lagrangian L(x; v), we de¯ne the \energy" (7.44)E(x; v) =Lv(x; v)¢v¡L(x; v) =X kLvk(x; v)vk¡L(x; v): Show that if u(t) solves the Lagrange equation (7.7), then (7.45)d dtE(u(t); u0(t))´0: This is energy conservation,onservation of energy in this setting. 2. Suppose (7.46) L(x; v) =m 2v¢G(x)v¡V(x); 7. Variational problems and the stationary action principle 297 Assume G(x)2M(n;R) is symmetric and invertible, and de¯ne E(x; v) as in (7.44). Show that (7.47) E(x; v) =m 2v¢G(x)v+V(x): 3. Let L(x; v) be given by (7.46). Show that the Lagrange equation (7.7) is (7.48) md dth G(u(t))u0(t)i ¡m 2u0(t)¢ rG(t)u0(t) =¡rV(u(t)); where the second term is evaluated as in (7.42){(7.43). Show in turn that this yields the ¯rst order system duk dt=vk mX jGkj(u(t))dvj dt+mX i;jvi(t)h@Gkj @xi¡1 2@Gij @xki vj(t) =¡@V @xk(u(t)): Produce a variant by symmetrizing the term in brackets in the second sum, with respect to iandj. 4. Consider the setting of constrained motion on M½­, as in (7.21){ (7.24), and consider the following generalization of (7.26): (7.49) L(x; v) =m 2kvk2¡V(x): Establish the following replacement for (7.32): (7.50) mu00(t) +mu0(t)¢³d dtn(u(t))´ n(u(t)) =¡PM(u(t))rV(u(t)); where, for x2M; w 2Rn, (7.51) PM(x)w=w¡³ n(x)¢w´ n(x): This describes motion of a particle in a force ¯eld F(x) =¡rV(x), constrained to move on M. 5. Motion of a spherical pendulum in R3, in the presence of Earth's grav- itational ¯eld, is described as in Exercise 4 with (7.52) M=fx2R3:kxk=`g; 298 4. Nonlinear Systems of Di®erential Equations andL(x; v) as in (7.49), with V(x) =mg(x¢k), where k= (0;0;1)t. Show that in this case, (7.50) produces, for (7.53) u(t) =`!(t); the system (7.54) !00(t) +k!0(t)k2!(t) =¡g `k+g `(!(t)¢k)!(t): 6. Results of Exercise 5 are also valid in the setting where R3is replaced byR2. Show that, in this setting, with (7.55) !(t) = (sin µ(t);¡cosµ(t))t; k = (0;1)t; the equation (7.54) leads to the (planar) pendulum equation (7.56) µ00(t) +g `sinµ(t) = 0 : 7. Let us return to the setting of Exercise 2, and set (7.57) p=Lv(x; v) =mG(x)v: Also set (7.58) E(x; p) =E(x; v) =E(x; G(x)¡1p=m): Show that (7.59) E(x; p) =1 2mp¢G(x)¡1p+V(x): Show that the Lagrange equation (7.48) for u(t) =x(t) is equivalent to the following Hamiltonian system: (7.60)dxk dt=@E @pk;dpk dt=¡@E @xk: Hint. To get started on (7.60), note that if (7.59) holds, then (7.61)@E @p=1 mG(x)¡1p=v; 7. Variational problems and the stationary action principle 299 and that the Lagrange equation implies (7.62)dpk dt=Lxk(x; v) =m 2v¢@G @xk(x)v¡@V @xk(x): Furthermore, as in (8.13) of Chapter 3, (7.63)@ @xkG(x)¡1=¡G(x)¡1@G @xk(x)G(x)¡1: Remark. More general cases in which the change of variable p=Lv(x; v) converts Lagrange's equation to Hamiltonian form are discussed in [ AM], [Ar], and Chapter 1 of [ T]. Exercises 8{11 study sufaces of revolution that are surfaces of \least area." To set this up, let u: [0;1]!(0;1) be smooth, and rotate the graph of y=u(x) about the x-axis in ( x; y; z )-space. Elementary calculus gives the formula (7.64) A(u) = 2 ¼Z1 0u(t)p 1 +u0(t)2dt for the area of the resulting surface of revolution. The problem is to ¯ndufor which the area is minimal, given constraints (7.65) u(0) = ®; u (1) = ¯; ®; ¯ > 0: 8. In (7.64), L(x; v) =xp 1 +v2. Show that the \energy" E(x; v) in (7.44) is given by (7.66) E(x; v) =¡xp 1 +v2: 9. Using (7.45), show that if u(t) solves the Lagrange equation (7.7) in this setting, then there is a constant asuch that (7.67)u(t)p 1 +u0(t)2=a; hence (7.68)du dt=§p b2u2¡1; b =1 a: 300 4. Nonlinear Systems of Di®erential Equations 10. Separate variables in (7.68) and use the substitution bu= cosh vto evaluate the u-integral and conclude that (7.69) u(t) =1 bcosh( bt+c); for some constant c. Equation (7.69) is the equation of a catenary, seen before in (3.24) of Chapter 1, for the hanging cable. 11. Consider the problem of ¯nding bandcin (7.69) such that the con- straints (7.65) are satis¯ed. Show that sometimes no solutions exist, and sometimes two solutions exist, but one gives a smaller area than the other. Exercises 12{15 take another look at the hanging cable problem men- tioned in Exercise 10. Here we state it as the problem of minimizing the potential energy, which is mgtimes (7.70) V(u) =ZA ¡Au(t)p 1 +u0(t)2dt; subject to the boundary conditions (7.71) u(¡A) =u(A) = 0 ; and the constraint that the curve y=u(x);¡A·x·A, have length L, (7.72) `(u) =ZA ¡Ap 1 +u0(t)2dt=L: Such a curve describes a cable, of length L, hanging from the two points (¡A;0) and ( A;0), under the force of gravity. To deal with the con- straint (7.72), we bring in the Lagrange multiplier method. That is, we set (7.73) I¸(u) =V(u) +¸`(u); ¯nd the stationary path for (7.73) (subject to (7.71)) as a function of ¸, and then ¯nd for which ¸the constraint (7.72) holds. Note that I¸(u) has the form (7.1) with (7.74) L¸(x; v) = (x+¸)p 1 +v2: 7. Variational problems and the stationary action principle 301 12. Show that the \energy" E¸(x; v) in (7.44) is given by (7.75) E¸(x; v) =x+¸p 1 +v2: 13. Using (7.45), show that if u(t) solves the Lagrange equation (7.7) in this setting, then there exists a constant a(maybe depending on ¸) such that (7.76)u(t) +¸p 1 +u0(t)2=a; hence (7.77)du dt=§p b2(u+¸)2¡1; b =1 a: 14. Separate variables in (7.77) and use the substitution b(u+¸) = cosh v to evaluate the u-integral and obtain u(t) =¡¸+1 bcosh( bt+c); for some constant c. Show that (7.71) forces c= 0, so (7.78) u(t) =¡¸+1 bcoshbt: 15. Calculate the length of the curve y=u(x);¡A·x·A, when uis given by (7.78), and show that the constraints (7.71){(7.72) yield the equations (7.79) sinh bA=bL 2; ¸ =1 bcoshbA: Note that the ¯rst equation has a unique solution b2(0;1) if and only ifL >2A. 16. Recall the planar pendulum problem illustrated in Fig. 7.1. Instead of assuming all the mass is at the end of the rod, assume the rod has a mass 302 4. Nonlinear Systems of Di®erential Equations Figure 8.1 distribution m(s)ds;0·s·`, so the total mass is m=R` 0m(s)ds. Show that for the potential energy Vyou replace (7.17) by (7.80) V=¡mag`cosµ; m a=Z` 0m(s)s `ds; and for the kinetic energy T, you replace (7.16) by (7.81) T=mb`2 2µ0(t)2; m b=Z` 0m(s)³s `´2 ds: Write down the replacement for the pendulum equation (7.20) in this setting. Specialize the calculation to the case (7.82) m(s) =m `;0·s·`; which represents a rod with uniform mass distribution. 8. The brachistochrone problem The early masters of calculus enjoyed posing challenging problems to each other. The most famous of these is called the brachistrochrone problem . It was posed by Johann Bernoulli in 1696, and solved by him, by his brother Jakob, and also by Newton and by Leibniz. The problem is to ¯nd the curve along which a particle will slide without friction in the minimum time, from one given point pin the ( x; y)-plane to another, q, starting at rest at p. Say p= (0;0) and q= (a; b). We assume a >0 and b <0; see Fig. 8.1. The force of gravity acts in the direction of the negative y-axis, with acceleration g. 8. The brachistochrone problem 303 Our approach to this problem will involve two applications of the vari- ational method developed in x7. (In fact, this problem helped spark the creation of the variational method.) First, let ': [0; a]!Rwith '(0) = 0; '(a) =b, and consider the constrained motion of a particle, (8.1) u: [0; t0]¡!M=f(x; '(x)) : 0·x·ag; under the force of gravity. Thus, in place of (7.27), we look for stationary paths for (8.2) I(u) =Za 0hm 2ku0(t)k2¡V(u(t))i dt; subject to the constraint (8.1), and with (8.3) V(x; y) =mgy: We can convert this to an unconstrained variational problem as was done in (7.35){(7.42), now with a nonzero V, and with lower dimension. We have (8.4) u(t) =¡ x(t); '(x(t))¢ ; and (8.5) ku0(t)k2=¡ 1 +'0(x(t))2¢ x0(t)2; so the problem of ¯nding a constrained stationary path u(t) for (8.2) is equivalent to the problem of ¯nding an unconstrained stationary path x(t) for (8.6) J(x) =Za 0L(x(t); x0(t))dt; with (8.7) L(x; v) =m 2¡ 1 +'0(x)2¢ v2¡mg'(x): The path x(t) is governed by the di®erential equation (8.8)d dtLv(x(t); x0(t))¡Lx(x(t); x0(t)) = 0 : We need not write this more explicitly, since by now our experience tells us that to describe solutions to such a single equation, all we need is conserva- tion of energy: (8.9) E(x; v) =m 2¡ 1 +'0(x)2¢ v2+mg'(x); 304 4. Nonlinear Systems of Di®erential Equations that is, for a solution to (8.8), (8.10)m 2¡ 1 +'0(x(t))2¢ x0(t)2+mg'(x(t)) =E is constant. In the current set-up, x(0) = 0 and x0(0) = 0, so E= 0. We get (8.11)dx dt=§s ¡2g'(x) 1 +'0(x)2; which separates to (8.12)1p2gZa 0s 1 +'0(x)2 ¡'(x)dx=Zt0 0dt: In other words, the elapsed time for the particle to move from p= (0;0) to q= (a; b) along the path y='(x) is given by the left side of (8.12). Hence the brachistochrone problem is reduced to the problem of ¯nding ': [0; a]¡!R, minimizing (8.13) K(') =Za 0L('(x); '0(x))dx; subject to the condition (8.14) '(0) = 0 ; '(a) =b; where (8.15) L('; Ã) =s 1 +Ã2 ¡': Stationary paths for (8.13) satisfy the Lagrange equation (8.16)d dtLÃ('(x); '0(x))¡ L'('(x); '0(x)) = 0 : Note that (8.17) LÃ('; Ã) =Ãp ¡'(1 +Ã2); L '('; Ã) =¡1 2p ¡'(1 +Ã2) '2: Solutions to (8.16) have the property that (8.18) E('(t); '0(t)) =E 8. The brachistochrone problem 305 is constant, where (parallel to (7.44)) (8.19) E('; Ã) =LÃ('; Ã)¡ L('; Ã): Using (8.15) and (8.17), we have (8.20)E('; Ã) =Ã2 p ¡'(1 +Ã2)¡s 1 +Ã2 ¡' =¡1p ¡'(1 +Ã2): Thus, if '(x) satis¯es (8.16), then (8.21) '(x)¡ 1 +'0(x)2¢ =¡k2;const. ; where we have written the constant as ¡k2to enforce the condition that '(x)<0 for 0 < x·a. For notational convenience, we make the change of variable (8.22) y(x) =¡'(x); so (8.21) becomes (8.23) y(x)¡ 1 +y0(x)2¢ =k2; giving (8.24)dy dx=s k2 y¡1: The equation (8.24) separates to (8.25)Zdyq k2 y¡1=Z dx: The left integral has the form of (5.15) in Chapter 1, with E0=¡1; Km = k2. Rather then recall the formulas (5.16){(5.22) of Chapter 1, we implement the method previewed in Exercise 3 of that section. We use the change of variable (8.26) y=k2sin2¿;2¿=µ: Then (8.27) dy= 2k2sin¿cos¿ d¿;s k2 y¡1 =cos¿ sin¿; 306 4. Nonlinear Systems of Di®erential Equations Figure 8.2 so (8.28)Zdyq k2 y¡1= 2k2Z sin2¿ d¿ =k2 2Z (1¡cosµ)dµ =k2 2(µ¡sinµ); the second identity because sin2¿= (1¡cos 2¿)=2. Thus the curve ( x; y(x)); x2 [0; a], is parametrized by (8.29)x=x(µ) =k2 2(µ¡sinµ); y=y(µ) =k2 2(1¡cosµ): The choice of k2>0 is dictated by the implication (8.30) 0 < µ < ¼k2;k2 2(µ¡sinµ) =a=)k2 2(1¡cosµ) =jbj: This solves the brachistochrone problem. The curve de¯ned by (8.29) is known as a cycloid . See Fig. 8.2. Here ½=k2=2. Remark. Note that y0(0) = + 1, so the optimal path starts directly down. 9. The double pendulum 307 Exercises 1. Show that for each a;jbj 2(0;1), there is a unique k2>0 such that (a;jbj)2R2 +lies on the curve (8.29), for some µ2(0; ¼k2). Hint. Consult Fig. 8.2. 2. In the setting of Exercise 1, show that if jbj=a < 2=¼, then µ > ¼k2=2, and the optimal path dips below bbefore reaching the endpoint q= (a; b). 3. With x(µ) and y(µ) as in (8.29), set '(µ) =¡y(µ). Let (8.31) µ1=k2 2¼; µ 02[0; µ1): Show that the time it takes a particle starting at rest at ( x(µ0); '(µ0)) to slide down the curve ( x(µ); '(µ)); µ0·µ·µ1, to the point ( x(µ1); '(µ1)) (the bottom of the cycloid) is independent of µ0. One says the cycloid also solves the tautochrone problem . 9. The double pendulum Here we study the motion of a double pendulum, such as illustrated in Fig. 9.1. We have a pair of rigid rods, of lengths `1and`2, of negligible mass except for objects of mass m1andm2attached to one end of each rod. The other end of rod 1 is attached to a ¯xed point, and the end of rod 2 not containing mass 2 is attached to rod 1 at mass 1. The rods are assumed free to swing back and forth in a plane. Thus the con¯guration at time tis described by the angles µ1(t) and µ2(t), that the rods make with the vertical. Gravity acts on the masses mj, with a downward force of mjg. We identify the plane mentioned above with the complex plane, with rod 1 attached to the origin and the real axis pointing down. Thus the position of mass 1 is (9.1) z1(t) =`1eiµ1(t); and the position of mass 2 is (9.2) z2(t) =z1(t) +`2eiµ2(t): Their velocities are (9.3)z0 1=i`1µ0 1eiµ1; z0 2=i`1µ0 1eiµ1+i`2µ0 2eiµ2; 308 4. Nonlinear Systems of Di®erential Equations Figure 9.1 with square norms (9.4)jz0 1j2=`2 1(µ0 1)2; jz0 2j2= (`1µ0 1eiµ1+`2µ0 2eiµ2)(`1µ0 1e¡iµ1+`2µ0 2e¡iµ2) =`2 1(µ0 1)2+`2 2(µ0 2)2+ 2`1`2µ0 1µ0 2cos(µ1¡µ2): The potential energy of this system is given by (9.5)V=¡m1gRez1(t)¡m2gRez2(t) =¡m1g`1cosµ1¡m2g(`1cosµ1+`2cosµ2); and the kinetic energy by (9.6) T=m1 2jz0 1(t)j2+m2 2jz0 2(t)j2: If we write (9.7) µ=µµ1 µ2¶ ; à =µÃ1 Ã2¶ =µµ0 1 µ0 2¶ ; then (9.4) gives (9.8) T=1 2âG(µ)Ã; 9. The double pendulum 309 with (9.9) G(µ) =µ(m1+m2)`2 1 m2`1`2cos(µ1¡µ2) m2`1`2cos(µ1¡µ2) m2`2 2¶ : Thus the Lagrangian L=T¡Vis given by (9.10) L(µ; Ã) =1 2âG(µ)áV(µ); with V(µ) as in (9.5), and the equation of motion for the double pendulum is (9.11)d dtLÃ(µ; µ0)¡Lµ(µ; µ0) = 0 : As in (7.48), this expands out to the 2 by 2 system (9.12)d dtX jGkj(µ(t))µ0 j(t)¡1 2X i;jµ0 i(t)@Gij @µkµ0 j(t) =¡@V @µk(µ(t)); fork= 1;2. Making explicit use of (9.5) and (9.9), we have (9.13)LÃ1(µ; Ã) = (m1+m2)`2 1Ã1+m2`1`2Ã2cos(µ1¡µ2); LÃ2(µ; Ã) =m2`2 2Ã2+m2`1`2Ã1cos(µ1¡µ2); and (9.14)Lµ1(µ; Ã) =¡m2`1`2Ã1Ã2sin(µ1¡µ2)¡(m1+m2)g`1sinµ1; Lµ2(µ; Ã) =m2`1`2Ã1Ã2sin(µ1¡µ2)¡m2g`2sinµ2: Thus the explicit version of (9.11){(9.12) is the pair of equations (9.15)(m1+m2)`2 1µ00 1+m2`1`2d dth µ0 2cos(µ1¡µ2)i =¡m2`1`2µ0 1µ0 2sin(µ1¡µ2)¡(m1+m2)g`1sinµ1; and (9.16)`2 2µ00 2+`1`2d dth µ0 1cos(µ1¡µ2)i =`1`2µ0 1µ0 2sin(µ1¡µ2)¡g`2sinµ2: Note that the masses m1andm2do not appear in (9.16); m1does not appear in either term of ( d=dt)LÃ2¡Lµ2, and m2factors out. 310 4. Nonlinear Systems of Di®erential Equations As in (7.44){(7.47), we have the energy (9.17) E(µ; Ã) =1 2âG(µ)Ã+V(µ); and if µ(t) solves (9.11), or equivalently (9.15){(9.16), then (9.18)d dtE(µ(t); µ0(t)) = 0 : By (9.5) and (9.9), the explicit form of the energy is (9.19)E(µ; Ã) =1 2(m1+m2)`2 1Ã2 1+m2`1`2Ã1Ã2cos(µ1¡µ2) +1 2m2`2 2Ã2 2¡m1g`1cosµ1¡m2g(`1cosµ1+`2cosµ2): As in (7.57){(7.60), we can convert the equations of motion to Hamil- tonian form, by setting (9.20) p=G(µ)Ã: The energy (9.17) becomes (9.21)E(µ; p) =E(µ; G(µ)¡1p) =1 2p¢G(µ)¡1p+V(µ); and (9.11) is equivalent to (9.22)dµk dt=@E @pk;dpk dt=¡@E @µk: Note that, for G(µ) given by (9.9), (9.23) G(µ)¡1=1 detG(µ)µm2`2 2 ¡m2`1`2cos(µ1¡µ2) ¡m2`1`2cos(µ1¡µ2) ( m1+m2)`2 1¶ ; and (9.24) det G(µ) =m1m2`2 1`2 2+m2 2`2 1`2 2sin2(µ1¡µ2): For notational simplicity we write (9.25) E(µ; p) =1 2p¢H(µ)p+V(µ); H (µ) =G(µ)¡1: 9. The double pendulum 311 Solutions to (9.22) are orbits of the °ow generated by the Hamiltonian vector ¯eld (9.26)XE(µ; p) =¡Jrµ;pE(µ; p) =µ0I ¡I0¶µrµE rpE¶ =µrpE ¡rµE¶ : Here I2M(2;R) is the identity matrix and J2M(4;R) is de¯ned by the second identity in (9.26). From this formula we see that the critical points ofXEcoincide with the critical points of E. Note that (9.27) rpE(µ; p) =H(µ)p; andH(µ) is invertible for all µ, so if Ehas a critical point at ( µ; p); p= 0. Now (9.28) rµE(µ;0) =rV(µ); so we deduce that ( µ; p) is a critical point of XEif and only if p= 0 and rV(µ) = 0. Rewriting (9.5) as (9.29) V(µ) =¡(m1+m2)g`1cosµ1¡m2g`2cosµ2; we see that (9.30) rV(µ) =µ(m1+m2)g`1sinµ1 m2g`2sinµ2¶ ; so the critical points of Vconsist of µ1=j¼; µ 2=k¼; j; k 2Z. In summary, the critical points of XEconsist of (9.31) ( µ1; µ2; p1; p2) = (j¼; k¼; 0;0); j; k 2Z: Towards the goal of understanding the behavior of XEnear these critical points, we examine its derivative. We have (9.32) DXE(µ;0) =µ0 H(µ) ¡D2V(µ) 0¶ : The matrix H(µ) is positive de¯nite for all µ, and in particular, since sin j¼= 0 and cos j¼= (¡1)j, (9.33) H(j¼; k¼ ) =1 m1m2`2 1`2 2µm2`2 2 (¡1)j¡k+1m2`1`2 (¡1)j¡k+1m2`1`2 (m1+m2)`2 1¶ : Also, (9.34) D2V(j¼; k¼ ) =µ(¡1)j(m1+m2)g`1 0 0 ( ¡1)km2g`2¶ : We are set up to examine the linearization of the °ow generated by XEat the critical points. This will be pursued, in a more general setting, in the next section. 312 4. Nonlinear Systems of Di®erential Equations Exercises 1. Pass to the limit m2!0 in the double pendulum system (9.15){(9.16) and derive the limiting system µ00 1+g `1sinµ1= 0;(9.35) µ00 2+`1 `2d dth µ0 1cos(µ1¡µ2)i =`1 `2µ0 1µ0 2sin(µ1¡µ2)¡g `2sinµ2:(9.36) 2. Recall the spherical pendulum, introduced in Exercise 5 of x7. Derive equations of motion for a double spherical pendulum. 3. Instead of assuming all the mass of rods 1 and 2 is concentrated at an end, assume that rod jhas mass distribution mj(s)ds;0·s·`j, so the total mass of rod jismj=R`j 0mj(s)ds; j = 1;2. Obtain formulas for the potential and kinetic energy, replacing (9.5) and (9.6), and then obtain equations of motion, replacing (9.15){(9.16). Note. See Exercise 16 in x7 to get started. 10. Momentum-quadratic Hamiltonian systems Most of the Lagrangians arising in the last three sections have been of the form (10.1) L(x; v) =1 2v¢G(x)v¡V(x); forx2­½Rn; v2Rn, where G(x)2M(n;R) is symmetric and invertible, in fact positive de¯nite, but for awhile we will work in this more general setting. As exercises in x7 have revealed, making the change of variables (x; v)7!(x; p) with p=G(x)v, one can convert the Lagrange system of di®erential equations to Hamiltonian form, (10.2)dxk dt=@E @pk;dpk dt=¡@E @xk; where (10.3) E(x; p) =1 2p¢H(x)p+V(x); H (x) =G(x)¡1: 10. Momentum-quadratic Hamiltonian systems 313 We call such systems momentum-quadratic Hamiltonian systems. Note that H(x) is also symmetric and invertible, and furthermore positive de¯nite if G(x) is. Solutions of (10.2) are orbits of the °ow generated by the Hamil- tonian vector ¯eld (10.4)XE(x; p) =¡Jrx;pE(x; p) =µ0I ¡I0¶µrxE rpE¶ =µrpE ¡rxE¶ : Here, I2M(n;R) is the identity matrix, and J2M(2n;R) is de¯ned by the second identity in (10.4). We record some general results about the critical points of such ¯elds, and their linearizations. To begin, the critical points of XEcoincide with the critical points of E. Note that (10.5) rpE(x; p) =H(x)p; so, since H(x) is invertible, we see that if Ehas a critical point at ( x; p), then p= 0. Now (10.6) rxE(x;0) =rV(x); so we deduce that the critical points of XEconsist of (10.7) f(x;0) :rV(x) = 0g: We next look at the linearization (cf. (3.32)) of XEat a critical point (x0;0), given by (10.8) DXE(x0;0) =µ0 H(x0) ¡D2V(x0) 0¶ : From here on, we assume H(x0) is positive de¯nite. For notational simplic- ity, we set (10.9) H=H(x0); W =D2V(x0); L =µ0H ¡W 0¶ : Then the linearization of (10.2) at ( x0;0) is (10.10)dx dt=Hp;dp dt=¡Wx: 314 4. Nonlinear Systems of Di®erential Equations To analyze the structure of solutions to (10.10), it is convenient to di- rectly tackle the second order system (10.11)d2x dt2=¡HWx; and to do this we bring in the following. Lemma 10.1. Given that H2M(n;R)is positive de¯nite, there exists a positive de¯nite A2M(n;R)such that (10.12) H=A2: Proof. From Chapter 2 we know that Rnhas an orthonormal basis fvjgof eigenvectors of H, soHvj=¸jvj;1·j·n. Each ¸jis positive, so we can de¯ne AbyAvj=p ¸jvj;1·j·n. If we make the change of variable (10.13) x=Ay; then (10.11) is converted to (10.14) y00+AWAy = 0: Note that W2M(n;R) is symmetric and so is AWA . Also AWA is invert- ible if and only if Wis. This invertibility is equivalent to the assertion that (x0;0) is a nondegenerate critical point of XE. We restrict attention to such cases. The following result will be useful. Lemma 10.2. LetW2M(n;R)be a symmetric matrix, and assume (10.15) Whaskpositive and n¡knegative eigenvalues. Then so does AWA , when A2M(n;R)is positive de¯nite. Proof. WriteRn=W+©W¡, where W+is the linear span of the eigenvec- tors of Wwith positive eigenvalue, W¡the linear span of the eigenvectors ofWwith negative eigenvalue. Similarly, write Rn=fW+©fW¡, with W replaced by AWA . The image AfW+offW+under Ais a linear subspace of Rn, and (10.16) v=Aw2AfW+=)v¢Wv=w¢AWAw ¸0 =)v2 W +: 10. Momentum-quadratic Hamiltonian systems 315 Thus (10.17) A:fW+¡! W +;injectively ; so (10.18) dim fW+·dimW+: A similar argument gives (10.19) dim fW¡·dimW¡; and ¯nishes the proof. To continue, under the hypotheses of Lemma 10.2, we have an orthonor- mal basis fu1; : : : ; u ngofRnsuch that, with ¹j2(0;1), (10.20)AWAu j=¹2 juj; j ·k; AWAu j=¡¹2 juj; j > k: in such a case, the general solution to (10.14) is (10.21)y(t) =X j·k(ajsin¹jt+bjcos¹jt)uj +X j>k(aje¹jt+bje¡¹jt)uj: Such y(t) leads to (10.22)µAy(t) A¡1y0(t)¶ =µx(t) p(t)¶ =etLµv0 v1¶ ; for general v0; v12Rn. As a result, we have the following. Proposition 10.3. Under the hypotheses of Lemma 10.2, L, given by (10.9), is diagonalizable, and its eigenvalues are (10.23)§i¹jfor j·k; §¹jfor j > k: Proof. The eigenvalues of Lare what appear in the exponents in the matrix coe±cients of etL. IfLwere not diagonalizable, some matrix coe±cients would also contain terms of the form t`et¸; `¸1, where ¹=§i¹jor§¹j in (10.23), depending on j. A critical point of XEis said to be hyperbolic if all of the eigenvalues of DXEhave nonzero real part. From the analysis above, we have the following. 316 4. Nonlinear Systems of Di®erential Equations Proposition 10.4. A critical point (x0;0)ofXEis hyperbolic if and only if (10.24) D2V(x0)is negative de¯nite. If (10.24) holds, DXE(x0;0)hasnpositive eigenvalues and nnegative eigen- values. Whenever a vector ¯eld X(Hamiltonian or not) has a hyperbolic critical point, say at z0, the phase portrait near z0for the °ow generated by Xhas a similar appearance to that for the °ow generated by its linearization at z0. This is a generalization of the two dimensional result mentioned below (3.69). See Appendix C for further discussion. The opposite extreme can also be read o® from (10.23). Proposition 10.5. At a critical point (x0;0)ofXE, all the eigenvalues of DXEare purely imaginary if and only if (10.25) D2V(x0)is positive de¯nite. Recalling that E(x; p) is given by (10.3), we see that (10.25) is equivalent to (10.26) D2E(x0;0)2M(2n;R) is positive de¯nite, in which case Ehas a local minimum at ( x0;0). In case (10.25) holds, we can deduce from (10.21){(10.22), with k=n, that the orbits of etLall lie in n-dimensional tori. As for the °ow generated byXEitself, we know that its orbits all lie on level surfaces of E. Near (x; p) = ( x0;0), these level sets look like (2 n¡1)-dimensional spheres in Rn. In case n= 1, these are closed curves in R2, and indeed the phase portrait for the °ow generated by XEnear ( x0;0) looks like that for the °ow generated by its linearization. In such a case, ( x0;0) is a center , discussed inx3. In case n > 1, the orbits of the °ow generated by XEnear ( x0;0) do not necessarily lie on n-dimensional tori. The analysis of this behavior is much more subtle than in the case of hyperbolic critical points. There will be n-dimensional invariant tori that are invariant under the °ow, arising rather densely near ( x0;0), but the °ow generated by XEoften has chaotic behavior on the complement of these tori. Study of this situation is part of the deep Kolmogorov-Arnold-Moser (KAM) theory. Discussion of this, and references to further work, can be found in [ AM], Chapter 8, and [ Ar], Appendices 7{8. Forn¸2, there can be cases intermediate between those covered by Proposition 10.4 and those covered by Proposition 10.5. 10. Momentum-quadratic Hamiltonian systems 317 Proposition 10.6. If(x0;0)is a critical point for XEand (10.27) D2V(x0)haskpositive eigenvalues and n¡knegative eigenvalues ; then (10.28)DXE(x0;0)has2kimaginary eigenvalues, and n¡kpositive, and n¡knegative eigenvalues. In such cases, with k¸1 and n¸2, the phase portrait for the °ow generated by XEnear ( x0;0) will generally di®er from that of its linearization in important details, with some exceptions, arising when XEis \integrable." We refer to the sources cited above for more on this. Let us specialize these results to the case of the double pendulum, dis- cussed in x9. There Vwas given by (9.29), and the critical points by (9.31), i.e., ( j¼; k¼; 0;0), and D2V(j¼; k¼ ) by (9.34). We have (10.29)jandkeven = )D2V(j¼; k¼ ) positive de¯nite ; jandkodd = )D2V(j¼; k¼ ) negative de¯nite ; jandkof opposite parity = )D2V(j¼; k¼ ) inde¯nite : In the ¯rst case Proposition 10.5 applies, in the second case Proposition 10.4 applies, and in the third case Proposition 10.6 applies, with k= 1 and n¡k= 1. Exercises 1. Establish analogues of Propositions 10.3, 10.5, and 10.6 in case His allowed to be inde¯nite (nondegenerate), and we assume (10.30) D2V(x0) is either positive de¯nite or negative de¯nite. Exercises 2{6 deal with the 2 £2 system (10.31)d2 dt2µx y¶ =¡rx;yV(x; y); for various functions V. The associated energy function, as in (10.3), is (10.32) E(x; y; p; q ) =1 2(p2+q2) +V(x; y): 318 4. Nonlinear Systems of Di®erential Equations In each case, do the following. (a) Find all the critical points of E. (b) Determine the type of each critical point of E. (c) Determine the behavior of the eigenvalues of DXEat each such critical point (via Proposition 10.6). 2. Take V(x; y) = (cos x)(cos y): 3. Take V(x; y) =x2+xy+y4: 4. Take V(x; y) =x4+xy+y4: 5. Take V(x; y) =x4¡xy+y4: 6. Take V(x; y) =x4¡x2y+y4: 7. Do analogues of Exercises 2{6 with (10.32) replaced by (10.33) E(x; y; p; q ) =1 2(p2¡q2) +V(x; y): Now Proposition 10.6 will not apply, but Exercise 1 might (or might not). 11. Numerical study { di®erence schemes We describe some ways of numerically approximating the solution to a sys- tem of di®erential equations (11.1)dx dt=F(x); x(t0) =x0: 11. Numerical study { di®erence schemes 319 Higher order systems can be transformed to ¯rst order systems and treated by these methods, which are known as di®erence schemes. To start, we pick a time step hand attempt an approximation to the solution to (11.1) at times t0+nh: (11.2) xn¼x(t0+nh): Noting that a smooth solution to (11.1) satis¯es (11.3)x(t+h) =x(t) +hx0(t) +O(h2) =x(t) +hF(x(t)) +O(h2); we have the following crude di®erence scheme: (11.4) xn+1=xn+hF(xn): This is said to be ¯rst order accurate, meaning that over an interval of unit length one carries out 1 =hsuch operations, each with error O(h2), giving an accumulated error O(h), i.e., on the order of hto the ¯rst power. This method of approximating the solution x(t) is often called the Euler method, though considering what a great master of computation Euler was, it is hard to believe he actually took it seriously. Shortly we will present a fourth order accurate method, which is generally satisfactory, after describing some second order accurate methods. These better di®erence schemes will be suggested by higher order ac- curate methods of numerical integration. The connection between the two comes from rewriting (11.1) as (11.5) x(t+h) =x(t) +Zh 0F(x(t+s))ds: Consider methods of approximating (11.6)Zh 0g(s)ds better than hg(0) + O(h2), for smooth g. Two simple improvements are (11.7)h 2h g(0) + g(h)i +O(h3); the trapezoidal method, and (11.8) hg³h 2´ +O(h3); 320 4. Nonlinear Systems of Di®erential Equations the midpoint method. These lead respectively to (11.9) x(t+h) =x(t) +h 2h F(x(t)) +F(x(t+h))i +O(h3) and (11.10) x(t+h) =x(t) +hF³ x³ t+h 2´´ +O(h3): Neither of them immediately converts to an explicit di®erence scheme, but in (11.9) we can substitute F(X(t+h)) =F¡ X(t)+hF(X(t))¢ +O(h2) and in (11.10) we can substitute F¡ X(t+h=2)¢ =F¡ X(t)+(h=2)F(X(t))¢ +O(h2); to obtain the second order accurate di®erence schemes (11.11) xn+1=xn+h 2h F(xn) +F¡ xn+hF(xn)¢i and (11.12) xn+1=xn+hF³ xn+h 2F(xn)´ : Often (11.11) is called Heun's method and (11.12) a modi¯ed Euler method. We now come to the heart of the matter for this section. The Runge- Kutta scheme for (11.1) is speci¯ed as follows. The approximation xnto x(t0+nh) is given recursively by (11.13) xn+1=xn+h 6³ Kn1+ 2Kn2+ 2Kn3+Kn4´ ; where (11.14)Kn1=F(xn); Kn2=F³ xn+1 2hKn1´ ; Kn3=F³ xn+1 2hKn2´ ; Kn4=F(xn+hKn3): This scheme is 4th order accurate. It is one of the most popular and impor- tant di®erence schemes used for numerical studies of systems of di®erential equations. We make some comments about its derivation. We will consider a method of deriving 4th order accurate di®erence schemes, based on Simpson's formula (11.15)Zh 0g(s)ds=h 6³ g¡ 0¢ + 4g¡h 2¢ +g¡ h¢´ +O(h5): 11. Numerical study { di®erence schemes 321 This formula is derived by producing a quadratic polynomial p(s) such that p(s) =g(s) ats= 0; h=2;andh;and then exactly integrating p(s):The formula can be veri¯ed by rewriting it as (11.16)Zh ¡hG(s)ds=h 3h G(¡h) + 4G(0) + G(h)i +O(h5): The main part on the right is exact for all odd G(s), and it is also exact for G(s) = 1 and G(s) =s2;so it is exact when G(s) is a polynomial of degree ·3:Making a power series expansion G(s) =P3 j=0ajsj+O(s4) then yields (11.16). Now, write the equation (11.1) as the integral equation (11.5). By (11.15), (11.17)Zh 0F(X(t+s))ds=h 6h F(X(t)) + 4 F³ X³ t+h 2´´ +F(X(t+h))i +O(h5): We then have as an immediate consequence the following result on producing accurate di®erence schemes. Proposition 11.1. Suppose the approximation (11.18) x(t+h)¼x(t) + ©( x(t); h) =X(x(t); h) produces a jth order accurate di®erence scheme for the solution to (11.1). Ifj·3;then a di®erence scheme accurate of order j+ 1is given by (11.19) xn+1=xn+h 6h F(xn) + 4F³ X³ xn;h 2´´ +F(X(xn; h))i : Furthermore, if x(t+h)¼ X `(x(t); h)both work in (11.18), `= 0;1;then you can use (11.20) xn+1=xn+h 6h F(xn) + 4F³ X0³ xn;h 2´´ +F(X1(xn; h))i : We apply this to two second order methods derived before: (11.21) X0(xn; h) =xn+h 2h F(xn) +F(xn+hF(xn))i ;Heun, and (11.22) X1(xn; h) =xn+hF³ xn+h 2F(xn)´ ;modi¯ed Euler. 322 4. Nonlinear Systems of Di®erential Equations Thus a third order accurate scheme is produced. The last term in (11.19) becomes (11.23) h 6h F(xn)+4F³ xn+h 4h F(xn)+F³ xn+h 2F´i´ +F³ xn+hF³ xn+h 2F´´i ; where F=F(xn):In terms of Kn1; Kn2as de¯ned in (11.14), we have (11.24)h 6h Kn1+ 4F³ xn+h 4[Kn1+Kn2]´ +F(xn+hKn2)i : This could be used in a 3rd order accurate scheme, but some simpli¯cation of the middle term is desirable. Note that, for smooth H; (11.25) H³ x+1 2´´ =1 2H(x) +1 2H(x+´) +O(j´j2): Consequently, as jKn1¡Kn2j=O(h);by (11.14), (11.26) F³ xn+h 4[Kn1+Kn2]´ =1 2F³ xn+h 2Kn1´ +1 2F³ xn+h 2Kn2´ +O(h4): Therefore we have the following. Proposition 11.2. A third order accurate di®erence scheme for (11.1) is given by (11.27) xn+1=xn+h 6[Kn1+ 2Kn2+ 2Kn3+Ln4] where Kn1; Kn2; Kn3are given by (11.14) and (11.28) Ln4=F(xn+hKn2): We can now produce a 4th order accurate di®erence scheme by applying Proposition 11.1 with X(xn; h) de¯ned by (11.27). Thus we obtain the di®erence scheme. (11.29)xn+1=xn+h 6n Kn1+ 4F³ xn+h 12[Kn1+ 2kn2+ 2kn3+`n4]´ +F³ xn+h 6[Kn1+ 2Kn2+ 2Kn3+Ln4]´o ; where Knj; Ln4are as above and (11.30)kn2=F³ xn+h 4Kn1´ ; kn3=F³ xn+h 4kn2´ ; `n4=F³ xn+h 2kn2´ : 11. Numerical study { di®erence schemes 323 This formula is more complicated than the Runge-Kutta formula (11.13). We say no more about how to obtain (11.13), which represents a masterpiece of insight. We have dealt speci¯cally with autonomous systems in (11.1), but a non-autonomous system (11.31)dx dt=G(t; x); x(t0) =x0; can be treated similarly, as one can see by writing its autonomous analogue (11.32)d dtµx y¶ =µG(y; x) 1¶ ;µx(t0) y(t0)¶ =µx0 t0¶ ; and applying the formulas just derived to (11.32). We move brie°y to another class of di®erence schemes, based on power series. It derives from the expansion (11.33) x(t+h) =x(t) +hx0(t) +h2 2x00(t) +¢¢¢+hk k!x(k)(t) +O(hk+1): To begin, di®erentiate (11.1), producing (11.38) x00(t) =F2(x; x0); F 2(x; x0) =DF(x)x0: Continue di®erentiating, getting (11.35) x(j)(t) =Fj(x; x0; : : : ; x(j¡1)); j·k: Then one obtains a di®erence scheme for an approximation xntox(t0+nh), of the form (11.36) xn+1=xn+hx0 n+h2 2x00 n+¢¢¢+hk k!x(k) n; where (11.37) x0 n=F(xn); x00 n=F2(xn; x0 n); and, inductively, (11.38) x(j) n=Fj(xn; x0 n; : : : ; x(j¡1) n): This di®erence scheme is kth order accurate. In practice, this is not usually a good method, because the formulas for Fjtend to become rapidly more 324 4. Nonlinear Systems of Di®erential Equations complex. However, in some cases the functions Fjhappen not to become very complex, and then this is a good method. To mention a couple of examples, ¯rst consider the central force problem (11.39)x0=v; y0=w; v0=¡x(x2+y2)¡3=2; w0=¡y(x2+y2)¡3=2: Here, the power series method is not nearly as convenient as the Runge- Kutta method. On the other hand, for the pendulum problem, which for g=`= 1 we can write as (11.40) µ0=Ã; Ã0=¡sinµ; we have (11.41)µ00=Ã0; Ã00=¡Ãcosµ; µ(3)=Ã00; Ã(3)=¡Ã0cosµ+Ã2sinµ; µ(4)=Ã(3); Ã(4)=¡Ã00cosµ+ 3Ã0Ãsinµ+Ã3cosµ; from which one can get a workable fourth order di®erence scheme of the form (11.36){(11.38). There are other classes of di®erence schemes, such as \predictor-corrector" methods, which we will not discuss here. More about this can be found in numerical analysis texts, such as [ At] and [ Sh]. Readers with a working knowledge of a general purpose computer pro- gramming language, such as FORTRAN or C, will ¯nd it interesting to implement the Runge-Kutta method on a variety of systems of di®erential equations, including (11.39) and (11.40). Be sure to use double precision arithmetic, which makes computations to 16 digits of accuracy. Alterna- tively, specialized programming tools such as MATLAB and Mathematica can be used. These tools have built-in graphics capability, with which one can produce phase portraits, and they also have built-in di®erential equa- tion solvers, whose output one can compare with the output from one's own program. Useful literature on these latter tools for the study of di®erential equations can be found in [ PA] and [ GMP ]. When running such programs, pay attention to the way solutions behave when the step size his changed. As a rule of thumb, if the solution does not change appreciably when the step size is halved, the solution is accurate. To be sure, there is frequently more to obtaining accurate solutions than just choosing a small step size. For more on this, we recommend numerical analysis texts, such as cited above, and of course we also recommend lots of practice on various systems of di®erential equations. 11. Numerical study { di®erence schemes 325 Exercises The following exercises are for readers who can use a programming language. 1. Write a program to apply the Runge-Kutta method to the pendulum problem (11.40). 2. Write a program to apply the power series method described in (11.33){ (11.38) to (11.40). Produce a fourth order accurate method. 3. Consider applying the Runge-Kutta scheme to the problem of motion in a planar force ¯eld, (11.42) x00=f(x; y); y00=g(x; y); which can be written as the ¯rst order system (11.43)x0=v; v0=f(x; y); y0=w; w0=g(x; y): Show that (11.13){(11.14) in this context become (11.44)x7!x+h 6(v+ 2v2+ 2v3+v4); y7!y+h 6(w+ 2w2+ 2w3+w4); v7!v+h 6(a1+ 2a2+ 2a3+a4); w7!w+h 6(b1+ 2b2+ 2b3+b4); where aj; bj; vj, and wjare computed as follows. First, (11.45) a1=f(x; y); b 1=g(x; y); then (11.46)x2=x+h 2v; y 2=y+h 2w; v2=v+h 2a1; w 2=w+h 2b1; 326 4. Nonlinear Systems of Di®erential Equations and (11.47) a2=f(x2; y2); b 2=g(x2; y2); then (11.48)x3=x+h 2v2; y 3=y+h 2w2; v3=v+h 2a2; w 3=w+h 2b2; and (11.49) a3=f(x3; y3); b 3=g(x3; y3); then (11.50)x4=x+hv3; y 4=y+hw3; v4=v+ha3; w 4=w+hb3; and ¯nally, (11.51) a4=f(x4; y4); b 4=g(x4; y4): Write a program to implement this di®erence scheme. Test it for various functions f(x; y) and g(x; y). Consider particularly (11.52) f(x; y) =¡x (x2+y2)3=2; g(x; y) =¡y (x2+y2)3=2; arising in the Kepler problem, (11.39). 4. Extend the scope of Exercise 3 to treat x00=f(x; y; x0; y0); y00=g(x; y; x0; y0): 5. Write a program to apply the Runge-Kutta method to the double pen- dulum problem (9.15){(9.16). 12. Limit sets and periodic orbits LetFbe a C1vector ¯eld on an open set O ½Rn, generating the °ow ©t. Take x2 O. If ©t(x) is well de¯ned for all t¸0, we de¯ne the !-limit setL!(x) to consist of all points y2 O such that there exist tk%+1 12. Limit sets and periodic orbits 327 Figure 12.1 Figure 12.2 with ©tk(x)!y. Similarly, if ©t(x) is well de¯ned for all t·0, we de¯ne the®-limit set L®(x) to consist of all points y2 O such that there exist tk& ¡1 with ©tk(x)!y. Sinks are !-limit sets for all nearby points. Other examples of !-limit sets are pictured in Figs. 12.1{12.3. In Fig. 12.1, L!(x) is a periodic orbit, i.e., for some T2(0;1);©T(y) =y. In Fig. 12.2, L!(x) is a ¯gure eight, containing a hyperbolic critical point of the vector ¯eld. In Fig. 12.3, L!(x) contains several critical points. The following result, characterizing !-limit sets in the plane without critical points (under a few additional hypotheses), is called the Poincar¶ e- Bendixson theorem. Theorem 12.1. LetObe a planar domain, and let Fgenerate a °ow ©t onO. Assume there is a set K½ O that is a closed, bounded subset of R2 and satis¯es ©t(K)½Kfor all t >0. Take x2K. IfL!(x)contains no critical point of F, then it is a periodic orbit of ©. 328 4. Nonlinear Systems of Di®erential Equations Figure 12.3 An important ingredient in the proof of the Poincar¶ e-Bendixson theorem is the following classical result about closed curves in the plane. Jordan Curve Theorem. LetCbe a simple closed curve in R2, i.e., a continuous, one-to-one image of the unit circle. Then R2nCconsists of two connected pieces. Any curve from a point in one of these pieces to a point in the other must cross C. We will not present a proof of the Jordan curve theorem. Proofs can be found in [ GrH ],x18, and in [ Mun ]. We do mention that actually we will need this result only for piecewise smooth simple closed curves, where a simpler proof exists; see [ Sto], pp. 34{40, or [ T], Chapter 1, x19. The ability of a simple closed curve to separate Rnfails for n¸3, which makes the Poincar¶ e-Bendixson theorem an essentially two-dimensional result. Ex- amples discussed in x15 illustrate how much more complex matters can be in higher dimension. To tackle Theorem 12.1, ¯rst note that the hypotheses imply L!(x) is a nonempty subset of K. Let y2L!(x), and say (12.1) yk= ©tk(x); t k%+1; y k!y: We have F(y)6= 0. Let ¡ be a smooth curve segment in O, containing y, such that the tangent to ¡ at yis linearly independent of F(y). Shrinking ¡ if necessary, we can assume that for each z2¡, the tangent to ¡ at zis linearly independent of F(z). We say Fis transverse to ¡; cf. Fig. 12.4. With ykas in (12.1), we can assume all ykare su±ciently close to yto lie in orbits through ¡, and adjusting each tkas needed, we can take (12.2) yk2¡;8k: 12. Limit sets and periodic orbits 329 Figure 12.4 At this point, is is useful to revise the list ftkgslightly. Let t12R+; y1= ©t1(x) be as above. Now let tk%+1denote all the successive times when ©t(x) intersects ¡, so we may be adding times to the set denoted tkin (12.1). Shortly we will show that (12.1) continues to hold for this expanded set of points yk= ©tk(x). First, we make the following useful observation. Lemma 12.2. With tj< tj+1< tj+2as above, (12.3) yj+1lies between yjand yj+2on¡: Proof. Consider the curve Cjstarting at yj, running to yj+1along ©t(x); tj· t·tj+1, and returning to yjalong ¡. Cf. Fig. 12.5. This is a simple closed curve, and the Jordan curve theorem applies. Now for sand¾small and positive, and z2¡, not on the opposite side of yj+1from yj, we have ©s(yj+1) = ©tj+s(x) and ©¡¾(z) in the two di®erent connected components of R2nCj. Since f©s(yj+1) :s¸0gcannot cross Cjat any point but a point in ¡, we must have ©¡¾(yj+2) = ©tj+2¡¾(x) in the opposite component of R2nCjfrom that containing such ©¡¾(z), so yj+2must be on the opposite side of yj+1from yjin ¡. Having Lemma 12.2, we see that the expanded set of points fykg ½¡ interlaces the original set, so (12.1) continues to hold. We see that the convergence of yktoyis monotone on ¡. If by chance some yj=y, then all yk=y. Otherwise, all the points yklie on the same side of y, i.e., on the same connected component of ¡ n fyg. The main thing we need to establish to prove Theorem 12.1 is that the orbit through yis periodic. The next result takes us closer to that goal. Lemma 12.3. Suppose s >0and©s(y)2¡. Then ©s(y) =y. 330 4. Nonlinear Systems of Di®erential Equations Figure 12.5 Proof. We have (12.4) sup 0·t·s+1k©t(yk)¡©t(y)k="k!0;ask! 1 : It follows that there exist ±k!0 such that ©s+±k(yk)2¡, and hence (12.5) ©s+±k(yk) =yk+`(k);for some `(k)2 f1;2;3; : : :g: Thus (12.6) ©s(y) = lim k!1©s+±k(yk) = lim k!1yk+`(k)=y; as asserted. We are ready for the endgame in the proof of Theorem 12.1. Let sj% +1and consider zj= ©sj(y). We have each zj2K, and passing to a subsequence, we can assume (12.7) zj= ©sj(y)¡!z2K: We have F(z)6= 0, so there is a curve segment e¡ through z, transverse to F. Adjusting sj, we can arrange (12.8) zj2e¡j: We need only two such points in such a curve e¡; say, upon relabeling, (12.9) z1= ©s1(y); z2= ©s2(y) = ©s2¡s1(z1)2e¡: 12. Limit sets and periodic orbits 331 Figure 12.6 See Fig. 12.6. Note that (12.10) ©tk+s1(x)¡!z1; so we can use the previous results, with tkreplaced by tk+s1andybyz1, and ¡ by e¡. In this case, the analogue of the hypothesis in Lemma 12.3 applies: (12.11) s2¡s1>0;©s2¡s1(z1)2e¡: The conclusion of Lemma 12.3 is (12.12) ©s2¡s1(z1) =z1; i.e., actually z2=z1. Thus the orbit of © through yis periodic, of period s2¡s1. Since y2L!(x), it follows that this periodic orbit is contained in L!(x). It is also readily seen that no other point in Ocan belong to L!(x), so Theorem 12.1 is proved. The following equation, known as the van der Pol equation, illustrates the workings of Theorem 12.1. The equation is (12.13) x00¡¹(1¡x2)x0+x= 0: Here ¹is a positive parameter. This models the current in a nonlinear circuit that ampli¯es a weak current ( jxj<1) and damps a strong current (jxj>1). See the exercises for more on this. The equation (12.13) converts to the ¯rst order system (12.14) x0=y; y0=¡x+¹(1¡x2)y: 332 4. Nonlinear Systems of Di®erential Equations Figure 12.7 Fig. 12.7 is a phase portrait for the case ¹= 1. The vector ¯eld Fassociated with (12.14) has one critical point, at the origin. The linearization of (12.14) at the origin is (12.15)d dtµ» ´¶ =µ0 1 ¡1¹¶µ» ´¶ ; and the eigenvalues of this matrix are (12.16)¹ 2§1 2p ¹2¡4: Thus the origin is a source whenever ¹ >0. It is a spiral source provided also¹ <2. Note that when ( x(t); y(t)) solves (12.14), (12.17)d dt(x2+y2) = 2 ¹(1¡x2)y2; which is ¸0 forjxj ·1, and in particular is ¸0 near the origin. An examination of Fig. 12.7 indicates the presence of a periodic orbit, attracting all the other orbits. Let us see how this ¯ts into the set-up of Theorem 12.1. To do this, we need to describe a closed bounded set K½R2 such that ©t(K)½Kfor all t >0, where ©tis the °ow generated by F, and such that Fhas no critical points in K. We construct Kas follows. Look at the orbit of Fstarting at the point Aon the positive y-axis, shown in Fig. 12.7 and again in Fig. 12.8. A numerical integration of (12.14) (using the Runge-Kutta scheme) shows 12. Limit sets and periodic orbits 333 Figure 12.8 that (12.18)©t(A) winds clockwise about the origin, and again hits the positive y-axis at the point B, lying below A. To this path from AtoB, one adds the line segment (on the y-axis) from BtoA, producing a simple closed curve C. It follows readily from (12.14) that on this line segment the vector ¯eld Fpoints to the right. Thus the closed region eKbounded by this curve has the invariance property (12.19) ©t(eK)½eK;8t¸0: We then pick " >0 small enough (in particular <1), and set (12.20) K=eKn f(x; y) :x2+y2< "2g: The fact that (12.21) ©t(K)½K;8t¸0 follows from (12.19) and (12.17). We have removed the only critical point ofF, soKcontains no critical points, and Theorem 12.1 applies. It must be said that the validity of the argument just given relies on the accuracy of the statement (12.18) about the orbit through A. Here we have relied on a numerical approximation to that orbit. We applied the Runge- Kutta scheme, described in x11, with step sizes h= 10¡2;10¡3, and 10¡4, 334 4. Nonlinear Systems of Di®erential Equations using double precision (16 digit) variables, and got consistent results in all three cases. The last case involves quite a small step size, and if one were to use 8 digit arithmetic, there could be a danger of accumulating truncation errors. In any case, with today's computers there is no point in using 8 digit arithmetic. Theorem 12.1 is a special case of the following result. Bendixson's Theorem. LetFbe aC1vector ¯eld on O ½R2, generating a °ow ©t. Assume there is a set K½ O that is a closed, bounded subset of R2and satis¯es ©t(K)½Kfor all t >0. Assume Fhas at most ¯nitely many critical points in K. Then if x2K; L !(x)is one of the following: (a) a critical point, (b) a periodic orbit, (c) a cyclic graph consisting of critical points joined by orbits. A proof can be found in [ CL], Chapter 16, or in [ Lef], Chapter 10. Note that alternative (c) is illustrated in Figs. 12.2 and 12.3. We emphasize that both this result and Theorem 12.1 are results for planar vetor ¯elds. In higher dimension, matters are completely di®erent, as we will discuss in x15. We recall a device already used to deal with alternative (a), and develop it a little further. Suppose Fis aC1vector ¯eld on O ½Rn, and there is a function V2C1(O). Assume Vhas a unique minimum, at p2K. If x(t) = ©t(x0), then, by the chain rule, (12.22)d dtV(x(t)) =rV(x(t))¢F(x(t)): If also Vhas the property (12.23) rV(y)¢F(y)<0;8y2 O n p; we say Vis a strong Lyapunov function for F. In such a case, (12.24)d dtV(x(t))<0;whenever x(t)6=p: If we replace (12.23) by the weaker property (12.25) rV(y)¢F(y)·0;8y2 O; we say Vis a Lyapunov function for F. In such a case, (12.26)d dtV(x(t))·0;8t¸0: Thus, as t%+1; V(x(t)) monotonically approaches a limit, V0, which must be ¸V(p), and furthermore, (12.27) lim t!+1d dtV(x(t)) = 0 : This has the following immediate consequence. 12. Limit sets and periodic orbits 335 Proposition 12.4. LetFbe a C1vector ¯eld on O ½Rn, generating a °ow©t. Assume there is a set K½ O that is a closed, bounded subset of Rn and satis¯es ©t(K)½Kfor all t >0. Take x02K. Assume V2C1(O)is a Lyapunov function for F. Then (12.28) L!(x0)½ fy2 O:rV(y)¢F(y) = 0g: IfVis a strong Lyapunov function, then (12.29) L!(x0) =fpg: Exercises 1. Let O ½Rnbe open and ­½ O a closed bounded set with smooth boundary @­, with outward pointing normal n. Let Fbe a C1vector ¯eld on O, generating the °ow ©t. Assume (12.30) F¢n·0 on @­: Show that (12.31) ©t(­)½­;8t¸0: 2. In the setting of Exercise 1, show that (12.32) ©t(­)½©s(­) for 0 < s < t: Set (12.33) B=\ t2R+©t(­) =\ k2Z+©k(­): Show that (12.34) ©t(B) =B;8t¸0: Remark. It can be shown from material in Appendix B that Bis nonempty, closed, and bounded. 336 4. Nonlinear Systems of Di®erential Equations 3. In the setting of Exercise 2, show that (12.35) 8x2­; L !(x)½ B: 4. In the setting of Exercise 2, show that (12.36) div F <0 on ­ =)Vol(B) = 0 : 5. In the setting of Exercise 4, assume that n= 2, and that Fhas no critical points in ­, so by Theorem 12.1 there is a periodic orbit of © in ­. Show that, due to (12.36), there can be only oneperiodic orbit of © in­. Hint. Feel free to use the Jordan Curve Theorem. Exercises 6{8 deal with a nonlinear RLC circuit, as pictured in Fig. 12.9. The setup is as in x13 of Chapter 1 (see also Chapter 3, x5), except that Ohm's law is modi¯ed. The voltage drop across the \resistor" is given by (12.37) V=f(I); where fcan be nonlinear, and not necessarily monotonic. As an exam- ple, one could have (12.38) f(I) =¹³1 3I3¡I´ : Vacuum tubes and transistors can behave as such circuit elements. The voltage drop across the capacitor and the inductor are, as before, given respectively by (12.39) V=LdI dt; V =Q C: Units of current, etc., are as in x13 of Chapter 1. 6. Modify the computations done in (14.1){(14.7) of Chapter 1 and show that the current I(t) satis¯es the di®erential equation (12.40)d2I dt2+f0(I) LdI dt+1 LCI=E0(t) L: 12. Limit sets and periodic orbits 337 Figure 12.9 Show that rescaling Iandtleads to (12.14), when f(I) is given by (12.38) and E´0. More generally, rescale (12.40) to (12.41)d2x dt2+f0(x)dx dt+x=g(t): 7. Assume g´0 in (12.41). Parallel to (12.14), one can convert this equation to the ¯rst order system x0=y; y0=¡x¡f0(x)y: Show that you can also convert it to the ¯rst order system (12.43)dx dt=y¡f(x); dy dt=¡x: This is called a Lienard equation. 8. Show that if ( x(t); y(t)) solves (12.43), then (12.44)d dt(x2+y2) =¡2xf(x): 338 4. Nonlinear Systems of Di®erential Equations Figure 13.1 13. Predator-prey equations Here and in the following section we consider di®erential equations that model population densities. We start with one species. The simplest model is the exponential growth model: (13.1)dx dt=ax: Here x(t) denotes the population of the species (or rather, an approximation to what would be an integer valued function). The model simply states that the rate of growth of the population is proportional to the population itself. The solution to (13.1) is our old friend x(t) =eatx(0). This unbounded increase in population is predicated on the existence of limitless resources to nourish the species. An alternative to (13.1) posits that the resources can support a population no greater than K. The following is called the logistic equation: (13.2)dx dt=ax(1¡bx); where b= 1=K. In this model, (13.1) is a good approximation for small x, but the rate of growth slows down to 0 as xapproaches its upper limit K. The equation (13.2) can be solved by separation of variables: (13.3)dx x(1¡bx)=a dt: The reader can perform the integration as an exercise. The function F(x) =ax(1¡bx) on the right side of (13.2) is a one- dimensional vector ¯eld, with critical points at x= 0 and x= 1=b. The intervals ( ¡1;0);(0;1=b), and (1 =b;1) are all invariant under the °ow generated by F, although only the interval (0 ;1=b) has biological relevance. See Fig. 13.1 for the \phase portrait." We turn to a class of 2 £2 systems called \predator-prey" equations. For this, we set (13.4)x(t) = population of predators ; y(t) = population of prey ; z(t) = rate at which each predator eats prey : 13. Predator-prey equations 339 Figure 13.2 Depending on the choice of the exponential growth model or the logistic model for the species of prey in the absence of predators, the following systems arise to model these populations: (13.5)dx dt=¡ax+bzx; dy dt=ry¡zx; or (13.6)dx dt=¡ax+bzx; dy dt=ry(1¡cy)¡zx: Here, a; b; c , and rare positive constants. As for the rate of feeding z, we assume (13.7) z=³(y): Clearly if y= 0 then z= 0. One possibility that is used is (13.8) ³(y) =·y; for some positive constant ·. This posits that the rate of feeding of a preda- tor is proportional to the rate of close encounters of that predator with members of the other species, which in turn is proportional to the popula- tiony. This seems intuitively reasonable if yis not large, but most creatures stop eating once they are full, so a more reasonable candidate for ³(y) might be as pictured in Fig. 13.2, representing a feeding rate bounded by ¯. A class of functions of this sort is given by (13.9) ³(y) =·y 1 +°y;· °=¯: 340 4. Nonlinear Systems of Di®erential Equations Another class is (13.10) ³(y) =¯(1¡e¡°y); ¯° =·: Let us examine various cases in more detail. Volterra-Lotka equations The case (13.5) with zgiven by (13.8) produces systems called Volterra- Lotka equations: (13.11)dx dt=¡ax+¾xy; ¾ =b·; dy dt=ry¡·xy: Note that the x-axis and y-axis are invariant under the °ow de¯ned by this system. We have x0=¡axon the x-axis and y0=ryon the y-axis. It follows that the ¯rst quadrant, where x¸0 and y¸0, is invariant under the °ow. This is the region in the ( x; y)-plane of biological signi¯cance. The vector ¯eld V(x; y) = (¡ax+¾xy; ry ¡·xy)thas two critical points. One is the origin. Note that (13.12) DV(0;0) =µ¡a0 0r¶ ; so the origin is a saddle. The other critical point is (13.13) ( x0; y0) =³r ·;a ¾´ : Note that (13.14) DV(x0; y0) =µ0 ¾x0 ¡·y00¶ ; with purely imaginary eigenvalues, so we have a center for the linearization ofVat (x0; y0). In fact, ( x0; y0) is a center for V, as we now show. From (13.11) we get (13.15)dy dx=y(r¡·x) x(¾y¡a); which separates to (13.16)³ ¾¡a y´ dy=³r x¡·´ dx: 13. Predator-prey equations 341 Figure 13.3 Integrating yields (13.17) ¾y¡alogy=rlogx¡·x+C: We deduce that the following smooth function on the region x; y > 0, (13.18) H(x; y) =¾y¡alogy+·x¡rlogx; is constant on orbits of (13.11), i.e., these orbits lie on level curves of H. Note that (13.19) rH(x; y) =µ·¡r x ¾¡a y¶ ; D2H(x; y) =µr x20 0a y2¶ ; hence, with ( x0; y0) as in (13.13), (13.20) rH(x0; y0) = 0 ; D2H(x0; y0) =µr x2 00 0a y2 0¶ ; the latter matrix being positive de¯nite, so Hhas a minimum at ( x0; y0), which implies that ( x0; y0) is a center for V. The phase portrait for orbits of (13.11) is pictured in Fig. 13.3. The system (13.11) was studied independently by Lotka and Volterra around 1925, by Lotka as a model of some chemical reactions and by Volterra as a predator-prey model, speci¯cally for sharks preying on another species of ¯sh. Volterra made the following further observation. Bring in another type of predator, ¯shermen. Assume the ¯shermen keep everything they 342 4. Nonlinear Systems of Di®erential Equations catch and that the probability of getting caught in their nets is the same for sharks and their prey. Then the system (13.11) gets revised to (13.21)dx dt=¡ax+¾xy¡ex; dy dt=ry¡·xy¡ey: Now (13.21) has the same form as (13.11), with areplaced by a+eand with rreplaced by r¡e, all these constants remaining positive as long as (13.22) 0 < e < r: Then the previous analysis applies. The system (13.21) has a stable critical point at (13.23) ( x1; y1) =³r¡e ·;a+e ¾´ : Note that at this critical point there are fewer sharks and more prey, com- pared to (13.13). Of course, this depends on the hypothesis (13.22). If e > r , things are catastrophically di®erent. First modi¯cation We turn from Volterra-Lotka equations to predator-prey models given by (13.6), still keeping (13.8). Then we have the following system: (13.24)dx dt=¡ax+¾xy; ¾ =b·; dy dt=ry(1¡cy)¡·xy: As with (13.11), the x-axis and y-axis are invariant under the °ow de¯ned by this system. We have x0=¡axon the x-axis and y0=ry(1¡cy) on the y-axis. Again, the ¯rst quadrant ( x¸0; y¸0) is invariant under the °ow. Note furthermore that, for (13.25) V(x; y) = (¡ax+¾xy; ry (1¡cy)¡·xy)t; we have (13.26) V³ x;1 c´ =³³¾ c¡a´ x;¡· cx´t ; which points downward for x >0. It follows that (13.27) R=n (x; y) :x¸0;0·y·1 co 13. Predator-prey equations 343 is invariant under this °ow. It is this region in the ( x; y)-plane that is of biological signi¯cance. To proceed, we ¯nd the critical points of V(x; y), given by (13.25). Two of these are (13.28) (0 ;0) and³ 0;1 c´ : DV(0;0) is again given by (13.12), so (0 ;0) is a saddle. Also, (13.29) DV³ 0;1 c´ =µ¡a+¾ c0 ¡· c¡r¶ : Vhas a third critical point, at (13.30) y0=a ¾; x 0=r ·³ 1¡ca ¾´ =rc ·¾³¾ c¡a´ : Note how this point is shifted to the left from the point (13.13). There are three cases to consider. Case I. ¾=c¡a <0: In this case, the critical point (13.30) is not in the ¯rst quadrant, so V has only the critical points (13.28) in R. In this case (13.29) has two neg- ative eigenvalues, so the critical point (0 ;1=c) is a sink. Note that the x-component of V(x; y) is (13.31) x(¾y¡a)·x³¾ c¡a´ ;forx¸0; y·1 c; soVpoints to the left everywhere in Rexcept the left edge. Consequently, the population of predators is driven to extinction as t!+1, whatever the initial condition. Case II. ¾=c¡a >0. In this case the third critical point ( x0; y0) is in the ¯rst quadrant. In fact, y0=a=¾ < 1=c, so ( x0; y0)2 R. Now (13.29) has one positive and one negative eigenvalue, so the critical point (0 ;1=c) is a saddle. As for the nature of ( x0; y0), we have (13.32)DV(x0; y0) =µ¡a+¾y0 ¾x0 ¡·y0 r(1¡2cy0)¡·x0¶ =µ0rc ·(¾ c¡a) ¡·a ¾¡rca ¾¶ : 344 4. Nonlinear Systems of Di®erential Equations Figure 13.4 Note that (13.33)detDV(x0; y0) =rca ¾³¾ c¡a´ >0; TrDV(x0; y0) =¡rca ¾<0: It follows that the eigenvalues of DV(x0; y0) are either both negative or have negative real part. Hence ( x0; y0) is a sink. We claim that the orbit through each point in Rnot on the xory- axis approaches ( x0; y0) ast!+1. To see this, we construct a Liapunov function. We do this by modifying H(x; y) in (13.18), which has a minimum at the point (13.13), to one that has a minimum at the point (13.30). We take (13.34) eH(x; y) =¾y¡alogy+·x¡r³ 1¡ca ¾´ logx: If (x(t); y(t)) solves (13.24), a computation gives (13.35)d dteH(x; y) =¡rc ¾(¾y¡a)2: By Proposition 12.4, if we take any point p2 R, with positive xandy- coordinates (so it is in the domain of eH), the !-limit set of psatis¯es (13.36) L!(p)½n (x; y)2 R:y=a ¾o : The right side is a horizontal line to which Vis clearly transverse except at the critical point ( x0; y0), so indeed L!(p) = (x0; y0). See Fig. 13.4 for a phase portrait treating Case II. 13. Predator-prey equations 345 Case III. ¾=c¡a= 0. In this case ( x0; y0) = (0 ;1=c). In (13.29) the eigenvalues are 0 and ¡r, so (0;1=c) is a degenerate critical point. In place of (13.31) we have that the x-component of V(x; y) is (13.37) x(¾y¡a)·0;forx¸0; y·1 c; and it is strictly negative for x > 0; y < 1=c. Hence, as in Case I, the population of predators is driven to extinction as t!+1. Second modi¯cation We now move to the next level of sophistication, using the system (13.6) with z=³(y), described as in Fig. 13.2. Thus, we look at systems of the form (13.37)dx dt=¡ax+bx³(y); dy dt=ry(1¡cy)¡x³(y): As before, a; b; c , and rare all positive constants. To be precise about what we mean when we say ³(y) behaves as in Fig. 13.2, we make the following hypotheses: (13.38)(a) ³: [0;1)![0;1) is smooth ; (b) ³(0) = 0 ; (c) ³0(y)>0;8y¸0; (d) sup ³(y) =¯ <1; (e) ³00(y)·0: All these conditions are satis¯ed by the examples (13.9) and (13.10). Hy- pothesis (c) implies ³is strictly monotone increasing, and hypothesis (e) implies ³is concave. In this case, the vector ¯eld is (13.39) V(x; y) = (x(b³(y)¡a); ry(1¡cy)¡x³(y))t: Parallel to (13.26), (13.40) V³ x;1 c´ =³³ b³³1 c´ ¡a´ x;¡³³1 c´ x´t ; 346 4. Nonlinear Systems of Di®erential Equations which points downward for x >0, and again it follows that the region R, given by (13.27), is invariant under the °ow ©tgenerated by V, for t¸0, and this is the region in the ( x; y)-plane that is of biological signi¯cance. Next, we ¯nd the critical points of V(x; y). Again, two of them are (0;0) and³ 0;1 c´ ; and again DV(0;0) is given by (13.12), so (0 ;0) is a saddle. This time, (13.41) DV³ 0;1 c´ =µb³(1 c)¡a0 ¡³(1 c)¡r¶ : Also, a critical point would occur at ( x0; y0) if these coordinates satisfy (13.42) ³(y0) =a b; x 0=b ary0(1¡cy0): Under the hypotheses (13.38), the ¯rst equation in (13.42) has a (unique) solution if and only if (13.43)a b< ¯: From here on we will assume (13.43) holds, and leave it to the reader to consider the behavior of the °ow when (13.43) fails. Given (13.43), x0and y0are well de¯ned by (13.42). Parallel to the study of (13.30), again we have three cases. Case I. 1¡cy0<0, Case II. 1¡cy0>0, Case III. 1¡cy0= 0. In Case I, ( x0; y0) is not in the ¯rst quadrant, and in Case III, ( x0; y0) = (0;1=c). Again we leave these cases to the reader to think about. We concentrate on Case II. In Case II, x0>0 and 0 < y0<1=c, so (13.44) ( x0; y0)2 R: Given ³(y0) =a=band the hypotheses (13.38) on ³, we have (13.45) ³³1 c´ >a b()1 c> y0()1¡cy0>0; 13. Predator-prey equations 347 and hence in Case II, DV(0;1=c) has one positive eigenvalue and one nega- tive eigenvalue, so (13.46)³ 0;1 c´ is a saddle. (In Case I, the eigenvalues of DV(0;1=c) are both negative, so (0 ;1=c) is a sink, and in Case III these eigenvalues are 0 and ¡r.) Next, a computation gives the following analogue of (13.32): (13.47)DV(x0; y0) =µb³(y0)¡a b³0(y0)x0 ¡³(y0)r(1¡2cy0)¡x0³0(y0)¶ =µ0 b³0(y0)x0 ¡a br(1¡2cy0)¡x0³0(y0)¶ ; and parallel to (13.33) we have (13.48)detDV(x0; y0) =ax0³0(y0)>0; TrDV(x0; y0) =r(1¡2cy0)¡x0³0(y0) =rh ¡cy0+ (1¡cy0)n 1¡³0(y0)y0 ³(y0)oi : Let us set (13.49) Z0= 1¡³0(y0)y0 ³(y0): Given ³, this is a function of a=b, but it is independent of candr. Note that, since ³(0) = 0, (13.50)³(y0) y0=³0(~y);for some ~ y2(0; y0); by the mean value theorem, so the hypotheses on ³in (13.38) imply (13.51) 0 < Z 0<1: (Note that in the context of the previous model, with ³(y) given by (13.8), Z0= 0.) We have (13.52) Tr DV(x0; y0) =r£ Z0(1¡cy0)¡cy0¤ : This gives rise to three cases. Case IIA. Z0< cy 0=(1¡cy0). Then Tr DV(x0; y0)<0, so, by (13.48), (13.53) ( x0; y0) is a sink. 348 4. Nonlinear Systems of Di®erential Equations Figure 13.5 Case IIB. Z0> cy 0=(1¡cy0). Then Tr DV(x0; y0)>0, so, by (13.48), (13.54) ( x0; y0) is a source. Case IIC. Z0=cy0=(1¡cy0). Then Tr DV(x0; y0) = 0, so, by (13.48), the eigenvalues of DV(x0; y0) are (nonzero) purely imaginary numbers. In this case, ( x0; y0) is a center for the linearization of V. We will concentrate on Cases IIA and IIB. Before pursuing these cases further, we want to describe a family of bounded domains in Rthat are invariant under the °ow ©tfort¸0. Namely, consider the triangle T¹with vertices at (0 ;1=c);(0;0), and ( ¹;0), as pictured in Fig. 13.5. Claim. If¹ >0 is large enough, the triangle T¹is invariant under ©t, for t¸0. Proof. Note that Vis vertical on the left edge of T¹, with critical points at the endpoints of this line segment. Also Vpoints horizontally to the left on the bottom edge of T¹. It remains to show that Vpoints into T¹along the line segment from (0 ;1=c) to ( ¹;0), provided ¹is su±ciently large. This line segment is given by (13.55) x=¹(1¡cy);0·y·1 c; and the vector (13.56) N¹=µ1 ¹c¶ 13. Predator-prey equations 349 is normal to this segment, and points away from T¹. We want to show that V¢N¹·0 along this line segment, for ¹large. Indeed, from (13.39), (13.57)V(¹(1¡cy); y)¢N¹= (1¡cy)£ ¹(b³(y)¡a) +¹cry¡¹2c³(y)¤ =¹(1¡cy)£ ¡a+cry¡(¹c¡b)³(y)¤ ; and under the hypotheses (13.38) on ³, this is (13.58) ·0;8y2h 0;1 ci ; if¹is su±ciently large, say ¹¸¹0. A similar computation shows that, if ¹1> ¹ 0, then, for each p2 R, ©t(p)2 T¹1for all su±ciently large t. Back to Cases IIA and IIB, as we have seen, in Case IIA ( x0; y0) is a sink. It is possible to show that (13.59) in Case IIA, ©t(p)¡!(x0; y0);ast!+1; for all pin the interior of R, so the phase portrait has qualitative features similar to Fig. 13.4. On the other hand, in Case IIB, ( x0; y0) is a source. Hence there is an open set Ucontaining ( x0; y0) such that (13.60) T¹0nUis invariant under ©t, fort¸0: This region does contain the two critical points (0 ;0) and (0 ;1=c), on its boundary, but since they are saddles, the argument used to establish the Poincar¶ e-Bendixson theorem, Theorem 12.1, shows that (13.61) in Case IIB, L!(p) is a periodic orbit, for all p6= (x0; y0) in the interior of R. The phase portrait is depicted in Fig. 13.6. Exercises Exercises 1{5 deal with the system (13.37), i.e., (13.62)x0=¡ax+bx³(y); y0=ry(1¡cy)¡x³(y); 350 4. Nonlinear Systems of Di®erential Equations Figure 13.6 where ³(y) is given by (13.9), i.e., (13.63) ³(y) =·y 1 +°y;· °=¯: As usual, a; b; c; ·; °; r 2(0;1). The exercises deal with when Cases I{III, speci¯ed below (13.43), hold. Recall these cases apply if and only if there is a critical point ( x0; y0) given by (13.42), i.e., of and only if (13.64)a b< ¯=· °: We will assume this holds. 1. Show that the critical point ( x0; y0) is given by (13.65) y0=a b·¡a°; x 0=b ary0(1¡cy0): 2. Show thatCase I ()ac > b· ¡a°; Case II ()ac < b· ¡a°; Case III ()ac=b·¡a°: 3. Let Z0be given by (13.49), i.e., (13.66) Z0= 1¡³0(y0)y0 ³(y0): 13. Predator-prey equations 351 Show that (13.67) Z0=a° b·: 4. In Case II, recall Cases IIA{IIC, speci¯ed below (13.52). Show that Case IIA ()° ·<c b·¡a°¡ac; Case IIB ()° ·>c b·¡a°¡ac; Case IIC ()° ·=c b·¡a°¡ac; 5. Let us take (13.68) a= 1; b = 2; · = 1; ° = 1: Note that (13.64) holds. Show that Case I ()c >1; Case II ()c <1; Case III ()c= 1: In Case II, show that Case IIA ()c >1 3; Case IIB ()c <1 3; Case IIC ()c=1 3: Exercises 6{10 deal with the system (13.62), where ³(y) is given by (13.10), i.e., (13.69) ³(y) =¯(1¡e¡°y); ¯° =·: Again there is a critical point ( x0; y0), given by (13.42), if and only if (13.64) holds. We assume this holds, so b¯ > a . 352 4. Nonlinear Systems of Di®erential Equations 6. Show that the critical point ( x0; y0) is given by (13.70) y0=1 °logb¯ b¯¡a; x 0=b ary0(1¡cy0): 7. For Z0, de¯ned by (13.66), show that (13.71) Z0= 1¡b¯¡a alogb¯ b¯¡a: 8. Parallel to Exercise 2, study when Cases I{III hold. 9. Parallel to Exercise 4, study when Cases IIA{IIC hold. 10. Take a; b; ·; , and °as in (13.68). Work out a parallel to Exercise 5. For Exercises 11{12, consider the following system, for xpredators and yprey, presented in [ Tau], p. 376: (13.72)x0=ax³ b¡x y´ ; y0=ry(1¡cy)¡x³(y): Here the equation for yis as in (13.62), modeling the population of prey in terms of the logistic equation, modi¯ed by how fast the prey is eaten. The equation for xhas a di®erent basis, a sort of logistic equation in which the population ydetermines the population limit of x, at any given time. 11. Work out an analysis of the system (13.72) as parallel as possible to the analysis done in this section for (13.62). 12. Take ³(y) as in (13.63) and work out results parallel to those of Exercises 1{5. Exercises 13{15 are for readers who can use a programming language, with graphics capabilities. 13. Predator-prey equations 353 13. The following system is known as the basic model of virus dynamics (cf. [NM], p. 100, [ W], p. 26): (13.73)dx dt=¸¡dx¡¯xv; dy dt=¯xv¡ay; dv dt=ky¡uv: Here, xrepresents the uninfected cell population, ythe infected cell pop- ulation, and vthe virus population. The positive parameters ¸; d; ¯; a; k , anduare taken to be constant. The ratio (13.74) R0=¸¯k adu is called the basic reproducive ratio. Graph solution curves for (13.73), with various choices of parameters. Account for the assertion that if R0<1 the virus cannot maintain an infection, but if R0>1 the system converges to an equilibrium, in which v >0. 14. The simplifying assumption that the virus population is proportional to the infected cell population (say ¯v=by) leads to the system (13.75)dx dt=¸¡dx¡bxy; dy dt=¡ay+bxy: Study this system, with an eye to comparison with the Volterra-Lotka system (13.11). Here, replace (13.74) by (13.76) R0=b¸ ad: 15. The following system modi¯es (13.75) by introducing z(t), the popula- tion of \killer T cells," which kill o® infected cells, thereby negatively a®ecting y: (13.77)dx dt=¸¡dx¡bxy; dy dt=bxy¡ay¡pyz; dz dt=cyz¡bz; 354 4. Nonlinear Systems of Di®erential Equations now with positive parameters ¸; d; b; a; p , and c. Continue to de¯ne R0 by (13.76). Consider particularly cases where (13.78) R0>1; c³¸ a¡d b´ > b: Account for the assertion that in this case the virus population ¯rst grows, stimulating the production of killer T cells, which in turn ¯ght the infection and lead to an equilibrium. For more on these models, see [ NM] and [ W], and references therein. 14. Competing species equations The following system models the populations x(t) and y(t) of two competing species: (14.1)dx dt=ax(1¡bx)¡cxy; dy dt=®y(1¡¯y)¡°xy: In this model, each population is governed by a logistic equation in the absence of the other species. The presence of the other species reduces the population of its opponent, at a rate proportional to xy. Setting X=bxand Y=¯yproduces an equation like (14.1), but with X(1¡X) and Y(1¡Y) in place of x(1¡bx) and y(1¡¯y), and with di®erent factors. A change of notation gives the system (14.2)dx dt=ax(1¡x)¡cxy; dy dt=®y(1¡y)¡°xy: which we will consider henceforth. We take a; c; ®; ° 2(0;1). Associated to this system is the vector ¯eld (14.3) V=µax(1¡x)¡cxy ®y(1¡y)¡°xy¶ : Note that V(x;0) = ( ax(1¡x);0)tandV(0; y) = (0 ; ®y(1¡y))t, so the x-axis and y-axis are invariant under the °ow ©tgenerated by V. Hence the quadrant fx¸0; y¸0g, which is the region of biological signi¯cance, is invariant under ©t. Note also that (14.4) V(x;1) =µax(1¡x)¡cx ¡°x¶ ; V (1; y) =µ¡cy ®y(1¡y)¡°y¶ ; 14. Competing species equations 355 so ©tleaves invariant the region (14.5) B=f(x; y) : 0·x; y·1g; fort¸0. The vector ¯eld Vhas the following critical points, (14.6) (0 ;0);(0;1);(1;0); and a fourth critical point ( x0; y0), satisfying (14.7) cy0=a(1¡x0); °x 0=®(1¡y0): A calculation gives (14.8) x0=®a¡c a®¡c°; y 0=a®¡° a®¡c°: The point ( x0; y0) may or may not lie in the ¯rst quadrant. We investigate this further below. We have (14.9) DV(0;0) =µa0 0®¶ ; so (0;0) is a source. Also, (14.10) DV(0;1) =µa¡c0 ¡°¡®¶ ; DV (1;0) =µ¡a¡c 0®¡°¶ ; and each of these might be a saddle or a sink, depending on the signs of a¡cand®¡°. Next, (14.11)DV(x0; y0) =µa(1¡2x0)¡cy0 ¡cx0 ¡°y0 ®(1¡2y0)¡°x0¶ =µ¡ax0¡cx0 ¡°y0¡®y0¶ ; the second identity by (14.7). Hence (14.12)detDV(x0; y0) = (a®¡c°)x0y0; TrDV(x0; y0) =¡ax0¡®y0: At this point, it is natural to consider the following cases: 356 4. Nonlinear Systems of Di®erential Equations Figure 14.1 Case I. a > c and® > ° . Case II. a < c and® < ° . Case III. a > c and® < ° . Case IV. a < c and® > ° . In Case I, we see from (14.10) that (14.13) (0 ;1) and (1 ;0) are saddles. In this case, a® > c° , so, by (14.8), (14.14) x0>0; y 0>0; and the critical point ( x0; y0) is in the ¯rst quadrant. Then we see from (14.12) that (14.15) det DV(x0; y0)>0;TrDV(x0; y0)<0; so (14.16) ( x0; y0) is a sink. We have (14.17) ©t(x; y)¡!(x0; y0) as t!+1; whenever x >0 and y >0. The two competing species tend to an equilib- rium of coexistence. The phase portrait for this case, with a= 2; ®= 2; c= 1; °= 1, is illustrated in Fig. 14.1. 14. Competing species equations 357 Figure 14.2 In Case II, we see from (14.10) that (14.18) (0 ;1) and (1 ;0) are sinks. In this case, a® < c° , so, by (14.8), again (14.14) holds, and the critical point ( x0; y0) is in the ¯rst quadrant. We see from (14.12) that (14.19) det DV(x0; y0)<0; so (14.20) ( x0; y0) is a saddle. The phase portrait for this case, with a= 1; ®= 1; c= 2; °= 2, is illustrated in Fig. 14.2. For almost all initial data ( x; y) in the ¯rst quadrant, ©t(x; y) tends to either (0 ;1) or (1 ;0) as t!+1. One species or the other tends toward extinction, depending on the initial conditions. In Case III, we see from (14.10) that (14.21) (0 ;1) is a saddle and (1 ;0) is a sink. From here two sub-cases arise, depending on the relative size of a®andc°. Case IIIA. a® > c° . This time, by (14.8), (14.22) x0>0; y 0<0; 358 4. Nonlinear Systems of Di®erential Equations Figure 14.3 so the critical point ( x0; y0) is not in the ¯rst quadrant. We see from (14.12) that (14.23) det DV(x0; y0)<0; so (14.24) ( x0; y0) is a saddle. The phase portrait for this case, with a= 2; ®= 1; c= 1=4; °= 2, is illustrated in Fig. 14.3. We have (14.25) ©t(x; y)¡!(1;0) as t!+1; whenever x >0 and y >0. Species ytends to extinction. Case IIIB. a® < c° . This time, by (14.8), (14.26) x0<0; y 0>0; and again the critical point ( x0; y0) is not in the ¯rst quadrant. We see from (14.22) that (14.27) det DV(x0; y0)>0: Thus (14.28) ( x0; y0) is a source or a sink, 14. Competing species equations 359 Figure 14.4 depending on the sign of Tr DV(x0; y0). The phase portrait for this case, with a= 2; ®= 1=2; c= 1; °= 2, is illustrated in Fig. 14.4. (In this example, ( x0; y0) is a sink.) Again (14.25) holds whenever x >0 and y >0. To summarize Case III, the °ows in the ¯rst quadrant have the same qualitative features in the two sub-cases; (14.25) holds. The features di®er outside the ¯rst quadrant. As for Case IV, this reduces to Case III by switching the roles of xand y. Exercises 1. Note that if xandysolve (14.2), then d dt(x+y) =¡ax2¡®y2¡(c+°)xy+ax+®y: Show that there exists R2(0;1) such that x; y¸0; x2+y2¸R2=)d dt(x+y)·0: Deduce global existence of solutions to (14.2), for t¸0, given ( x(0); y(0)) in the ¯rst quadrant. 360 4. Nonlinear Systems of Di®erential Equations 2. In the setting of Exercise 1, show that whenever x(0)>0 and y(0)>0, we have ( x(t); y(t))2 B, given by (14.5), for t >0 su±ciently large. 3. Consider the system dx dt=x(1¡x)¡xy; dy dt=y(1¡y)¡°xy; with °2(0;1). Specify when Cases I{IV hold. Record the possible outcomes, as regards coexistence/extinction. 4. Consider the system dx dt=1 2x(1¡x)¡cxy; dy dt=y(1¡y)¡2xy; with c2(0;1). Specify when Cases I{IV hold. Record the possible outcomes, as regards coexistence/extinction. 5. Consider the system dx dt=ax(1¡x)¡xy; dy dt= 2y(1¡y)¡xy; with a2(0;1). Specify when Cases I{IV hold. Record the possible outcomes, as regards coexistence/extinction. 15. Chaos in multidimensional systems As previewed in the introduction to this chapter, two phenomena conspire to limit the complexity of °ows generated by autonomous planar vector ¯elds. One is that orbits cannot cross each other, due to uniqueness (this holds in any number of dimensions). The other is that a directed curve (with nonzero velocity) in the plane divides a neighborhood of each of its points into two parts, the left and the right. This latter fact played an important role in x12. In dimension 3 and higher, this breaks down completely, and allows for far more complex °ows. Newtonian motion in a force ¯eld in the plane is described by a second order 2 £2 system of di®erential equations, which is converted to a 4 £4 ¯rst 15. Chaos in multidimensional systems 361 order system. Energy conservation con¯nes the motion to a 3-dimensional constant energy surface. If the force is a central force, there is also conser- vation of angular momentum. These two conservation laws make for regular motion, as seen in xx5{6. These are \integrable" systems. Such integrability is special. Most systems from physics and other sources do not possess it. For example, the double pendulum equation, derived in x9, does not have this property. (We do not prove this here.) Flows generated by vector ¯elds on n-dimensional domains with n¸3 are thus sometimes regular, but often they lack regularity to such a degree that they are deemed \chaotic." Signatures of chaos include the inability to predict the long time behavior of orbits. This inability arises not only from the lack of a formula for the solution in terms of elementary functions. In addition, numerical approximations to the orbits of these °ows reveal a \sen- sitive dependence" on initial conditions and other parameters. Furthermore, phase portraits of these orbits lookcomplex. Research into these chaotic °ows takes the study of di®erential equations to the next level, beyond this introduction. We end this chapter with a discussion of two special cases of 3 £3 systems, to give a °avor of the complexities that lie beyond, and we provide pointers to literature that addresses the deep questions raised by e®orts to understand such systems. Lorenz equations The ¯rst example is the following system, produced by E. Lorenz in 1963 to model some aspects of °uid turbulence: (15.1)x0=¾(y¡x); y0=rx¡y¡xz; z0=xy¡bz: An alternative presentation is (15.2)d dt0 @x y z1 A=0 @¡¾ ¾ 0 r¡1 0 0 0 ¡b1 A0 @x y z1 A+0 @0 0 0 0 0 ¡x 0x01 A0 @x y z1 A: Denoting the right side of (15.2) by V(x; y; z ), we see that the ¯rst matrix on the right side is DV(0;0;0). One assumes the parameters ¾; b, and rare all positive. Lorenz took (15.3) ¾= 10; b =8 3; and considered various values of r, with emphasis on (15.4) r= 28: 362 4. Nonlinear Systems of Di®erential Equations Phase portraits of some orbits for (15.1), with ¾andbgiven by (15.3) and with various values of rare given in Fig. 15.1. Each of the six portraits depicts the forward orbits through the three points (15.5) x=k 100; y = 0; z = 5; k =¡1;0;1: The portraits start out simple, execute a sequence of changes, as rincreases, reaching substantial apparent complexity at r= 28. We discuss some as- pects of this. First, some global results. Global forward solvability of (15.1) can be established with the help of the remarkable function (15.6) f(x; y; z ) =rx2+¾y2+¾(z¡2r)2: A calculation shows that if ( x(t); y(t); z(t)) solves (15.1), then (15.7)d dtf(x; y; z ) =¡2¾(rx2+y2+bz2¡2brz): Clearly there exists K2(0;1) such that (15.8) B=f(x; y; z )2R3:f(x; y; z )·Kg is a closed, bounded subset of R3and the right side of (15.7) is <0 on the complement of B. Hence (15.9) ©t(B)½B;8t >0; where ©tis the °ow generated by V(x; y; z ). Moreover, for each ( x; y; z )2 R3, (15.10) ©t(x; y; z )2B; for all su±ciently large t >0: Note that (15.9) plus the identity ©t= ©s±©t¡simplies (15.11) ©t(B)½©s(B) for 0 < s < t; soB(t) = ©t(B) is a family of closed, bounded sets that is decreasing as t%+1. Now set (15.12) B=\ t2R+B(t) =\ k2Z+B(k): 15. Chaos in multidimensional systems 363 Figure 15.1 364 4. Nonlinear Systems of Di®erential Equations The set Bis called the attractor for (15.1). We have (15.13) ©t(B) =B;8t¸0: Note that (15.14) div V=¡¾¡1¡b <0; so results of x3 imply (15.15) Vol B= 0: This attractor has a simple description for small r, but becomes very com- plex for larger r. To proceed with the analysis, consider the critical points. The origin is a critical point of Vfor all ¾; b; r 2(0;1). Since DV(0) is the ¯rst matrix on the right side of (15.2), we see its eigenvalues are (15.16) ¸§=¡¾+ 1 2§1 2p (¾+ 1)2+ 4¾(r¡1); ¸ 3=¡b; with eigenvectors (15.17) v§=0 @¾ ¸§+¾ 01 A; v 3=0 @0 0 11 A: It follows from (15.16) that (15.18) 0< r < 1 =)DV(0) has 3 negative eigenvalues ; r >1 =)DV(0) has 2 negative and one positive eigenvalue. Forr > 1, the positive eigenvalue is ¸+and its associated eigenvector isv+. There is a parallel to the results in (3.31) describing saddles. It is shown in [ Hart ] that there is a smooth 2-dimensional surface through the origin consisting of points psuch that ©t(p)!0 as t!+1and a smooth 1-dimensional curve through the origin consisting of points psuch that ©t(p)!0 ast! ¡1 . In general, a smooth k-dimensional surface in Rnis called a k-dimensional manifold. The sets described above are called a \stable manifold" and an \unstable manifold," respectively. See also Ap- pendic C for further discussion. Forr >1,Vhas two additional critical points, satisfying (15.19) x=y;(r¡1¡z)x= 0; bz =x2; 15. Chaos in multidimensional systems 365 i.e., (15.20) C§= (§p b(r¡1);§p b(r¡1); r¡1): We have (15.21) DV(C§) =0 @¡¾ ¾ 0 1¡1§» §»§»¡b1 A; » =p b(r¡1): Note that DV(C+) and DV(C¡) are conjugate by the action of (15.22)0 @¡1 ¡1 11 A; so they have the same eigenvalues. This mirrors the fact that (15.1) is invari- ant under the transformation ( x; y; z )7!(¡x;¡y; z). Further calculations give the following results, when ¾andbare given by (15.3): (15.23) DV(C§) has 3 negative eigenvalues for 1 < r < 1:346¢¢¢ 1 negative and 2 with negative real part for 1 :346¢¢¢< r < 24:74¢¢¢ 1 negative and 2 with positive real part for r >24:74¢¢¢: In the ¯rst two cases in (15.23), Proposition 3.4 applies, and for all points p su±ciently close to C+;©t(p)!C+ast!+1, and similarly for C¡. The third case in (15.23) is like the second case in (15.18), except the numbers are reversed. In such a case, there are a 2-dimensional unstable manifold and a 1-dimensional stable manifold through C+, and similarly for C¡, in the language introduced below (15.18). With these calculations in hand, let's take a closer look at the six phase portraits depicted in Fig. 15.1, orbits with initial data given by (15.5). In all cases there is a vertical line segment from (0 ;0;5) to (0 ;0;0), and we see from (15.1) that the z-axis is invariant under the °ow for all values of the parameters. Furthermore, on the z-axis, z0=¡bz. Now the initial points ( §0:01;0;5) are close by, but for all r-values depicted, DV(0) has one positive eigenvalue, and the orbits push away from the origin, in a direction close to §v+, where v+is given by (15.17). The orbit from (+0 :01;0;5) spirals into the critical point C+, and the orbit from ( ¡0:01;0;5) spirals intoC¡, in the ¯rst two portraits, where r= 4:667 and 9 :333. Around r¼14, something new happens. These orbits pass close to the origin. 366 4. Nonlinear Systems of Di®erential Equations Fig. 15.2 shows such a transition in more detail. Here we have two orbits with slightly di®erent initial conditions, namely (15.24) p§=§"v+; with "chosen small, to capture the unstable manifold of the origin more accurately. At a certain critical value rh¼13:926, the unstable manifold is actually a pair of homoclinic orbits, approaching the origin both as t! ¡1 and as t!+1. For larger values of r, the orbit from p+(and that from (+0:01;0;5)) crosses over and spirals into C¡, while the orbit from p¡(and that from ( ¡0:01;0;5)) spirals into C+, as depicted in the fourth portrait in Fig. 15.1. This spiraling into C§does not endure as rincreases. As stated in (15.23), there is a critical rc¼24:74 past which DV(C§) has two eigenvalues with positive rather than negative real part. At r= 23:333, this spiraling has slowed. In fact, the ¯fth portrait in Fig. 15.1 does not reveal spiraling all the way in. The orbits pictured there are of the form ©t(pj) for t2[0;40]. Iftis taken somewhat larger, one sees asymptotic approaches to C§, with r= 23:333. In the sixth phase portrait of Fig. 15.1, we have r= 28 > rc. The or- bits approach the unstable manifolds of C§and then spiral out from these critical points. After some spiraling out from C¡, the orbit starting from (+0:01;0;5) makes a jump to the vicinity of C+, approaches its unstable manifold, and starts spiraling out from C+. After a while, the orbit jumps back to the vicinity of C¡, and this spiraling and jumping is endlessly re- peated. The six phase portraits in Fig. 15.3 show (15.25) ©t(0:01;0;5);80j < t < 80(j+ 1);0·j·5: The portraits di®er in ¯ne detail from each other, but they are fairly similar, and seem to reveal what is called a strange attractor. Figures 15.1{3 were produced by numerically integrating (15.1), using a fourth-order Runge-Kutta scheme, described in x11, with step size (15.26) h= 0:0005: Use of the step size h= 0:001 produced apparently identical phase portraits in Fig. 15.2, and in all but the last portrait in Fig. 15.1. There were notice- able di®erences in the last phase portrait of Fig. 15.1 and in the portraits of Fig. 15.3. This phenomenon gives evidence of sensitive dependence of the orbits on initial conditions, and leads to unpredictibility of orbits, which is part of the signature of chaos. 15. Chaos in multidimensional systems 367 Figure 15.2 368 4. Nonlinear Systems of Di®erential Equations Figure 15.3 15. Chaos in multidimensional systems 369 We make one further comment about Figs. 15.1{15.3. Of course, the orbits depicted are curves ( x(t); y(t); z(t)) inR3. What is shown in these ¯gures are 2-dimensional projections, namely ( u(t); v(t)), with u(t) =x(t)+ y(t)=2; v(t) =z(t)¡y(t)=2. Periodically forced Du±ng equation Our second example arises from motion in 1 dimension, in a nonlinear background ¯eld, with a periodic forcing term added: (15.27)d2x dt2=f(x) +rcost: Here ris a parameter. When converted to a ¯rst order system and put in autonomous form, this becomes (15.28)dx dt=y; dy dt=f(x) +rcosz; dz dt= 1: We take (15.29) f(x) =x¡x3: Then (15.27) is called a periodically forced Du±ng equation if r6= 0. For r= 0 it is called Du±ng's equation, and it reduces to a 2 £2 system, whose phase portrait is given in Fig. 15.4. There are two homoclinic orbits, each tending to the origin as t! §1 . All the other orbits are closed, and lie on level curves of (15.30) E(x; y) =y2 2¡x2 2+x4 4: Fig. 15.5 shows six individual orbits of (15.28) with r= 0, orbits through the six points (15.31) x=p 2 +3k 10;¡4·k·1; y = 0: The orbits for (15.28) are curves ( x(t); y(t); z(t)), but for r= 0 we simply plot ( x(t); y(t)) in Fig. 15.5. Forr6= 0, matters are more complicated, since zis coupled to ( x; y) in (15.28). We need a di®erent way to portray the orbits ( x(t); y(t); z(t)). In 370 4. Nonlinear Systems of Di®erential Equations Figure 15.4 this case, unlike for the Lorenz system, a linear projection of ( x; y; z ) space onto ( u; v) space is not the best way to proceed. Taking into account the periodicity of the right side of (15.28) in z, we treat z=tas an angular variable, and transfer ( x; y; z ) space to (~ x;~y;~z) space, with ~x= (x+ 2) cos t;~y=y;~z= (x+ 2) sin t: This corresponds to taking the ( x; y) plane pictured in Fig. 15.5 and rotating it about the vertical axis x=¡2. We follow this with the linear map to the ( u; v) plane, u= ~z¡~x=2,v= ~y¡~x=2. Consequently, to produce Figs. 15.6{15.8, we draw curves ( u(t); v(t)), with (15.32) u(t) = (x(t) + 2)³ sint¡cost 2´ ; v(t) =y(t)¡(x(t) + 2)cost 2: For initial data, we take xandyas in (15.31) and z= 0. We use a fourth order Runge-Kutta scheme. Fig. 15.6 draws such curves when ( x; y; z ) solve (15.28) with r= 0. Note that in all but the ¯fth portrait, the orbits lie on smooth donut-shaped surfaces (called tori). The ¯fth portrait depicts the homoclinic orbit, which spends most of its time near the origin in ( x; y)-space. It lies on a surface that is smooth except along a curve, where it has a corner. Figure 15.7 gives this representation of orbits of (15.28), with (15.33) r= 0:1: Two of the six orbits seem to lie on smooth tori (one very thin, the other somewhat deformed). The other four are all apparently a mess, and also, 15. Chaos in multidimensional systems 371 Figure 15.5 372 4. Nonlinear Systems of Di®erential Equations Figure 15.6 15. Chaos in multidimensional systems 373 Figure 15.7 374 4. Nonlinear Systems of Di®erential Equations Figure 15.8 15. Chaos in multidimensional systems 375 Figure 15.9 apparently, about the same mess. In Fig. 15.8 we present such orbits with initial data (15.34) x=p 2 +k 20;0·k·5; y = 0; interpolating from the ¯fth orbit of Fig. 15.7 halfway to the sixth orbit. Here the ¯rst three orbits appear chaotic and the last three appear to lie on smooth surfaces. An alternative to depicting orbits of the system (15.28){(15.29) is to depict orbits of the associated Poincar¶ e map , characterized as follows. Take an initial point p= (x0; y0;0). Solve (15.28) with this initial data, and then set q= (x(2¼); y(2¼);2¼). The nature of the mapping on the third coordinate is trivial in this case, so we just consider (15.35) ( x(0); y(0))7!(x(2¼); y(2¼)): This is the Poincar¶ e map associated to the system (15.28). The Poincar¶ e map is de¯ned in a more general context. Let Xbe a smooth vector ¯eld on ­ ½Rnand let Sbe an ( n¡1)-dimensional surface transversal to X, i.e., Xis nowhere tangent to S. Under certain circum- stances, one has a Poincar¶ e map (15.36) P:O ¡! S; de¯ned on an open subset O ½ S, where p2 O andP(p) =qis the point ©t X(p) with smallest t >0 such that ©t X(p)2S. See Fig. 15.9. 376 4. Nonlinear Systems of Di®erential Equations In the setting of (15.28), (15.35), the orbit for Poincar¶ e map ( x(0); y(0)) = (x;0), with xas in (15.34), is presented in Fig. 15.10, which can be appre- ciated in light of Fig. 15.8. Each picture in Fig. 15.10 is made of 9000 points in the orbit of the Poincar¶ e map (or rather an approximation via a Runge-Kutta di®erence scheme). The ¯rst three pictures seem to show or- bits spread out in a 2-dimensional region, while the last three seem to show orbits lying on smooth curves. Going further, each of the last three pictures in Fig. 15.10 suggest the following: Assertion. There is a region ­½R2, smoothly equivalent to the disk (15.37) D=fx2R2:kxk ·1g; that is to say, there is a smooth one-to-one map ':­!Dwith smooth inverse '¡1:D!­, and the Poincar¶ e map takes ­ into itself, i.e., (15.38) P:­¡!­: Granted this, we can make use of the following result, known as Brouwer's ¯xed-point theorem. Theorem. Each smooth map (15.39) Ã:D¡!D has a ¯xed point, i.e., there exists p2Dsuch that Ã(p) =p. See Appendix F for a proof of this result. Given the assertion above, we can take Ã='±P±'¡1and conclude that (15.40) P(q) =q; q ='¡1(p): Such ¯xed points of the Poincar¶ e map give rise to periodic solutions to the associated systems of di®erential equations (in this case, (15.27)). Estab- lishing the existence of periodic solutions is one of many uses for Poincar¶ e maps. We refer to references cited in the next paragraph for discussions of other uses. Understanding how the chaotic looking orbits for the Lorenz and Du±ng systems and other systems arechaotic has engendered a lot of work. For more material on this, we particularly recommend the Introduction to Chaos in Chapter 2 of [ GH], which treats four examples, including the Lorenz system and the forced Du±ng system. Other material on chaotic systems can be found in [ AS], [AP], [HK], [HSD ], [J], and [ LL]. A detailed study of the Lorenz system is given in [ Sp]. 15. Chaos in multidimensional systems 377 Figure 15.10 378 4. Nonlinear Systems of Di®erential Equations Exercises 1. Consider the double pendulum system, in the limit m2= 0, given by (9.35){(9.36). Substitute (15.41) µ1(t)¼rcos!t; ! =rg `1 into (9.36), expand in powers of r, and throw away terms containing second and higher powers of r. Show that you get (15.42) µ00 2(t) +g `2sinµ2(t) =r!2`1 `2cosµ2(t) cos !t: Exercises 2{9 are for readers who can use a programming language, with graphics capabilities. 2. Write a program to exhibit solution curves of (15.42), in a fashion anal- ogous to the treatment of (15.27), involving an analogue of (15.28). Try various values of r; g=` 2, etc., and see when the behavior is more chaotic or less chaotic. 3. Write a program to exhibit solutions to the full double pendulum system (9.15){(9.16). Take, e.g., m1=m2= 1; `1=`2= 1, and variants. 4. Examine orbits and Poincar¶ e maps for the periodically forced Du±ng equation for other values of r, such as r= 0:2;0:05;10¡2;10¡3, etc. Exercises 5{8 deal with systems of the form (15.43)d2 dt2µx y¶ =¡rV(x; y): These are 2 £2 second order systems, which convert to 4 £4 ¯rst order systems. Energy conservation leads to °ows on 3-dimensional constant energy surfaces. In each case, write a program to exhibit solution curves (x(t); y(t)). See whether the displayed solutions seem to be regular or chaotic. 5. Take V(x; y) =x2+axy+y4: A. The derivative in several variables 379 Try various a2[0;10]. 6. Take V(x; y) =x4+axy+y4; a2[¡2;2]: 7. Take V(x; y) =x4+ax2y+y4; a2[¡1;1]: 8. Take V(x; y) =1 2(x2+y2) +a(x4¡x2y+y4); a2[0;1]: 9. Taking o® from models in xx13{14, see if you can construct models of interactions of 3 species that exhibit chaotic behavior. A. The derivative in several variables To start this section o®, we de¯ne the derivative and discuss some of its basic properties. Let Obe an open subset of Rn;andF:O !Rma continuous function. We say Fis di®erentiable at a point x2 O;with derivative L;if L:Rn!Rmis a linear transformation such that, for y2Rn;small, (A.1) F(x+y) =F(x) +Ly+R(x; y) withkR(x; y)k=o(kyk), i.e., (A.2)kR(x; y)k kyk!0 as y!0: We denote the derivative at xbyDF(x) =L:With respect to the standard bases ofRnandRm; DF (x) is simply the matrix of partial derivatives, (A.3) DF(x) =µ@Fj @xk¶ ; so that, if v= (v1; : : : ; v n)t;(regarded as a column vector) then (A.4) DF(x)v=³X k@F1 @xkvk; : : : ;X k@Fm @xkvk´t : 380 4. Nonlinear Systems of Di®erential Equations It will be shown below that Fis di®erentiable whenever all the partial derivatives exist and are continuous onO:In such a case we say Fis aC1 function on O:More generally, Fis said to be Ckif all its partial derivatives of order ·kexist and are continuous. If FisCkfor all k, we say FisC1. Sometimes one might want to di®erentiate an Rm-valued function F(x; t) only with respect to x. In that case, if F(x+y; t) =F(x; t) +Ly+R(x; y; t ); withkR(x; y; t )k=o(kyk), we write DxF(x; t) =L. We now derive the chain rule for the derivative. Let F:O !Rmbe di®erentiable at x2 O;as above, let Ube a neighborhood of z=F(x) in Rm;and let G:U!Rkbe di®erentiable at z:Consider H=G±F:We have (A.5)H(x+y) =G(F(x+y)) =G¡ F(x) +DF(x)y+R(x; y)¢ =G(z) +DG(z)¡ DF(x)y+R(x; y)¢ +R1(x; y) =G(z) +DG(z)DF(x)y+R2(x; y) with kR2(x; y)k kyk!0 as y!0: Thus G±Fis di®erentiable at x;and (A.6) D(G±F)(x) =DG(F(x))¢DF(x): In case k= 1, so G:U!R, we can rewrite (A.6) as (A.7) D(G±F)(x) =rG(F(x))tDF(x); where rG(y)t= (@G=@y 1; : : : ; @G=@y m). If in addition n= 1, so Fis a function of one variable x2 O ½R, with values in Rm, this in turn leads to (A.8)d dxG(F(x)) =rG(F(x))¢F0(x): This leads to such formulas as (3.10). Another useful remark is that, by the Fundamental Theorem of Calculus, applied to '(t) =F(x+ty); (A.9) F(x+y) =F(x) +Z1 0DF(x+ty)y dt; A. The derivative in several variables 381 provided DFis continuous. A closely related application of the Fundamental Theorem of Calculus is that, if we assume F:O !Rmis di®erentiable in each variable separately, and that each @F=@x jis continuous on O;then (A.10) F(x+y) =F(x) +nX j=1£ F(x+zj)¡F(x+zj¡1)¤ =F(x) +nX j=1Aj(x; y)yj; Aj(x; y) =Z1 0@F @xj¡ x+zj¡1+tyjej¢ dt; where z0= 0; zj= (y1; : : : ; y j;0; : : : ; 0);andfejgis the standard basis of Rn:Now (A.10) implies Fis di®erentiable on O;as we stated below (A.4). Thus we have established the following. Proposition A.1. IfOis an open subset of RnandF:O !Rmis of class C1, then Fis di®erentiable at each point x2 O. For the study of higher order derivatives of a function, the following result is fundamental. Proposition A.2. Assume F:O !Rmis of class C2, with Oopen inRn. Then, for each x2 O;1·j; k·n, (A.11)@ @xj@F @xk(x) =@ @xk@F @xj(x): To prove Proposition A.2, it su±ces to treat real valued functions, so consider f:O !R. For 1 ·j·n, we set (A.12) ¢ j;hf(x) =1 h¡ f(x+hej)¡f(x)¢ ; where fe1; : : : ; e ngis the standard basis of Rn. The mean value theorem (for functions of xjalone) implies that if @jf=@f=@x jexists on O, then, for x2 O; h > 0 su±ciently small, (A.13) ¢ j;hf(x) =@jf(x+®jhej); for some ®j2(0;1), depending on xandh. Iterating this, if @j(@kf) exists onO, then, for x2­; h > 0 su±ciently small, (A.14)¢k;h¢j;hf(x) =@k(¢j;hf)(x+®khek) = ¢ j;h(@kf)(x+®khek) =@j@kf(x+®khek+®jhej); with ®j; ®k2(0;1). Here we have used the elementary result (A.15) @k¢j;hf= ¢ j;h(@kf): We deduce the following. 382 4. Nonlinear Systems of Di®erential Equations Proposition A.3. If@kfand@j@kfexist on Oand@j@kfis continuous atx02 O, then (A.16) @j@kf(x0) = lim h!0¢k;h¢j;hf(x0): Clearly (A.17) ¢ k;h¢j;hf= ¢ j;h¢k;hf; so we have the following, which easily implies Proposition A.2. Corollary A.4. In the setting of Proposition A.3, if also @jfand@k@jf exist on Oand@k@jfis continuous at x0, then (A.18) @j@kf(x0) =@k@jf(x0): IfUandVare open subsets of RnandF:U!Vis aC1map, we sayFis a di®eomorphism of Uonto Vprovided Fmaps Uone-to-one and onto V, and its inverse G=F¡1is aC1map. If Fis a di®eomorphism, it follows from the chain rule that DF(x) is invertible for each x2U. We now state a partial converse of this, the Inverse Function Theorem, which is a fundamental result in multivariable calculus. Theorem A.5. LetFbe aCkmap from an open neighborhood ­ofp02Rn toRn;with q0=F(p0):Assume k¸1. Suppose the derivative DF(p0)is invertible. Then there is a neighborhood Uofp0and a neighborhood Vof q0such that F:U!Vis one-to-one and onto, and F¡1:V!Uis aCk map. (So F:U!Vis a di®eomorphism.) Proofs of Theorem A.5 can be found in a number of texts, including [ LS] and Chapter 1 of [ T]. B. Convergence, compactness, and continuity We discuss a number of notions and results related to convergence in Rn, of use in this chapter. First, a sequence of points ( pj) inRnconverges to a limit p2Rn(we write pj!p) if and only if (B.1) kpj¡pk ¡! 0: Herek ¢ k is the norm on Rnarising in x10 of Chapter 2, and the meaning of (B.1) is that for every " >0 there exists Nsuch that (B.2) j¸N=) kpj¡pk< ": B. Convergence, compactness, and continuity 383 A set S½Rnis said to be closed if and only if (B.3) pj2S; p j!p=)p2S: The complement RnnSof a closed set Sisopen. Alternatively, ­ ½Rnis open if and only if, given q2­, there exists " >0 such that B"(q)½­, where (B.4) B"(q) =fp2Rn:kp¡qk< "g; soqcannot be a limit of a sequence of points in Rnn­. An important property of Rniscompleteness , a property de¯ned as follows. A sequence ( pj) of points in Rnis called a Cauchy sequence if and only if (B.5) kpj¡pkk ¡! 0;asj; k! 1 : It is easy to see that if pj!pfor some p2Rn, then (B.5) holds. The completeness property is the converse. Theorem B.1. If(pj)is a Cauchy sequence in Rn, then it has a limit, i.e., (B.1) holds for some p2Rn. Since convergence pj!pinRnis equivalent to convergence in Rof each component, it is the fundamental property of completeness of Rthat is the issue. This is discussed in [ BS], from an axiomatic viewpoint, and in [ Kr], and also [ T2], from a more constructive viewpoint. Completeness provides a path to the following key notion of compactness . A set K½Rnis compact if and only if the following property holds. (B.6)Each in¯nite sequence ( pj) inKhas a subsequence that converges to a point in K. It is clear that if Kis compact, then it must be closed. It must also be bounded, i.e., there exists R <1such that K½BR(0). Indeed, if Kis not bounded, there exist pj2Ksuch that kpj+1k ¸ k pjk+ 1. In such a case, kpj¡pkk ¸1 whenever j6=k, so (pj) cannot have a convergent subsequence. The following converse statement is a key result. Theorem B.2. IfK½Rnis closed and bounded, then it is compact. We start with a special case. Proposition B.3. Each closed bounded interval I= [a; b]½Ris compact. 384 4. Nonlinear Systems of Di®erential Equations Proof. Let (pj) be an in¯nite sequence in [ a; b]; j2Z+. Divide Iinto two halves, I0= [a;(a+b)=2]; I1= [(a+b)=2; b]. Ifpj2I0for in¯nitely many j, pick some pj02I0, and set a1= 0. Otherwise, pick some pj02I1, and seta1= 1. Set q0=pj0. Now divide Ia1into two equal intervals, Ia10andIa11. Ifpj2Ia10for in¯nitely many j, pick pj12Ia10; j1> j0. Otherwise, pick pj12Ia11; j1> j0. Set q1=pj1. Continue. One gets ( qj), a subsequence of ( pj), with the property that (B.7) jqj¡qj+kj ·2¡jjb¡aj;8k¸0: Thus ( qj) is a Cauchy sequence, so by the completeness of R, it converges, to the desired limit p2[a; b]. From Proposition B.3 it is easy enough to show that any closed, bounded box (B.8) B=f(x1; : : : ; x n)2Rn:aj·xj·bj;8jg is compact. If K½Rnis closed and bounded, it is a subset of such a box, and clearly every closed subset of a compact set is compact, so we have Theorem B.2. We next discuss continuity. If S½Rn, a function (B.9) f:S¡!Rm is said to be continuous at p2Sprovided (B.10) pj2S; p j!p=)f(pj)!f(p): Iffis continuous at each p2S, we say fis continuous on S. The following two results give important connections between continuity and compactness. Proposition B.4. IfK½Rnis compact and f:K!Rmis continuous, then f(K)is compact. Proof. If (qk) is an in¯nite sequence of points in f(K), pick pk2Ksuch thatf(pk) =qk. IfKis compact, we have a subsequence pkº!pinK, and then qkº!f(p) inRm. This leads to the second connection. B. Convergence, compactness, and continuity 385 Proposition B.5. IfK½Rnis compact and f:K!Rmis continuous, then there exists p2Ksuch that (B.11) kf(p)k= max x2Kkf(x)k; and there exists q2Ksuch that (B.12) kf(q)k= min x2Kkf(x)k: The meaning of (B.11) is that kf(p)k ¸ k f(x)kfor all x2K, and the meaning of (B.12) is similar. For the proof, consider (B.13) g:K¡!R; g(p) =kf(p)k: This is continuous, so g(K) is compact. Hence g(K) is bounded; say g(K)½ I= [a; b]. Repeatedly subdividing Iinto equal halves, as in the proof of Proposition B.3, at each stage throwing out subintervals that do not intersect g(K) and keeping only the leftmost and rightmost amongst those remaining, we obtain ®2g(K) and ¯2g(K) such that g(K)½[®; ¯]. Then ®=f(q) and¯=f(p) for some pandq2Ksatisfying (B.11){(B.12). A variant of Proposition B.5, with a very similar proof, is that if K½Rn is compact and f:K!Ris continuous, then there exist p; q2Ksuch that (B.14) f(p) = max x2Kf(x); f(q) = min x2Kf(x): We next de¯ne the closure Sof a set S½Rn, to consist of all points p2Rnsuch that B"(p)\S6=;for all " >0. Equivalently, p2Sif and only if there exists an in¯nite sequence ( pj) of points in Ssuch that pj!p. Now we de¯ne sup Sand inf S. First, let S½Rbe nonempty and bounded from above, i.e., there exists R <1such that x·Rfor all x2S. Hence x·Rfor all x2S. In such a case, there exists an interval [ R¡k; R] whose intersection with Sis nonempty, hence compact. We set (B.15) sup S= max S\[R¡k;R]x; the right side well de¯ned by (B.14), with f(x) =x. There is a similar de¯nition of (B.16) inf S; when Sis bounded from below. We establish some further properties of compact sets K½Rn, leading to the important result, Proposition B.9 below. 386 4. Nonlinear Systems of Di®erential Equations Proposition B.6. LetK½Rnbe compact. Assume X1¾X2¾X3¾ ¢¢¢ form a decreasing sequence of closed subsets of K. If each Xm6=;, then \mXm6=;. Proof. Pick xm2Xm. If Kis compact, ( xm) has a convergent subse- quence, xmk!y. Since fxmk:k¸`g ½Xm`, which is closed, we have y2 \mXm. Corolary B.7. LetK½Rnbe compact. Assume U1½U2½U3½ ¢¢¢ form an increasing sequence of open sets in Rn. If[mUm¾K, then UM¾K for some M. Proof. Consider Xm=KnUm. Before getting to Proposition B.9, we bring in the following. Let Q denote the set of rational numbers, and let Qndenote the set of points in Rnall of whose components are rational. The set Qn½Rnhas the following \denseness" property: given p2Rnand" >0, there exists q2Qnsuch thatkp¡qk< ". Let (B.17) R=fBrj(qj) :qj2Qn; rj2Q\(0;1)g: Note thatQandQnarecountable , i.e., they can be put in one-to-one corre- spondence with N. Hence Ris a countable collection of balls. The following lemma is left as an exercise for the reader. Lemma B.8. Let­½Rnbe a nonempty open set. Then (B.18) ­ =[ fB:B2 R; B½­g: To state the next result, we say that a collection fU®:®2 Ag covers K ifK½ [ ®2AU®. If each U®½Rnis open, it is called an open cover of K. If B ½ A andK½ [ ¯2BU¯, we say fU¯:¯2 Bg is a subcover. Proposition B.9. IfK½Rnis compact, then it has the following property. (B.19) Every open cover fU®:®2 Ag ofKhas a ¯nite subcover. Proof. By Lemma B.8, it su±ces to prove the following. (B.20)Every countable cover fBj:j2NgofKby open balls has a ¯nite subcover. For this, we set (B.21) Um=B1[ ¢¢¢ [ Bm and apply Corollary B.7. C. Critical points that are saddles 387 C. Critical points that are saddles LetFbe aC3vector ¯eld on ­ ½Rn, with a critical point at p2­. We say pis a simple critical point if L=DF(p) has no eigenvalues that are purely imaginary (or zero). From here on we assume this condition holds. As seen in Chapter 2, we can write (C.1) Cn=W+©W¡; where W+is the direct sum of the generalized eigenspaces of Lassociated to eigenvalues with positive real part and W¡is the direct sum of the gen- eralized eigenspaces associated to eigenvalues with negative real part. Since L2M(n;R), non-real eigenvalues of Lmust occur in complex conjugate pairs, and (C.2) Rn=V+©V¡; V§=W§\Rn: We have (C.3) v2W§=)etLv!0 as t! ¨1 ; and a fortiori (C.4) v2V§=)etLv!0 as t! ¨1 : We say the critical point at pis a source if V¡= 0, a sink if V+= 0, and a saddle if V¡6= 0 and V+6= 0. The fact that (C.5) V¡=Rn=)©t F(x)!past!+1; forxsu±ciently close to p, where ©t Fis the °ow generated by F, was proven inx3 (cf. Proposition 3.4), and similarly we have (C.6) V+=Rn=)©t F(x)!past! ¡1 ; forxsu±ciently close to p. The purpose of this appendix is to discuss the saddle case, where n+= dim V+>0 and n¡= dim V¡>0. In such a case, as advertised in x3, there is a neighborhood Uofpand there are C1surfaces S§, of dimension n§, such that (C.7) fpg=S+\S¡; and (C.8) x2S§=)©t F(x)!past! ¨1 : 388 4. Nonlinear Systems of Di®erential Equations The surfaces S¡andS+are called, respectively, the stable and unstable manifolds of Fatp. They have the further property that if °is aC1 curve in S+(respectively, S¡), and °(0) = p, then °0(0)2V+(respectively, V¡). In addition, given " >0, there exists ± >0 such that if x2UnS¡ but dist( x; S¡)< ±; then for some t1>0;k©t1 F(x)¡pk< ", and for all t¸t1;dist(©t F(x); S+)< ", at least until ©t F(x) exits U. We want to demonstrate this result. For simplicity of presentation, we concentrate on the case n= 2 (and n+=n¡= 1). However, the argument we present can be modi¯ed to treat saddles in higher dimension. We make some preliminary constructions. Relabeling the coordinates, we can assume p= 0. Altering Foutside some neighborhood of p= 0 if necessary, we can assume Fis aC3vector ¯eld on Rnand there exists C <1such that (C.9) kF(x)k ·Ckxk;8x2Rn: Hence, as seen in x3 (Exercise 3), ©t F(x) is well de¯ned for all x2Rn; t2R. Applying the fundamental theorem of calculus twice gives (C.10) F(x) =Lx+X j;kxjxkGjk(x); where (C.11) L=DF(0); andGjkareC1vector ¯elds, given by (C.12) Gjk(x) =Z1 0Z1 0@2 @xk@xjF(stx)ds dt: We de¯ne the family of vector ¯elds F"by (C.13) F"(x) =1 "F("x); for" >0. By (C.10), (C.14) F"(x) =Lx+"G"(x); where (C.15) G"(x) =X j;kxjxkGjk("x): C. Critical points that are saddles 389 Passing to the limit "!0 gives F0(x) = Lx. Results of x2 yield the following. Lemma C.1. Given ± >0; T < 1, there exists "0="0(±; T; F )>0such that for all "2(0; "0], (C.16)kxk ·2;jtj ·T;k©s F"(x)k ·28s2[0; t] =) k©t F"(x)¡etLxk ·±: Specializing to n= 2, we can assume that (C.17) L=µa ¡b¶ ; a; b > 0: We take the box (C.18) O=f(x1; x2) :jx1j;jx2j ·1g; and set (C.19) Ok= 2¡kO: We de¯ne four families of maps (C.20) '"j; Ã"j:h ¡1 2;1 2i ¡![¡1;1]; j = 1;2;0·"·1; as follows. For j= 1, de¯ne t"(s) as the smallest positive number such that ©¡t"(s) F"³ s;1 2´ 2 f(¾;1) :¡1·¾·1g; and then set '"1(s) to be the x1-coordinate of ©¡t"(s) F"(s;1=2). To give an alternative description, we are mapping the top edge of O1(identi¯ed with [¡1=2;1=2]) to the top edge of O(identi¯ed with [ ¡1;1]) by the backward °ow generated by F". Similarly de¯ne '"2via the backward °ow map of the bottom edge of O1to the bottom edge of O, and de¯ne Ã"1andÃ"2via the forward °ow maps of the right and left edges of O1to the corresponding edges of O. See Fig. C.1. It is readily veri¯ed that these maps are contrac- tions for "= 0, where F0(x) =Lx, i.e., there exists A=A(a; b)<1 such that (C.21)j'"j(s)¡'"j(t)j ·Ajs¡tj; jÃ"j(s)¡Ã"j(t)j ·Ajs¡tj; 390 4. Nonlinear Systems of Di®erential Equations Figure C.1 for all s; t2[¡1=2;1=2]. Results of x2 then establish the following. Lemma C.2. There exists "1="1(F)>0andA=A(F)<1such that whenever 0·"·"1, the maps '"jandÃ"jin (C.20) are well de¯ned on [¡1=2;1=2]and (C.21) holds for all s; t2[¡1=2;1=2]. We make a further adjustment. Take "2·min("1(F); "0(1=10;10; F)). Further shrinking "2is nesessary, arrange that, whenever "2(0; "2], (C.22) kxk ·2 =) kF"(x)¡Lxk ·1 2kLxk; so that, if ( x1; x2)2 O,F"(x1; x2) points down if x22[1=2;1], up if x22 [¡1;¡1=2], left if x12[1=2;1], and right if x12[¡1;¡1=2]. Now replace F byF"2, denoting this scaled vector ¯eld by F. Then (C.16) holds with T= 10 and±= 1=10 for all "2(0;1] and (C.21) holds for all s; t2[¡1=2;1=2], with A <1, for all "2(0;1]. For notational simplicity, set (C.23) ©t k= ©t F"; " = 2¡k: Note that dilation by the factor 2ktakes the °ow ©t FonOkto the °ow ©t k onO. With these preliminaries done, we start in earnest our demonstration that the °ow generated by Fhas saddle-like behavior near the critical point p= 0. Denote by T;B;L, andRthe top, bottom, left, and right edges of O, and similarly denote by Tk;Bk;Lk, andRkthe top, bottom, left, and right C. Critical points that are saddles 391 Figure C.2 sides of Ok. Then the maps (C.20) can by slight abuse of terminology be labeled (C.24)'"1:T1! T; ' "2:B1! B; Ã"1:R1! R; à "2:L1! L: Pick two points p0`; p0r2 Tsuch that for some t0`; t0r2(0;1), (C.25) p¤ 0`= ©t0` F(p0`)2 L; p¤ 0r= ©t0r F(p0r)2 R; and pick two points q0`; q0r2 Bsuch that for some s0`; s0r2(0;1), (C.26) q¤ 0`= ©s0` F(q0`)2 L; q¤ 0r= ©s0r F(q0r)2 R: See Fig. C.2. The possibility to do this is guaranteed by Lemma C.1. Denote the orbits of ©t 0= ©t Fthrough p0`; p0rby°0`; °0rand those through q0`; q0r by¾0`; ¾0r. Let us call the construction just described Step 0. To continue, at Step 1, pick ~p1`;~p1r2 T1and ~q1`;~q1r2 B1such that the following holds. Note that 2~p1`;2~p1r2 Tand 2~ q1`;2~q1r2 B. We require that, for some t1`; t1r2(0; T1), (C.27) ©t1` 1(2~p1`)2 L;©t1r 1(2~p1r)2 R; and for some s1`; s1r2(0; T1), (C.28) ©s1` 1(2~q1`)2 L;©s1r 1(2~q1r)2 R: The conditions (C.27) and (C.28) are equivalent to (C.29) ~ p¤ 1`= ©t1` F(~p1`)2 L1;~p¤ 1r= ©t1r F(~p1r)2 R 1; 392 4. Nonlinear Systems of Di®erential Equations Figure C.3 and (C.30) ~ q¤ 1`= ©s1` F(~q1`)2 L1;~q¤ 1r= ©s1r F(~q1r)2 R 1: See Fig. C.3. Denote the orbits of Fthrough ~ p1`;~p1rby°1`; °1rand those through ~ q1`;~q1rby¾1`; ¾1r. When picking ~ p1`;~p1r;~q1`, and ~ q1r, one can and should enforce the following condition. If °0`intersects T1, ~p1`should be to the right of such an intersection, if °0rintersects T1, ~p1rshould be to the left of such an intersection, and similarly for cases when ¾0`or¾0rintersect B1. Also, we can take T1>1. (More on this below.) Now we continue the orbits °1`; °1r; ¾1`, and ¾1rforward and back- ward, until they intersect the boundary of O, at points p1`; p1r; q1`; q1rand p¤ 1`; p¤ 1r; q¤ 1`; q¤ 1r, as illustrated in Fig. C.4. That such an intersection must occur is guaranteed by (C.22). This, together with the fact that orbits of F cannot intersect, guarantees that (C.31) p0`< p1`< p1r< p0r; in the sense that p < p0means pis to the left of p0. In a similar sense, made clear in Fig. C.4, we have (C.32)q0`< q1`< q1r< q0r; p¤ 0r< p¤ 1r< q¤ 1r< q¤ 0r; p¤ 0`< p¤ 1`< q¤ 1`< q¤ 0`: C. Critical points that are saddles 393 Figure C.4 Furthermore, as a consequence of (C.21), we have (C.33)jp1`¡p1rj ·Aj~p1`¡~p1rj ·A; jq1`¡q1rj ·Aj~q1`¡~q1rj ·A; jp¤ 1r¡q¤ 1rj ·Aj~p¤ 1r¡~q¤ 1rj ·A; jp¤ 1`¡q¤ 1`j ·Aj~p¤ 1`¡~q¤ 1rj ·A: We proceed iteratively. At step k, pick ~ pk`;~pkr2 Tkand ~qk`;~qkr2 Bk such that the following holds. Note that 2k~pk`;2k~pkr2 Tand 2k~qk`;2k~qkr2 B. We require that, for some tk`; tkr2(0; Tk), (C.34) ©tk` k(2k~pk`)2 L;©tkr k(2k~pkr)2 R; and for some sk`; skr2(0; Tk), (C.35) ©sk` k(2k~qk`)2 L;©skr k(2k~qkr)2 R: The conditions (C.34) and (C.35) are equivalent to (C.36) ~ p¤ k`= ©tk` F(~pk`)2 Lk;~p¤ kr= ©tkr F(~pkr)2 R k; and (C.37) ~ q¤ k`= ©sk` F(~qk`)2 Lk;~q¤ kr= ©skr F(~qkr)2 R k: Denote the orbits of Fthrough ~ pk`;~pkrby°k`; °kr, and those through ~ qk`;~qkr by¾k`; ¾kr. When picking ~ pk`;~pkr;~qk`, and ~ qkr, one can and should enforce 394 4. Nonlinear Systems of Di®erential Equations Cigure C.5 the following condition. If °k¡1;`intersects Tk, ~pk`should lie to the right of such a point of intersection, if °k¡1;rintersects Tk, ~pkrshould lie to the left of such a point of intersection, and similarly for cases where ¾k¡1;`or¾k¡1;r intersect Bk. At this point it is useful to note that, by Lemma D.1, we can take (C.38) Tk! 1 ask! 1 ; and hence take (with z=`orr) (C.39) k2k~pkz¡(0;1)k ·´k;k2k~qkz¡(0;1)k ·´k; ´ k!0 ask! 1 : It then follows that (again with z=`orr) (C.40) k~p¤ kz¡(2¡k;0)k ·2¡k~´k;k~q¤ kz¡(2¡k;0)k ·2¡k~´k;~´k!0 ask! 1 : Now we continue the orbits °k`; °kr; ¾k`, and ¾krforward and back- ward, until they intersect the boundary of O, at points pk`; pkr; qk`; qkr, andp¤ k`; p¤ kr; q¤ k`; q¤ kr, as illustrated in Fig. C.5. That such intersections must occur is guaranteed by (C.22). As before, the fact that orbits of ©t Fcannot intersect guarantees that (C.41) p0`<¢¢¢< pk`< pkr<¢¢¢< p0r; in the sense speci¯ed in (C.31), and, as in (C.32), (C.42)q0`<¢¢¢< qk`< qkr<¢¢¢< q0r; p¤ 0r<¢¢¢< p¤ kr< q¤ kr<¢¢¢< q¤ 0r; p¤ 0`<¢¢¢< p¤ k`< q¤ k`<¢¢¢< q¤ 0`: C. Critical points that are saddles 395 Figure C.6 Furthermore, as a consequence of (C.21), we have (C.43)jpk`¡pkrj ·Akj~pk`¡~pkrj ·Ak2¡k´k; jqk`¡qkrj ·Akj~qk`¡~qkrj ·Ak2¡k´k; jp¤ k`¡p¤ krj ·Akj~p¤ k`¡~p¤ krj ·Ak2¡k~´k; jq¤ k`¡q¤ krj ·Akj~q¤ k`¡~q¤ krj ·Ak2¡k~´k: In particular, these distances are converging to 0 quite rapidly. We obtain limits (C.44)pk`; pkr!pt2 T; q k`; qkr!pb2 B; p¤ k`; q¤ kr!p¤ r2 R; p¤ k`; q¤ k`!p¤ `2 L: See Fig. C.6. We have (C.45) ©t F(pt);©t F(pb)!0 as t!+1; and (C.46) ©t F(p¤ `);©t F(p¤ r)!0 as t! ¡1 ; since the paths in (C.45) meet each Okfor large positive tand those in (C.46) meet each Okfor large negative t. Furthermore, by (C.39){(C.40), plus the fact that all these paths solve dx=dt =F(x), the curves in (C.45) ¯t together to form a C1curve tangent to the x2-axis at p= 0, and those in (C.46) ¯t together to form a C1curve tangent to the x1-axis at p= 0. 396 4. Nonlinear Systems of Di®erential Equations We sketch how to treat the case n= 3; n+= 2; n¡= 1. In place of (C.17), we can take (C.47) L=µA ¡b¶ ; A2M(2;R); b > 0; and, via Lemma 3.5, arrange that (C.48) Av¢v¸akvk2; a > 0;8v2R2: In place of (C.18), we use the cylinder (C.49) O=f(x1; x2; x3) :x2 1+x2 2·1;jx3j ·1g: with boundary (C.50) @O=T [ B [ L ; where TandB(the top and bottom) are disks and L(the side) is S1£[¡1;1]. We then take Ok= 2¡kO, with boundary Tk[ Bk[ Lk. Parallel to (C.24), we have (at least for small ") maps (C.51)'"1:T1! T; ' "2:B1! B; Ã":L1! L; with '"jde¯ned by backward °ow of ©t F"andÃ"de¯ned by forward °ow. Again the maps '"jare contractions for small ". The maps Ã"are not con- tractions, but composing them on the left with the projection S1£[¡1;1]! [¡1;1] produces a contraction, for small ", and this is what one needs. In place of a pair of initial data on Tand a pair on B, one takes a circle of initial data on Tand one on B. Applying ©t Fyields a pair of °ared tubes, as pictured in Fig. C.7. From here, an iteration produces nested families of such °ared tubes, converging in on the one-dimensional stable manifold S¡and the two-dimensional unstable manifold S+. The interested reader is invited to ¯ll in the details, and work out the higher dimensional cases. See also [CL] and [ Hart ] for other approaches to this result. D. Periodic solutions of x00+x="Ã(x) Equations of the form (D.1) x00+x="Ã(x) with \small" "arise in a number of cases, and it is of interest to analyze various features of these solutions. For example, as mentioned in x6, the D. Periodic solutions of x00+x="Ã(x) 397 Figure C.7 relativistic correction for planetary motion gives rise to the equation (6.53), which takes the form (D.1) for x=u¡A, with (D.2) Ã(x) = (x+A)2: Another example, (D.3) Ã(x) =¡x3; yields a special case of Du±ng's equation. As we mentioned in x6, solutions to (6.53) tend not to be periodic of period 2 ¼, and this leads to the phe- nomenon of precession of perihelia. It is of general interest to compute the period of a solution to (D.1), and we discuss this problem here. We assume Ãis smooth. We rewrite (D.1) as a ¯rst order system and also explicitly record the dependence on ": (D.4)x0 "(t) =y"(t); y0 "(t) =¡x"(t) +"Ã(x"(t)): We pick a2(0;1) and impose the initial conditions (D.5) x"(0) = a; y "(0) = 0 : Note that if we take (D.6) F"(x; y) =y2 2+x2 2¡"ª(x); 398 4. Nonlinear Systems of Di®erential Equations where ª0(x) =Ã(x), then ( d=dt)F"(x"(t); y"(t)) = 0 for solutions to (D.4), so orbits of (D.4) lie on level curves of F". For "su±ciently small with respect to a, the level curves of F"onf(x; y) :x2+y2·2a2gwill be close to those of F0, that is to say, such level curves of F"will be closed curves, close to circles, and the associated solutions to (D.4){(D.5) will be periodic. The period T(") will have the following two properties, at least for small ": (D.7) T(") = 2 ¼+O("); y "(T(")) = 0 : We will calculate a more precise approximation to T("), accurate for small ". The ¯rst order of business is to calculate accurate approximations to x"(t) and y"(t), valid uniformly for tin an interval containing [0 ;2¼]. It follows from x2 that x"(t) and y"(t) are smooth functions of ". Hence, for each N2N, we can write (D.8)x"(t) =acost+NX k=1Xk(t)"k+R1N(t; "); y"(t) =¡asint+NX k=1Yk(t)"k+R2N(t; "); where (D.9) jRjN(t; ")j ·CKN"N+1;8jtj ·K: We write RjN(t; ") =O("N+1). The coe±cients Xk(t) and Yk(t) satisfy di®erential equations, obtained as follows. We have from (D.8) (D.10) x00 "(t) +x"(t) =NX k=1£ X00 k(t) +Xk(t)¤ "k+O("N+1); while (D.11) "Ã(x"(t)) ="ó acost+NX k=1Xk(t)"k´ +O("N+1) ="h Ã(acost) +NX j=11 j!Ã(j)(acost)³NX k=1Xk(t)"k´ji +O("N+1): We match up the coe±cients of "kin (D.10) and (D.11) to obtain equations forXk(t). The case k= 1 gives (D.12) X00 1(t) +X1(t) =Ã(acost); D. Periodic solutions of x00+x="Ã(x) 399 and from (D.5) the initial conditions are seen to be (D.13) X1(0) = 0 ; X0 1(0) = 0 : The solution to (D.12){(D.13) is given by Duhamel's formula, cf. (4.10) of Chapter 3: (D.14) X1(t) =Zt 0sin(t¡s)Ã(acoss)ds: It is convenient to expand sin( t¡s) and rewrite (D.14) as (D.15) X1(t) = (sin t)Zt 0coss Ã(acoss)ds¡(cost)Zt 0sins Ã(acoss)ds: Regarding Yk(t), we have (D.16) Yk(t) =X0 k(t); for all k, and in particular (D.17) Y1(t) = (cos t)Zt 0coss Ã(acoss)ds+ (sin t)Zt 0sins Ã(acoss)ds: In case Ã(x) is given by (D.2), we have (D.18) Ã(acoss) =a2 2cos 2s+ 2aAcoss+³ A2+a2 2´ ; hence (D.19)Zt 0coss Ã(acoss)ds =³ A2+3a2 4´ sint+aA 2sin 2t+a2 12sin 3t+aAt; and (D.20)Zt 0sins Ã(acoss)ds =³ A2+aA 2+a2 2´ ¡³ A2+5a2 12´ cost¡aA 2cos 2t¡a2 12cos 3t: One can compute higher terms in (D.8). For example, matching up coe±cients of "2in (D.10) and (D.11) yields (D.21) X00 2(t) +X2(t) =Ã0(acost)X1(t): 400 4. Nonlinear Systems of Di®erential Equations Again X2(0) = X0 2(0) = 0, and, parallel to (D.14), we have (D.22) X2(t) =Zt 0sin(t¡s)Ã0(acoss)X1(s)ds: Y2(t) is given by (D.16). One can continue this, but we will leave o® at this point. We return to the problem of approximating the period T("), making use of (D.7). A very e®ective method for solving y"(T) = 0 with T¼2¼is Newton's method, which gives T(") as the limit of Tn("), de¯ned recursively by (D.23) T0(") = 2 ¼; T n+1(") =Tn(")¡y"(Tn(")) y0"(Tn(")): This sequence converges fast: (D.24) T(") =Tn(") +O("2n); provided one has y"(t) evaluated exactly. Given an approximation to y"(t), (D.25) y"(t) = ~y"(t) +O("N); y0 "(t) = ~y0 "(t) +O("N); we can work with eTn("), given by (D.26) eT0(") = 2 ¼;eTn+1(") =eTn(")¡~y"(eTn(")) ~y0(eTn(")); and we get (D.27) T(") =eTn(") +O("N);provided 2n¸N: In particular, taking (D.27) y"(t) =asint+Y1(t)"+O("2); we have (D.28) T(") =eT1(") +O("2); with (D.29)eT1(") = 2 ¼¡~y"(2¼) ~y0"(2¼) = 2¼+1 aY1(2¼)"; D. Periodic solutions of x00+x="Ã(x) 401 hence, by (D.17), (D.30) T(") = 2 ¼+" aZ2¼ 0coss Ã(acoss)ds+O("2): In case Ã(x) is given by (D.2), we have from (D.19) that (D.31)Z2¼ 0coss Ã(acoss)ds= 2¼aA; so in this case (D.32) T(") = 2 ¼(1 +A") +O("2): Given an approximation ~ y"(t) satisfying (D.25) with N= 3 or 4, we can iterate (D.26) once more, obtaining T2(") =T(") +O("N), and so on. We will not pursue the details. We now return to the problem of approximating the solution ( x"(t); y"(t)) of (D.4), and address a limitation of the approximations of the form (D.8). As follows from (D.15){(D.20), the ¯rst order approximation has the form (D.33)x"(t) =acost+X1(t)"+O("2); y"(t) =¡asint+Y1(t)"+O("2); and, in the case that Ã(x) is given by (D.2), (D.34)X1(t) =Xb 1(t) +aAtsint; Y1(t) =Yb 1(t)¡aAtcost; where Xb 1(t) and Yb 1(t) are periodic in t, of period 2 ¼, being sums of products of sin ktand cos kt(0·k·3). In (D.33), the notation O("2) means that, for any given bounded interval [ ¡K; K ], the remainder is bounded by CK"2, fort2[¡K; K ]. However, it is apparent from (D.34) that the accuracy of this approximation breaks down severely on intervals of length ¼1=". In fact, both x"(t) and y"(t) are uniformly bounded, being periodic of period T("). As far as the terms on the right side of (D.33) are concerned, acost+Xb 1(t)"and ¡asint+Yb 1(t)" are uniformly bounded, of period 2 ¼, but (D.35) aA"t sintand¡aA"t cost 402 4. Nonlinear Systems of Di®erential Equations are unbounded as jtj ! 1 . These terms are called secular terms , and it is desirable to have a replacement for (D.8), in which such secular terms do not appear. To get this, we proceed as follows. The functions (D.36) x# "(t) =x"³T(")t 2¼´ ; y# "(t) =y"³T(")t 2¼´ are periodic of period 2 ¼intand smooth in ". Hence we have expansions (D.37)x# "(t) =acost+NX k=1X# k(t)"k+O("N+1); y# "(t) =¡asint+NX k=1Y# k(t)"k+O("N+1): Note that (D.38)d dtx# "(t) =T(") 2¼y# "(t); which leads to a variant of (D.16). We have the following. Proposition D.1. The solution to (D.4){(D.5) has the expansion (D.39)x"(t) =acos2¼t T(")+NX k=1X# k³2¼t T(")´ "k+O("N+1); y"(t) =¡asin2¼t T(")+NX k=1Y# k³2¼t T(")´ "k+O("N+1): Each term in this series is periodic in tof period T("), and the remainders areO("N+1)uniformly for all t2R. It is natural and convenient to set (D.40) X0(t) =X# 0(t) =acost; Y 0(t) =Y# 0(t) =¡asint: It remains to compute X# k(t) and Y# k(t) for k¸1. To this end, set (D.41)T(") 2¼= 1 + °("); °(") ="X `¸0°`"`: D. Periodic solutions of x00+x="Ã(x) 403 If we compare the expressions for x"(t) in (D.8) and (D.39) and make the substitution s= 2¼t=T ("), we obtain (D.42)X k¸0X# k(s)"k=X i¸0Xi(s+°(")s)"i =X i¸0X j¸01 j!X(j) i(s)sj°(")j"i =X i¸0X j¸01 j!X(j) i(s)sj³X `¸0°`"`´j "i+j: We conclude that X# k(s) is equal to the coe±cient of "kin the last power series. For k= 0, we get (D.43) X# 0(s) =X0(s) =acoss; as already noted in (D.40). For k= 1, we get (D.44)X# 1(s) =X1(s) +°0sX0 0(s) =X1(s)¡°0assins: When Ã(x) is given by (D.2), we have from (D.34) that this is =Xb 1(s) +aAssins¡°0assins =Xb 1(s); the last identity by (D.41) and (D.32), which gives °0=Ain this case. Alternatively, since X# 1(s) and Xb 1(s) are periodic in sand the other terms are secular, these secular terms have to cancel. This holds for general Ã(x); X# 1(s) is obtained from X1(s) by striking out the secular terms. One can similarly characterize the higher order terms X# k(t) in (D.37). We forego the details. We end this appendix with an indication of how to extend the scope of (D.1). We treat the pendulum equation (D.45) u00+ sin u= 0; and seek information on small oscillations, solving (D.45) with initial data (D.46) u(0) = ap"; u0(0) = 0 : Thus we set (D.47) x(t) =p" u(t); 404 4. Nonlinear Systems of Di®erential Equations which solves (D.48) x00+sinp"xp"= 0; x(0) = a; x0(0) = 0 : If we set (D.49)sin¿ ¿= 1¡¿2F(¿); F (¿) =1 3!¡¿2 5!+¢¢¢; we get (D.50)x00+x="x3F(p"x) ="x3 3!¡"2x5 5!+¢¢¢: This has a form similar to (D.1), but generalized to (D.51) x00+x="Ã("; x); with Ãsmooth in ( "; x). Treatments of the solutions to (D.1) and their periods T(") extend to the case (D.51). The reader is invited to work out details. E. A dram of potential theory Newton's law of gravitation states that the force a particle of mass m1 located at p2R3exerts on a particle of mass m2located at x2R3is (E.1) F(x) =Gm1m2p¡x kp¡xk3: Here Gis the gravitational constant, given by (6.63). As indicated in Exer- cise 6 of x6, the force that a planet exerts on an external body is the same as what would be exerted if all the mass of the planet were concentrated at its center, in the Newtonian theory. In this appendix we explain why this is true, and in the course of doing so introduce an area of mathematical analysis known as potential theory. We establish this identity of force ¯elds under the hypothesis that the mass distribution of the planet is spherically symmetric about its center. That is to say, we assume the planet, centered atp, has mass density ½, and (E.2) ½(p+Ry) =½(p+y);8R2O(3); y2R3; where we recall from Chapter 2 that O(3) is the set of orthogonal transfor- mations ofR3. Say the planet has radius a, so (E.3) kyk> a=)½(p+y) = 0 : E. A dram of potential theory 405 The planet's mass is (E.4) m1=Z ½(y)dy: If a particle of mass m2is located at x2R3andkp¡xk> a, then the force the planet exerts on this particle is given by (E.5) G(x) =Gm2Zy¡x ky¡xk3½(y)dy: We will show that if (E.2){(E.4) hold and kp¡xk> a, then F(x) =G(x). For notational simplicity, we may as well take (E.6) p= 0; so (E.7) F(x) =¡Gm1m2x kxk3: Note that (E.8) F(x) =¡rV(x); G (x) =¡rW(x); with (E.9) V(x) =¡Gm1m2 kxk; W (x) =¡Gm2Z kyk·a1 kx¡yk½(y)dy; so it su±ces to prove that these potential energies coincide for kxk> a, i.e., (E.10) kxk> a=)V(x) =W(x): As a ¯rst step toward proving (E.10), note that clearly, for all R2O(3), (E.11) V(Rx) =V(x); and furthermore (E.12)W(Rx) =¡Gm1Z1 kRx¡yk½(y)dy =¡Gm1Z1 kRx¡Rzk½(Rz)dz =¡Gm1Z1 kx¡zk½(z)dz =W(x); 406 4. Nonlinear Systems of Di®erential Equations the second identity by change of variable and the third by (E.2). Conse- quently, we have (E.13) V(x) =v(r); W (x) =w(r); r =kxk; and it remains to show that (E.14) r > a =)v(r) =w(r): As another step toward showing this, we note that, given a2(0;1), there exists C <1such that (E.15) kyk ·a;kxk ¸a+ 1 =)¯¯¯1 kxk¡1 kx¡yk¯¯¯·C kxk2; and hence, by (E.4), (E.9) and (E.13), there exists C2<1such that (E.16)r=kxk ¸a+ 1 =) jV(x)¡W(x)j ·C2 kxk2 =) jv(r)¡w(r)j ·C2 r2: The next step toward establishing (E.14) involves the following har- monicity, (E.17) ¢ V(x) = 0 ;8x2R3n0; where ¢ is the Laplace operator, (E.18) ¢ f(x) =@2f @x2 1+@2f @x2 2+@2f @x2 3: To see this, recall from (A.11) of Chapter 1 that (on R3) (E.19) f(x) =g(r) =)¢f(x) =g00(r) +2 rg0(r); and by results on Euler equations from x15 of Chapter 1, (E.20) g00(r) +2 rg0(r) = 0()g(r) =c1 r+c2; Since (E.21) V(x) =v(r) =¡Gm1m2 r; F. Brouwer's ¯xed-point theorem 407 we have (E.17), and hence we also have (E.22) ¢³1 kx¡yk´ = 0 for x6=y; so a direct consequence of the integral formula (E.9) for W(x) is (E.23) ¢ W(x) = 0 for kxk> a: Hence, by (E.13), (E.19), and (E.20), (E.24)r > a =)w00(r) +2 rw0(r) = 0 =)w(r) =c1 r+c2; for some constants c1andc2. This identity together with (E.21) and (E.16) proves (E.14). Hence we have (E.10), so indeed, under the hypotheses (E.2){ (E.4) (and with p= 0), (E.25) kxk> a=)F(x) =G(x): We mention the following re¯nement of (E.23), (E.26) ¢ W= 4¼Gm 2½: This is not needed to establish (E.14), so we will not prove it here. A proof can be found in [ T], Chapter 3, x4. Further exploration of the relation be- tween the Laplace operator and the \potential" function W, through (E.9), leads to the subject of potential theory, addressed in Chapters 3{5 of [ T] and in other books on partial di®erential equations. The earth, the sun, and other planets and stars are approximately spher- ically symmetric, but not exactly so. This leads to further corrections in cal- culations in celestial mechanics. In addition, measurements of the strength of the earth's gravitational ¯eld give information on the inhomogeneities of the earth's composition, leading to the ¯eld of physical geodesy; cf. [ HM]. F. Brouwer's ¯xed-point theorem Here we prove the following ¯xed-point theorem of L. Brouwer, which arose inx15. Take (F.1) D=fx2R2:kxk ·1g: Theorem F.1. Each smooth map F:D!Dhas a ¯xed point. 408 4. Nonlinear Systems of Di®erential Equations The proof proceeds by contradiction. We are claiming that F(x) =x for some x2D. If not, then for each x2Dde¯ne '(x) to be the endpoint of the ray from F(x) tox, continued until it hits (F.2) @D=fx2R2:kxk= 1g: An explicit formula is (F.3)'(x) =x+t(x¡F(x)); t=p b2+ 4ac¡b 2a; a=kx¡F(x)k2; b = 2x¢(x¡F(x)); c = 1¡ kxk2: Here tis picked to solve the equation kx+t(x¡F(x))k2= 1. Note that ac¸0, so t¸0. It is clear that 'would have the following properties: (F.4) ':D!@Dsmoothly ; x2@D)'(x) =x: Such a map is called a smooth retraction. The contradiction that proves Theorem F.1 is provided by the following result, called Brouwer's no-retraction theorem. Theorem F.2. There is no smooth retraction ':D!@DofDonto its boundary. Proof. This proof, also by contradiction, brings in material developed in x4. Suppose we had such a retraction '. Consider the closed curve (F.5) °: [0;2¼]¡!@D; ° (t) = (cos t;sint); and form (F.6) °s(t) ='(s°(t));0·s·1: This would be a smooth family of maps (F.7) °s: [0;2¼]¡!@D; ° s(0) = °s(2¼); such that °1=°and°0(t) ='(0) for all t. The variant of Lemma 4.2 given in Exercise 13 of x4 implies (F.7)Z °sF(y)¢dyis independent of s2[0;1]; for each C1vector ¯eld Fde¯ned on a neighborhood of @Dand satisfying (4.4). Clearly the line integral (F.7) is 0 for s= 0, so we deduce that (F.8)Z °F(y)¢dy= 0 F. Brouwer's ¯xed-point theorem 409 for each such vector ¯eld. In particular, this would apply to the vector ¯eld given by (4.19){(4.20), i.e., (F.9) F(x) =1 kxk2µ¡x2 x1¶ ; which is smooth on R2n0 and satis¯es (4.4) (cf. (4.21)). On the other hand, we compute (F.10)Z °F(y)¢dy=Z2¼ 0(¡sint;cost)¢(¡sint;cost)dt = 2¼; contradicting (F.8) and hence contradicting the existence of such a retrac- tion. The ¯xed-point theorem is valid for all continuous F:D!D. In fact, an approximation argument, which we omit here, can be used to show that if such continuous Fhas no ¯xed point, there is a smooth approximation eF:D!Dthat would also have no ¯xed point. Furthermore, Theorem F.1 holds in ndimensions, i.e., when (F.11) D=fx2Rn:kxk ·1g: The reduction to Theorem F.2, in the setting of (F.11), is the same as above, but the proof of Theorem F.2 in the n-dimensional setting requires a further argument. Proofs using topology can be found in [ GrH ] and [ Mun ]. Proofs using di®erential forms can be found in [ Kan], [T], Chapter 1, and [ T3], Appendix G. We have no space to introduce di®erential forms here, but as shown in [ T], and also in [ AM] and [ Ar], they give rise to many important results in the study of di®erential equations, at the next level.