dynamicalSystems

Chapter 10 Phase Plane Methods Contents 10.1 Autonomous Planar Systems . . . . . . . . . 542 10.2 Planar Constant Linear Systems . . . . . . . 551 10.3 Biological Models . . . . . . . . . . . . . . . . 558 10.4 Mechanical Models . . . . . . . . . . . . . . . 566 Studied here are planar autonomous systems of differential equations. The topics: • Autonomous Planar Systems – Phase Portraits – Stability • Constant Linear Planar Systems – Classification of isolated equilibria – Almost linear systems – Phase diagrams – Nonlinear classifications of equilibria • Biological Models – Predator-prey models – Competition models – Survival of one species – Co-existence – Alligators, doomsday and extinction • Mechanical Models – Nonlinear spring-mass system – Soft and hard springs – Energy conservation – Phase plane and scenes 548 Phase Plane Methods 10.1 Autonomous Planar Systems A set of two scalar differential equations of the form x′(t) = F (x(t), y(t)), y′(t) = G(x(t), y(t)).(1) is called a planar autonomous system. The term autonomous means self-governing, justified by the absence of the time variable t in the functions F (x, y), G(x, y). To obtain the vector form, let x(t) = ( x(t) y(t) ) , f(x, y) = ( F (x, y) G(x, y) ) and write (1) as the first order vector-matrix system x′(t) = f(x(t)).(2) It is assumed that F , G are continuously differentiable in some region D in the xy-plane. This assumption makes f continuously differentiable in D and guarantees that Picard’s existence-uniqueness theorem for ini- tial value problems applies to the initial value problem x′(t) = f(x(t)), x(0) = x0. Accordingly, to each x0 = (x0, y0) in D there corresponds a unique solution x(t) = (x(t), y(t)), represented as a planar curve in the xy-plane, which passes through x0 at t = 0. Such a planar curve is called a trajectory of the system and its param- eter interval is some maximal interval of existence T1 < t < T2, where T1 and T2 might be infinite. The graphic of a trajectory drawn as a parametric curve in the xy-plane is called a phase portrait and the xy-plane in which it is drawn is called the phase plane. Trajectories don’t cross. Autonomy of the planar system plus uniqueness of initial value problems implies that trajectories (x1(t), y1(t)) and (x2(t), y2(t)) cannot touch or cross. Hand-drawn phase portraits are accordingly limited: you cannot draw a solution trajectory that touches another solution curve! Theorem 1 (Identical trajectories) Assume that Picard’s existence-uniqueness theorem applies to initial value problems in D for the planar system (1). Let (x1(t), y1(t)) and (x2(t), y2(t)) be two trajectories of system (1). If times t1, t2 exist such that x1(t1) = x2(t2), y1(t1) = y2(t2),(3) then for the value c = t1−t2 the equations x1(t+c) = x2(t) and y1(t+c) = y2(t) are valid for all allowed values of t. This means that the two trajectories are on one and the same planar curve, or in the contrapositive, two different trajectories cannot touch or cross in the phase plane. 10.1 Autonomous Planar Systems 549 Proof: Define x(t) = x1(t+ c), y(t) = y1(t+ c). By the chain rule, (x(t), y(t)) is a solution of the planar system, because x′(t) = x′1(t+c) = F (x1(t+c), y1(t+ c)) = F (x(t), y(t)), and similarly for the second differential equation. Further, (3) implies x(t2) = x2(t2) and y(t2) = y2(t2), therefore Picard’s uniqueness theorem implies that x(t) = x2(t) and y(t) = y2(t) for all allowed values of t. The proof is complete. Equilibria. A trajectory that reduces to a point, or a constant so- lution x(t) = x0, y(t) = y0, is called an equilibrium solution. The equilibrium solutions or equilibria are found by solving the nonlinear equations F (x0, y0) = 0, G(x0, y0) = 0. Each such (x0, y0) in D is a trajectory whose graphic in the phase plane is a single point, called an equilibrium point. In applied literature, it may be called a critical point, stationary point or rest point. Theorem 1 has the following geometrical interpretation. Assuming uniqueness, no other trajectory (x(t), y(t)) in the phase plane can touch an equilibrium point (x0, y0). Equilibria (x0, y0) are often found from linear equations ax0 + by0 = e, cx0 + dy0 = f, which are solved by linear algebra methods. They constitute an impor- tant subclass of algebraic equations which can be solved symbolically. In this special case, symbolic solutions exist for the equilibria. It is interesting to report that in a practical sense the equilibria may be reported incorrectly, due to the limitations of computer software, even in this case when exact symbolic solutions are available. An example is x′ = x+ y, y′ = �y− � for small � > 0. The root of the problem is trans- lation of � to a machine constant, which is zero for small enough �. The result is that computer software detects infinitely many equilibria when in fact there is exactly one equilibrium point. This example suggests that symbolic computation be used by default. Practical methods for computing equilibria. There is no sup- porting theory to find equilibria for all choices of F and G. However, there is a rich library of special methods for solving nonlinear algebraic equations, including celebrated numerical methods such as Newton’s method and the bisection method. Computer algebra systems like maple and mathematica offer convenient codes to solve the equations, when possible, including symbolic solutions. Applied mathematics relies heavily on the dynamically expanding library of special methods, which grows monthly due to new mathematical discoveries. 550 Phase Plane Methods Population biology. Planar autonomous systems have been applied to two-species populations like two species of trout, who compete for food from the same supply, and foxes and rabbits, who compete in a predator-prey situation. Trout system. Certain equilibria are significant, because they repre- sent the population sizes for cohabitation. A point in the phase space that is not an equilibrium point corresponds to population sizes that cannot coexist, they must change with time. Some equilibria are con- sequently observable or average population sizes while non-equilibria correspond to snapshot population sizes that are subject to flux. Biolo- gists expect population sizes of such two-species competition models to undergo change until they reach approximately the observable values. Rabbit-fox system. This is an example of a predator-prey sys- tem, in which the expected observable population sizes oscillate periodi- cally over time. Certain equilibria for these systems represent ideal co- habitation. Biological experiments suggest that initial population sizes close to the equilibrium values cause populations to stay near the initial sizes, even though the populations oscillate periodically. Observations by biologists of large population variations seem to verify that individual populations oscillate periodically around the ideal cohabitation sizes. Trout system. Consider a population of two species of trout who compete for the same food supply. A typical autonomous planar system for the species x and y is x′(t) = x(−2x− y + 180), y′(t) = y(−x− 2y + 120). Equilibria. The equilibrium solutions for this system are (0, 0), (90, 0), (0, 60), (80, 20). Only nonnegative population sizes are physically significant. Units for the population sizes might be in hundreds or thousands of fish. The equi- librium (0, 0) corresponds to extinction of both species, while (0, 60) and (90, 0) correspond to the unusual situation of extinction for one species. The last equilibrium (80, 20) corresponds to co-existence of the two trout species with observable population sizes of 80 and 20. Phase Portraits A graphic which contains all the equilibria and typical trajectories or orbits of a planar autonomous system (1) is called a phase portrait. 10.1 Autonomous Planar Systems 551 While graphing equilibria is not a challenge, graphing typical trajectories seems to imply that we are going to solve the differential system. This is not the case. The plan is this: Equilibria Plot in the xy-plane all equilibria of (1). Window Select an x-range and a y-range for the graph window which includes all significant equilibria (Figure 3). Grid Plot a uniform grid of N grid points (N ≈ 50 for hand work) within the graph window, to populate the graph- ical white space (Figure 4). The isocline method might also be used to select grid points. Field Draw at each grid point a short tangent vector, a re- placement curve for a solution curve through a grid point on a small time interval (Figure 5). Orbits Draw additional threaded trajectories on long time inter- vals into the remaining white space of the graphic (Figure 6). This is guesswork, based upon tangents to threaded trajectories matching nearby field tangents drawn in the previous step. See Figure 1 for matching details. C y x b a Figure 1. Badly threaded orbit. Threaded solution curve C correctly matches its tangent to the tangent at nearby grid point a, but it fails to match at grid point b. Why does a tangent ~T1 have to match a tangent ~T2 at a nearby grid point (see Figure 2)? A tangent vector is given by ~T = x′(t) = f(x(t)). Hence ~T1 = f(u1), ~T2 = f(u2). However, u1 ≈ u2 in the graphic, hence by continuity of f it follows that ~T1 ≈ ~T2. y C ~T1 ~T2 u2 u1 x Figure 2. Tangent matching. Threaded solution curve C matches its tangent ~T1 at u1 to direction field tangent ~T2 at nearby grid point u2. It is important to emphasize that solution curves starting at a grid point are defined for a small t-interval about t = 0, and therefore their graphics extend on both sides of the grid point. We intend to shorten these curves until they appear to be straight line segments, graphically identical to the tangent line. Adding an arrowhead pointing in the tangent vector direction is usual. After all this construction, the shaft of the arrow is graphically identical to a solution curve segment. In fact, if 50 grid points 552 Phase Plane Methods were used, then 50 solution curve segments have already been entered onto the graphic! Threaded orbits are added to show what happens to solutions that are plotted on longer and longer t-intervals. Phase portrait illustration. The method outlined above will be applied to the illustration x′(t) = x(t) + y(t), y′(t) = 1− x2(t).(4) The equilibria are (1,−1) and (−1, 1). The graph window is selected as |x| ≤ 2, |y| ≤ 2, in order to include both equilibria. The uniform grid will be 11× 11, although for hand work 5× 5 is normal. Tangents at the grid points are short line segments which do not touch each another – they are graphically the same as short solution curves. x (−1, 1) (1,−1) −2 2 −2 2 y Figure 3. Equilibria (1,−1), (−1, 1) for (4) and graph window. The equilibria (x, y) are calculated from equations 0 = x+y, 0 = 1−x2. The graph window |x| ≤ 2, |y| ≤ 2 is invented ini- tially, then updated until Figure 5 reveals sufficiently rich field details. x y −2 2 −2 2 Figure 4. Equilibria (1,−1), (−1, 1) for (4) with 11× 11 uniform grid. The equilibria (squares) happen to cover up two grid points (circles). The size 11× 11 is invented to fill the white space in the graphic. y x 1−1 −1 1 Figure 5. Equilibria for (4). Equilibria (1,−1), (−1, 1) with an 11× 11 uniform grid and direction field. An arrow shaft at a grid point represents a solution curve over a small time interval. Threaded solution curves on long time in- tervals have tangents matching nearby ar- row shaft directions. 10.1 Autonomous Planar Systems 553 y x −1 1 1 −1 Figure 6. Equilibria for (4). Equilibria (1,−1), (−1, 1) with an 11× 11 uniform grid, threaded solution curves and arrow shafts from some direction field arrows. Threaded solution curve tangents are to match nearby direction field arrow shafts. See Figure 1 for how to match tangents. y x −1 1 1 −1 Figure 7. Phase portrait for (4). Shown are typical solution curves and an 11× 11 grid. The direction field has been removed for clar- ity. Threaded solution curves do not actu- ally cross, even though graphical resolution might suggest otherwise. Phase plot by computer. Illustrated here is how to make the phase plot in Figure 8 with the computer algebra system maple. y 1 −1 −1 1 x Figure 8. Phase portrait for (4). The graphic shows typical solution curves and a direction field. Produced in maple us- ing a 13× 13 grid. Before the computer work begins, the differential equation is defined and the equilibria are computed. Defaults supplied by maple allow an initial phase portrait to be plotted, from which the graph window is selected. The initial plot code: with(DEtools): des:=diff(x(t),t)=x(t)+y(t),diff(y(t),t)=1-x(t)^2: wind:=x=-2..2,y=-2..2: DEplot({des},[x(t),y(t)],t=-20..20,wind); The initial plot suggests which initial conditions near the equilibria should be selected in order to create typical orbits on the graphic. The final code with initial data and options: with(DEtools): des:=diff(x(t),t)=x(t)+y(t),diff(y(t),t)=1-x(t)^2: wind:=x=-2..2,y=-2..2: opts:=stepsize=0.05,dirgrid=[13,13], axes=none,thickness=3,arrows=small: 554 Phase Plane Methods ics:=[[x(0)=-1,y(0)=1.1],[x(0)=-1,y(0)=1.5], [x(0)=-1,y(0)=.9],[x(0)=-1,y(0)=.6],[x(0)=-1,y(0)=.3], [x(0)=1,y(0)=-0.9],[x(0)=1,y(0)=-0.6],[x(0)=1,y(0)=-0.6], [x(0)=1,y(0)=-0.3],[x(0)=1,y(0)=-1.6],[x(0)=1,y(0)=-1.3], [x(0)=1,y(0)=-1.1]]: DEplot({des},[x(t),y(t)],t=-20..20,wind,ics,opts); Direction field by computer. While maple can produce direction fields with its DEplot tool, the basic code that produces a field can be written with minimal outside support, therefore it applies to other programming languages. The code below applies to the example x′ = x+ y, y′ = 1− x2 treated above. # 2D phase plane direction field with uniform nxm grid. # Tangent length is 9/10 the grid box width W0. a:=-2:b:=2:c:=-2:d:=2:n:=11:m:=11: H:=evalf((b-a)/(n+1)):K:=evalf((d-c)/(m+1)):W0:=min(H,K): X:=t->a+H*(t):Y:=t->c+K*(t):P:=[]: F1:=(x,y)->evalf(x+y):F2:=(x,y)->evalf(1-x^2): for i from 1 to n do for j from 1 to m do x:=X(i):y:=Y(j):M1:=F1(x,y): M2:=F2(x,y): if (M1 =0 and M2 =0) then # no tangent, make a box h:=W0/5:V:=plottools[rectangle]([x-h,y+h],[x+h,y-h]): else h:=evalf(((1/2)*9*W0/10)/sqrt(M1^2+M2^2)): p1:=x-h*M1:p2:=y-h*M2:q1:=x+h*M1:q2:=y+h*M2: V:=plottools[arrow]([p1,p2],[q1,q2],0.2*W0,0.5*W0,1/4): fi: if (P = []) then P:=V: else P:=P,V: fi: od:od: plots[display](P); Maple libraries plots and plottools are used. The routine rectangle requires two arguments ul, lr, which are the upper left (ul) and lower right (lr) vertices of the rectangle. The routine arrow requires five ar- guments P , Q, sw, aw, af : the two points P , Q which define the arrow shaft and direction, plus the shaft width sw, arrowhead width aw and arrowhead length fraction af (fraction of the shaft length). These prim- itives plot a polygon from its vertices. The rectangle computes four vertices and the arrow seven vertices, which are then passed on to the PLOT primitive to make the graphic. 10.1 Autonomous Planar Systems 555 Stability Consider an autonomous system x′(t) = f(x(t)) with f continuously differentiable in a region D in the plane. Stable equilibrium. An equilibrium point x0 in D is said to be stable provided for each � > 0 there corresponds δ > 0 such that (a) given x(0) in D with ‖x(0)−x0‖ < δ, then the solution x(t) exists on 0 ≤ t <∞ and (b) ‖x(t)− x0‖ < � for 0 ≤ t <∞. Unstable equilibrium. The equilibrium point x0 is called unstable provided it is not stable, which means (a) or (b) fails (or both). Asymptotically stable equilibrium. The equilibrium point x0 is said to be asymptotically stable provided (a) and (b) hold (it is stable), and additionally (c) limt→∞ ‖x(t)− x0‖ = 0 for ‖x(0)− x0‖ < δ. Applied accounts of stability tend to emphasize item (b). Careful appli- cation of stability theory requires attention to (a), which is the question of extension of solutions of initial value problems to the half-axis. Basic extension theory for solutions of autonomous equations says that (a) will be satisfied provided (b) holds for those values of t for which x(t) is already defined. Stability verifications in mathematical and applied literature often implicitly use extension theory, in order to present details compactly. The reader is advised to adopt the same predisposition as researchers, who assume the reader to be equally clever as they. Physical stability. In the model x′(t) = f(x(t)), physical stability ad- dresses changes in f as well as changes in x(0). The meaning is this: physical parameters of the model, e.g., the mass m > 0, damping con- stant c > 0 and Hooke’s constant k > 0 in a damped spring-mass system x′ = y, y′ = − c m y − k m x, may undergo small changes without significantly affecting the solution. In physical stability, stable equilibria correspond to physically ob- served data whereas other solutions correspond to transient obser- vations that disappear over time. A typical instance is the trout system x′(t) = x(−2x− y + 180), y′(t) = y(−x− 2y + 120).(5) 556 Phase Plane Methods Physically observed data in the trout system (5) corresponds to the car- rying capacity, represented by the stable equilibrium point (80, 20), whereas transient observations are snapshot population sizes that are subject to change over time. The strange extinction equilibria (90, 0) and (0, 60) are unstable equilibria, which disagrees with intuition about zero births for less than two individuals, but agrees with graphical repre- sentations of the trout system in Figure 9. Changing f for a trout system adjusts the physical constants which describe the birth and death rates, whereas changing x(0) alters the initial population sizes of the two trout species. Figure 9. Phase portrait for the trout system (5). Shown are typical solution curves and a direction field. Equilibrium (80, 20) is asymptotically stable (a square). Equilibria (0, 0), (90, 0), (0, 60) are unstable (circles). 10.2 Planar Constant Linear Systems 557 10.2 Planar Constant Linear Systems A constant linear planar system is a set of two scalar differential equa- tions of the form x′(t) = ax(t) + by(t)), y′(t) = cx(t) + dy(t)),(1) where a, b, c and d are constants. In matrix form, x′(t) = Ax(t), A = ( a b c d ) , x(t) = ( x(t) y(t) ) . Solutions drawn in phase portraits don’t cross, because the system is autonomous. The origin is always an equilibrium solution. There can be infinitely many equilibria, found by solving Ax = 0 for the constant vector x, when A is not invertible. Recipe. A recipe exists for solving system (1), which parallels the recipe for second order constant coefficient equations Ay′′+By′+Cy = 0. The reader should view the result as an advertisement for learning Putzer’s spectral method, page 618, which is used to derive the formulas. Theorem 2 (Planar Constant Linear System Recipe) Consider the real planar system x′(t) = Ax(t). Let λ1, λ2 be the roots of the characteristic equation det(A−λI) = 0. The real general solution x(t) is given by the formula x(t) = Φ(t)x(0) where the 2× 2 real invertible matrix Φ(t) is defined as follows. Real λ1 6= λ2 Φ(t) = eλ1t I + e λ2t − eλ1t λ2 − λ1 (A− λ1I). Real λ1 = λ2 Φ(t) = eλ1t I + teλ1t (A− λ1I). Complex λ1 = λ2, λ1 = a+ bi, b > 0 Φ(t) = eat ( cos(bt) I + (A− aI)sin(bt) b ) . Continuity and redundancy. The formulas are continuous in the sense that limiting λ1 → λ2 in the first formula or b → 0 in the last formula produces the middle formula for real double roots. The first formula is also valid for complex conjugate roots λ1, λ2 = λ1 and it reduces to the third when λ1 = a + ib, therefore the third formula is technically redundant, but nevertheless useful, because it contains no complex numbers. 558 Phase Plane Methods Illustrations. Typical cases are represented by the following 2×2 ma- tric