--- PAGE 1 ---


--- PAGE 2 ---
DIFFERENTIAL EQUATIONS, DYNAMICAL SYSTEMS, AND AN INTRODUCTION TO CHAOS Morris W. Hirsch University of California, Berkeley Stephen Smale University of California, Berkeley Robert L. Devaney Boston University AMSTERDAM • BOSTON • HEIDELBERG • LONDON NEW YORK • OXFORD • PARIS • SAN DIEGO SAN FRANCISCO • SINGAPORE • SYDNEY • TOKYO Academic Press is an imprint of Elsevier

--- PAGE 3 ---
Academic Press is an imprint of Elsevier 225 Wyman Street, Waltham, MA 02451, USA The Boulevard, Langford Lane, Kidlington, Oxford, OX5 1GB, UK c © 2013 Elsevier Inc. All rights reserved. No part of this publication may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher. Details on how to seek permission, further information about the Publisher’s permissions policies and our arrangements with organizations such as the Copyright Clearance Center and the Copyright Licensing Agency, can be found at our website: www.elsevier.com/permissions . This book and the individual contributions contained in it are protected under copyright by the Publisher (other than as may be noted herein). Notices Knowledge and best practice in this field are constantly changing. As new research and experience broaden our understanding, changes in research methods, professional practices, or medical treatment may become necessary. Practitioners and researchers must always rely on their own experience and knowledge in evaluating and using any information, methods, compounds, or experiments described herein. In using such information or methods they should be mindful of their own safety and the safety of others, including parties for whom they have a professional responsibility. To the fullest extent of the law, neither the Publisher nor the authors, contributors, or editors, assume any liability for any injury and/or damage to persons or property as a matter of products liability, negligence or otherwise, or from any use or operation of any methods, products, instructions, or ideas contained in the material herein. Library of Congress Cataloging-in-Publication Data Hirsch, Morris W., 1933- Differential equations, dynamical systems, and an introduction to chaos. — 3rd ed. / Morris W. Hirsch, Stephen Smale, Robert L. Devaney. p. cm. ISBN 978-0-12-382010-5 (hardback) 1. Differential equations. 2. Algebras, Linear. 3. Chaotic behavior in systems. I. Smale, Stephen, 1930– II. Devaney, Robert L., 1948– III. Title. QA372.H67 2013 515’.35–dc23 2012002951 British Library Cataloguing-in-Publication Data A catalogue record for this book is available from the British Library. For information on all Academic Press publications visit our Website at www.elsevierdirect.com Printed in the United States 12 13 14 15 16 10 9 8 7 6 5 4 3 2 1

--- PAGE 4 ---
Contents Preface to the Third Edition ix Preface xi CHAPTER 1 First-Order Equations 1 1.1 The Simplest Example 1 1.2 The Logistic Population Model 4 1.3 Constant Harvesting and Bifurcations 7 1.4 Periodic Harvesting and Periodic Solutions 10 1.5 Computing the Poincar ´ e Map 11 1.6 Exploration: A Two-Parameter Family 15 CHAPTER 2 Planar Linear Systems 21 2.1 Second-Order Differential Equations 23 2.2 Planar Systems 24 2.3 Preliminaries from Algebra 26 2.4 Planar Linear Systems 29 2.5 Eigenvalues and Eigenvectors 30 2.6 Solving Linear Systems 33 2.7 The Linearity Principle 36 iii

--- PAGE 5 ---
iv Contents CHAPTER 3 Phase Portraits for Planar Systems 39 3.1 Real Distinct Eigenvalues 39 3.2 Complex Eigenvalues 44 3.3 Repeated Eigenvalues 47 3.4 Changing Coordinates 49 CHAPTER 4 Classification of Planar Systems 61 4.1 The Trace–Determinant Plane 61 4.2 Dynamical Classification 64 4.3 Exploration: A 3D Parameter Space 71 CHAPTER 5 Higher-Dimensional Linear Algebra 73 5.1 Preliminaries from Linear Algebra 73 5.2 Eigenvalues and Eigenvectors 82 5.3 Complex Eigenvalues 85 5.4 Bases and Subspaces 88 5.5 Repeated Eigenvalues 93 5.6 Genericity 100 CHAPTER 6 Higher-Dimensional Linear Systems 107 6.1 Distinct Eigenvalues 107 6.2 Harmonic Oscillators 114 6.3 Repeated Eigenvalues 120 6.4 The Exponential of a Matrix 123 6.5 Nonautonomous Linear Systems 130 CHAPTER 7 Nonlinear Systems 139 7.1 Dynamical Systems 140 7.2 The Existence and Uniqueness Theorem 142 7.3 Continuous Dependence of Solutions 147 7.4 The Variational Equation 149 7.5 Exploration: Numerical Methods 153 7.6 Exploration: Numerical Methods and Chaos 156 CHAPTER 8 Equilibria in Nonlinear Systems 159 8.1 Some Illustrative Examples 159 8.2 Nonlinear Sinks and Sources 165

--- PAGE 6 ---
Contents v 8.3 Saddles 168 8.4 Stability 174 8.5 Bifurcations 175 8.6 Exploration: Complex Vector Fields 182 CHAPTER 9 Global Nonlinear Techniques 187 9.1 Nullclines 187 9.2 Stability of Equilibria 192 9.3 Gradient Systems 202 9.4 Hamiltonian Systems 206 9.5 Exploration: The Pendulum with Constant Forcing 209 CHAPTER 10 Closed Orbits and Limit Sets 213 10.1 Limit Sets 213 10.2 Local Sections and Flow Boxes 216 10.3 The Poincar ´ e Map 218 10.4 Monotone Sequences in Planar Dynamical Systems 220 10.5 The Poincar ´ e–Bendixson Theorem 222 10.6 Applications of Poincar ´ e–Bendixson 225 10.7 Exploration: Chemical Reactions that Oscillate 228 CHAPTER 11 Applications in Biology 233 11.1 Infectious Diseases 233 11.2 Predator–Prey Systems 237 11.3 Competitive Species 244 11.4 Exploration: Competition and Harvesting 250 11.5 Exploration: Adding Zombies to the SIR Model 251 CHAPTER 12 Applications in Circuit Theory 257 12.1 An RLC Circuit 257 12.2 The Li ´ enard Equation 261 12.3 The van der Pol Equation 263 12.4 A Hopf Bifurcation 270 12.5 Exploration: Neurodynamics 272

--- PAGE 7 ---
vi Contents CHAPTER 13 Applications in Mechanics 277 13.1 Newton’s Second Law 277 13.2 Conservative Systems 280 13.3 Central Force Fields 282 13.4 The Newtonian Central Force System 285 13.5 Kepler’s First Law 290 13.6 The Two-Body Problem 293 13.7 Blowing Up the Singularity 294 13.8 Exploration: Other Central Force Problems 298 13.9 Exploration: Classical Limits of Quantum Mechanical Systems 299 13.10 Exploration: Motion of a Glider 301 CHAPTER 14 The Lorenz System 305 14.1 Introduction 306 14.2 Elementary Properties of the Lorenz System 308 14.3 The Lorenz Attractor 312 14.4 A Model for the Lorenz Attractor 316 14.5 The Chaotic Attractor 321 14.6 Exploration: The R ¨ ossler Attractor 326 CHAPTER 15 Discrete Dynamical Systems 329 15.1 Introduction 329 15.2 Bifurcations 334 15.3 The Discrete Logistic Model 337 15.4 Chaos 340 15.5 Symbolic Dynamics 344 15.6 The Shift Map 349 15.7 The Cantor Middle–Thirds Set 351 15.8 Exploration: Cubic Chaos 354 15.9 Exploration: The Orbit Diagram 355 CHAPTER 16 Homoclinic Phenomena 361 16.1 The Shilnikov System 361 16.2 The Horseshoe Map 368 16.3 The Double Scroll Attractor 375

--- PAGE 8 ---
Contents vii 16.4 Homoclinic Bifurcations 377 16.5 Exploration: The Chua Circuit 381 CHAPTER 17 Existence and Uniqueness Revisited 385 17.1 The Existence and Uniqueness Theorem 385 17.2 Proof of Existence and Uniqueness 387 17.3 Continuous Dependence on Initial Conditions 394 17.4 Extending Solutions 397 17.5 Nonautonomous Systems 401 17.6 Differentiability of the Flow 404 Bibliography 411 Index 415

--- PAGE 9 ---
This page intentionally left blank

--- PAGE 10 ---
Preface to Third Edition The main new features in this edition consist of a number of additional explo- rations together with numerous proof simplifications and revisions. The new explorations include a sojourn into numerical methods that highlights how these methods sometimes fail, which in turn provides an early glimpse of chaotic behavior. Another new exploration involves the previously treated SIR model of infectious diseases, only now considered with zombies as the infected population. A third new exploration involves explaining the motion of a glider. This edition has benefited from numerous helpful comments from a variety of readers. Special thanks are due to Jamil Gomes de Abreu, Eric Adams, Adam Leighton, Tiennyu Ma, Lluis Fernand Mello, Bogdan Przeradzki, Charles Pugh, Hal Smith, and Richard Venti for their valuable insights and corrections. ix

--- PAGE 11 ---
This page intentionally left blank

--- PAGE 12 ---
Preface In the thirty years since the publication of the first edition of this book, much has changed in the field of mathematics known as dynamical systems . In the early 1970s, we had very little access to high-speed computers and computer graphics. The word chaos had never been used in a mathematical setting. Most of the interest in the theory of differential equations and dynamical systems was confined to a relatively small group of mathematicians. Things have changed dramatically in the ensuing three decades. Comput- ers are everywhere, and software packages that can be used to approximate solutions of differential equations and view the results graphically are widely available. As a consequence, the analysis of nonlinear systems of differential equations is much more accessible than it once was. The discovery of com- plicated dynamical systems, such as the horseshoe map, homoclinic tangles, the Lorenz system, and their mathematical analysis, convinced scientists that simple stable motions such as equilibria or periodic solutions were not always the most important behavior of solutions of differential equations. The beauty and relative accessibility of these chaotic phenomena motivated scientists and engineers in many disciplines to look more carefully at the important differen- tial equations in their own fields. In many cases, they found chaotic behavior in these systems as well. Now dynamical systems phenomena appear in virtually every area of sci- ence, from the oscillating Belousov–Zhabotinsky reaction in chemistry to the chaotic Chua circuit in electrical engineering, from complicated motions in celestial mechanics to the bifurcations arising in ecological systems. xi

--- PAGE 13 ---
xii Preface As a consequence, the audience for a text on differential equations and dynamical systems is considerably larger and more diverse than it was in the 1970s. We have accordingly made several major structural changes to this book, including: 1. The treatment of linear algebra has been scaled back. We have dispensed with the generalities involved with abstract vector spaces and normed lin- ear spaces. We no longer include a complete proof of the reduction of all n × n matrices to canonical form. Rather, we deal primarily with matrices no larger than 4 × 4. 2. We have included a detailed discussion of the chaotic behavior in the Lorenz attractor, the Shil’nikov system, and the double-scroll attractor. 3. Many new applications are included; previous applications have been updated. 4. There are now several chapters dealing with discrete dynamical systems. 5. We deal primarily with systems that are C ∞ , thereby simplifying many of the hypotheses of theorems. This book consists of three main parts. The first deals with linear systems of differential equations together with some first-order nonlinear equations. The second is the main part of the text: here we concentrate on nonlinear systems, primarily two-dimensional, as well as applications of these systems in a wide variety of fields. Part three deals with higher dimensional systems. Here we emphasize the types of chaotic behavior that do not occur in planar systems, as well as the principal means of studying such behavior—the reduction to a discrete dynamical system. Writing a book for a diverse audience whose backgrounds vary greatly poses a significant challenge. We view this one as a text for a second course in differ- ential equations that is aimed not only at mathematicians, but also at scientists and engineers who are seeking to develop sufficient mathematical skills to analyze the types of differential equations that arise in their disciplines. Many who come to this book will have strong backgrounds in linear algebra and real analysis, but others will have less exposure to these fields. To make this text accessible to both groups, we begin with a fairly gentle introduction to low-dimensional systems of differential equations. Much of this will be a review for readers with a more thorough background in differential equations, so we intersperse some new topics throughout the early part of the book for those readers. For example, the first chapter deals with first-order equations. We begin it with a discussion of linear differential equations and the logistic popula- tion model, topics that should be familiar to anyone who has a rudimentary acquaintance with differential equations. Beyond this review, we discuss the logistic model with harvesting, both constant and periodic. This allows us to introduce bifurcations at an early stage as well as to describe Poincar´ e maps

--- PAGE 14 ---
Preface xiii and periodic solutions. These are topics that are not usually found in elemen- tary differential equations courses, yet they are accessible to anyone with a background in multivariable calculus. Of course, readers with a limited back- ground may wish to skip these specialized topics at first and concentrate on the more elementary material. Chapters 2 through 6 deal with linear systems of differential equations. Again we begin slowly, with Chapters 2 and 3 dealing only with planar sys- tems of differential equations and two-dimensional linear algebra. Chapters 5 and 6 introduce higher dimensional linear systems; however, our emphasis remains on three- and four-dimensional systems rather than completely gen- eral n -dimensional systems, even though many of the techniques we describe extend easily to higher dimensions. The core of the book lies in the second part. Here, we turn our atten- tion to nonlinear systems. Unlike linear systems, nonlinear systems present some serious theoretical difficulties such as existence and uniqueness of solu- tions, dependence of solutions on initial conditions and parameters, and the like. Rather than plunge immediately into these difficult theoretical questions, which require a solid background in real analysis, we simply state the impor- tant results in Chapter 7 and present a collection of examples that illustrate what these theorems say (and do not say). Proofs of all of the results are included in the final chapter of the book. In the first few chapters in the nonlinear part of the book, we introduce important techniques such as linearization near equilibria, nullcline analysis, stability properties, limit sets, and bifurcation theory. In the latter half of this part, we apply these ideas to a variety of systems that arise in biology, electrical engineering, mechanics, and other fields. Many of the chapters conclude with a section called “Exploration.” These sections consist of a series of questions and numerical investigations dealing with a particular topic or application relevant to the preceding material. In each Exploration we give a brief introduction to the topic at hand and provide references for further reading about this subject. But, we leave it to the reader to tackle the behavior of the resulting system using the material presented ear- lier. We often provide a series of introductory problems as well as hints as to how to proceed, but in many cases, a full analysis of the system could become a major research project. You will not find “answers in the back of the book” for the questions; in many cases, nobody knows the complete answer. (Except, of course, you!) The final part of the book is devoted to the complicated nonlinear behav- ior of higher dimensional systems known as chaotic behavior . We introduce these ideas via the famous Lorenz system of differential equations. As is often the case in dimensions three and higher, we reduce the problem of com- prehending the complicated behavior of this differential equation to that of understanding the dynamics of a discrete dynamical system or iterated

--- PAGE 15 ---
xiv Preface function. So we then take a detour into the world of discrete systems, dis- cussing along the way how symbolic dynamics can be used to describe certain chaotic systems completely. We then return to nonlinear differential equations to apply these techniques to other chaotic systems, including those that arise when homoclinic orbits are present. We maintain a website at math.bu.edu/hsd devoted to issues regarding this text. Look here for errata, suggestions, and other topics of interest to teachers and students of differential equations. We welcome any contributions from readers at this site.

--- PAGE 16 ---
1 First-Order Equations The purpose of this chapter is to develop some elementary yet important examples of first-order differential equations. The examples here illustrate some of the basic ideas in the theory of ordinary differential equations in the simplest possible setting. We anticipate that the first few examples will be familiar to readers who have taken an introductory course in differential equations. Later examples, such as the logistic model with harvesting, are included to give the reader a taste of certain topics (e.g., bifurcations, periodic solutions, and Poincar´ e maps) that we will return to often throughout this book. In later chapters, our treatment of these topics will be much more systematic. 1.1 The Simplest Example The differential equation familiar to all calculus students, dx dt = ax , is the simplest. It is also one of the most important. First, what does it mean? Here x = x ( t ) is an unknown real-valued function of a real variable t and dx / dt is its derivative (we will also use x ′ or x ′ ( t ) for the derivative). In addi- tion, a is a parameter; for each value of a we have a different differential Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00001-4 c © 2013 Elsevier Inc. All rights reserved. 1

--- PAGE 17 ---
2 Chapter 1 First-Order Equations equation. The equation tells us that for every value of t the relationship x ′ ( t ) = ax ( t ) is true. The solutions of this equation are obtained from calculus: if k is any real number, then the function x ( t ) = ke at is a solution since x ′ ( t ) = ake at = ax ( t ) . Moreover, there are no other solutions . To see this, let u ( t ) be any solution and compute the derivative of u ( t ) e − at : d dt ( u ( t ) e − at ) = u ′ ( t ) e − at + u ( t )( − ae − at ) = au ( t ) e − at − au ( t ) e − at = 0. Therefore, u ( t ) e − at is a constant k , so u ( t ) = ke at . This proves our assertion. Thus, we have found all possible solutions of this differential equation. We call the collection of all solutions of a differential equation the general solution of the equation. The constant k appearing in this solution is completely determined if the value u 0 of a solution at a single point t 0 is specified. Suppose that a function x ( t ) satisfying the differential equation is also required to satisfy x ( t 0 ) = u 0 . Then we must have ke at 0 = u 0 , so that k = u 0 e − at 0 . Thus, we have determined k and this equation therefore has a unique solution satisfying the specified initial condition x ( t 0 ) = u 0 . For simplicity, we often take t 0 = 0; then k = u 0 . There is no loss of generality in taking t 0 = 0, for if u ( t ) is a solution with u ( 0 ) = u 0 , then the function v ( t ) = u ( t − t 0 ) is a solution with v ( t 0 ) = u 0 . It is common to restate this in the form of an initial value problem : x ′ = ax , x ( 0 ) = u 0 . A solution x ( t ) of an initial value problem must not only solve the differential equation, but must also take on the prescribed initial value u 0 at t = 0. Note that there is a special solution of this differential equation when k = 0. This is the constant solution x ( t ) ≡ 0. A constant solution like this is called an equilibrium solution or equilibrium point for the equation. Equilibria are often among the most important solutions of differential equations. The constant a in the equation x ′ = ax can be considered as a parameter. If a changes, the equation changes and so do the solutions. Can we describe qualitatively the way the solutions change? The sign of a is crucial here: 1. If a > 0, lim t →∞ ke at equals ∞ when k > 0, and equals −∞ when k < 0

--- PAGE 18 ---
1.1 The Simplest Example 3 2. If a = 0, ke at = constant 3. If a < 0, lim t →∞ ke at = 0 The qualitative behavior of solutions is vividly illustrated by sketching the graphs of solutions as in Figure 1.1. Note that the behavior of solutions is quite different when a is positive and negative. When a > 0, all nonzero solutions tend away from the equilibrium point at 0 as t increases, whereas when a < 0, solutions tend toward the equi- librium point. We say that the equilibrium point is a source when nearby solutions tend away from it. The equilibrium point is a sink when nearby solutions tend toward it. We also describe solutions by drawing them on the phase line . As the solu- tion x ( t ) is a function of time, we may view x ( t ) as a particle moving along the real line. At the equilibrium point, the particle remains at rest (indicated by a solid dot), while any other solution moves up or down the x -axis, as indicated by the arrows in Figure 1.2. The equation x ′ = ax is stable in a certain sense if a 6 = 0. More precisely, if a is replaced by another constant b with a sign that is the same as a , then x t Figure 1.1 The solution graphs and phase line for x ′ = ax for a > 0. Each graph represents a particular solution. t x Figure 1.2 The solution graphs and phase line for x ′ = ax for a < 0.

--- PAGE 19 ---
4 Chapter 1 First-Order Equations the qualitative behavior of the solutions does not change. But if a = 0, the slightest change in a leads to a radical change in the behavior of solutions. We therefore say that we have a bifurcation at a = 0 in the one-parameter fam- ily of equations x ′ = ax . The concept of a bifurcation is one that will arise over and over in subsequent chapters of this book. 1.2 The Logistic Population Model The differential equation x ′ = ax can be considered as a simplistic model of population growth when a > 0. The quantity x ( t ) measures the population of some species at time t . The assumption that leads to the differential equa- tion is that the rate of growth of the population (namely, dx / dt ) is directly proportional to the size of the population. Of course, this naive assumption omits many circumstances that govern actual population growth, including, for example, the fact that actual populations cannot increase without bound. To take this restriction into account, we can make the following further assumptions about the population model: 1. If the population is small, the growth rate remains directly proportional to the size of the population. 2. If the population grows too large, however, the growth rate becomes negative. One differential equation that satisfies these assumptions is the logistic popu- lation growth model . This differential equation is x ′ = ax ( 1 − x N ) . Here a and N are positive parameters: a gives the rate of population growth when x is small, while N represents a sort of “ideal” population or “carrying capacity.” Note that if x is small, the differential equation is essentially x ′ = ax (since the term 1 − ( x / N ) ≈ 1), but if x > N , then x ′ < 0. Thus, this simple equation satisfies the preceding assumptions. We should add here that there are many other differential equations that correspond to these assumptions; our choice is perhaps the simplest. Without loss of generality, we will assume that N = 1. That is, we will choose units so that the carrying capacity is exactly 1 unit of population and x ( t ) therefore represents the fraction of the ideal population present at time t . Therefore, the logistic equation reduces to x ′ = f a ( x ) = ax ( 1 − x ) .

--- PAGE 20 ---
1.2 The Logistic Population Model 5 This is an example of a first-order, autonomous, nonlinear differential equation. It is first order since only the first derivative of x appears in the equation. It is autonomous since the right side of the equation depends on x alone, not on time t . Plus, it is nonlinear since f a ( x ) is a nonlinear func- tion of x . The previous example, x ′ = ax , is a first-order, autonomous, linear differential equation. The solution of the logistic differential equation is easily found by the tried- and-true calculus method of separation and integration: ∫ dx x ( 1 − x ) = ∫ a dt . The method of partial fractions allows us to rewrite the left integral as ∫ ( 1 x + 1 1 − x ) dx . Integrating both sides and then solving for x yields x ( t ) = Ke at 1 + Ke at , where K is the arbitrary constant that arises from integration. Evaluating this expression at t = 0 and solving for K gives K = x ( 0 ) 1 − x ( 0 ) . Using this, we may rewrite this solution as x ( 0 ) e at 1 − x ( 0 ) + x ( 0 ) e at . So this solution is valid for any initial population x ( 0 ) . When x ( 0 ) = 1, we have an equilibrium solution, since x ( t ) reduces to x ( t ) ≡ 1. Similarly, x ( t ) ≡ 0 is an equilibrium solution. Thus, we have “existence” of solutions for the logistic differential equation. We have no guarantee that these are all of the solutions of this equation at this stage; we will return to this issue when we discuss the existence and uniqueness problem for differential equations in Chapter 7. To get a qualitative feeling for the behavior of solutions, we sketch the slope field for this equation. The right side of the differential equation determines the slope of the graph of any solution at each time t . Thus, we may plot little slope lines in the tx –plane as in Figure 1.3, with the slope of the line at ( t , x )

--- PAGE 21 ---
6 Chapter 1 First-Order Equations x t x = 1 x = 0 Figure 1.3 Slope field, solution graphs, and phase line for x ′ = ax (1 − x ). 0.8 0.5 1 Figure 1.4 The graph of the function f ( x ) = ax (1 − x ) with a = 3.2. given by the quantity ax ( 1 − x ) . Our solutions must therefore have graphs that are tangent to this slope field everywhere. From these graphs, we see immediately that, in agreement with our assumptions, all solutions for which x ( 0 ) > 0 tend to the ideal population x ( t ) ≡ 1. For x ( 0 ) < 0, solutions tend to −∞ , although these solutions are irrelevant in the context of a population model. Note that we can also read this behavior from the graph of the function f a ( x ) = ax ( 1 − x ) . This graph, displayed in Figure 1.4, crosses the x -axis at the two points x = 0 and x = 1, so these represent our equilibrium points. When 0 < x < 1, we have f ( x ) > 0. Therefore, slopes are positive at any ( t , x ) with 0 < x < 1, so solutions must increase in this region. When x < 0 or x > 1, we have f ( x ) < 0, so solutions must decrease, as we see in both the solution graphs and the phase lines in Figure 1.3. We may read off the fact that x = 0 is a source and x = 1 is a sink from the graph of f in similar fashion. Near 0, we have f ( x ) > 0 if x > 0, so slopes are positive and solutions increase, but if x < 0, then f ( x ) < 0, so slopes are negative and solutions decrease. Thus, nearby solutions move away from 0, so 0 is a source. Similarly, 1 is a sink.

--- PAGE 22 ---
1.3 Constant Harvesting and Bifurcations 7 x x = 1 x = − 1 x = 0 t Figure 1.5 Slope field, solution graphs, and phase line for x ′ = x − x 3 . We may also determine this information analytically. We have f ′ a ( x ) = a − 2 ax so that f ′ a ( 0 ) = a > 0 and f ′ a ( 1 ) = − a < 0. Since f ′ a ( 0 ) > 0, slopes must increase through the value 0 as x passes through 0. That is, slopes are negative below x = 0 and positive above x = 0. Thus, solutions must tend away from x = 0. In similar fashion, f ′ a ( 1 ) < 0 forces solutions to tend toward x = 1, making this equilibrium point a sink. We will encounter many such “derivative tests” like this that predict the qualitative behavior near equilibria in subsequent chapters. Example. As a further illustration of these qualitative ideas, consider the differential equation x ′ = g ( x ) = x − x 3 . There are three equilibrium points at x = 0, ± 1. Since g ′ ( x ) = 1 − 3 x 2 , we have g ′ ( 0 ) = 1, so the equilibrium point 0 is a source. Also, g ′ ( ± 1 ) = − 2, so the equilibrium points at ± 1 are both sinks. Between these equilibria, the sign of the slope field of this equation is nonzero. From this information we can immediately display the phase line, which is shown in Figure 1.5.  1.3 Constant Harvesting and Bifurcations Now let’s modify the logistic model to take into account harvesting of the pop- ulation. Suppose that the population obeys the logistic assumptions with the parameter a = 1, but it is also harvested at the constant rate h . The differential equation becomes x ′ = x ( 1 − x ) − h , where h ≥ 0 is a new parameter.

--- PAGE 23 ---
8 Chapter 1 First-Order Equations x 0.5 h < 1/4 h = 1/4 h > 1/4 f h ( x ) Figure 1.6 The graphs of the function f h ( x ) = x (1 − x ) − h . Rather than solving this equation explicitly (which can be done—see Exer- cise 6 of this chapter), we use the graph of the function f h ( x ) = x ( 1 − x ) − h to “read off ” the qualitative behavior of solutions. In Figure 1.6 we display the graph of f h in three different cases: 0 < h < 1 / 4, h = 1 / 4, and h > 1 / 4. It is straightforward to check that f h has two roots when 0 ≤ h < 1 / 4, one root when h = 1 / 4, and no roots if h > 1 / 4, as illustrated in the graphs. As a consequence, the differential equation has two equilibrium points, x ` and x r , with 0 ≤ x ` < x r when 0 < h < 1 / 4. It is also easy to check that f ′ h ( x ` ) > 0 so that x ` is a source, and f ′ h ( x r ) < 0 so that x r is a sink. As h passes through h = 1 / 4, we encounter another example of a bifurca- tion. The two equilibria, x ` and x r , coalesce as h increases through 1 / 4 and then disappear when h > 1 / 4. Moreover, when h > 1 / 4, we have f h ( x ) < 0 for all x . Mathematically, this means that all solutions of the differential equation decrease to −∞ as time goes on. We record this visually in the bifurcation diagram . In Figure 1.7, we plot the parameter h horizontally. Over each h -value we plot the corresponding phase line. The curve in this picture represents the equilibrium points for each value of h . This gives another view of the sink and source merging into a single equilibrium point and then disappearing as h passes through 1 / 4. Ecologically, this bifurcation corresponds to a disaster for the species under study. For rates of harvesting 1 / 4 or lower, the population persists, provided the initial population is sufficiently large ( x ( 0 ) ≥ x ` ) . But a very small change in the rate of harvesting when h = 1 / 4 leads to a major change in the fate of the population: at any rate of harvesting h > 1 / 4, the species becomes extinct. This phenomenon highlights the importance of detecting bifurcations in families of differential equations—a procedure that we will encounter many times in later chapters. We should also mention that, despite the simplicity of

--- PAGE 24 ---
1.3 Constant Harvesting and Bifurcations 9 x 1/4 h Figure 1.7 The bifurcation diagram for f h ( x ) = x (1 − x ) − h . x x = a a Figure 1.8 The bifurcation diagram for x ′ = x 2 − ax . this population model, the prediction that small changes in harvesting rates can lead to disastrous changes in population has been observed many times in real situations on earth. Example. As another example of a bifurcation, consider the family of differ- ential equations x ′ = g a ( x ) = x 2 − ax = x ( x − a ) , which depends on a parameter a . The equilibrium points are given by x = 0 and x = a . We compute that g ′ a ( 0 ) = − a , so 0 is a sink if a > 0 and a source if a < 0. Similarly, g ′ a ( a ) = a , so x = a is a sink if a < 0 and a source if a > 0. We have a bifurcation at a = 0 since there is only one equilibrium point when a = 0. Moreover, the equilibrium point at 0 changes from a source to a sink as a increases through 0. Similarly, the equilibrium at x = a changes from a sink to a source as a passes through 0. The bifurcation diagram for this family is shown in Figure 1.8. 

--- PAGE 25 ---
10 Chapter 1 First-Order Equations 1.4 Periodic Harvesting and Periodic Solutions Now let’s change our assumptions on the logistic model to reflect the fact that harvesting does not always occur at a constant rate. For example, populations of many species of fish are harvested at a higher rate in warmer months than in colder months. So, we assume that the population is harvested at a periodic rate. One such model is then x ′ = f ( t , x ) = ax ( 1 − x ) − h ( 1 + sin ( 2 π t )) , where again a and h are positive parameters. Thus, the harvesting reaches a maximum rate − 2 h at time t = 1 4 + n where n is an integer (representing the year), and the harvesting reaches its minimum value 0 when t = 3 4 + n , exactly one half year later. Note that this differential equation now depends explicitly on time; this is an example of a nonautonomous differential equation. As in the autonomous case, a solution x ( t ) of this equation must satisfy x ′ ( t ) = f ( t , x ( t )) for all t . Also, this differential equation is no longer separable, so we cannot generate an analytic formula for its solution using the usual methods from calculus. Thus, we are forced to take a more qualitative approach (see Figure 1.9). To describe the fate of the population in this case, we first note that the right side of the differential equation is periodic with period 1 in the time variable; that is, f ( t + 1, x ) = f ( t , x ) . This fact simplifies the problem of find- ing solutions somewhat. Suppose that we know the solution of all initial value problems, not for all times but only for 0 ≤ t ≤ 1. Then in fact we know the solutions for all time . For example, suppose x 1 ( t ) is the solution that is defined for 0 ≤ t ≤ 1 and satisfies x 1 ( 0 ) = x 0 . Suppose that x 2 ( t ) is the solution that satisfies x 2 ( 0 ) = Figure 1.9 The slope field for f ( x ) = x (1 − x ) − h ( 1 + sin (2 π t )).

--- PAGE 26 ---
1.5 Computing the Poincar ´ e Map 11 x 1 ( 1 ) . Then we can extend the solution x 1 by defining x 1 ( t + 1 ) = x 2 ( t ) for 0 ≤ t ≤ 1. The extended function is a solution since we have x ′ 1 ( t + 1 ) = x ′ 2 ( t ) = f ( t , x 2 ( t )) = f ( t + 1, x 1 ( t + 1 )) . Thus, if we know the behavior of all solutions in the interval 0 ≤ t ≤ 1, then we can extrapolate in similar fashion to all time intervals and thereby know the behavior of solutions for all time. Second, suppose that we know the value at time t = 1 of the solution satisfy- ing any initial condition x ( 0 ) = x 0 . Then, to each such initial condition x 0 , we can associate the value x ( 1 ) of the solution x ( t ) that satisfies x ( 0 ) = x 0 . This gives us a function p ( x 0 ) = x ( 1 ) . If we compose this function with itself, we derive the value of the solution through x 0 at time 2; that is, p ( p ( x 0 )) = x ( 2 ) . If we compose this function with itself n times, then we can compute the value of the solution curve at time n and hence we know the fate of the solution curve. The function p is called a Poincar´ e map for this differential equation. Having such a function allows us to move from the realm of continuous dynami- cal systems (differential equations) to the often easier-to-understand realm of discrete dynamical systems (iterated functions). For example, suppose that we know that p ( x 0 ) = x 0 for some initial condition x 0 ; that is, x 0 is a fixed point for the function p . Then, from our previous observations, we know that x ( n ) = x 0 for each integer n . Moreover, for each time t with 0 < t < 1, we also have x ( t ) = x ( t + 1 ) and thus x ( t + n ) = x ( t ) for each integer n . That is, the solution satisfying the initial condition x ( 0 ) = x 0 is a periodic function of t with period 1. Such solutions are called periodic solutions of the differential equation. In Figure 1.10, we have displayed several solutions of the logistic equation with periodic harvesting. Note that the solution satisfying the initial condi- tion, x ( 0 ) = x 0 , is a periodic solution, and we have x 0 = p ( x 0 ) = p ( p ( x 0 )). . . . Similarly, the solution satisfying the initial condition, x ( 0 ) = ˆ x 0 , also appears to be a periodic solution, so we should have p ( ˆ x 0 ) = ˆ x 0 . Unfortunately, it is usually the case that computing a Poincar´ e map for a differential equation is impossible, but for the logistic equation with periodic harvesting we get lucky. 1.5 Computing the Poincar ´ e Map Before computing the Poincar´ e map for this equation, we need to introduce some important terminology. To emphasize the dependence of a solution on

--- PAGE 27 ---
12 Chapter 1 First-Order Equations t = 0 x 0 x (1) = p ( x 0 ) x (2) = p ( p ( x 0 )) x 0 t = 1 t = 2 ∧ Figure 1.10 The Poincar ´ e map for x ′ = 5 x (1 − x ) − 0.8(1 + sin (2 π t )). the initial value x 0 , we will denote the corresponding solution by φ( t , x 0 ) . This function, φ : R × R → R , is called the flow associated with the differential equation. If we hold the variable x 0 fixed, then the function t → φ( t , x 0 ) is just an alternative expression for the solution of the differential equation satisfying the initial condition x 0 . Sometimes we write this function as φ t ( x 0 ) . Example. For our first example, x ′ = ax , the flow is given by φ( t , x 0 ) = x 0 e at . For the logistic equation (without harvesting), the flow is φ( t , x 0 ) = x ( 0 ) e at 1 − x ( 0 ) + x ( 0 ) e at . Now we return to the logistic differential equation with periodic harvesting, x ′ = f ( t , x ) = ax ( 1 − x ) − h ( 1 + sin ( 2 π t )) . The solution that satisfies the initial condition, x ( 0 ) = x 0 , is given by t → φ( t , x 0 ) . Although we do not have a formula for this expression, we do

--- PAGE 28 ---
1.5 Computing the Poincar ´ e Map 13 know that, by the Fundamental Theorem of Calculus, this solution satisfies φ( t , x 0 ) = x 0 + t ∫ 0 f ( s , φ( s , x 0 )) ds since ∂φ ∂ t ( t , x 0 ) = f ( t , φ( t , x 0 )) and φ( 0, x 0 ) = x 0 . If we differentiate this solution with respect to x 0 , using the Chain Rule, we obtain: ∂φ ∂ x 0 ( t , x 0 ) = 1 + t ∫ 0 ∂ f ∂ x 0 ( s , φ( s , x 0 )) · ∂φ ∂ x 0 ( s , x 0 ) ds . Now let z ( t ) = ∂φ ∂ x 0 ( t , x 0 ) . Note that z ( 0 ) = ∂φ ∂ x 0 ( 0, x 0 ) = 1. Differentiating z with respect to t , we find z ′ ( t ) = ∂ f ∂ x 0 ( t , φ( t , x 0 )) · ∂φ ∂ x 0 ( t , x 0 ) = ∂ f ∂ x 0 ( t , φ( t , x 0 )) · z ( t ) . Again, we do not know φ( t , x 0 ) explicitly, but this equation does tell us that z ( t ) solves the differential equation z ′ ( t ) = ∂ f ∂ x 0 ( t , φ( t , x 0 )) z ( t )

--- PAGE 29 ---
14 Chapter 1 First-Order Equations with z ( 0 ) = 1. Consequently, via separation of variables, we may compute that the solution of this equation is z ( t ) = exp t ∫ 0 ∂ f ∂ x 0 ( s , φ( s , x 0 )) ds , and so we find ∂φ ∂ x 0 ( 1, x 0 ) = exp 1 ∫ 0 ∂ f ∂ x 0 ( s , φ( s , x 0 )) ds . Since p ( x 0 ) = φ( 1, x 0 ) , we have determined the derivative p ′ ( x 0 ) of the Poin- car´ e map; note that p ′ ( x 0 ) > 0. Therefore, p is an increasing function. Differentiating once more, we find p ′′ ( x 0 ) = p ′ ( x 0 )   1 ∫ 0 ∂ 2 f ∂ x 0 ∂ x 0 ( s , φ( s , x 0 )) · exp   s ∫ 0 ∂ f ∂ x 0 ( u , φ( u , x 0 )) du   ds   , which looks pretty intimidating. However, since f ( t , x 0 ) = ax 0 ( 1 − x 0 ) − h ( 1 + sin ( 2 π t )) , we have ∂ 2 f ∂ x 0 ∂ x 0 ≡ − 2 a . Thus, we know in addition that p ′′ ( x 0 ) < 0. Consequently, the graph of the Poincar´ e map is concave down. This implies that the graph of p can cross the diagonal line y = x at most two times; that is, there can be at most two values of x for which p ( x ) = x . Therefore, the Poincar´ e map has at most two fixed points. These fixed points yield periodic solutions of the original differential equation. These are solutions that satisfy x ( t + 1 ) = x ( t ) for all t . Another way to say this is that the flow, φ( t , x 0 ) , is a periodic function in t with period 1 when the initial condition x 0 is one of the fixed points. We saw these two solutions in the particular case when h = 0.8 in Figure 1.10. In Figure 1.11, we again see two solutions that appear to be periodic. Note that one of these appears to attract all nearby solutions, while the other appears to repel them. We’ll return to these concepts often and make them more precise later in the book.

--- PAGE 30 ---
1.6 Exploration: A Two-Parameter Family 15 1 1 2 3 4 5 Figure 1.11 Several solutions of x ′ = 5 x (1 − x ) − 0.8(1 + sin (2 π t )). Recall that the differential equation also depends on the harvesting param- eter h . For small values of h , there will be two fixed points such as shown in Figure 1.11. Differentiating f with respect to h , we find ∂ f ∂ h ( t , x 0 ) = − ( 1 + sin 2 π t ) . Thus, ∂ f /∂ h < 0 (except when t = 3 / 4). This implies that the slopes of the slope field lines at each point ( t , x 0 ) decrease as h increases. As a consequence, the values of the Poincar´ e map also decrease as h increases. There is a unique value h ∗ , therefore, for which the Poincar´ e map has exactly one fixed point. For h > h ∗ , there are no fixed points for p , so p ( x 0 ) < x 0 for all initial values. It then follows that the population again dies out.  1.6 Exploration: A Two-Parameter Family Consider the family of differential equations x ′ = f a , b ( x ) = ax − x 3 − b , which depends on two parameters, a and b . The goal of this exploration is to combine all of the ideas in this chapter to put together a complete picture of the two-dimensional parameter plane (the ab –plane) for this differential equation. Feel free to use a computer to experiment with this differential

--- PAGE 31 ---
16 Chapter 1 First-Order Equations equation at first, but then try the following to verify your observations rigorously: 1. First fix a = 1. Use the graph of f 1, b to construct the bifurcation diagram for this family of differential equations depending on b . 2. Repeat the previous question for a = 0 and then for a = − 1. 3. What does the bifurcation diagram look like for other values of a ? 4. Now fix b and use the graph to construct the bifurcation diagram for this family, which this time depends on a . 5. In the ab –plane, sketch the regions where the corresponding differential equation has different numbers of equilibrium points, including a sketch of the boundary between these regions. 6. Describe, using phase lines and the graph of f a , b ( x ) , the bifurcations that occur as the parameters pass through this boundary. 7. Describe in detail the bifurcations that occur at a = b = 0 as a and/or b vary. 8. Consider the differential equation x ′ = x − x 3 − b sin ( 2 π t ) , where | b | is small. What can you say about solutions of this equation? Are there any periodic solutions? 9. Experimentally, what happens as | b | increases? Do you observe any bifurcations? Explain what you observe. E X E R C I S E S 1. Find the general solution of the differential equation x ′ = ax + 3 where a is a parameter. What are the equilibrium points for this equation? For which values of a are the equilibria sinks? For which are they sources? 2. For each of the following differential equations, find all equilibrium solu- tions and determine whether they are sinks, sources, or neither. Also sketch the phase line. (a) x ′ = x 3 − 3 x (b) x ′ = x 4 − x 2 (c) x ′ = cos x (d) x ′ = sin 2 x (e) x ′ = | 1 − x 2 | 3. Each of the following families of differential equations depends on a parameter a . Sketch the corresponding bifurcation diagrams. (a) x ′ = x 2 − ax (b) x ′ = x 3 − ax (c) x ′ = x 3 − x + a

--- PAGE 32 ---
Exercises 17 x b f ( x ) Figure 1.12 Graph of the function f . 4. Consider the function f ( x ) with a graph that is displayed in Figure 1.12. (a) Sketch the phase line corresponding to the differential equation x ′ = f ( x ) . (b) Let g a ( x ) = f ( x ) + a . Sketch the bifurcation diagram corresponding to the family of differential equations x ′ = g a ( x ) . (c) Describe the different bifurcations that occur in this family. 5. Consider the family of differential equations x ′ = ax + sin x , where a is a parameter. (a) Sketch the phase line when a = 0. (b) Use the graphs of ax and sin x to determine the qualitative behavior of all of the bifurcations that occur as a increases from − 1 to 1. (c) Sketch the bifurcation diagram for this family of differential equations. 6. Find the general solution of the logistic differential equation with con- stant harvesting, x ′ = x ( 1 − x ) − h , for all values of the parameter h > 0. 7. Consider the nonautonomous differential equation x ′ = { x − 4 if t < 5, 2 − x if t ≥ 5. (a) Find a solution of this equation satisfying x ( 0 ) = 4. Describe the qualitative behavior of this solution.

--- PAGE 33 ---
18 Chapter 1 First-Order Equations (b) Find a solution of this equation satisfying x ( 0 ) = 3. Describe the qualitative behavior of this solution. (c) Describe the qualitative behavior of any solution of this system as t → ∞ . 8. Consider a first-order linear equation of the form x ′ = ax + f ( t ) , where a ∈ R . Let y ( t ) be any solution of this equation. Prove that the general solution is y ( t ) + c exp ( at ) where c ∈ R is arbitrary. 9. Consider a first-order, linear, nonautonomous equation of the form x ′ ( t ) = a ( t ) x . (a) Find a formula involving integrals for the solution of this system. (b) Prove that your formula gives the general solution of this system. 10. Consider the differential equation x ′ = x + cos t . (a) Find the general solution of this equation. (b) Prove that there is a unique periodic solution for this equation. (c) Compute the Poincar´ e map p : { t = 0 } → { t = 2 π } for this equation and use this to verify again that there is a unique periodic solution. 11. First-order differential equations need not have solutions that are defined for all time. (a) Find the general solution of the equation x ′ = x 2 . (b) Discuss the domains over which each solution is defined. (c) Give an example of a differential equation for which the solution satisfying x ( 0 ) = 0 is defined only for − 1 < t < 1. 12. First-order differential equations need not have unique solutions satisfy- ing a given initial condition. (a) Prove that there are infinitely many different solutions of the differ- ential equations x ′ = x 1 / 3 satisfying x ( 0 ) = 0. (b) Discuss the corresponding situation that occurs for x ′ = x / t , x ( 0 ) = x 0 . (c) Discuss the situation that occurs for x ′ = x / t 2 , x ( 0 ) = 0. 13. Let x ′ = f ( x ) be an autonomous first-order differential equation with an equilibrium point at x 0 . (a) Suppose f ′ ( x 0 ) = 0. What can you say about the behavior of solu- tions near x 0 ? Give examples. (b) Suppose f ′ ( x 0 ) = 0 and f ′′ ( x 0 ) 6 = 0. What can you say now? (c) Suppose f ′ ( x 0 ) = f ′′ ( x 0 ) = 0 but f ′′′ ( x 0 ) 6 = 0. What can you say now?

--- PAGE 34 ---
Exercises 19 14. Consider the first-order nonautonomous equation x ′ = p ( t ) x , where p ( t ) is differentiable and periodic with period T . Prove that all solutions of this equation are periodic with period T if and only if T ∫ 0 p ( s ) ds = 0. 15. Consider the differential equation x ′ = f ( t , x ) , where f ( t , x ) is continu- ously differentiable in t and x . Suppose that f ( t + T , x ) = f ( t , x ) for all t . Suppose there are constants p , q such that f ( t , p ) > 0, f ( t , q ) < 0 for all t . Prove that there is a periodic solution x ( t ) for this equation with p < x ( 0 ) < q . 16. Consider the differential equation x ′ = x 2 − 1 − cos ( t ) . What can be said about the existence of periodic solutions for this equation?

--- PAGE 35 ---
This page intentionally left blank

--- PAGE 36 ---
2 Planar Linear Systems In this chapter we begin the study of systems of differential equations . A system of differential equations is a collection of n interrelated differential equations of the form x ′ 1 = f 1 ( t , x 1 , x 2 , . . . , x n ) x ′ 2 = f 2 ( t , x 1 , x 2 , . . . , x n ) . . . x ′ n = f n ( t , x 1 , x 2 , . . . , x n ) . Here the functions f j are real-valued functions of the n + 1 variables x 1 , x 2 , . . . , x n , and t . Unless otherwise specified, we will always assume that the f j are C ∞ functions. This means that the partial derivatives of all orders of the f j exist and are continuous. To simplify notation, we will use vector notation: X =    x 1 . . . x n    . We often write the vector X as ( x 1 , . . . , x n ) to save space. Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00002-6 c © 2013 Elsevier Inc. All rights reserved. 21

--- PAGE 37 ---
22 Chapter 2 Planar Linear Systems Our system may then be written more concisely as X ′ = F ( t , X ) , where F ( t , X ) =    f 1 ( t , x 1 , . . . , x n ) . . . f n ( t , x 1 , . . . , x n )    . A solution of this system is therefore a function of the form X ( t ) = ( x 1 ( t ) , . . . , x n ( t )) that satisfies the equation, so that X ′ ( t ) = F ( t , X ( t )) , where X ′ ( t ) = ( x ′ 1 ( t ) , . . . , x ′ n ( t )) . Of course, at this stage, we have no guarantee that there is such a solution, but we will begin to discuss this complicated question in Section 2.7. The system of equations is called autonomous if none of the f j depends on t , so the system becomes X ′ = F ( X ) . For most of the rest of this book we will be concerned with autonomous systems. In analogy with first-order differential equations, a vector X 0 for which F ( X 0 ) = 0 is called an equilibrium point for the system. An equilibrium point corresponds to a constant solution X ( t ) ≡ X 0 of the system as before. Just to set some notation once and for all, we will always denote real variables by lowercase letters such as x , y , x 1 , x 2 , t , and so forth. Real-valued functions will also be written in lowercase such as f ( x , y ) or f 1 ( x 1 , . . . , x n , t ) . We will reserve capital letters for vectors, such as X = ( x 1 , . . . , x n ) , or for vector-valued functions such as F ( x , y ) = ( f ( x , y ) , g ( x , y )) or H ( x 1 , . . . , x n ) =    h 1 ( x 1 , . . . , x n ) . . . h n ( x 1 , . . . , x n )    . We will denote n -dimensional Euclidean space by R n , so that R n consists of all vectors of the form X = ( x 1 , . . . , x n ) .

--- PAGE 38 ---
2.1 Second-Order Differential Equations 23 2.1 Second-Order Differential Equations Many of the most important differential equations encountered in science and engineering are second-order differential equations. These are differential equations of the form x ′′ = f ( t , x , x ′ ) . Important examples of second-order equations include Newton’s equation, mx ′′ = f ( x ) , the equation for an RLC circuit in electrical engineering, LCx ′′ + RCx ′ + x = v ( t ) , and the mainstay of most elementary differential equations courses, the forced harmonic oscillator, mx ′′ + bx ′ + kx = f ( t ) . We discuss these and more complicated relatives of these equations at length as we go along. First, however, we note that these equations are a special sub- class of two-dimensional systems of differential equations that are defined by simply introducing a second variable y = x ′ . For example, consider a second-order constant coefficient equation of the form x ′′ + ax ′ + bx = 0. If we let y = x ′ , then we may rewrite this equation as a system of first-order equations: x ′ = y y ′ = − bx − ay . Any second-order equation may be handled similarly. Thus, for the remainder of this book, we will deal primarily with systems of equations.

--- PAGE 39 ---
24 Chapter 2 Planar Linear Systems 2.2 Planar Systems In this chapter we will deal with autonomous systems in R 2 , which we will write in the form x ′ = f ( x , y ) y ′ = g ( x , y ) , thus eliminating the annoying subscripts on the functions and variables. As before, we often use the abbreviated notation X ′ = F ( X ) , where X = ( x , y ) and F ( X ) = F ( x , y ) = ( f ( x , y ) , g ( x , y )) . In analogy with the slope fields of Chapter 1, we regard the right side of this equation as defining a vector field on R 2 . That is, we think of F ( x , y ) as representing a vector with x - and y -components that are f ( x , y ) and g ( x , y ) , respectively. We visualize this vector as being based at the point ( x , y ) . For example, the vector field associated with the system, x ′ = y y ′ = − x , is displayed in Figure 2.1. Note that, in this case, many of the vectors overlap, making the pattern difficult to visualize. For this reason, we always draw a direction field instead, which consists of scaled versions of the vectors. A solution of this system should now be thought of as a parametrized curve in the plane of the form ( x ( t ) , y ( t )) such that, for each t , the tangent vector at the point ( x ( t ) , y ( t )) is F ( x ( t ) , y ( t )) . That is, the solution curve ( x ( t ) , y ( t )) winds its way through the plane always tangent to the given vector F ( x ( t ) , y ( t )) based at ( x ( t ) , y ( t )) . Figure 2.1 Vector field, direction field, and several solutions for the system x ′ = y , y ′ = − x .

--- PAGE 40 ---
2.2 Planar Systems 25 Example. The curve ( x ( t ) y ( t ) ) = ( a sin t a cos t ) for any a ∈ R is a solution of the system x ′ = y y ′ = − x since x ′ ( t ) = a cos t = y ( t ) y ′ ( t ) = − a sin t = − x ( t ) , as required by the differential equation. These curves define circles of radius | a | in the plane, which are traversed in the clockwise direction as t increases. When a = 0, the solutions are the constant functions x ( t ) ≡ 0 ≡ y ( t ) .  Note that this example is equivalent to the second-order differential equa- tion x ′′ = − x by simply introducing the second variable y = x ′ . This is an example of a linear second-order differential equation, which, in more general form, may be written a ( t ) x ′′ + b ( t ) x ′ + c ( t ) x = f ( t ) . An important special case of this is the linear, constant coefficient equation ax ′′ + bx ′ + cx = f ( t ) , which we write as a system as x ′ = y y ′ = − c a x − b a y + f ( t ) a . An even more special case is the homogeneous equation in which f ( t ) ≡ 0. Example. One of the simplest yet most important second-order, linear, constant-coefficient differential equations is the equation for a harmonic oscil- lator . This equation models the motion of a mass attached to a spring. The spring is attached to a vertical wall and the mass is allowed to slide along a

--- PAGE 41 ---
26 Chapter 2 Planar Linear Systems horizontal track. We let x denote the displacement of the mass from its natu- ral resting place (with x > 0 if the spring is stretched and x < 0 if the spring is compressed). Therefore the velocity of the moving mass is x ′ ( t ) and the acceleration is x ′′ ( t ) . The spring exerts a restorative force proportional to x ( t ) . In addition, there is a frictional force proportional to x ′ ( t ) in the direction opposite to that of the motion. There are three parameters for this system: m denotes the mass of the oscillator, b ≥ 0 is the damping constant , and k > 0 is the spring constant . New- ton’s law states that the force acting on the oscillator is equal to mass times acceleration. Therefore the differential equation for the damped harmonic oscillator is mx ′′ + bx ′ + kx = 0. If b = 0, the oscillator is said to be undamped ; otherwise, we have a damped harmonic oscillator. This is an example of a second-order, linear, constant coefficient, homogeneous differential equation. As a system, the harmonic oscillator equation becomes x ′ = y y ′ = − k m x − b m y . More generally, the motion of the mass-spring system can be subjected to an external force (such as moving the vertical wall back and forth periodically). Such an external force usually depends only on time, not position, so we have a more general forced harmonic oscillator system, mx ′′ + bx ′ + kx = f ( t ) , where f ( t ) represents the external force. This is now a nonautonomous, second-order, linear equation.  2.3 Preliminaries from Algebra Before proceeding further with systems of differential equations, we need to recall some elementary facts regarding systems of algebraic equations. We will often encounter simultaneous equations of the form ax + by = α cx + dy = β ,

--- PAGE 42 ---
2.3 Preliminaries from Algebra 27 where the values of a , b , c , and d as well as α and β are given. In matrix form, we may write this equation as ( a b c d ) ( x y ) = ( α β ) . We denote by A the 2 × 2 coefficient matrix A = ( a b c d ) . This system of equations is easy to solve, assuming that there is a solution. There is a unique solution of these equations if and only if the determinant of A is nonzero. Recall that this determinant is the quantity given by det A = ad − bc . If det A = 0, we may or may not have solutions, but if there is a solution, then in fact there must be infinitely many solutions. In the special case where α = β = 0, we always have infinitely many solutions of A ( x y ) = ( 0 0 ) when det A = 0. Indeed, if the coefficient a of A is nonzero, we have x = − ( b / a ) y and so − c ( b a ) y + dy = 0. Thus, ( ad − bc ) y = 0. Since det A = 0, the solutions of the equation assume the form ( − ( b / a ) y , y ) , where y is arbitrary. This says that every solution lies on a straight line through the origin in the plane. A similar line of solutions occurs as long as at least one of the entries of A is nonzero. We will not worry too much about the case where all entries of A are 0; in fact, we will completely ignore it. Let V and W be vectors in the plane. We say that V and W are linearly independent if V and W do not lie along the same straight line through the origin. The vectors V and W are linearly dependent if either V or W is the zero vector or both lie on the same line through the origin. A geometric criterion for two vectors in the plane to be linearly indepen- dent is that they do not point in the same or opposite directions. That is, two

--- PAGE 43 ---
28 Chapter 2 Planar Linear Systems nonzero vectors V and W are linearly independent if and only if V 6 = λ W for any real number λ . An equivalent algebraic criterion for linear independence is given in the following proposition. Proposition. Suppose V = ( v 1 , v 2 ) and W = ( w 1 , w 2 ) . Then V and W are linearly independent if and only if det ( v 1 w 1 v 2 w 2 ) 6 = 0. For a proof, see Exercise 11 of this chapter.  Whenever we have a pair of linearly independent vectors V and W , we may always write any vector Z ∈ R 2 in a unique way as a linear combination of V and W . That is, we may always find a pair of real numbers α and β such that Z = α V + β W . Moreover, α and β are unique. To see this, suppose Z = ( z 1 , z 2 ) . Then we must solve the equations z 1 = α v 1 + β w 1 z 2 = α v 2 + β w 2 , where v i , w i , and z i are known. But this system has a unique solution (α , β) since det ( v 1 w 1 v 2 w 2 ) 6 = 0. The linearly independent vectors V and W are said to define a basis for R 2 . Any vector Z has unique “coordinates” relative to V and W . These coordinates are the pair (α , β) for which Z = α V + β W . Example. The unit vectors E 1 = ( 1, 0 ) and E 2 = ( 0, 1 ) obviously form a basis called the standard basis of R 2 . The coordinates of Z in this basis are just the “usual” Cartesian coordinates ( x , y ) of Z .  Example. The vectors V 1 = ( 1, 1 ) and V 2 = ( − 1, 1 ) also form a basis of R 2 . Relative to this basis, the coordinates of E 1 are ( 1 / 2, − 1 / 2 ) and those of

--- PAGE 44 ---
2.4 Planar Linear Systems 29 E 2 are ( 1 / 2, 1 / 2 ) because ( 1 0 ) = 1 2 ( 1 1 ) − 1 2 ( − 1 1 ) ( 0 1 ) = 1 2 ( 1 1 ) + 1 2 ( − 1 1 ) These “changes of coordinates” will become important later.  Example. The vectors V 1 = ( 1, 1 ) and V 2 = ( − 1, − 1 ) do not form a basis of R 2 since these vectors are collinear. Any linear combination of these vectors is of the form α V 1 + β V 2 = ( α − β α − β ) , which yields only vectors on the straight line through the origin, that is, V 1 and V 2 .  2.4 Planar Linear Systems We now further restrict our attention to the most important class of planar systems of differential equations, namely linear systems. In the autonomous case, these systems assume the simple form x ′ = ax + by y ′ = cx + dy , where a , b , c , and d are constants. We may abbreviate this system by using the coefficient matrix A , where A = ( a b c d ) . Then the linear system may be written as X ′ = AX .

--- PAGE 45 ---
30 Chapter 2 Planar Linear Systems Note that the origin is always an equilibrium point for a linear system. To find other equilibria, we must solve the linear system of algebraic equations ax + by = 0 cx + dy = 0. This system has a nonzero solution if and only if det A = 0. As we saw in the preceding, if det A = 0, then there is a straight line through the origin on which each point is an equilibrium. Thus we have Proposition. The planar linear system X ′ = AX has 1. A unique equilibrium point ( 0, 0 ) if det A 6 = 0 2. A straight line of equilibrium points if det A = 0 (and A is not the 0 -matrix)  2.5 Eigenvalues and Eigenvectors Now we turn to the question of finding nonequilibrium solutions of the linear system X ′ = AX . The key observation here is this: suppose V 0 is a nonzero vector for which we have AV 0 = λ V 0 , where λ ∈ R . Then the function X ( t ) = e λ t V 0 is a solution of the system. To see this, we compute X ′ ( t ) = λ e λ t V 0 = e λ t (λ V 0 ) = e λ t ( AV 0 ) = A ( e λ t V 0 ) = AX ( t ) , so X ( t ) does indeed solve the system of equations. Such a vector V 0 and its associated scalar have names as follows. Definition A nonzero vector V 0 is called an eigenvector of A if AV 0 = λ V 0 for some λ . The constant λ is called an eigenvalue of A .

--- PAGE 46 ---
2.5 Eigenvalues and Eigenvectors 31 As we observed, there is an important relationship between eigenvalues, eigenvectors, and solutions of systems of differential equations: Theorem. Suppose that V 0 is an eigenvector for the matrix A with associ- ated eigenvalue λ . Then the function X ( t ) = e λ t V 0 is a solution of the system X ′ = AX.  Note that if V 0 is an eigenvector for A with eigenvalue λ , then any nonzero scalar multiple of V 0 is also an eigenvector for A with eigenvalue λ . Indeed, if AV 0 = λ V 0 , then A (α V 0 ) = α AV 0 = λ(α V 0 ) for any nonzero constant α . Example. Consider A = ( 1 3 1 − 1 ) . Then A has an eigenvector V 0 = ( 3, 1 ) with associated eigenvalue λ = 2 since ( 1 3 1 − 1 ) ( 3 1 ) = ( 6 2 ) = 2 ( 3 1 ) . Similarly, V 1 = ( 1, − 1 ) is an eigenvector with associated eigenvalue λ = − 2.  Thus, for the system X ′ = ( 1 3 1 − 1 ) X we now know three solutions: the equilibrium solution at the origin together with X 1 ( t ) = e 2 t ( 3 1 ) and X 2 ( t ) = e − 2 t ( 1 − 1 ) . We will see that we can use these solutions to generate all solutions of this sys- tem in a moment, but first we address the question of how to find eigenvectors and eigenvalues.

--- PAGE 47 ---
32 Chapter 2 Planar Linear Systems To produce an eigenvector V = ( x , y ) , we must find a nonzero solution ( x , y ) of the equation A ( x y ) = λ ( x y ) . Note that there are three unknowns in this system of equations: the two components of V as well as λ . Let I denote the 2 × 2 identity matrix I = ( 1 0 0 1 ) . Then we may rewrite the equation in the form ( A − λ I ) V = 0, where 0 denotes the vector ( 0, 0 ) . Now A − λ I is just a 2 × 2 matrix (having entries involving the variable λ ), so this linear system of equations has nonzero solutions if and only if det ( A − λ I ) = 0, as we saw previously. But this equation is just a quadratic equation in λ , and so its roots are easy to find. This equation will appear over and over in the sequel; it is called the characteristic equation . As a function of λ , we call det ( A − λ I ) the characteristic polynomial . Thus the strategy to generate eigenvectors is first to find the roots of the characteristic equation. This yields the eigenvalues. Then we use each of these eigenvalues to generate in turn an associated eigenvector. Example. We return to the matrix A = ( 1 3 1 − 1 ) . We have A − λ I = ( 1 − λ 3 1 − 1 − λ ) . So the characteristic equation is det ( A − λ I ) = ( 1 − λ)( − 1 − λ) − 3 = 0. Simplifying, we find λ 2 − 4 = 0,

--- PAGE 48 ---
2.6 Solving Linear Systems 33 which yields the two eigenvalues λ = ± 2. Then, for λ = 2, we next solve the equation ( A − 2 I ) ( x y ) = ( 0 0 ) . In component form, this reduces to the system of equations ( 1 − 2 ) x + 3 y = 0 x + ( − 1 − 2 ) y = 0, or − x + 3 y = 0, as these equations are redundant. Thus any vector of the form ( 3 y , y ) with y 6 = 0 is an eigenvector associated with λ = 2. In similar fashion, any vector of the form ( y , − y ) with y 6 = 0 is an eigenvector associated with λ = − 2.  Of course, the astute reader will notice that there is more to the story of eigenvalues, eigenvectors, and solutions of differential equations than what we have described previously. For example, the roots of the characteristic equa- tion may be complex or they may be repeated real numbers. We will handle all of these cases shortly, but first we return to the problem of solving linear systems. 2.6 Solving Linear Systems As we saw in the example in the previous section, if we find two real roots λ 1 and λ 2 (with λ 1 6 = λ 2 ) of the characteristic equation, then we may generate a pair of solutions of the system of differential equations of the form X i ( t ) = e λ i t V i , where V i is the eigenvector associated with λ i . Note that each of these solutions is a straight-line solution . Indeed, we have X i ( 0 ) = V i , which is a nonzero point in the plane. For each t , e λ i t V i is a scalar multiple of V i and so lies on the straight ray emanating from the origin and passing through V i . Note that, if λ i > 0, then lim t →∞ | X i ( t ) | = ∞ and lim t →−∞ X i ( t ) = ( 0, 0 ) . The magnitude of the solution X i ( t ) increases monotonically to ∞ along the ray through V i as t increases, and X i ( t ) tends to the origin along this ray in

--- PAGE 49 ---
34 Chapter 2 Planar Linear Systems backward time. The exact opposite situation occurs if λ i < 0, whereas, if λ i = 0, the solution X i ( t ) is the constant solution X i ( t ) = V i for all t . So how do we find all solutions of the system given this pair of special solu- tions? The answer is now easy and important. Suppose we have two distinct real eigenvalues λ 1 and λ 2 with eigenvectors V 1 and V 2 . Then V 1 and V 2 are linearly independent, as is easily checked (see Exercise 14 of this chapter). Thus V 1 and V 2 form a basis of R 2 , so, given any point Z 0 ∈ R 2 , we may find a unique pair of real numbers α and β for which α V 1 + β V 2 = Z 0 . Now consider the function Z ( t ) = α X 1 ( t ) + β X 2 ( t ) , where the X i ( t ) are the preceding straight-line solutions. We claim that Z ( t ) is a solution of X ′ = AX . To see this we compute Z ′ ( t ) = α X ′ 1 ( t ) + β X ′ 2 ( t ) = α AX 1 ( t ) + β AX 2 ( t ) = A (α X 1 ( t ) + β X 2 ( t )) . This last step follows from the linearity of matrix multiplication (see Exercise 13 of this chapter). Thus, we have shown that Z ′ ( t ) = AZ ( t ) , so Z ( t ) is a solu- tion. Moreover, Z ( t ) is a solution that satisfies Z ( 0 ) = Z 0 . Finally, we claim that Z ( t ) is the unique solution of X ′ = AX that satisfies Z ( 0 ) = Z 0 . Just as in Chapter 1, we suppose that Y ( t ) is another such solution with Y ( 0 ) = Z 0 . Then we may write Y ( t ) = ζ ( t ) V 1 + μ( t ) V 2 , with ζ( 0 ) = α , μ( 0 ) = β . Thus, AY ( t ) = Y ′ ( t ) = ζ ′ ( t ) V 1 + μ ′ ( t ) V 2 . But AY ( t ) = ζ ( t ) AV 1 + μ( t ) AV 2 = λ 1 ζ ( t ) V 1 + λ 2 μ( t ) V 2 . Therefore, we have ζ ′ ( t ) = λ 1 ζ ( t ) μ ′ ( t ) = λ 2 μ( t ) ,

--- PAGE 50 ---
2.6 Solving Linear Systems 35 with ζ( 0 ) = α , μ( 0 ) = β . As we saw in Chapter 1, it follows that ζ ( t ) = α e λ 1 t , μ( t ) = β e λ 2 t , so that Y ( t ) is indeed equal to Z ( t ) . As a consequence, we have now found the unique solution to the system X ′ = AX that satisfies X ( 0 ) = Z 0 for any Z 0 ∈ R 2 . The collection of all such solutions is called the general solution of X ′ = AX . That is, the general solution is the collection of solutions of X ′ = AX that features a unique solution of the initial value problem X ( 0 ) = Z 0 for each Z 0 ∈ R 2 . We therefore have shown the theorem that follows. Theorem. Suppose A has a pair of real eigenvalues λ 1 6 = λ 2 and associated eigenvectors V 1 and V 2 . Then the general solution of the linear system X ′ = AX is given by X ( t ) = α e λ 1 t V 1 + β e λ 2 t V 2 .  Example. Consider the second-order differential equation: x ′′ + 3 x ′ + 2 x = 0. This is a specific case of the damped harmonic oscillator discussed earlier, where the mass is 1, the spring constant is 2, and the damping constant is 3. As a system, this equation may be rewritten: X ′ = ( 0 1 − 2 − 3 ) X = AX . The characteristic equation is λ 2 + 3 λ + 2 = (λ + 2 )(λ + 1 ) = 0, so the system has eigenvalues − 1 and − 2. The eigenvector corresponding to the eigenvalue − 1 is given by solving the equation: ( A + I ) ( x y ) = ( 0 0 ) . In component form this equation becomes x + y = 0 − 2 x − 2 y = 0.

--- PAGE 51 ---
36 Chapter 2 Planar Linear Systems Thus, one eigenvector associated with the eigenvalue − 1 is ( 1, − 1 ) . In similar fashion we compute that an eigenvector associated with the eigenvalue − 2 is ( 1, − 2 ) . Note that these two eigenvectors are linearly independent. Therefore, by the previous theorem, the general solution of this system is X ( t ) = α e − t ( 1 − 1 ) + β e − 2 t ( 1 − 2 ) . That is, the position of the mass is given by the first component of the solution, x ( t ) = α e − t + β e − 2 t , and the velocity is given by the second component, y ( t ) = x ′ ( t ) = − α e − t − 2 β e − 2 t .  2.7 The Linearity Principle The theorem discussed in the previous section is a very special case of the fundamental theorem for n -dimensional linear systems, which we shall prove in Chapter 6, Section 6.1, “Distinct Eigenvalues.” For the two-dimensional version of this result, note that if X ′ = AX is a planar linear system for which Y 1 ( t ) and Y 2 ( t ) are both solutions, then, just as before, the function α Y 1 ( t ) + β Y 2 ( t ) is also a solution of this system. We do not need real and distinct eigenvalues to prove this. This fact is known as the Linearity Principle. More important, if the initial conditions Y 1 ( 0 ) and Y 2 ( 0 ) are linearly inde- pendent vectors, then these vectors form a basis of R 2 . Thus, given any vector X 0 ∈ R 2 , we may determine constants α and β such that X 0 = α Y 1 ( 0 ) + β Y 2 ( 0 ) . Then the Linearity Principle tells us that the solution X ( t ) satisfying the initial condition X ( 0 ) = X 0 is given by X ( t ) = α Y 1 ( t ) + β Y 2 ( t ) . We have therefore produced a solution of the system that solves any given initial value problem. The Existence and Uniqueness Theorem for linear systems in Chap- ter 6 will show that this solution is also unique. This important result may then be summarized: Theorem. Let X ′ = AX be a planar system. Suppose that Y 1 ( t ) and Y 2 ( t ) are solutions of this system, and that the vectors Y 1 ( 0 ) and Y 2 ( 0 ) are linearly independent. Then X ( t ) = α Y 1 ( t ) + β Y 2 ( t ) is the unique solution of this system that satisfies X ( 0 ) = α Y 1 ( 0 ) + β Y 2 ( 0 ) . 

--- PAGE 52 ---
Exercises 37 E X E R C I S E S 1. Find the eigenvalues and eigenvectors of each of the following 2 × 2 matrices: ( a ) ( 3 1 1 3 ) ( b ) ( 2 1 1 1 ) ( c ) ( a b 0 c ) ( d ) ( 1 3 √ 2 3 √ 2 ) 2. Find the general solution of each of the following linear systems: ( a ) X ′ = ( 1 2 0 3 ) X ( b ) X ′ = ( 1 2 3 6 ) X ( c ) X ′ = ( 1 2 1 0 ) X ( d ) X ′ = ( 1 2 3 − 3 ) X 3. In Figure 2.2, you see four direction fields. Match each of these direction fields with one of the systems in the previous exercise. 4. Find the general solution of the system X ′ = ( a b c a ) X , where bc > 0. 1. 2. 3. 4. Figure 2.2 Match these direction fields with the systems in Exercise 2.

--- PAGE 53 ---
38 Chapter 2 Planar Linear Systems 5. Find the general solution of the system X ′ = ( 0 0 0 0 ) X . 6. For the harmonic oscillator system x ′′ + bx ′ + kx = 0, find all values of b and k for which this system has real, distinct eigenvalues. Find the general solution of this system in these cases. Find the solution of the system that satisfies the initial condition ( 0, 1 ) . Describe the motion of the mass in this particular case. 7. Consider the 2 × 2 matrix A = ( a 1 0 1 ) . Find the value a 0 of the parameter a for which A has repeated real eigenvalues. What happens to the eigenvectors of this matrix as a approaches a 0 ? 8. Describe all possible 2 × 2 matrices with eigenvalues of 0 and 1. 9. Give an example of a linear system for which ( e − t , α) is a solution for every constant α . Sketch the direction field for this system. What is the general solution of this system? 10. Give an example of a system of differential equations for which ( t , 1 ) is a solution. Sketch the direction field for this system. What is the general solution of this system? 11. Prove that two vectors V = ( v 1 , v 2 ) and W = ( w 1 , w 2 ) are linearly inde- pendent if and only if det ( v 1 w 1 v 2 w 2 ) 6 = 0. 12. Prove that if λ , μ are real eigenvalues of a 2 × 2 matrix, then any nonzero column of the matrix A − λ I is an eigenvector for μ . 13. Let A be a 2 × 2 matrix and let V 1 and V 2 vectors in R 2 . Prove that A (α V 1 + β V 2 ) = α AV 1 + β AV 2 . 14. Prove that the eigenvectors of a 2 × 2 matrix corresponding to distinct real eigenvalues are always linearly independent.

--- PAGE 54 ---
3 Phase Portraits for Planar Systems Given the Linearity Principle from the previous chapter, we may now com- pute the general solution of any planar system. There is a seemingly endless number of distinct cases, but we will see that these represent in the simplest possible form nearly all of the types of solutions we will encounter in the higher-dimensional case. 3.1 Real Distinct Eigenvalues Consider X ′ = AX and suppose that A has two real eigenvalues λ 1 < λ 2 . Assuming for the moment that λ i 6 = 0, there are three cases to consider: 1. λ 1 < 0 < λ 2 2. λ 1 < λ 2 < 0 3. 0 < λ 1 < λ 2 We give a specific example of each case; any system that falls into any one of these three categories may be handled similarly, as we show later. Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00003-8 c © 2013 Elsevier Inc. All rights reserved. 39

--- PAGE 55 ---
40 Chapter 3 Phase Portraits for Planar Systems Example. (Saddle) First consider the simple system X ′ = AX , where A = ( λ 1 0 0 λ 2 ) with λ 1 < 0 < λ 2 . This can be solved immediately since the system decouples into two unrelated first-order equations: x ′ = λ 1 x y ′ = λ 2 y . We already know how to solve these equations, but, having in mind what comes later, let’s find the eigenvalues and eigenvectors. The characteristic equation is (λ − λ 1 )(λ − λ 2 ) = 0, so λ 1 and λ 2 are the eigenvalues. An eigenvector corresponding to λ 1 is ( 1, 0 ) and to λ 2 is ( 0, 1 ) . Thus, we find the general solution X ( t ) = α e λ 1 t ( 1 0 ) + β e λ 2 t ( 0 1 ) . Since λ 1 < 0, the straight-line solutions of the form α e λ 1 t ( 1, 0 ) lie on the x -axis and tend to ( 0, 0 ) as t → ∞ . This axis is called the stable line . Since λ 2 > 0, the solutions β e λ 2 t ( 0, 1 ) lie on the y -axis and tend away from ( 0, 0 ) as t → ∞ ; this axis is the unstable line . All other solutions (with α , β 6 = 0) tend to ∞ in the direction of the unstable line, as t → ∞ , since X ( t ) comes closer and closer to ( 0, β e λ 2 t ) as t increases. In backward time, these solutions tend to ∞ in the direction of the stable line.  In Figure 3.1 we have plotted the phase portrait of this system. The phase portrait is a picture of a collection of representative solution curves of the system in R 2 , which we call the phase plane . The equilibrium point of a system of this type (eigenvalues satisfying λ 1 < 0 < λ 2 ) is called a saddle . For a slightly more complicated example of this type, consider X ′ = AX , where A = ( 1 3 1 − 1 ) . As we saw in Chapter 2, the eigenvalues of A are ± 2. The eigenvector associ- ated with λ = 2 is the vector ( 3, 1 ) ; the eigenvector associated with λ = − 2 is

--- PAGE 56 ---
3.1 Real Distinct Eigenvalues 41 Figure 3.1 Saddle phase portrait for x ′ = − x , y ′ = y . ( 1, − 1 ) . Thus, we have an unstable line that contains straight-line solutions of the form X 1 ( t ) = α e 2 t ( 3 1 ) , each of which tends away from the origin as t → ∞ . The stable line contains the straight-line solutions X 2 ( t ) = β e − 2 t ( 1 − 1 ) , which tend toward the origin as t → ∞ . By the Linearity Principle, any other solution assumes the form X ( t ) = α e 2 t ( 3 1 ) + β e − 2 t ( 1 − 1 ) for some α , β . Note that, if α 6 = 0, as t → ∞ , we have X ( t ) ∼ α e 2 t ( 3 1 ) = X 1 ( t ) , whereas, if β 6 = 0, as t → −∞ , X ( t ) ∼ β e − 2 t ( 1 − 1 ) = X 2 ( t ) . Thus, as time increases, the typical solution approaches X 1 ( t ) while, as time decreases, this solution tends toward X 2 ( t ) , just as in the previous case. Figure 3.2 displays this phase portrait.

--- PAGE 57 ---
42 Chapter 3 Phase Portraits for Planar Systems Figure 3.2 Saddle phase portrait for x ′ = x + 3 y , y ′ = x − y . In the general case where A has a positive and negative eigenvalue, we always find a similar stable and unstable line on which solutions tend toward or away from the origin. All other solutions approach the unstable line as t → ∞ , and tend toward the stable line as t → −∞ . Example. (Sink) Now consider the case X ′ = AX where A = ( λ 1 0 0 λ 2 ) but λ 1 < λ 2 < 0. As before, we find two straight-line solutions and then the general solution X ( t ) = α e λ 1 t ( 1 0 ) + β e λ 2 t ( 0 1 ) . Unlike the saddle case, now all solutions tend to ( 0, 0 ) as t → ∞ . The question is this: How do they approach the origin? To answer this, we compute the slope dy / dx of a solution with β 6 = 0. We write x ( t ) = α e λ 1 t y ( t ) = β e λ 2 t and compute dy dx = dy / dt dx / dt = λ 2 β e λ 2 t λ 1 α e λ 1 t = λ 2 β λ 1 α e (λ 2 − λ 1 ) t . Since λ 2 − λ 1 > 0, it follows that these slopes approach ±∞ (provided β 6 = 0). Thus these solutions tend to the origin tangentially to the y -axis. 

--- PAGE 58 ---
3.1 Real Distinct Eigenvalues 43 (a) (b) Figure 3.3 Phase portraits for a sink and a source. Since λ 1 < λ 2 < 0, we call λ 1 the stronger eigenvalue and λ 2 the weaker eigenvalue. The reason for this in this particular case is that the x -coordinates of solutions tend to 0 much more quickly than the y -coordinates. This accounts for why solutions (except those on the line corresponding to λ 1 - eigenvector) tend to “hug” the straight-line solution corresponding to the weaker eigenvalue as they approach the origin. The phase portrait for this system is displayed in Figure 3.3a. In this case the equilibrium point is called a sink . More generally, if the system has eigenvalues λ 1 < λ 2 < 0 with eigenvectors ( u 1 , u 2 ) and ( v 1 , v 2 ) respectively, then the general solution is α e λ 1 t ( u 1 u 2 ) + β e λ 2 t ( v 1 v 2 ) . The slope of this solution is given by dy dx = λ 1 α e λ 1 t u 2 + λ 2 β e λ 2 t v 2 λ 1 α e λ 1 t u 1 + λ 2 β e λ 2 t v 1 = ( λ 1 α e λ 1 t u 2 + λ 2 β e λ 2 t v 2 λ 1 α e λ 1 t u 1 + λ 2 β e λ 2 t v 1 ) e − λ 2 t e − λ 2 t = λ 1 α e (λ 1 − λ 2 ) t u 2 + λ 2 β v 2 λ 1 α e (λ 1 − λ 2 ) t u 1 + λ 2 β v 1 , which tends to the slope v 2 / v 1 of the λ 2 -eigenvector, unless we have β = 0. If β = 0, our solution is the straight-line solution corresponding to the eigen- value λ 1 . Thus, in this case as well, all solutions (except those on the straight line corresponding to the stronger eigenvalue) tend to the origin tangentially to the straight-line solution corresponding to the weaker eigenvalue.

--- PAGE 59 ---
44 Chapter 3 Phase Portraits for Planar Systems Example. (Source) When the matrix A = ( λ 1 0 0 λ 2 ) satisfies 0 < λ 2 < λ 1 , our vector field may be regarded as the negative of the previous example. The general solution and phase portrait remain the same, except that all solutions now tend away from ( 0, 0 ) along the same paths. See Figure 3.3b.  One may argue that we are presenting examples here that are much too simple. Although this is true, we will soon see that any system of differential equations with a matrix that has real distinct eigenvalues can be manipulated into the preceding special forms by changing coordinates. Finally, a special case occurs if one of the eigenvalues is equal to 0. As we have seen, there is a straight line of equilibrium points in this case. If the other eigenvalue λ is nonzero, then the sign of λ determines whether the other solu- tions tend toward or away from these equilibria (see Exercises 10 and 11 of this chapter). 3.2 Complex Eigenvalues It may happen that the roots of the characteristic polynomial are complex numbers. In analogy with the real case, we call these roots complex eigenvalues . When the matrix A has complex eigenvalues, we no longer have straight-line solutions. However, we can still derive the general solution as before by using a few tricks involving complex numbers and functions. The following examples indicate the general procedure. Example. (Center) Consider X ′ = AX with A = ( 0 β − β 0 ) and β 6 = 0. The characteristic polynomial is λ 2 + β 2 = 0, so the eigenvalues are now the imaginary numbers ± i β . Without worrying about the resulting complex vectors, we react just as before to find the eigenvector corresponding to λ = i β . We therefore solve ( − i β β − β − i β ) ( x y ) = ( 0 0 ) ,

--- PAGE 60 ---
3.2 Complex Eigenvalues 45 or i β x = β y , since the second equation is redundant. Thus we find a complex eigenvector ( 1, i ) , and so the function X ( t ) = e i β t ( 1 i ) is a complex solution of X ′ = AX . Now in general it is not polite to hand someone a complex solution to a real system of differential equations, but we can remedy this with the help of Euler’s formula: e i β t = cos β t + i sin β t . Using this fact, we rewrite the solution as X ( t ) = ( cos β t + i sin β t i ( cos β t + i sin β t ) ) = ( cos β t + i sin β t − sin β t + i cos β t ) . Better yet, by breaking X ( t ) into its real and imaginary parts, we have X ( t ) = X re ( t ) + iX im ( t ) , where X re ( t ) = ( cos β t − sin β t ) , X im ( t ) = ( sin β t cos β t ) . But now we see that both X re ( t ) and X im ( t ) are (real!) solutions of the original system. To see this, we simply check X ′ re ( t ) + iX ′ im ( t ) = X ′ ( t ) = AX ( t ) = A ( X re ( t ) + iX im ( t )) = AX re + iAX im ( t ) . Equating the real and imaginary parts of this equation yields X ′ re = AX re and X ′ im = AX im , which shows that both are indeed solutions. Moreover, since X re ( 0 ) = ( 1 0 ) , X im ( 0 ) = ( 0 1 ) , the linear combination of these solutions, X ( t ) = c 1 X re ( t ) + c 2 X im ( t ) ,

--- PAGE 61 ---
46 Chapter 3 Phase Portraits for Planar Systems Figure 3.4 Phase portrait for a center. where c 1 and c 2 are arbitrary constants, provides a solution to any initial value problem. We claim that this is the general solution of this equation. To prove this, we need to show that these are the only solutions of this equation. So suppose that this is not the case. Let Y ( t ) = ( u ( t ) v ( t ) ) be another solution. Consider the complex function f ( t ) = ( u ( t ) + iv ( t )) e i β t . Differentiating this expression and using the fact that Y ( t ) is a solution of the equation yields f ′ ( t ) = 0. Thus, u ( t ) + iv ( t ) is a complex constant times e − i β t . From this it follows directly that Y ( t ) is a linear combination of X re ( t ) and X im ( t ) . Note that each of these solutions is a periodic function with period 2 π/β . Indeed, the phase portrait shows that all solutions lie on circles centered at the origin. These circles are traversed in the clockwise direction if β > 0, counterclockwise if β < 0. See Figure 3.4. This type of system is called a center .  Example. (Spiral Sink, Spiral Source) More generally, consider X ′ = AX , where A = ( α β − β α ) and α , β 6 = 0. The characteristic equation is now λ 2 − 2 αλ + α 2 + β 2 , so the eigenvalues are λ = α ± i β . An eigenvector associated with α + i β is determined by the equation (α − (α + i β)) x + β y = 0.

--- PAGE 62 ---
3.3 Repeated Eigenvalues 47 Figure 3.5 Phase portraits for a spiral sink and a spiral source. Thus ( 1, i ) is again an eigenvector, and so we have complex solutions of the form X ( t ) = e (α + i β) t ( 1 i ) = e α t ( cos β t − sin β t ) + ie α t ( sin β t cos β t ) = X re ( t ) + iX im ( t ) . As before, both X re ( t ) and X im ( t ) yield real solutions of the system with initial conditions that are linearly independent. Thus we find the general solution, X ( t ) = c 1 e α t ( cos β t − sin β t ) + c 2 e α t ( sin β t cos β t ) . Without the term e α t , these solutions would wind periodically around circles centered at the origin. The e α t term converts solutions into spirals that either spiral into the origin (when α < 0) or away from the origin (α > 0 ) . In these cases the equilibrium point is called a spiral sink or spiral source respectively. See Figure 3.5.  3.3 Repeated Eigenvalues The only remaining cases occur when A has repeated real eigenvalues. One simple case occurs when A is a diagonal matrix of the form A = ( λ 0 0 λ ) .

--- PAGE 63 ---
48 Chapter 3 Phase Portraits for Planar Systems The eigenvalues of A are both equal to λ . In this case every nonzero vector is an eigenvector since AV = λ V for any V ∈ R 2 . Thus, solutions are of the form X ( t ) = α e λ t V . Each such solution lies on a straight line through ( 0, 0 ) and either tends to ( 0, 0 ) (if λ < 0) or away from ( 0, 0 ) (if λ > 0). So this is an easy case. A more interesting case occurs when A = ( λ 1 0 λ ) . Again, both eigenvalues are equal to λ , but now there is only one linearly inde- pendent eigenvector that is given by ( 1, 0 ) . Thus, we have one straight-line solution X 1 ( t ) = α e λ t ( 1 0 ) . To find other solutions note that the system may be written x ′ = λ x + y y ′ = λ y . Thus, if y 6 = 0, we must have y ( t ) = β e λ t . Therefore, the differential equation for x ( t ) reads x ′ = λ x + β e λ t . This is a nonautonomous, first-order differential equation for x ( t ) . One might first expect solutions of the form e λ t , but the nonautonomous term is also in this form. As you perhaps saw in calculus, the best option is to guess a solution of the form x ( t ) = α e λ t + μ te λ t for some constants α and μ . This technique is often called the method of undetermined coefficients. Inserting this guess into the differential equation

--- PAGE 64 ---
3.4 Changing Coordinates 49 Figure 3.6 Phase portrait for a system with repeated negative eigenvalues. shows that μ = β while α is arbitrary. Thus, the solution of the system may be written α e λ t ( 1 0 ) + β e λ t ( t 1 ) . This is in fact the general solution (see Exercise 12 of this chapter). Note that, if λ < 0, each term in this solution tends to 0 as t → ∞ . This is clear for the α e λ t and β e λ t terms. For the term β te λ t , this is an immediate con- sequence of l’Hˆ opital’s rule. Thus, all solutions tend to ( 0, 0 ) as t → ∞ . When λ > 0, all solutions tend away from ( 0, 0 ) . See Figure 3.6. In fact, solutions tend toward or away from the origin in a direction tangent to the eigenvector ( 1, 0 ) (see Exercise 7 at the end of this chapter). 3.4 Changing Coordinates Despite differences in the associated phase portraits, we really have dealt with only three type of matrices in these past four sections: ( λ 0 0 μ ) , ( α β − β α ) , ( λ 1 0 λ ) . Any 2 × 2 matrix that is in one of these three forms is said to be in canonical form . Systems in this form may seem rather special, but they are not. Given any linear system X ′ = AX , we can always “change coordinates” so that the new system’s coefficient matrix is in canonical form and so is easily solved. Here is how to do this.

--- PAGE 65 ---
50 Chapter 3 Phase Portraits for Planar Systems A linear map (or linear transformation ) on R 2 is a function T : R 2 → R 2 of the form T ( x y ) = ( ax + by cx + dy ) . That is, T simply multiplies any vector by the 2 × 2 matrix ( a b c d ) . We will thus think of the linear map and its matrix as being interchangeable, so that we also write T = ( a b c d ) . Hopefully no confusion will result from this slight imprecision. Now suppose that T is invertible . This means that the matrix T has an inverse matrix S that satisfies TS = ST = I where I is the 2 × 2 identity matrix. It is traditional to denote the inverse of a matrix T by T − 1 . As is easily checked, the matrix S = 1 det T ( d − b − c a ) serves as T − 1 if det T 6 = 0. If det T = 0, we know from Chapter 2 that there are infinitely many vectors ( x , y ) for which T ( x y ) = ( 0 0 ) . Thus, there is no inverse matrix in this case, for we would need ( x y ) = T − 1 T ( x y ) = T − 1 ( 0 0 ) for each such vector. We have shown this. Proposition. T he 2 × 2 matrix T is invertible if and only if det T 6 = 0 . 

--- PAGE 66 ---
3.4 Changing Coordinates 51 Now, instead of considering a linear system X ′ = AX , suppose we consider a different system, Y ′ = ( T − 1 AT ) Y , for some invertible linear map T . Note that if Y ( t ) is a solution of this new system, then X ( t ) = TY ( t ) solves X ′ = AX . Indeed, we have ( TY ( t )) ′ = TY ′ ( t ) = T ( T − 1 AT ) Y ( t ) = A ( TY ( t )) , as required. That is, the linear map T converts solutions of Y ′ = ( T − 1 AT ) Y to solutions of X ′ = AX . Alternatively, T − 1 takes solutions of X ′ = AX to solutions of Y ′ = ( T − 1 AT ) Y . We therefore think of T as a change of coordinates that converts a given linear system into one with a different coefficient matrix. What we hope to be able to do is find a linear map T that converts the given system into a sys- tem of the form Y ′ = ( T − 1 AT ) Y that is easily solved. And, as you may have guessed, we can always do this by finding a linear map that converts a given linear system to one in canonical form. Example. (Real Eigenvalues) Suppose the matrix A has two real, distinct eigenvalues λ 1 and λ 2 with associated eigenvectors V 1 and V 2 . Let T be the matrix with columns V 1 and V 2 . Thus, TE j = V j for j = 1, 2 where the E j form the standard basis of R 2 . Also, T − 1 V j = E j . Therefore, we have ( T − 1 AT ) E j = T − 1 AV j = T − 1 (λ j V j ) = λ j T − 1 V j = λ j E j . Thus the matrix T − 1 AT assumes the canonical form T − 1 AT = ( λ 1 0 0 λ 2 ) and the corresponding system is easy to solve.  Example. As a further specific example, suppose A = ( − 1 0 1 − 2 ) .

--- PAGE 67 ---
52 Chapter 3 Phase Portraits for Planar Systems The characteristic equation is λ 2 + 3 λ + 2, which yields eigenvalues λ = − 1 and λ = − 2. An eigenvector corresponding to λ = − 1 is given by solving ( A + I ) ( x y ) = ( 0 0 1 − 1 ) ( x y ) = ( 0 0 ) , which yields an eigenvector ( 1, 1 ) . Similarly an eigenvector associated with λ = − 2 is given by ( 0, 1 ) . We therefore have a pair of straight-line solutions, each tending to the origin as t → ∞ . The straight-line solution corresponding to the weaker eigen- value lies along the line y = x ; the straight-line solution corresponding to the stronger eigenvalue lies on the y -axis. All other solutions tend to the origin tangentially to the line y = x . To put this sytem in canonical form, we choose T to be the matrix with columns that are these eigenvectors: T = ( 1 0 1 1 ) , so that T − 1 = ( 1 0 − 1 1 ) . Finally, we compute T − 1 AT = ( − 1 0 0 − 2 ) , so T − 1 AT is in canonical form. The general solution of the system Y ′ = ( T − 1 AT ) Y is Y ( t ) = α e − t ( 1 0 ) + β e − 2 t ( 0 1 ) , so the general solution of X ′ = AX is TY ( t ) = ( 1 0 1 1 ) ( α e − t ( 1 0 ) + β e − 2 t ( 0 1 )) = α e − t ( 1 1 ) + β e − 2 t ( 0 1 ) .

--- PAGE 68 ---
3.4 Changing Coordinates 53 T Figure 3.7 Change of variables T in the case of a (real) sink. Thus the linear map T converts the phase portrait for the system, Y ′ = ( − 1 0 0 − 2 ) Y , to that of X ′ = AX as shown in Figure 3.7.  Note that we really do not have to go through the step of converting a specific system to one in canonical form; once we have the eigenvalues and eigenvectors, we can simply write down the general solution. We take this extra step because, when we attempt to classify all possible linear systems, the canonical form of the system will greatly simplify this process. Example. (Complex Eigenvalues) Now suppose that the matrix A has complex eigenvalues α ± i β with β 6 = 0. Then we may find a complex eigen- vector V 1 + iV 2 corresponding to α + i β , where both V 1 and V 2 are real vectors. We claim that V 1 and V 2 are linearly independent vectors in R 2 . If this were not the case, then we would have V 1 = cV 2 for some c ∈ R . But then we have A ( V 1 + iV 2 ) = (α + i β)( V 1 + iV 2 ) = (α + i β)( c + i ) V 2 . But we also have A ( V 1 + iV 2 ) = ( c + i ) AV 2 . So we conclude that AV 2 = (α + i β) V 2 . This is a contradiction since the left side is a real vector while the right is complex. Since V 1 + iV 2 is an eigenvector associated with α + i β , we have A ( V 1 + iV 2 ) = (α + i β)( V 1 + iV 2 ) .

--- PAGE 69 ---
54 Chapter 3 Phase Portraits for Planar Systems Equating the real and imaginary components of this vector equation, we find AV 1 = α V 1 − β V 2 AV 2 = β V 1 + α V 2 . Let T be the matrix with columns V 1 and V 2 . Thus TE j = V j for j = 1, 2. Now consider T − 1 AT . We have ( T − 1 AT ) E 1 = T − 1 (α V 1 − β V 2 ) = α E 1 − β E 2 and similarly ( T − 1 AT ) E 2 = β E 1 + α E 2 . Thus the matrix T − 1 AT is in the canonical form T − 1 AT = ( α β − β α ) . We saw that the system Y ′ = ( T − 1 AT ) Y has phase portrait corresponding to a spiral sink, center, or spiral source depending on whether α < 0, α = 0, or α > 0. Therefore, the phase portrait of X ′ = AX is equivalent to one of these after changing coordinates using T .  Example. (Another Harmonic Oscillator) Consider the second-order equa- tion x ′′ + 4 x = 0. This corresponds to an undamped harmonic oscillator with mass 1 and spring constant 4. As a system, we have X ′ = ( 0 1 − 4 0 ) X = AX . The characteristic equation is λ 2 + 4 = 0,

--- PAGE 70 ---
3.4 Changing Coordinates 55 so that the eigenvalues are ± 2 i . A complex eigenvector associated with λ = 2 i is a solution of the system − 2 ix + y = 0 − 4 x − 2 iy = 0. One such solution is the vector ( 1, 2 i ) . So we have a complex solution of the form e 2 it ( 1 2 i ) . Breaking this solution into its real and imaginary parts, we find the general solution X ( t ) = c 1 ( cos 2 t − 2 sin 2 t ) + c 2 ( sin 2 t 2 cos 2 t ) . Thus the position of this oscillator is given by x ( t ) = c 1 cos 2 t + c 2 sin 2 t , which is a periodic function of period π . Now, let T be the matrix with columns that are the real and imaginary parts of the eigenvector ( 1, 2 i ) ; that is T = ( 1 0 0 2 ) . Then we compute easily that T − 1 AT = ( 0 2 − 2 0 ) , which is in canonical form. The phase portraits of these systems are shown in Figure 3.8. Note that T maps the circular solutions of the system Y ′ = ( T − 1 AT ) Y to elliptic solutions of X ′ = AX .  Example. (Repeated Eigenvalues) Suppose A has a single real eigenvalue λ . If there exists a pair of linearly independent eigenvectors, then in fact A must be in the form A = ( λ 0 0 λ ) , so the system X ′ = AX is easily solved (see Exercise 15 of this chapter).

--- PAGE 71 ---
56 Chapter 3 Phase Portraits for Planar Systems T Figure 3.8 Change of variables T in the case of a center. For the more complicated case, let’s assume that V is an eigenvector and that every other eigenvector is a multiple of V . Let W be any vector for which V and W are linearly independent. Then we have AW = μ V + ν W for some constants μ , ν ∈ R . Note that μ 6 = 0, for otherwise we would have a second linearly independent eigenvector W with eigenvalue ν . We claim that ν = λ . If ν − λ 6 = 0, a computation shows that A ( W + ( μ ν − λ ) V ) = ν ( W + ( μ ν − λ ) V ) . This says that ν is a second eigenvalue different from λ . Thus, we must have ν = λ . Finally, let U = ( 1 /μ) W . Then AU = V + λ μ W = V + λ U . Thus if we define TE 1 = V , TE 2 = U , we get T − 1 AT = ( λ 1 0 λ ) , as required. X ′ = AX is therefore again in canonical form after this change of coordinates. 

--- PAGE 72 ---
Exercises 57 E X E R C I S E S 1. In Figure 3.9, you see six phase portraits. Match each of these phase portraits with one of the following linear systems: ( a ) ( 3 5 − 2 − 2 ) ( b ) ( − 3 − 2 5 2 ) ( c ) ( 3 − 2 5 − 2 ) ( d ) ( − 3 5 − 2 3 ) ( e ) ( 3 5 − 2 − 3 ) ( f ) ( − 3 5 − 2 2 ) 2. For each of the following systems of the form X ′ = AX (a) Find the eigenvalues and eigenvectors of A . (b) Find the matrix T that puts A in canonical form. (c) Find the general solution of both X ′ = AX and Y ′ = ( T − 1 AT ) Y . (d) Sketch the phase portraits of both systems. ( i ) A = ( 0 1 1 0 ) ( ii ) A = ( 1 1 1 0 ) 1. 4. 2. 5. 3. 6. Figure 3.9 Match these phase portraits with the systems in Exercise 1.

--- PAGE 73 ---
58 Chapter 3 Phase Portraits for Planar Systems ( iii ) A = ( 1 1 − 1 0 ) ( iv ) A = ( 1 1 − 1 3 ) ( v ) A = ( 1 1 − 1 − 3 ) ( vi ) A = ( 1 1 1 − 1 ) 3. Find the general solution of the following harmonic oscillator equations: (a) x ′′ + x ′ + x = 0 (b) x ′′ + 2 x ′ + x = 0 4. Consider the harmonic oscillator system X ′ = ( 0 1 − k − b ) X , where b ≥ 0, k > 0, and the mass m = 1. (a) For which values of k , b does this system have complex eigenvalues? Repeated eigenvalues? Real and distinct eigenvalues? (b) Find the general solution of this system in each case. (c) Describe the motion of the mass when the mass is released from the initial position x = 1 with zero velocity in each of the cases in part (a). 5. Sketch the phase portrait of X ′ = AX where A = ( a 1 2 a 2 ) . For which values of a do you find a bifurcation? Describe the phase portrait for a -values above and below the bifurcation point. 6. Consider the system X ′ = ( 2 a b b 0 ) X . Sketch the regions in the ab -plane where this system has different types of canonical forms. 7. Consider the system X ′ = ( λ 1 0 λ ) X with λ 6 = 0. Show that all solutions tend to (respectively, away from) the origin tangentially to the eigenvector ( 1, 0 ) when λ < 0 (respectively, λ > 0).

--- PAGE 74 ---
Exercises 59 8. Find all 2 × 2 matrices that have pure imaginary eigenvalues. That is, determine conditions on the entries of a matrix that guarantee the matrix has pure imaginary eigenvalues. 9. Determine a computable condition that guarantees that, if a matrix A has complex eigenvalues with nonzero imaginary parts, then solutions of X ′ = AX travel around the origin in the counterclockwise direction. 10. Consider the system X ′ = ( a b c d ) X , where a + d 6 = 0 but ad − bc = 0. Find the general solution of this system and sketch the phase portrait. 11. Find the general solution and describe completely the phase portrait for X ′ = ( 0 1 0 0 ) X . 12. Prove that α e λ t ( 1 0 ) + β e λ t ( t 1 ) is the general solution of X ′ = ( λ 1 0 λ ) X . 13. Prove that a 2 × 2 matrix A always satisfies its own characteristic equa- tion. That is, if λ 2 + αλ + β = 0 is the characteristic equation associated with A , then the matrix A 2 + α A + β I is the 0-matrix. 14. Suppose the 2 × 2 matrix A has repeated eigenvalues λ . Let V ∈ R 2 . Using the previous problem, show that either V is an eigenvector for A or else ( A − λ I ) V is an eigenvector for A . 15. Suppose the matrix A has repeated real eigenvalues λ and there exist, a pair of linearly independent eigenvectors associated with A . Prove that A = ( λ 0 0 λ ) . 16. Consider the (nonlinear) system x ′ = | y | y ′ = − x . Use the methods of this chapter to describe the phase portrait.

--- PAGE 75 ---
This page intentionally left blank

--- PAGE 76 ---
4 Classification of Planar Systems In this chapter, we summarize what we have accomplished so far using a dynamical systems point of view. Among other things, this means that we would like to have a complete “dictionary” of all possible behaviors of 2 × 2 linear systems. One of the dictionaries we present here is geometric: the trace– determinant plane. The other dictionary is more dynamic: This involves the notion of conjugate systems. 4.1 The Trace–Determinant Plane For a matrix A = ( a b c d ) , we know that the eigenvalues are the roots of the characteristic equation, which may be written λ 2 − ( a + d )λ + ( ad − bc ) = 0. Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00004-X c © 2013 Elsevier Inc. All rights reserved. 61

--- PAGE 77 ---
62 Chapter 4 Classification of Planar Systems The constant term in this equation is det A . The coefficient of λ also has a name: The quantity a + d is called the trace of A and is denoted by tr A . Thus the eigenvalues satisfy λ 2 − ( tr A )λ + det A = 0 and are given by λ ± = 1 2 ( tr A ± √ ( tr A ) 2 − 4 det A ) . Note that λ + + λ − = tr A and λ + λ − = det A , so the trace is the sum of the eigenvalues of A while the determinant is the product of the eigenvalues of A . We will also write T = tr A and D = det A . Knowing T and D tells us the eigenvalues of A and therefore virtually everything about the geometry of solutions of X ′ = AX . For example, the values of T and D tell us whether solutions spiral into or away from the origin, whether we have a center, and so forth. We may display this classification visually by painting a picture in the trace– determinant plane . In this picture a matrix with trace T and determinant D corresponds to the point with coordinates ( T , D ) . The location of this point in the TD -plane then determines the geometry of the phase portrait as before. For example, the sign of T 2 − 4 D tells us that the eigenvalues are 1. Complex with nonzero imaginary part if T 2 − 4 D < 0 2. Real and distinct if T 2 − 4 D > 0 3. Real and repeated if T 2 − 4 D = 0 Thus the location of ( T , D ) relative to the parabola T 2 − 4 D = 0 in the TD - plane tells us all we need to know about the eigenvalues of A from an algebraic point of view. In terms of phase portraits, however, we can say more. If T 2 − 4 D < 0, then the real part of the eigenvalues is T / 2, and so we have a 1. Spiral sink if T < 0 2. Spiral source if T > 0 3. Center if T = 0 If T 2 − 4 D > 0, we have a similar breakdown into cases. In this region, both eigenvalues are real. If D < 0, then we have a saddle. This follows since D is the product of the eigenvalues, one of which must be positive, the other negative. Equivalently, if D < 0, we compute T 2 < T 2 − 4 D

--- PAGE 78 ---
4.1 The Trace–Determinant Plane 63 so that ± T < √ T 2 − 4 D . Thus we have T + √ T 2 − 4 D > 0 T − √ T 2 − 4 D < 0, so the eigenvalues are real and have different signs. If D > 0 and T < 0, then both T ± √ T 2 − 4 D < 0, so we have a (real) sink. Similarly, T > 0 and D > 0 lead to a (real) source. When D = 0 and T 6 = 0 we have one zero eigenvalue, while both eigenvalues vanish if D = T = 0. Plotting all of this verbal information in the TD -plane gives us a visual sum- mary of all of the different types of linear systems. The preceding equations partition the TD -plane into various regions in which systems of a particular type reside. See Figure 4.1. This yields a geometric classification of 2 × 2 linear systems. A couple of remarks are in order. First, the trace–determinant plane is a two- dimensional representation of what really is a four-dimensional space, since 2 × 2 matrices are determined by four parameters, the entries of the matrix. Thus there are infinitely many different matrices corresponding to each point in the TD -plane. Although all of these matrices share the same eigenvalue con- figuration, there may be subtle differences in the phase portraits, such as the direction of rotation for centers and spiral sinks and sources, or the possibility of one or two independent eigenvectors in the repeated eigenvalue case. We also think of the trace–determinant plane as the analogue of the bifur- cation diagram for planar linear systems. A one-parameter family of linear systems corresponds to a curve in the TD -plane. When this curve crosses the T -axis, the positive D -axis, or the parabola T 2 − 4 D = 0, the phase portrait of the linear system undergoes a bifurcation: there is a major change in the geometry of the phase portrait. Finally, note that we may obtain quite a bit of information about the system from D and T without ever computing the eigenvalues. For example, if D < 0, we know that we have a saddle at the origin. Similarly, if both D and T are positive, then we have a source at the origin.

--- PAGE 79 ---
64 Chapter 4 Classification of Planar Systems Det T r T 2 = 4 D Figure 4.1 The trace–determinant plane. Any resemblance to any of the authors’ faces is purely coincidental. 4.2 Dynamical Classification In this section we give a different, more dynamical classification of planar lin- ear systems. From a dynamical systems point of view, we are usually interested primarily in the long-term behavior of solutions of differential equations. Thus two systems are equivalent if their solutions share the same fate. To make this precise we recall some terminology introduced in Chapter 1, Section 1.5. To emphasize the dependence of solutions on both time and the initial con- ditions X 0 , we let φ t ( X 0 ) denote the solution that satisfies the initial condition X 0 . That is, φ 0 ( X 0 ) = X 0 . The function φ( t , X 0 ) = φ t ( X 0 ) is called the flow of the differential equation, while φ t is called the time t map of the flow. For example, let X ′ = ( 2 0 0 3 ) X . Then the time t map is given by φ t ( x 0 , y 0 ) = ( x 0 e 2 t , y 0 e 3 t ) .

--- PAGE 80 ---
4.2 Dynamical Classification 65 Thus the flow is a function that depends on both time and initial values. We will consider two systems to be dynamically equivalent if there is a func- tion h that takes one flow to the other. We require that this function be a homeomorphism ; that is, h is a one-to-one, onto, and continuous function with an inverse that is also continuous. Definition Suppose X ′ = AX and X ′ = BX have flows φ A and φ B . These two systems are (topologically) conjugate if there exists a homeomorphism h : R 2 → R 2 that satisfies φ B ( t , h ( X 0 )) = h (φ A ( t , X 0 )) . The homeomorphism h is called a conjugacy . Thus a conjugacy takes the solution curves of X ′ = AX to those of X ′ = BX . Example. For the one-dimensional linear differential equations x ′ = λ 1 x and x ′ = λ 2 x , we have the flows φ j ( t , x 0 ) = x 0 e λ j t for j = 1, 2. Suppose that λ 1 and λ 2 are nonzero and have the same sign. Then let h ( x ) = { x λ 2 /λ 1 if x ≥ 0 −| x | λ 2 /λ 1 if x < 0, where we recall that x λ 2 /λ 1 = exp ( λ 2 λ 1 log ( x ) ) . Note that h is a homeomorphism of the real line. We claim that h is a conjugacy between x ′ = λ 1 x and x ′ = λ 2 x . To see this, we check that when x 0 > 0, h (φ 1 ( t , x 0 )) = ( x 0 e λ 1 t ) λ 2 /λ 1 = x λ 2 /λ 1 0 e λ 2 t = φ 2 ( t , h ( x 0 )) , as required. A similar computation works when x 0 < 0. 

--- PAGE 81 ---
66 Chapter 4 Classification of Planar Systems There are several things to note here. First, λ 1 and λ 2 must have the same sign, for otherwise we have | h ( 0 ) | = ∞ , in which case h is not a homeo- morphism. This agrees with our notion of dynamical equivalence: If λ 1 and λ 2 have the same sign, then their solutions behave similarly as either both tend to the origin or both tend away from the origin. Also, note that if λ 2 < λ 1 , then h is not differentiable at the origin, whereas if λ 2 > λ 1 , then h − 1 ( x ) = x λ 1 /λ 2 is not differentiable at the origin. This is the reason we require h to be only a homeomorphism and not a diffeomor- phism (a differentiable homeomorphism with a differentiable inverse): If we assume differentiability, then we must have λ 1 = λ 2 , which does not yield a very interesting notion of “equivalence.” This gives a classification of (autonomous) linear, a first-order differen- tial equations which agrees with our qualitative observations in Chapter 1. There are three conjugacy “classes”: the sinks, the sources, and the special “in-between” case, x ′ = 0, where all solutions are constants. Now we move to the planar version of this scenario. We first note that we only need to decide on conjugacies among systems with matrices in canonical form. For, as we saw in Chapter 3, if the linear map T : R 2 → R 2 puts A in canonical form, then T takes the time t map of the flow of Y ′ = ( T − 1 AT ) Y to the time t map for X ′ = AX . Our classification of planar linear systems now proceeds just as in the one- dimensional case. We will stay away from the case where the system has eigenvalues with real part equal to 0, but you will tackle this case in the Exercises at the end of this chapter. Definition A matrix A is hyperbolic if none of its eigenvalues has real part 0. We also say that the system X ′ = AX is hyperbolic. Theorem. Suppose that the 2 × 2 matrices A 1 and A 2 are hyperbolic. Then the linear systems X ′ = A i X are conjugate if and only if each matrix has the same number of eigenvalues with negative real part.  Thus, two matrices yield conjugate linear systems if both sets of eigenvalues fall into the same category: 1. One eigenvalue is positive and the other is negative. 2. Both eigenvalues have negative real parts. 3. Both eigenvalues have positive real parts. Before proving this, note that this theorem implies that a system with a spiral sink is conjugate to a system with a (real) sink. Of course! Even though their phase portraits look very different, it is nevertheless the case that all solutions of both systems share the same fate: They tend to the origin as t → ∞ .

--- PAGE 82 ---
4.2 Dynamical Classification 67 Proof: Recall from before that we may assume that all of the systems are in canonical form. Then the proof divides into three distinct cases. Case 1 Suppose we have two linear systems X ′ = A i X for i = 1, 2 such that each A i has eigenvalues λ i < 0 < μ i . Thus each system has a saddle at the origin. This is the easy case. As we saw earlier, the real differential equations x ′ = λ i x have conjugate flows via the homeomorphism h 1 ( x ) = { x λ 2 /λ 1 if x ≥ 0 −| x | λ 2 /λ 1 if x < 0 . Similarly, the equations y ′ = μ i y have conjugate flows via an analogous function h 2 . Now define H ( x , y ) = ( h 1 ( x ) , h 2 ( y )) . Then one checks immediately that H provides a conjugacy between these two systems. Case 2 Consider the system X ′ = AX where A is in canonical form with eigenvalues that have negative real parts. We further assume that the matrix A is not in the form ( λ 1 0 λ ) with λ < 0. Thus, in canonical form, A assumes one of the two forms ( a ) ( α β − β α ) ( b ) ( λ 0 0 μ ) with α , λ , μ < 0. We will show that, in either (a) or (b), the system is conjugate to X ′ = BX where B = ( − 1 0 0 − 1 ) . It then follows that any two systems of this form are conjugate. Consider the unit circle in the plane parametrized by the curve X (θ) = ( cos θ , sin θ) , 0 ≤ θ ≤ 2 π . We denote this circle by S 1 . We first claim that the

--- PAGE 83 ---
68 Chapter 4 Classification of Planar Systems vector field determined by a matrix in the preceding form must point inside S 1 . In case 2(a), we have that the vector field on S 1 is given by AX (θ) = ( α cos θ + β sin θ − β cos θ + α sin θ ) . The outward-pointing normal vector to S 1 at X (θ) is N (θ) = ( cos θ sin θ ) . The dot product of these two vectors satisfies AX (θ) · N (θ) = α( cos 2 θ + sin 2 θ) < 0 since α < 0. This shows that AX (θ) does indeed point inside S 1 . Case 2(b) is even easier. As a consequence, each nonzero solution of X ′ = AX crosses S 1 exactly once. Let φ A t denote the time t map for this system, and let τ = τ ( x , y ) denote the time at which φ A t ( x , y ) meets S 1 . Thus ∣ ∣ ∣ φ A τ ( x , y ) ( x , y ) ∣ ∣ ∣ = 1. Let φ B t denote the time t map for the system X ′ = BX . Clearly, φ B t ( x , y ) = ( e − t x , e − t y ) . We now define a conjugacy H between these two systems. If ( x , y ) 6 = ( 0, 0 ) , let H ( x , y ) = φ B − τ ( x , y ) φ A τ ( x , y ) ( x , y ) and set H ( 0, 0 ) = ( 0, 0 ) . Geometrically, the value of H ( x , y ) is given by fol- lowing the solution curve of X ′ = AX exactly τ ( x , y ) time units (forward or backward) until the solution reaches S 1 , and then following the solution of X ′ = BX starting at that point on S 1 and proceeding in the opposite time direction exactly τ time units. See Figure 4.2. To see that H gives a conjugacy, note first that τ ( φ A s ( x , y ) ) = τ ( x , y ) − s

--- PAGE 84 ---
4.2 Dynamical Classification 69 φ A τ ( x , y ) φ B −τ ( φ A τ ( x , y )) = H ( x , y ) ( x , y ) S 1 Figure 4.2 The definition of τ ( x , y ) . since φ A τ − s φ A s ( x , y ) = φ A τ ( x , y ) ∈ S 1 . Therefore, we have H ( φ A s ( x , y ) ) = φ B − τ + s φ A τ − s ( φ A s ( x , y ) ) = φ B s φ B − τ φ A τ ( x , y ) = φ B s ( H ( x , y ) ) . So H is a conjugacy. Now we show that H is a homeomorphism. We can construct an inverse for H by simply reversing the process defining H . That is, let G ( x , y ) = φ A − τ 1 ( x , y ) φ B τ 1 ( x , y ) ( x , y ) and set G ( 0, 0 ) = ( 0, 0 ) . Here τ 1 ( x , y ) is the time for the solution of X ′ = BX through ( x , y ) to reach S 1 . An easy computation shows that τ 1 ( x , y ) = log r where r 2 = x 2 + y 2 . Clearly, G = H − 1 , so H is one-to-one and onto. Also, G is continuous at ( x , y ) 6 = ( 0, 0 ) since G may be written G ( x , y ) = φ A − log r ( x r , y r ) , which is a composition of continuous functions. For continuity of G at the origin, suppose that ( x , y ) is close to the origin, so that r is small. Observe that as r → 0, − log r → ∞ . Now ( x / r , y / r ) is a point on S 1 and for r sufficiently small, φ A − log r maps the unit circle very close to ( 0, 0 ) . This shows that G is continuous at ( 0, 0 ) .

--- PAGE 85 ---
70 Chapter 4 Classification of Planar Systems We thus need only show continuity of H . For this, we need to show that τ ( x , y ) is continuous. But τ is determined by the equation ∣ ∣ ∣ φ A t ( x , y ) ∣ ∣ ∣ = 1. We write φ A t ( x , y ) = ( x ( t ) , y ( t )) . Taking the partial derivative of | φ A t ( x , y ) | with respect to t , we find ∂ ∂ t ∣ ∣ ∣ φ A t ( x , y ) ∣ ∣ ∣ = ∂ ∂ t √ ( x ( t )) 2 + ( y ( t )) 2 = 1 √ ( x ( t )) 2 + ( y ( t )) 2 ( x ( t ) x ′ ( t ) + y ( t ) y ′ ( t ) ) = 1 ∣ ∣ φ A t ( x , y ) ∣ ∣ (( x ( t ) y ( t ) ) · ( x ′ ( t ) y ′ ( t ) )) . But the latter dot product is nonzero when t = τ ( x , y ) since the vector field given by ( x ′ ( t ) , y ′ ( t )) points inside S 1 . So ∂ ∂ t ∣ ∣ ∣ φ A t ( x , y ) ∣ ∣ ∣ 6 = 0 at (τ ( x , y ) , x , y ) . Thus we may apply the Implicit Function Theorem to show that τ is differentiable at ( x , y ) and thus continuous. Continuity of H at the origin follows as in the case of G = H − 1 . Thus H is a homeomorphism and we have a conjugacy between X ′ = AX and X ′ = BX . Note that this proof works equally well if the eigenvalues have positive real parts. Case 3 Finally, suppose that A = ( λ 1 0 λ ) with λ < 0. The associated vector field need not point inside the unit circle in this case. However, if we let T = ( 1 0 0  ) , then the vector field given by Y ′ = ( T − 1 AT ) Y

--- PAGE 86 ---
Exercises 71 now does have this property, provided  > 0 is sufficiently small. Indeed, T − 1 AT = ( λ  0 λ ) , so ( T − 1 AT ( cos θ sin θ )) · ( cos θ sin θ ) = λ +  sin θ cos θ . Thus if we choose  < − λ , this dot product is negative. Therefore, the change of variables T puts us into the situation where the same proof as shown in Case 2 applies. This completes the proof in one direction. The “only if ” part of the proof follows immediately.  4.3 Exploration: A 3D Parameter Space Consider the three-parameter family of linear systems given by X ′ = ( a b c 0 ) X , where a , b , and c are parameters. 1. First fix a > 0. Describe the analogue of the trace–determinant plane in the bc -plane. That is, identify the bc -values in this plane where the corre- sponding system has saddles, centers, spiral sinks, and so on. Sketch these regions in the bc -plane. 2. Repeat the previous task when a < 0 and when a = 0. 3. Describe the bifurcations that occur as a changes from positive to negative. 4. Now put all of the previous pieces of information together and give a description of the full three-dimensional parameter space for this system. You could build a 3D model of this space, create a flip-book animation of the changes as, say, a varies, or use a computer model to visualize this image. In any event, your model should accurately capture all of the distinct regions in this space. E X E R C I S E S 1. Consider the one-parameter family of linear systems given by X ′ = ( a √ 2 + ( a / 2 ) √ 2 − ( a / 2 ) 0 ) X .

--- PAGE 87 ---
72 Chapter 4 Classification of Planar Systems (a) Sketch the path traced out by this family of linear systems in the trace–determinant plane as a varies. (b) Discuss any bifurcations that occur along this path and compute the corresponding values of a . 2. Sketch the analogue of the trace–determinant plane for the two- parameter family of systems, X ′ = ( a b b a ) X , in the ab -plane. That is, identify the regions in the ab -plane where this system has similar phase portraits. 3. Consider the harmonic oscillator equation (with m = 1), x ′′ + bx ′ + kx = 0, where b ≥ 0 and k > 0. Identify the regions in the relevant portion of the bk -plane where the corresponding system has similar phase portraits. 4. Prove that H ( x , y ) = ( x , − y ) provides a conjugacy between X ′ = ( 1 1 − 1 1 ) X and Y ′ = ( 1 − 1 1 1 ) Y . 5. For each of the following systems, find an explicit conjugacy between their flows. (a) X ′ = ( − 1 1 0 2 ) X and Y ′ = ( 1 0 1 − 2 ) Y . (b) X ′ = ( 0 1 − 4 0 ) X and Y ′ = ( 0 2 − 2 0 ) Y . 6. Prove that any two linear systems with the same eigenvalues ± i β , β 6 = 0 are conjugate. What happens if the systems have eigenvalues ± i β and ± i γ with β 6 = γ ? What if γ = − β ? 7. Consider all linear systems with exactly one eigenvalue equal to 0. Which of these systems are conjugate? Prove this. 8. Consider all linear systems with two zero eigenvalues. Which of these systems are conjugate? Prove this. 9. Provide a complete description of the conjugacy classes for 2 × 2 systems in the nonhyperbolic case.

--- PAGE 88 ---
5 Higher-Dimensional Linear Algebra As in Chapter 2, we need to make another detour into the world of linear algebra before proceeding to the solution of higher-dimensional linear sys- tems of differential equations. There are many different canonical forms for matrices in higher dimensions, but most of the algebraic ideas involved in changing coordinates to put matrices into these forms are already present in the 2 × 2 case. In particular, the case of matrices with distinct (real or complex) eigenvalues can be handled with minimal additional algebraic com- plications, so we deal with this case first. This is the “generic case,” as we show in Section 5.6. Matrices with repeated eigenvalues demand more sophisticated concepts from linear algebra; we provide this background in Section 5.4. We assume throughout this chapter that the reader is familiar with solving systems of linear algebraic equations by putting the associated matrix in (reduced) row echelon form. 5.1 Preliminaries from Linear Algebra In this section we generalize many of the algebraic notions of Section 2.3 to higher dimensions. We denote a vector X ∈ R n in coordinate form as X =    x 1 . . . x n    . Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00005-1 c © 2013 Elsevier Inc. All rights reserved. 73

--- PAGE 89 ---
74 Chapter 5 Higher-Dimensional Linear Algebra In the plane, we called a pair of vectors V and W linearly independent if they were not collinear. Equivalently, V and W were linearly independent if there were no (nonzero) real numbers α and β such that α V + β W is the zero vector. More generally, in R n , a collection of vectors V 1 , . . . , V k in R n is said to be linearly independent if, whenever α 1 V 1 + · · · + α k V k = 0 with α j ∈ R , it follows that each α j = 0. If we can find such α 1 , . . . , α k , not all of which are 0, then the vectors are linearly dependent . Note that if V 1 , . . . , V k are linearly independent and W is the linear combination, W = β 1 V 1 + · · · + β k V k , then the β j are unique. This follows since, if we could also write W = γ 1 V 1 + · · · + γ k V k , then we would have 0 = W − W = (β 1 − γ 1 ) V 1 + · · · (β k − γ k ) V k , which forces β j = γ j for each j , by linear independence of the V j . Example. The vectors ( 1, 0, 0 ) , ( 0, 1, 0 ) , and ( 0, 0, 1 ) are clearly linearly independent in R 3 . More generally, let E j be the vector in R n where the j th component is 1 and all other components are 0. Then the vectors E 1 , . . . , E n are linearly independent in R n . The collection of vectors E 1 , . . . , E n is called the standard basis of R n . We will discuss the concept of a basis in Section 5.4.  Example. The vectors ( 1, 0, 0 ) , ( 1, 1, 0 ) , and ( 1, 1, 1 ) in R 3 are also linearly independent, for if we have α 1   1 0 0   + α 2   1 1 0   + α 3   1 1 1   =   α 1 + α 2 + α 3 α 2 + α 3 α 3   =   0 0 0   , then the third component says that α 3 = 0. The fact that α 3 = 0 in the second component then says that α 2 = 0, and finally the first component similarly

--- PAGE 90 ---
5.1 Preliminaries from Linear Algebra 75 tells us that α 1 = 0. On the other hand, the vectors ( 1, 1, 1 ) , ( 1, 2, 3 ) , and ( 2, 3, 4 ) are linearly dependent, for we have 1   1 1 1   + 1   1 2 3   − 1   2 3 4   =   0 0 0   .  When solving linear systems of differential equations, we often encounter special subsets of R n called subspaces . A subspace of R n is a collection of all possible linear combinations of a given (nonempty) set of vectors. More precisely, given V 1 , . . . , V k ∈ R n , the set S = { α 1 V 1 + · · · + α k V k | α j ∈ R } is a subspace of R n . In this case we say that S is spanned by V 1 , . . . , V k . Equiv- alently, it can be shown (see Exercise 12 at the end of this chapter) that a subspace S is a nonempty subset of R n having the following two properties: 1. If X , Y ∈ S , then X + Y ∈ S 2. If X ∈ S and α ∈ R , then α X ∈ S Note that the zero vector lies in every subspace of R n and that any linear combination of vectors in a subspace S also lies in S . Example. Any straight line through the origin in R n is a subspace of R n , since this line may be written as { tV | t ∈ R } for some nonzero V ∈ R n . The single vector V spans this subspace. The plane P defined by x + y + z = 0 in R 3 is a subspace of R 3 . Indeed, any vector V in P may be written in the form ( x , y , − x − y ) or V = x   1 0 − 1   + y   0 1 − 1   , which shows that the vectors ( 1, 0, − 1 ) and ( 0, 1, − 1 ) span P .  In linear algebra, one often encounters rectangular n × m matrices, but in differential equations, most often these matrices are square ( n × n ). Conse- quently, we will assume that all matrices in this chapter are n × n . We write

--- PAGE 91 ---
76 Chapter 5 Higher-Dimensional Linear Algebra such a matrix, A =      a 11 a 12 · · · a 1 n a 21 a 22 · · · a 2 n . . . a n 1 a n 2 · · · a nn      , more compactly as A = [ a ij ]. For X = ( x 1 , . . . , x n ) ∈ R n , we define the product AX to be the vector AX =    ∑ n j = 1 a 1 j x j . . . ∑ n j = 1 a nj x j    , so that the i th entry in this vector is the dot product of the i th row of A with the vector X . Matrix sums are defined in the obvious way. If A = [ a ij ] and B = [ b ij ] are n × n matrices, then we define A + B = C , where C = [ a ij + b ij ]. Matrix arithmetric has some obvious linearity properties: 1. A ( k 1 X 1 + k 2 X 2 ) = k 1 AX 1 + k 2 AX 2 , where k j ∈ R , X j ∈ R n 2. A + B = B + A 3. ( A + B ) + C = A + ( B + C ) The product of the n × n matrices A and B is defined to be the n × n matrix AB = [ c ij ], where c ij = n ∑ k = 1 a ik b kj , so that c ij is the dot product of the i th row of A with the j th column of B . One checks easily that, if A , B , and C are n × n matrices, then 1. ( AB ) C = A ( BC ) 2. A ( B + C ) = AB + AC 3. ( A + B ) C = AC + BC 4. k ( AB ) = ( kA ) B = A ( kB ) for any k ∈ R All of the preceding properties of matrix arithmetic are easily checked by writing out the ij -entries of the corresponding matrices. It is important to remember that matrix multiplication is not commutative, so that AB 6 = BA in general. For example, ( 1 0 1 1 ) ( 1 1 0 1 ) = ( 1 1 1 2 ) ,

--- PAGE 92 ---
5.1 Preliminaries from Linear Algebra 77 whereas ( 1 1 0 1 ) ( 1 0 1 1 ) = ( 2 1 1 1 ) . Also, matrix cancellation is usually forbidden; if AB = AC , then we do not necessarily have B = C as in ( 1 1 1 1 ) ( 1 0 0 0 ) = ( 1 0 1 0 ) = ( 1 1 1 1 ) ( 1 / 2 1 / 2 1 / 2 − 1 / 2 ) . In particular, if AB is the zero matrix, it does not follow that one of A or B is also the zero matrix. The n × n matrix A is invertible if there exists an n × n matrix C for which AC = CA = I , where I is the n × n identity matrix that has 1s along the diag- onal and 0s elsewhere. The matrix C is called the inverse of A . Note that if A has an inverse, then this inverse is unique. For if AB = BA = I as well, then C = CI = C ( AB ) = ( CA ) B = IB = B . The inverse of A is denoted by A − 1 . If A is invertible, then the vector equation AX = V has a unique solution for any V ∈ R n . Indeed, A − 1 V is one solution. Moreover, it is the only one, for if Y is another solution, then we have Y = ( A − 1 A ) Y = A − 1 ( AY ) = A − 1 V . For the converse of this statement, recall that the equation AX = V has unique solutions if and only if the reduced row echelon form of the matrix A is the identity matrix. The reduced row echelon form of A is obtained by applying to A a sequence of elementary row operations of the form 1. Add k times row i of A to row j 2. Interchange row i and j 3. Multiply row i by k 6 = 0 Note that these elementary row operations correspond exactly to the opera- tions that are used to solve linear systems of algebraic equations: 1. Add k times equation i to equation j 2. Interchange equations i and j 3. Multiply equation i by k 6 = 0

--- PAGE 93 ---
78 Chapter 5 Higher-Dimensional Linear Algebra Each of these elementary row operations may be represented by multiplying A by an elementary matrix. For example, if L = [ ` ij ] is the matrix that has 1s along the diagonal, ` ji = k for some choice of i and j , i 6 = j , and all other entries are 0, then LA is the matrix obtained by performing row operation 1 on A . Similarly, if L has 1s along the diagonal with the exception that ` ii = ` jj = 0, but ` ij = ` ji = 1, and all other entries are 0, then LA is the matrix that results after performing row operation 2 on A . Finally, if L is the identity matrix with a k instead of 1 in the ii position, then LA is the matrix obtained by performing row operation 3. A matrix L in one of these three forms is called an elementary matrix. Each elementary matrix is invertible, since its inverse is given by the matrix that simply “undoes” the corresponding row operation. As a con- sequence, any product of elementary matrices is invertible. Therefore, if L 1 , . . . , L n are the elementary matrices that correspond to the row operations that put A into the reduced row echelon form that is the identity matrix, then ( L n · · · · · L 1 ) = A − 1 . That is, if the vector equation AX = V has unique solu- tions for any V ∈ R n , then A is invertible. Thus we have our first important result. Proposition. Let A be an n × n matrix. Then the system of algebraic equa- tions AX = V has a unique solution for any V ∈ R n if and only if A is invertible.  Thus the natural question now is: How do we tell if A is invertible? One answer is provided by the following result. Proposition. The matrix A is invertible if and only if the columns of A form a linearly independent set of vectors. Proof: Suppose first that A is invertible and has columns V 1 , . . . , V n . We have AE j = V j where the E j form the standard basis of R n . If the V j are not lin- early independent, we may find real numbers α 1 , . . . , α n , not all zero, such that ∑ j α j V j = 0. But then 0 = n ∑ j = 1 α j AE j = A   n ∑ j = 1 α j E j   . Thus the equation AX = 0 has two solutions, the nonzero vector (α 1 , . . . , α n ) and the 0 vector. This contradicts the previous proposition. Conversely, suppose that the V j are linearly independent. If A is not invertible, then we may find a pair of vectors X 1 and X 2 with X 1 6 = X 2 and AX 1 = AX 2 . Therefore, the nonzero vector Z = X 1 − X 2 satisfies AZ = 0. Let

--- PAGE 94 ---
5.1 Preliminaries from Linear Algebra 79 Z = (α 1 , . . . , α n ) . Then we have 0 = AZ = n ∑ j = 1 α j V j , so that the V j are not linearly independent. This contradiction establishes the result.  A more computable criterion for determining whether or not a matrix is invertible, as in the 2 × 2 case, is given by the determinant of A . Given the n × n matrix A , we will denote by A ij the ( n − 1 ) × ( n − 1 ) matrix obtained by deleting the i th row and j th column of A . Definition The determinant of A = [ a ij ] is defined inductively by det A = n ∑ k = 1 ( − 1 ) 1 + k a 1 k det A 1 k . Note that we know the determinant of a 2 × 2 matrix, so this induction makes sense for k > 2. Example. From the definition we compute det   1 2 3 4 5 6 7 8 9   = 1 det ( 5 6 8 9 ) − 2 det ( 4 6 7 9 ) + 3 det ( 4 5 7 8 ) = − 3 + 12 − 9 = 0. We remark that the definition of det A just given involves “expanding along the first row” of A . One can equally well expand along the j th row so that det A = n ∑ k = 1 ( − 1 ) j + k a jk det A jk . We will not prove this fact; the proof is an entirely straightforward though tedious calculation. Similarly, det A can be calculated by expanding along a given column (see Exercise 1 of this chapter). 

--- PAGE 95 ---
80 Chapter 5 Higher-Dimensional Linear Algebra Example. Expanding the matrix in the previous example along the second and third rows yields the same result: det   1 2 3 4 5 6 7 8 9   = − 4 det ( 2 3 8 9 ) + 5 det ( 1 3 7 9 ) − 6 det ( 1 2 7 8 ) = 24 − 60 + 36 = 0 = 7 det ( 2 3 5 6 ) − 8 det ( 1 3 4 6 ) + 9 det ( 1 2 4 5 ) = − 21 + 48 − 27 = 0. Incidentally, note that this matrix is not invertible, since   1 2 3 4 5 6 7 8 9     1 − 2 1   =   0 0 0   .  The determinant of certain types of matrices is easy to compute. A matrix [ a ij ] is called upper triangular if all entries below the main diagonal are 0. That is, a ij = 0 if i > j . Lower triangular matrices are defined similarly. We have the following proposition. Proposition. If A is an upper or lower triangular n × n matrix, then det A is the product of the entries along the diagonal. That is, det[ a ij ] = a 11 . . . a nn .  The proof is a straightforward application of induction. The following proposition describes the effects that elementary row operations have on the determinant of a matrix. Proposition. Let A and B be n × n matrices. 1. Suppose matrix B is obtained by adding a multiple of one row of A to another row of A. Then det B = det A. 2. Suppose B is obtained by interchanging two rows of A. Then det B = − det A. 3. Suppose B is obtained by multiplying each element of a row of A by k. Then det B = k det A. Proof: The proof of the proposition is straightforward when A is a 2 × 2 matrix, so we use induction. Suppose A is k × k with k > 2. To compute det B , we expand along a row that is left untouched by the row operation. By induc- tion on k , we see that det B is a sum of determinants of size ( k − 1 ) × ( k − 1 ) .

--- PAGE 96 ---
5.1 Preliminaries from Linear Algebra 81 Each of these subdeterminants has precisely the same row operation per- formed on them as in the case of the full matrix. By induction, it follows that each of these subdeterminants is multiplied by 1, − 1, or k in 1 to 3 respectively. Thus det B has the same property.  In particular, we note that if L is an elementary matrix, then det ( LA ) = ( det L )( det A ) . Indeed, det L = 1, − 1, or k in cases 1 through 3 (see Exercise 7 at the end of this chapter). The preceding proposition now yields a criterion for A to be invertible: Corollary. (Invertibility Criterion) The matrix A is invertible if and only if det A 6 = 0 . Proof: By elementary row operations, we can manipulate any matrix A into an upper triangular matrix. Then A is invertible if and only if all diagonal entries of this row-reduced matrix are nonzero. In particular, the determinant of this matrix is nonzero. Now, by the preceding observation, row operations multiply det A by nonzero numbers, so we see that all of the diagonal entries are nonzero if and only if det A is also nonzero. This concludes the proof.  This section concludes with a further important property of determinants. Proposition. det ( AB ) = ( det A )( det B ) . Proof: If either A or B is noninvertible, then AB is also noninvertible (see Exercise 11 of this chapter). Thus the proposition is true since both sides of the equation are zero. If A is invertible, then we can write A = L 1 . . . L n · I , where each L j is an elementary matrix. Thus det ( AB ) = det ( L 1 · · · L n B ) = det ( L 1 ) det ( L 2 · · · L n B ) = det ( L 1 )( det L 2 ) · · · ( det L n )( det B ) = det ( L 1 · · · L n ) det ( B ) = det ( A ) det ( B ) . 

--- PAGE 97 ---
82 Chapter 5 Higher-Dimensional Linear Algebra 5.2 Eigenvalues and Eigenvectors As we saw in Chapter 3, eigenvalues and eigenvectors play a central role in the process of solving linear systems of differential equations. Definition A vector V is an eigenvector of an n × n matrix A if V is a nonzero solution to the system of linear equations ( A − λ I ) V = 0 . The quan- tity λ is called an eigenvalue of A , and V is an eigenvector associated with λ . Just as in Chapter 2, the eigenvalues of a matrix A may be real or complex and the associated eigenvectors may have complex entries. By the Invertibility Criterion of the previous section, it follows that λ is an eigenvalue of A if and only if λ is a root of the characteristic equation det ( A − λ I ) = 0. Since A is n × n , this is a polynomial equation of degree n , which therefore has exactly n roots (counted with multiplicity). As we saw in R 2 , there are many different types of solutions of systems of differential equations, and these types depend on the configuration of the eigenvalues of A and the resulting canonical forms. There are many, many more types of canonical forms in higher dimensions. We will describe these types in this and the following sections, but we will relegate some of the more specialized proofs of these facts to the exercises at the end of this chapter. Suppose first that λ 1 , . . . , λ ` are real and distinct eigenvalues of A with associated eigenvectors V 1 , . . . , V ` . Here “distinct” means that no two of the eigenvalues are equal. Thus AV k = λ k V k for each k . We claim that the V k are linearly independent. If not, we may choose a maximal subset of the V i that are linearly independent, say V 1 , . . . , V j . Then any other eigenvector may be written in a unique way as a linear combination of V 1 , . . . , V j . Say V j + 1 is one such eigenvector. Then we may find α i , not all 0, such that V j + 1 = α 1 V 1 + · · · + α j V j . Multiplying both sides of this equation by A , we find λ j + 1 V j + 1 = α 1 AV 1 + · · · + α j AV j = α 1 λ 1 V 1 + · · · + α j λ j V j .

--- PAGE 98 ---
5.2 Eigenvalues and Eigenvectors 83 Now λ j + 1 6 = 0, for otherwise we would have α 1 λ 1 V 1 + · · · + α j λ j V j = 0, with each λ i 6 = 0. This contradicts the fact that V 1 , . . . , V j are linearly indepen- dent. Thus we have V j + 1 = α 1 λ 1 λ j + 1 V 1 + · · · + α j λ j λ j + 1 V j . Since the λ i are distinct, we have now written V j + 1 in two different ways as a linear combination of V 1 , . . . , V j . This contradicts the fact that this set of vectors is linearly independent. We have proved the following proposition. Proposition. Suppose λ 1 , . . . , λ ` are real and distinct eigenvalues for A with associated eigenvectors V 1 , . . . , V ` . Then the V j are linearly independent.  Of primary importance when we return to differential equations is the following corollary. Corollary. Suppose A is an n × n matrix with real, distinct eigenvalues. Then there is a matrix T such that T − 1 AT =    λ 1 . . . λ n    , where all of the entries off the diagonal are 0 . Proof: Let V j be an eigenvector associated to λ j . Consider the linear map T for which TE j = V j , where the E j form the standard basis of R n . That is, T is the matrix with columns that are V 1 , . . . , V n . Since the V j are linearly independent, T is invertible and we have ( T − 1 AT ) E j = T − 1 AV j = λ j T − 1 V j = λ j E j . That is, the j th column of T − 1 AT is just the vector λ j E j , as required. 

--- PAGE 99 ---
84 Chapter 5 Higher-Dimensional Linear Algebra Example. Let A =   1 2 − 1 0 3 − 2 0 2 − 2   . Expanding det ( A − λ I ) along the first column, we find that the characteristic equation of A is det ( A − λ I ) = ( 1 − λ) det ( 3 − λ − 2 2 − 2 − λ ) = ( 1 − λ)(( 3 − λ)( − 2 − λ) + 4 ) = ( 1 − λ)(λ − 2 )(λ + 1 ) , so the eigenvalues are 2, 1, and − 1. The eigenvector corresponding to λ = 2 is given by solving the equations ( A − 2 I ) X = 0, which yields − x + 2 y − z = 0 y − 2 z = 0 2 y − 4 z = 0. These equations reduce to x − 3 z = 0 y − 2 z = 0. Thus V 1 = ( 3, 2, 1 ) is an eigenvector associated to λ = 2. In similar fashion we find that ( 1, 0, 0 ) is an eigenvector associated to λ = 1, while ( 0, 1, 2 ) is an eigenvector associated to λ = − 1. Then we set T =   3 1 0 2 0 1 1 0 2   . A simple calculation shows that AT = T   2 0 0 0 1 0 0 0 − 1   .

--- PAGE 100 ---
5.3 Complex Eigenvalues 85 Since det T = − 3, T is invertible and we have T − 1 AT =   2 0 0 0 1 0 0 0 − 1   .  5.3 Complex Eigenvalues Now we treat the case where A has nonreal (complex) eigenvalues. Suppose α + i β is an eigenvalue of A with β 6 = 0. Since the characteristic equation for A has real coefficients, it follows that if α + i β is an eigenvalue, then so is its complex conjugate α + i β = α − i β . Another way to see this is the following. Let V be an eigenvector associated to α + i β . Then the equation AV = (α + i β) V shows that V is a vector with complex entries. We write V =    x 1 + iy 1 . . . x n + iy n    . Let V denote the complex conjugate of V : V =    x 1 − iy 1 . . . x n − iy n    . Then we have AV = AV = (α + i β) V = (α − i β) V , which shows that V is an eigenvector associated to the eigenvalue α − i β . Notice that we have (temporarily) stepped out of the “real” world of R n and into the world C n of complex vectors. This is not really a problem, since all of the previous linear algebraic results hold equally well for complex vectors. Now suppose that A is a 2 n × 2 n matrix with distinct nonreal eigenval- ues α j ± i β j for j = 1, . . . , n . Let V j and V j denote the associated eigenvectors.

--- PAGE 101 ---
86 Chapter 5 Higher-Dimensional Linear Algebra Then, just as in the previous proposition, this collection of eigenvectors is linearly independent. That is, if we have n ∑ j = 1 ( c j V j + d j V j ) = 0, where the c j and d j are now complex numbers, then we must have c j = d j = 0 for each j . Now we change coordinates to put A into canonical form. Let W 2 j − 1 = 1 2 ( V j + V j ) W 2 j = − i 2 ( V j − V j ) . Note that W 2 j − 1 and W 2 j are both real vectors. Indeed, W 2 j − 1 is just the real part of V j while W 2 j is its imaginary part. So working with the W j brings us back home to R n . Proposition. The vectors W 1 , . . . , W 2 n are linearly independent. Proof: Suppose not. Then we can find real numbers c j and d j for j = 1, . . . , n such that n ∑ j = 1 ( c j W 2 j − 1 + d j W 2 j ) = 0, but not all of the c j and d j are zero. So we have 1 2 n ∑ j = 1 ( c j ( V j + V j ) − id j ( V j − V j ) ) = 0, from which we find n ∑ j = 1 ( ( c j − id j ) V j + ( c j + id j ) V j ) = 0. Since the V j and the V j are linearly independent, we must have c j ± id j = 0, from which we conclude c j = d j = 0 for all j . This contradiction establishes the result. 

--- PAGE 102 ---
5.3 Complex Eigenvalues 87 Note that we have AW 2 j − 1 = 1 2 ( AV j + AV j ) = 1 2 ( (α + i β) V j + (α − i β) V j ) = α 2 ( V j + V j ) + i β 2 ( V j − V j ) = α W 2 j − 1 − β W 2 j . Similarly, we compute AW 2 j = β W 2 j − 1 + α W 2 j . Now consider the linear map T for which TE j = W j for j = 1, . . . , 2 n . That is, the matrix associated to T has columns W 1 , . . . , W 2 n . Note that this matrix has real entries. Since the W j are linearly independent, it follows from Section 5.1 that T is invertible. Now consider the matrix T − 1 AT . We have ( T − 1 AT ) E 2 j − 1 = T − 1 AW 2 j − 1 = T − 1 (α W 2 j − 1 − β W 2 j ) = α E 2 j − 1 − β E 2 j and similarly ( T − 1 AT ) E 2 j = β E 2 j − 1 + α E 2 j . Therefore, the matrix associated to T − 1 AT is T − 1 AT =    D 1 . . . D n    , where each D j is a 2 × 2 matrix of the form D j = ( α j β j − β j α j ) . This is our canonical form for matrices with distinct nonreal eigenvalues. Combining the results of this and the previous section, we have the follow- ing theorem.

--- PAGE 103 ---
88 Chapter 5 Higher-Dimensional Linear Algebra Theorem. Suppose that the n × n matrix A has distinct eigenvalues. Then we may choose a linear map T so that T − 1 AT =           λ 1 . . . λ k D 1 . . . D `           , where the D j are 2 × 2 matrices in the form D j = ( α j β j − β j α j ) .  5.4 Bases and Subspaces To deal with the case of a matrix with repeated eigenvalues, we need some further algebraic concepts. Recall that the collection of all linear combinations of a given finite set of vectors is called a subspace of R n . More precisely, given V 1 , . . . , V k ∈ R n , the set S = { α 1 V 1 + · · · + α k V k | α j ∈ R } is a subspace of R n . In this case we say that S is spanned by V 1 , . . . , V k . Definition Let S be a subspace of R n . A collection of vectors V 1 , . . . , V k is a basis of S if the V j are linearly independent and span S . Note that a subspace always has a basis, for if S is spanned by V 1 , . . . , V k , we can always throw away certain of the V j to reach a linearly independent subset of these vectors that spans S . More precisely, if the V j are not linearly independent, then we may find one of these vectors, say V k , for which V k = β 1 V 1 + · · · + β k − 1 V k − 1 .

--- PAGE 104 ---
5.4 Bases and Subspaces 89 Thus we can write any vector in S as a linear combination of the V 1 , . . . , V k − 1 alone; the vector V k is extraneous. Continuing in this fashion, we eventually reach a linearly independent subset of the V j that spans S . More important for our purposes is the following proposition. Proposition. Every basis of a subspace S ⊂ R n has the same number of elements. Proof: We first observe that the system of k linear equations in k + ` unknowns given by a 11 x 1 + · · · + a 1 k + ` x k + ` = 0 . . . a k 1 x 1 + · · · + a k k + ` x k + ` = 0 always has a nonzero solution. Indeed, using row reduction, we may first solve for one unknown in terms of the others, and then we may eliminate this unknown to obtain a system of k − 1 equations in k + ` − 1 unknowns. Thus we are finished by induction (the first case, k = 1, being obvious). Now suppose that V 1 , . . . , V k is a basis for the subspace S . Suppose that W 1 , . . . , W k + ` is also a basis of S , with ` > 0. Then each W j is a linear combination of the V i , so we have constants a ij such that W j = k ∑ i = 1 a ij V i , for j = 1, . . . , k + ` . By the preceding observation, the system of k equations k + ` ∑ j = 1 a ij x j = 0, for i = 1, . . . , k has a nonzero solution ( c 1 , . . . , c k + ` ) . Then k + ` ∑ j = 1 c j W j = k + ` ∑ j = 1 c j ( k ∑ i = 1 a ij V i ) = k ∑ i = 1 ( k + ` ∑ j = 1 a ij c j ) V i = 0, so that the W j are linearly dependent. This contradiction completes the proof. 

--- PAGE 105 ---
90 Chapter 5 Higher-Dimensional Linear Algebra As a consequence of this result, we may define the dimension of a subspace S as the number of vectors that form any basis for S . In particular, R n is a subspace of itself, and its dimension is clearly n . Furthermore, any other subspace of R n must have dimension less than n , for otherwise we would have a collection of more than n vectors in R n that are linearly independent. This cannot happen by the previous proposition. The set consisting of only the 0 vector is also a subspace, and we define its dimension to be zero. We write dim S for the dimension of the subspace S . Example. A straight line through the origin in R n forms a one-dimensional subspace of R n , since any vector on this line may be written uniquely as tV where V ∈ R n is a fixed nonzero vector lying on the line and t ∈ R is arbitrary. Clearly, the single vector V forms a basis for this subspace.  Example. The plane P in R 3 defined by x + y + z = 0 is a two-dimensional subspace of R 3 . The vectors ( 1, 0, − 1 ) and ( 0, 1, − 1 ) both lie in P and are linearly independent. If W ∈ P , we may write W =   x y − y − x   = x   1 0 − 1   + y   0 1 − 1   , so these vectors also span P .  As in the planar case, we say that a function T : R n → R n is linear if T ( X ) = AX for some n × n matrix A . T is called a linear map or linear transformation. Using the properties of matrices discussed in Section 5.1, we have T (α X + β Y ) = α T ( X ) + β T ( Y ) for any α , β ∈ R and X , Y ∈ R n . We say that the linear map T is invertible if the matrix A associated to T has an inverse. For the study of linear systems of differential equations, the most important types of subspaces are the kernels and ranges of linear maps. We define the kernel of T , denoted Ker T , to be the set of vectors mapped to 0 by T . The range of T consists of all vectors W for which there exists a vector V for which TV = W . This, of course, is a familiar concept from calculus. The difference here is that the range of T is always a subspace of R n .

--- PAGE 106 ---
5.4 Bases and Subspaces 91 Example. Consider the linear map T ( X ) =   0 1 0 0 0 1 0 0 0   X . If X = ( x , y , z ) , then T ( X ) =   y z 0   . Thus Ker T consists of all vectors of the form (α , 0, 0 ) while Range T is the set of vectors of the form (β , γ , 0 ) , where α , β , γ ∈ R . Both sets are clearly subspaces.  Example. Let T ( X ) = AX =   1 2 3 4 5 6 7 8 9   X . For Ker T , we seek vectors X that satisfy AX = 0. Using row reduction, we find that the reduced row echelon form of A is the matrix   1 0 − 1 0 1 2 0 0 0   . Thus the solutions X = ( x , y , z ) of AX = 0 satisfy x = z , y = − 2 z . Therefore, any vector in Ker T is of the form ( z , − 2 z , z ) , so Ker T has dimension 1. For Range T , note that the columns of A are vectors in Range T , since they are the images of ( 1, 0, 0 ) , ( 0, 1, 0 ) , and ( 0, 0, 1 ) respectively. These vectors are not linearly independent since − 1   1 4 7   + 2   2 5 8   =   3 6 9   . However, ( 1, 4, 7 ) and ( 2, 5, 8 ) are linearly independent, so these two vectors give a basis of Range T . 

--- PAGE 107 ---
92 Chapter 5 Higher-Dimensional Linear Algebra Proposition. Let T : R n → R n be a linear map. Then Ker T and Range T are both subspaces of R n . Moreover, dim Ker T + dim Range T = n . Proof: First suppose that Ker T = { 0 } . Let E 1 , . . . , E n be the standard basis of R n . Then we claim that TE 1 , . . . , TE n are linearly independent. If this is not the case, then we may find α 1 , . . . , α n , not all 0, such that n ∑ j = 1 α j TE j = 0. But then we have T   n ∑ j = 1 α j E j   = 0, which implies that ∑ α j E j ∈ Ker T , so that ∑ α j E j = 0. Thus each α j = 0, which is a contradiction, so the vectors TE j are linearly independent. But then, given V ∈ R n , we may write V = n ∑ j = 1 β j TE j for some β 1 , . . . , β n . Thus V = T   n ∑ j = 1 β j E j   , which shows that Range T = R n . Thus both Ker T and Range T are subspaces of R n and we have dim Ker T = 0 and dim Range T = n . If Ker T 6 = { 0 } , we may find a nonzero vector V 1 ∈ Ker T . Clearly, T (α V 1 ) = 0 for any α ∈ R , so all vectors of the form α V 1 lie in Ker T . If Ker T contains additional vectors, choose one and call it V 2 . Then Ker T contains all linear combinations of V 1 and V 2 , since T (α 1 V 1 + α 2 V 2 ) = α 1 TV 1 + α 2 TV 2 = 0. Continuing in this fashion we obtain a set of linearly independent vectors that span Ker T , thus showing that Ker T is a subspace. Note that this process must

--- PAGE 108 ---
5.5 Repeated Eigenvalues 93 end, since every collection of more than n vectors in R n is linearly dependent. A similar argument works to show that Range T is a subspace. Now suppose that V 1 , . . . , V k form a basis of Ker T where 0 < k < n (the case where k = n being obvious). Choose vectors W k + 1 , . . . , W n so that V 1 , . . . , V k , W k + 1 , . . . , W n form a basis of R n . Let Z j = TW j for each j . Then the vectors Z j are linearly independent, for if we had α k + 1 Z k + 1 + · · · + α n Z n = 0, then we would also have T (α k + 1 W k + 1 + · · · + α n W n ) = 0. This implies that α k + 1 W k + 1 + · · · + α n W n ∈ Ker T . But this is impossible, since we cannot write any W j (and thus any linear com- bination of the W j ) as a linear combination of the V i . This proves that the sum of the dimensions of Ker T and Range T is n .  We remark that it is easy to find a set of vectors that spans Range T ; simply take the set of vectors that comprise the columns of the matrix associated to T . This works since the i th column vector of this matrix is the image of the stan- dard basis vector E i under T . In particular, if these column vectors are linearly independent, then Ker T = { 0 } and there is a unique solution to the equation T ( X ) = V for every V ∈ R n . Thus we have this corollary. Corollary 1. If T : R n → R n is a linear map with dim Ker T = 0 , then T is invertible.  5.5 Repeated Eigenvalues In this section we describe the canonical forms that arise when a matrix has repeated eigenvalues. Rather than spending an inordinate amount of time developing the general theory in this case, we will give the details only for 3 × 3 and 4 × 4 matrices with repeated eigenvalues. More general cases are relegated to the exercises of this chapter. We justify this omission in the next section where we show that the “typical” matrix has distinct eigenvalues, and thus can be handled as in the previous sec- tion. (If you happen to meet a random matrix while walking down the street,

--- PAGE 109 ---
94 Chapter 5 Higher-Dimensional Linear Algebra the chances are very good that this matrix will have distinct eigenvalues!) The most general result regarding matrices with repeated eigenvalues is given by the following proposition. Proposition. Let A be an n × n matrix. Then there is a change of coordinates T for which T − 1 AT =    B 1 . . . B k    , where each of the B j s is a square matrix (and all other entries are zero) of one of the following forms ( i )         λ 1 λ 1 . . . . . . . . . 1 λ         ( ii )         C 2 I 2 C 2 I 2 . . . . . . . . . I 2 C 2         , where C 2 = ( α β − β α ) , I 2 = ( 1 0 0 1 ) , and where α , β , λ ∈ R with β 6 = 0 . The special cases where B j = (λ) or B j = ( α β − β α ) are, of course, allowed.  We first consider the case of R 3 . If A has repeated eigenvalues in R 3 , then all eigenvalues must be real. There are then two cases. Either there are two distinct eigenvalues, one of which is repeated, or else all eigenvlaues are the same. The former case can be handled by a process similar to that described in Chapter 3, so we restrict our attention here to the case where A has a single eigenvalue λ of multiplicity 3.

--- PAGE 110 ---
5.5 Repeated Eigenvalues 95 Proposition. Suppose A is a 3 × 3 matrix for which λ is the only eigenvalue. Then we may find a change of coordinates T such that T − 1 AT assumes one of the following three forms: ( i )   λ 0 0 0 λ 0 0 0 λ   ( ii )   λ 1 0 0 λ 0 0 0 λ   ( iii )   λ 1 0 0 λ 1 0 0 λ   . Proof: Let K be the kernel of A − λ I . Any vector in K is an eigenvector of A . There are then three subcases depending on whether the dimension of K is 1, 2, or 3. If the dimension of K is 3, then ( A − λ I ) V = 0 for any V ∈ R 3 . Thus A = λ I . This yields matrix (i). Suppose the dimension of K is 2. Let R be the range of A − λ I . Then R has dimension 1 since dim K + dim R = 3, as we saw in the previous section. We claim that R ⊂ K . If this is not the case, let V ∈ R be a nonzero vector. Since ( A − λ I ) V ∈ R and R is one-dimensional, we must have ( A − λ I ) V = μ V for some μ 6 = 0. But then AV = (λ + μ) V , so we have found a new eigenvalue λ + μ . This contradicts our assumption, so we must have R ⊂ K . Now let V 1 ∈ R be nonzero. Since V 1 ∈ K , V 1 is an eigenvector and so ( A − λ I ) V 1 = 0. Since V 1 also lies in R , we may find V 2 ∈ R 3 − K with ( A − λ I ) V 2 = V 1 . Since K is two-dimensional we may choose a second vec- tor V 3 ∈ K such that V 1 and V 3 are linearly independent. Note that V 3 is also an eigenvector. If we now choose the change of coordinates TE j = V j for j = 1, 2, 3, then it follows easily that T − 1 AT assumes the form of case (ii). Finally, suppose that K has dimension 1. Thus R has dimension 2. We claim that, in this case, K ⊂ R . If this is not the case, then ( A − λ I ) R = R and so A − λ I is invertible on R . Thus, if V ∈ R , there is a unique W ∈ R for which ( A − λ I ) W = V . In particular, we have AV = A ( A − λ I ) W = ( A 2 − λ A ) W = ( A − λ I )( AW ) . This shows that, if V ∈ R , then so too is AV . Thus A also preserves the subspace R . It then follows immediately that A must have an eigenvector in R , but this then says that K ⊂ R and we have a contradiction. Next we claim that ( A − λ I ) R = K . To see this, note that ( A − λ I ) R is one- dimensional, since K ⊂ R . If ( A − λ I ) R 6 = K , there is a nonzero vector V 6 ∈ K for which ( A − λ I ) R = { tV } , where t ∈ R . But then ( A − λ I ) V = tV for some

--- PAGE 111 ---
96 Chapter 5 Higher-Dimensional Linear Algebra t ∈ R , t 6 = 0, and so AV = ( t + λ) V yields another new eigenvalue. Thus we must in fact have ( A − λ I ) R = K . Now let V 1 ∈ K be an eigenvector for A . As before, there exists V 2 ∈ R such that ( A − λ I ) V 2 = V 1 . Since V 2 ∈ R there exists V 3 such that ( A − λ I ) V 3 = V 2 . Note that ( A − λ I ) 2 V 3 = V 1 . The V j are easily seen to be linearly inde- pendent. Moreover, the linear map defined by TE j = V j finally puts A into canonical form (iii). This completes the proof.  Example. Suppose A =   2 0 − 1 0 2 1 − 1 − 1 2   . Expanding along the first row, we find det ( A − λ I ) = ( 2 − λ) [ ( 2 − λ) 2 + 1] − ( 2 − λ) = ( 2 − λ) 3 , so the only eigenvalue is 2. Solving ( A − 2 I ) V = 0 yields only one indepen- dent eigenvector V 1 = ( 1, − 1, 0 ) , so we are in case (iii) of the proposition. We compute ( A − 2 I ) 2 =   1 1 0 − 1 − 1 0 0 0 0   , so that the vector V 3 = ( 1, 0, 0 ) solves ( A − 2 I ) 2 V 3 = V 1 . We also have ( A − 2 I ) V 3 = V 2 = ( 0, 0, − 1 ) . As in the preceding, we let TE j = V j for j = 1, 2, 3, so that T =   1 0 1 − 1 0 0 0 − 1 0   . Then T − 1 AT assumes the canonical form T − 1 AT =   2 1 0 0 2 1 0 0 2   . 

--- PAGE 112 ---
5.5 Repeated Eigenvalues 97 Example. Now suppose A =   1 1 0 − 1 3 0 − 1 1 2   . Again expanding along the first row, we find det ( A − λ I ) = ( 1 − λ) [ ( 3 − λ)( 2 − λ) ] + ( 2 − λ) = ( 2 − λ) 3 , so again the only eigenvalue is 2. This time, however, we have A − 2 I =   − 1 1 0 − 1 1 0 − 1 1 0   , so that we have two linearly independent eigenvectors ( x , y , z ) for which we must have x = y while z is arbitrary. Note that ( A − 2 I ) 2 is the zero matrix, so we may choose any vector that is not an eigenvector as V 2 , say V 2 = ( 1, 0, 0 ) . Then ( A − 2 I ) V 2 = V 1 = ( − 1, − 1, − 1 ) is an eigenvector. A second linearly independent eigenvector is then V 3 = ( 0, 0, 1 ) , for example. Defining TE j = V j as usual then yields the canonical form T − 1 AT =   2 1 0 0 2 0 0 0 2   .  Now we turn to the 4 × 4 case. The case of all real eigenvalues is similar to the 3 × 3 case (though a little more complicated algebraically) and is left as an exercise at the end of this chapter. Thus we assume that A has repeated complex eigenvalues α ± i β with β 6 = 0. There are just two cases; either we can find a pair of linearly independent eigenvectors corresponding to α + i β , or we can find only one such eigenvec- tor. In the former case, let V 1 and V 2 be the independent eigenvectors. The V 1 and V 2 are linearly independent eigenvectors for α − i β . As before, choose the real vectors W 1 = ( V 1 + V 1 )/ 2 W 2 = − i ( V 1 − V 1 )/ 2 W 3 = ( V 2 + V 2 )/ 2 W 4 = − i ( V 2 − V 2 )/ 2.

--- PAGE 113 ---
98 Chapter 5 Higher-Dimensional Linear Algebra If we set TE j = W j , then changing coordinates via T puts A in canonical form, T − 1 AT =     α β 0 0 − β α 0 0 0 0 α β 0 0 − β α     . If we find only one eigenvector V 1 for α + i β , then we solve the system of equations ( A − (α + i β) I ) X = V 1 as in the case of repeated real eigenvalues. The proof of the previous proposition shows that we can always find a nonzero solution V 2 of these equations. Then choose the W j as before and set TE j = W j . Then T puts A into the canonical form T − 1 AT =     α β 1 0 − β α 0 1 0 0 α β 0 0 − β α     . For example, we compute ( T − 1 AT ) E 3 = T − 1 AW 3 = T − 1 A ( V 2 + V 2 )/ 2 = T − 1 ( ( V 1 + (α + i β) V 2 )/ 2 + ( V 1 + (α − i β) V 2 )/ 2 ) = T − 1 ( ( V 1 + V 1 )/ 2 + α ( V 2 + V 2 ) / 2 + i β ( V 2 − V 2 ) / 2 ) = E 1 + α E 3 − β E 4 . Example. Let A =     1 − 1 0 1 2 − 1 1 0 0 0 − 1 2 0 0 − 1 1     . The characteristic equation, after a little computation, is (λ 2 + 1 ) 2 = 0. Thus A has eigenvalues ± i , each repeated twice. Solving the system ( A − iI ) X = 0 yields one linearly independent com- plex eigenvector V 1 = ( 1, 1 − i , 0, 0 ) associated to i . Then V 1 is an eigenvector associated to the eigenvalue − i .

--- PAGE 114 ---
5.5 Repeated Eigenvalues 99 Next we solve the system ( A − iI ) X = V 1 to find V 2 = ( 0, 0, 1 − i , 1 ) . Then V 2 solves the system ( A − iI ) X = V 1 . Finally, choose W 1 = ( V 1 + V 1 ) / 2 = Re V 1 W 2 = − i ( V 1 − V 1 ) / 2 = Im V 1 W 3 = ( V 2 + V 2 ) / 2 = Re V 2 W 4 = − i ( V 2 − V 2 ) / 2 = Im V 2 and let TE j = W j for j = 1, . . . , 4. We have T =     1 0 0 0 1 − 1 0 0 0 0 1 − 1 0 0 1 0     , T − 1 =     1 0 0 0 1 − 1 0 0 0 0 0 1 0 0 − 1 1     , and we find the canonical form T − 1 AT =     0 1 1 0 − 1 0 0 1 0 0 0 1 0 0 − 1 0     .  Example. Let A =     2 0 1 0 0 2 0 1 0 0 2 0 0 − 1 0 2     . The characteristic equation for A is ( 2 − λ) 2 (( 2 − λ) 2 + 1 ) = 0, so the eigenvalues are 2 ± i and 2 (with multiplicity 2). Solving the equations ( A − ( 2 + i ) I ) X = 0 yields an eigenvector V = ( 0, − i , 0, 1 ) for 2 + i . Let W 1 = ( 0, 0, 0, 1 ) and W 2 = ( 0, − 1, 0, 0 ) be the real and imaginary parts of V . Solving the equations ( A − 2 I ) X = 0 yields only one eigenvector associated to 2, namely W 3 = ( 1, 0, 0, 0 ) . Then we solve ( A − 2 I ) X = W 3 to find W 4 =

--- PAGE 115 ---
100 Chapter 5 Higher-Dimensional Linear Algebra ( 0, 0, 1, 0 ) . Setting TE j = W j as usual puts A into the canonical form T − 1 AT =     2 1 0 0 − 1 2 0 0 0 0 2 1 0 0 0 2     , as is easily checked.  5.6 Genericity We have mentioned several times that “most” matrices have distinct eigenval- ues. Our goal in this section is to make this precise. Recall that a set U ⊂ R n is open if whenever X ∈ U there is an open ball about X contained in U ; that is, for some a > 0 (depending on X ) the open ball about X of radius a , { Y ∈ R n ∣ ∣ | Y − X | < a } , is contained in U . Using geometrical language we say that if X belongs to an open set U , any point sufficiently near to X also belongs to U . Another kind of subset of R n is a dense set: U ⊂ R n is dense if there are points in U arbitrarily close to each point in R n . More precisely, if X ∈ R n , then for every  > 0 there exists some Y ∈ U with | X − Y | <  . Equivalently, U is dense in R n if V ∩ U is nonempty for every nonempty open set V ⊂ R n . For example, the rational numbers form a dense subset of R , as do the irrational numbers. Similarly, { ( x , y ) ∈ R 2 | both x and y are rational } is a dense subset of the plane. An interesting kind of subset of R n is a set that is both open and dense. Such a set U is characterized by the following properties: Every point in the com- plement of U can be approximated arbitrarily closely by points of U (since U is dense), but no point in U can be approximated arbitrarily closely by points in the complement (because U is open). Here is a simple example of an open and dense subset of R 2 : V = { ( x , y ) ∈ R 2 | xy 6 = 1 } .

--- PAGE 116 ---
5.6 Genericity 101 This, of course, is the complement in R 2 of the hyperbola defined by xy = 1. Suppose ( x 0 , y 0 ) ∈ V . Then x 0 y 0 6 = 1 and if | x − x 0 | , | y − y 0 | are small enough, then xy 6 = 1; this proves that V is open. Given any ( x 0 , y 0 ) ∈ R 2 , we can find ( x , y ) as close as we like to ( x 0 , y 0 ) with xy 6 = 1; this proves that V is dense. An open and dense set is a very fat set, as the following proposition shows. Proposition. Let V 1 , . . . , V m be open and dense subsets of R n . Then V = V 1 ∩ . . . ∩ V m is also open and dense. Proof: It can be easily shown that the intersection of a finite number of open sets is open, so V is open. To prove that V is dense let U ⊂ R n be a nonempty open set. Then U ∩ V 1 is nonempty since V 1 is dense. Because U and V 1 are open, U ∩ V 1 is also open. Since U ∩ V 1 is open and nonempty, ( U ∩ V 1 ) ∩ V 2 is nonempty because V 2 is dense. Since V 2 is open, U ∩ V 1 ∩ V 2 is open. Thus ( U ∩ V 1 ∩ V 2 ) ∩ V 3 is nonempty, and so on. So U ∩ V is nonempty, which proves that V is dense in R n .  We therefore think of a subset of R n as being large if this set contains an open and dense subset. To make precise what we mean by “most” matrices, we need to transfer the notion of an open and dense set to the set of all matrices. Let L ( R n ) denote the set of n × n matrices, or, equivalently, the set of linear maps of R n . In order to discuss open and dense sets in L ( R n ) , we need to have a notion of how far apart two given matrices in L ( R n ) are. But we can do this by simply writing all of the entries of a matrix as one long vector (in a specified order) and thereby thinking of L ( R n ) as R n 2 . Theorem. The set M of matrices in L ( R n ) that have n distinct eigenvalues is open and dense in L ( R n ) . Proof: We first prove that M is dense. Let A ∈ L ( R n ) . Suppose that A has some repeated eigenvalues. The proposition from the previous section states that we can find a matrix T such that T − 1 AT assumes one of two forms. Either we have a canonical form with blocks along the diagonal of the form ( i )         λ 1 λ 1 . . . . . . . . . 1 λ         or ( ii )         C 2 I 2 C 2 I 2 . . . . . . . . . I 2 C 2         ,

--- PAGE 117 ---
102 Chapter 5 Higher-Dimensional Linear Algebra where α , β , λ ∈ R with β 6 = 0 and C 2 = ( α β − β α ) , I 2 = ( 1 0 0 1 ) , or else we have a pair of separate diagonal blocks (λ) or C 2 . Either case can be handled as follows. Choose distinct values λ j such that | λ − λ j | is as small as desired, and replace the preceding block (i) with         λ 1 1 λ 2 1 . . . . . . . . . 1 λ j         . This new block now has distinct eigenvalues. In block (ii) we may similarly replace each 2 × 2 block, ( α β − β α ) , with distinct α i s. The new matrix thus has distinct eigenvalues α i ± β . In this fashion, we find a new matrix B arbitrarily close to T − 1 AT with distinct eigen- values. Then the matrix TBT − 1 also has distinct eigenvalues, and, moreover, this matrix is arbitrarily close to A . Indeed, the funtion F : L ( R n ) → L ( R n ) given by F ( M ) = TMT − 1 where T is a fixed invertible matrix is a continuous function on L ( R n ) and thus takes matrices close to T − 1 AT to new matrices close to A . This shows that M is dense. To prove that M is open, consider the characteristic polynomial of a matrix A ∈ L ( R n ) . If we vary the entries of A slightly, then the characteristic polyno- mial’s coefficients vary only slightly. Therefore, the roots of this polynomial in C move only slightly as well. Thus, if we begin with a matrix that has distinct eigenvalues, nearby matrices have this property as well. This proves that M is open.  A property P of matrices is a generic property if the set of matrices hav- ing property P contains an open and dense set in L ( R n ) . Thus a property is generic if it is shared by some open and dense set of matrices (and perhaps other matrices as well). Intuitively speaking, a generic property is one that “almost all” matrices have. Thus, having all distinct eigenvalues is a generic property of n × n matrices.

--- PAGE 118 ---
Exercises 103 E X E R C I S E S 1. Prove that the determinant of a 3 × 3 matrix can be computed by expanding along any row or column. 2. Find the eigenvalues and eigenvectors of the following matrices: ( a )   0 0 1 0 1 0 1 0 0   ( b )   0 0 1 0 2 0 3 0 0   ( c )   1 1 1 1 1 1 1 1 1   ( d )   0 0 2 0 2 0 − 2 0 0   ( e )     3 0 0 1 0 1 2 2 1 − 2 − 1 − 4 − 1 0 0 3     3. Describe the regions in a , b , c -space where the matrix   0 0 a 0 b 0 c 0 0   has real, complex, and repeated eigenvalues. 4. Describe the regions in a , b , c -space where the matrix     a 0 0 a 0 a b 0 0 c a 0 a 0 0 a     has real, complex, and repeated eigenvalues. 5. Put the following matrices in canonical form: ( a )   0 0 1 0 1 0 1 0 0   ( b )   1 0 1 0 1 0 0 0 1   ( c )   0 1 0 − 1 0 0 1 1 1   ( d )   0 1 0 1 0 0 1 1 1   ( e )   1 0 1 0 1 0 1 0 1   ( f )   1 1 0 1 1 1 0 1 1  

--- PAGE 119 ---
104 Chapter 5 Higher-Dimensional Linear Algebra ( g )   1 0 − 1 − 1 1 − 1 0 0 1   ( h )     1 0 0 1 0 1 1 0 0 0 1 0 1 0 0 0     6. Suppose that a 5 × 5 matrix has eigenvalues 2 and 1 ± i . List all possible canonical forms for a matrix of this type. 7. Let L be the elementary matrix that interchanges the i th and j th rows of a given matrix. That is, L has 1s along the diagonal, with the exception that ` ii = ` jj = 0 but ` ij = ` ji = 1. Prove that det L = − 1. 8. Find a basis for both Ker T and Range T when T is the matrix ( a ) ( 1 2 2 4 ) ( b )   1 1 1 1 1 1 1 1 1   ( c )   1 9 6 1 4 1 2 7 1   9. Suppose A is a 4 × 4 matrix that has a single real eigenvalue λ and only one independent eigenvector. Prove that A may be put in canonical form:     λ 1 0 0 0 λ 1 0 0 0 λ 1 0 0 0 λ     . 10. Suppose A is a 4 × 4 matrix with a single real eigenvalue and two linearly independent eigenvectors. Describe the possible canonical forms for A and show that A may indeed be transformed into one of these canonical forms. Describe explicitly the conditions under which A is transformed into a particular form. 11. Show that if A and/or B are noninvertible matrices, then AB is also non- invertible. 12. Suppose that S is a subset of R n having the following properties: (a) If X , Y ∈ S , then X + Y ∈ S (b) If X ∈ S and α ∈ R , then α X ∈ S Prove that S may be written as the collection of all possible linear combinations of a finite set of vectors. 13. Which of the following subsets of R n are open and/or dense? Give a brief reason in each case. (a) U 1 = { ( x , y ) | y > 0 } (b) U 2 = { ( x , y ) | x 2 + y 2 6 = 1 } (c) U 3 = { ( x , y ) | x is irrational }

--- PAGE 120 ---
Exercises 105 (d) U 4 = { ( x , y ) | x and y are not integers } (e) U 5 is the complement of a set C 1 where C 1 is closed and not dense (f) U 6 is the complement of a set C 2 that contains exactly 6 billion and 2 distinct points 14. Each of the following properties defines a subset of real n × n matrices. Which of these sets are open and/or dense in the L ( R n ) ? Give a brief reason in each case. (a) Det A 6 = 0 (b) Trace A is rational (c) Entries of A are not integers (d) 3 ≤ det A < 4 (e) − 1 < | λ | < 1 for every eigenvalue λ (f) A has no real eigenvalues (g) Each real eigenvalue of A has multiplicity 1 15. Which of the following properties of linear maps on R n are generic? (a) | λ | 6 = 1 for every eigenvalue λ (b) n = 2; one eigenvalue is not real (c) n = 3; one eigenvalue is not real (d) No solution of X ′ = AX is periodic (except the zero solution) (e) There are n distinct eigenvalues, each with distinct imaginary parts (f) AX 6 = X and AX 6 = − X for all X 6 = 0

--- PAGE 121 ---
This page intentionally left blank

--- PAGE 122 ---
6 Higher-Dimensional Linear Systems After our little sojourn into the world of linear algebra, it’s time to return to differential equations and, in particular, to the task of solving higher- dimensional linear systems with constant coefficients. As in the linear algebra chapter, we have to deal with a number of different cases. 6.1 Distinct Eigenvalues Consider first a linear system X ′ = AX where the n × n matrix A has n dis- tinct, real eigenvalues λ 1 , . . . , λ n . By the results in Chapter 5, there is a change of coordinates T so that the new system Y ′ = ( T − 1 AT ) Y assumes the particularly simple form y ′ 1 = λ 1 y 1 . . . y ′ n = λ n y n . Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00006-3 c © 2013 Elsevier Inc. All rights reserved. 107

--- PAGE 123 ---
108 Chapter 6 Higher-Dimensional Linear Systems The linear map T is the map that takes the standard basis vector E j to the eigenvector V j associated with λ j . Clearly, a function of the form Y ( t ) =    c 1 e λ 1 t . . . c n e λ n t    is a solution of Y ′ = ( T − 1 AT ) Y that satisfies the initial condition Y ( 0 ) = ( c 1 , . . . , c n ) . As in Chapter 3, this is the only such solution, for if W ( t ) =    w 1 ( t ) . . . w n ( t )    is another solution, then differentiating each expression w j ( t ) exp ( − λ j t ) , we find d dt w j ( t ) e − λ j t = ( w ′ j − λ j w j ) e − λ j t = 0. Thus w j ( t ) = c j exp (λ j t ) for each j . Therefore, the collection of solutions Y ( t ) yields the general solution of Y ′ = ( T − 1 AT ) Y . It then follows that X ( t ) = TY ( t ) is the general solution of X ′ = AX , so this general solution may be written in the form X ( t ) = n ∑ j = 1 c j e λ j t V j . Now suppose that the eigenvalues λ 1 , . . . , λ k of A are negative, while the eigenvalues λ k + 1 , . . . , λ n are positive. Since there are no zero eigenvalues, the system is hyperbolic. Then any solution that starts in the subspace spanned by the vectors V 1 , . . . , V k must first of all stay in that subspace for all time since c k + 1 = . . . = c n = 0. Second, each such solution tends to the origin as t → ∞ . In analogy with the terminology introduced for planar systems, we call this subspace the stable subspace . Similarly, the subspace spanned by V k + 1 , . . . , V n contains solutions that move away from the origin. This subspace is the unstable subspace . All other solutions tend toward the stable subspace as time goes backward and toward the unstable subspace as time increases. Therefore, this system is a higher-dimensional analogue of a saddle .

--- PAGE 124 ---
6.1 Distinct Eigenvalues 109 Example. Consider X ′ =   1 2 − 1 0 3 − 2 0 2 − 2   X . In Chapter 5, Section 5.2, we showed that this matrix has eigenvalues 2, 1, and − 1 with associated eigenvectors ( 3, 2, 1 ) , ( 1, 0, 0 ) , and ( 0, 1, 2 ) respectively. Therefore, the matrix T =   3 1 0 2 0 1 1 0 2   converts X ′ = AX to Y ′ = ( T − 1 AT ) Y =   2 0 0 0 1 0 0 0 − 1   Y , which we can solve immediately. Multiplying the solution by T then yields the general solution X ( t ) = c 1 e 2 t   3 2 1   + c 2 e t   1 0 0   + c 3 e − t   0 1 2   of X ′ = AX . The line through the origin and ( 0, 1, 2 ) is the stable line, while the plane spanned by ( 3, 2, 1 ) and ( 1, 0, 0 ) is the unstable plane. A collection of solutions of this system as well as the system Y ′ = ( T − 1 AT ) Y is displayed in Figure 6.1.  Example. If the 3 × 3 matrix A has three real, distinct eigenvalues that are negative, then we may find a change of coordinates so that the system assumes the form Y ′ = ( T − 1 AT ) Y =   λ 1 0 0 0 λ 2 0 0 0 λ 3   Y , where λ 3 < λ 2 < λ 1 < 0. All solutions therefore tend to the origin and so we have a higher-dimensional sink . See Figure 6.2. For an initial condition ( x 0 , y 0 , z 0 ) with all three coordinates nonzero, the corresponding solution tends to the origin tangentially to the x -axis (see Exercise 2 at the end of the chapter). 

--- PAGE 125 ---
110 Chapter 6 Higher-Dimensional Linear Systems z x y T (0, 1, 2) Figure 6.1 Stable and unstable subspaces of a saddle in dimension 3. On the left, the system is in canonical form. z y x Figure 6.2 A sink in three dimensions. Now suppose that the n × n matrix A has n distinct eigenvalues, of which k 1 are real and k 2 are nonreal, so that n = k 1 + 2 k 2 . Then, as in Chapter 5, we may change coordinates so that the system assumes the form x ′ j = λ j x j u ′ ` = α ` u ` + β ` v ` v ′ ` = − β ` u ` + α ` v ` for j = 1, . . . , k 1 and ` = 1, . . . , k 2 . As in Chapter 3, we therefore have solutions of the form x j ( t ) = c j e λ j t u ` ( t ) = p ` e α ` t cos β ` t + q ` e α ` t sin β ` t v ` ( t ) = − p ` e α ` t sin β ` t + q ` e α ` t cos β ` t . As before, it is straightforward to check that this is the general solution. We have therefore shown the following theorem.

--- PAGE 126 ---
6.1 Distinct Eigenvalues 111 Theorem. Consider the system X ′ = AX where A has distinct eigenvalues λ 1 , . . . , λ k 1 ∈ R and α 1 + i β 1 , . . . , α k 2 + i β k 2 ∈ C . Let T be the matrix that puts A in the canonical form T − 1 AT =           λ 1 . . . λ k 1 B 1 . . . B k 2           , where B j = ( α j β j − β j α j ) . Then the general solution of X ′ = AX is TY ( t ) , where Y ( t ) =               c 1 e λ 1 t . . . c k 1 e λ k 1 t a 1 e α 1 t cos β 1 t + b 1 e α 1 t sin β 1 t − a 1 e α 1 t sin β 1 t + b 1 e α 1 t cos β 1 t . . . a k 2 e α k 2 t cos β k 2 t + b k 2 e α k 2 t sin β k 2 t − a k 2 e α k 2 t sin β k 2 t + b k 2 e α k 2 t cos β k 2 t               .  As usual, the columns of the matrix T in this theorem are the eigenvectors (or the real and imaginary parts of the eigenvectors) corresponding to each eigenvalue. Also, as before, the subspace spanned by the eigenvectors corres- ponding to eigenvalues with negative (resp., positive) real parts is the stable (resp., unstable) subspace. Example. Consider the system X ′ =   0 1 0 − 1 0 0 0 0 − 1   X

--- PAGE 127 ---
112 Chapter 6 Higher-Dimensional Linear Systems x 2 + y 2 = a 2 Figure 6.3 Phase portrait for a spiral center. with a matrix that is already in canonical form. The eigenvalues are ± i , − 1. The solution satisfying the initial condition ( x 0 , y 0 , z 0 ) is given by Y ( t ) = x 0   cos t − sin t 0   + y 0   sin t cos t 0   + z 0 e − t   0 0 1   , so this is the general solution. The phase portrait for this system is displayed in Figure 6.3. The stable line lies along the z -axis, whereas all solutions in the xy -plane travel around circles centered at the origin. In fact, each solution that does not lie on the stable line actually lies on a cylinder in R 3 given by x 2 + y 2 = constant. These solutions spiral toward the periodic solution in the xy -plane if z 0 6 = 0.  Example. Now consider X ′ = AX where A =   − 0.1 0 1 − 1 1 − 1.1 − 1 0 − 0.1   . The characteristic equation is − λ 3 + 0.8 λ 2 − 0.81 λ + 1.01 = 0,

--- PAGE 128 ---
6.1 Distinct Eigenvalues 113 which we have kindly factored for you into ( 1 − λ)(λ 2 + 0.2 λ + 1.01 ) = 0. Therefore, the eigenvalues are the roots of this equation, which are 1 and − 0.1 ± i . Solving ( A − ( − 0.1 + i ) I ) X = 0 yields the eigenvector ( − i , 1, 1 ) asso- ciated with − 0.1 + i . Let V 1 = Re ( − i , 1, 1 ) = ( 0, 1, 1 ) and V 2 = Im ( − i , 1, 1 ) = ( − 1, 0, 0 ) . Solving ( A − I ) X = 0 yields V 3 = ( 0, 1, 0 ) as an eigenvector corre- sponding to λ = 1. Then the matrix with columns that are the V i , T =   0 − 1 0 1 0 1 1 0 0   , converts X ′ = AX into Y ′ =   − 0.1 1 0 − 1 − 0.1 0 0 0 1   Y . This system has an unstable line along the z -axis, while the xy -plane is the stable plane. Note that solutions spiral into 0 in the stable plane. We call this system a spiral saddle. See Figure 6.4. Typical solutions off the stable plane spi- ral toward the z -axis while the z -coordinate meanwhile increases or decreases. See Figure 6.5.  Figure 6.4 A spiral saddle in canonical form.

--- PAGE 129 ---
114 Chapter 6 Higher-Dimensional Linear Systems Figure 6.5 Typical spiral saddle solutions tend to spiral toward the unstable line. 6.2 Harmonic Oscillators Consider a pair of undamped harmonic oscillators with equations x ′′ 1 = − ω 2 1 x 1 x ′′ 2 = − ω 2 2 x 2 . We can almost solve these equations by inspection as visions of sin ω t and cos ω t pass through our minds. But let’s push on a bit, first to illustrate the theorem in the previous section in the case of nonreal eigenvalues, but more importantly to introduce some interesting geometry. We first introduce the new variables y j = x ′ j for j = 1, 2 so that the equations may be written as a system: x ′ j = y j y ′ j = − ω 2 j x j . In matrix form, this system is X ′ = AX , where X = ( x 1 , y 1 , x 2 , y 2 ) and A =     0 1 − ω 2 1 0 0 1 − ω 2 2 0     . This system has eigenvalues ± i ω 1 and ± i ω 2 . An eigenvector corresponding to i ω 1 is V 1 = ( 1, i ω 1 , 0, 0 ) , while V 2 = ( 0, 0, 1, i ω 2 ) is associated with i ω 2 . Let W 1

--- PAGE 130 ---
6.2 Harmonic Oscillators 115 and W 2 be the real and imaginary parts of V 1 , and let W 3 and W 4 be the same for V 2 . Then, as usual, we let TE j = W j and the linear map T puts this system into canonical form with the matrix T − 1 AT =     0 ω 1 − ω 1 0 0 ω 2 − ω 2 0     . We then see that the general solution of Y ′ = T − 1 AT · Y is Y ( t ) =     x 1 ( t ) y 1 ( t ) x 2 ( t ) y 2 ( t )     =     a 1 cos ω 1 t + b 1 sin ω 1 t − a 1 sin ω 1 t + b 1 cos ω 1 t a 2 cos ω 2 t + b 2 sin ω 2 t − a 2 sin ω 2 t + b 2 cos ω 2 t     , just as we originally expected. We could say that this is the end of the story and stop here since we have the formulas for the solution. However, let’s push on a bit more. Each pair of solutions ( x j ( t ) , y j ( t )) for j = 1, 2 is clearly a periodic solu- tion of the equation with period 2 π/ω j , but this does not mean that the full four-dimensional solution is a periodic function. Indeed, the full solution is a periodic function with period τ if and only if there exist integers m and n such that ω 1 τ = m · 2 π and ω 2 τ = n · 2 π . Thus, for periodicity, we must have τ = 2 π m ω 1 = 2 π n ω 2 or, equivalently, ω 2 ω 1 = n m . That is, the ratio of the two frequencies of the oscillators must be a rational number. In Figure 6.6 we have plotted ( x 1 ( t ) , x 2 ( t )) for the particular solution of this system when the ratio of the frequencies is 5 / 2. When the ratio of the frequencies is irrational, something very different happens. To understand this, we make another (and much more familiar) change of coordinates. In canonical form, our system currently is x ′ j = ω j y j y ′ j = − ω j x j .

--- PAGE 131 ---
116 Chapter 6 Higher-Dimensional Linear Systems x 1 x 2 Figure 6.6 A solution with frequency ratio 5 / 2 projected into the x 1 x 2 -plane. Note that x 2 ( t ) oscillates five times and x 1 ( t ) only twice before returning to the initial position. Let’s now introduce polar coordinates ( r j , θ j ) in place of the x j and y j variables. Differentiating r 2 j = x 2 j + y 2 j , we find 2 r j r ′ j = 2 x j x ′ j + 2 y j y ′ j = 2 x j y j ω j − 2 x j y j ω j = 0. Therefore, r ′ j = 0 for each j . Also, differentiating the equation tan θ j = y j x j yields ( sec 2 θ j )θ ′ j = y ′ j x j − y j x ′ j x 2 j = − ω j r 2 j r 2 j cos 2 θ j , from which we find θ ′ j = − ω j .

--- PAGE 132 ---
6.2 Harmonic Oscillators 117 So, in polar coordinates, these equations really are quite simple: r ′ j = 0 θ ′ j = − ω j . The first equation tells us that both r 1 and r 2 remain constant along any solution. Then, no matter what we pick for our initial r 1 and r 2 values, the θ j equations remain the same. Thus we may as well restrict our attention to r 1 = r 2 = 1. The resulting set of points in R 4 is a torus —the surface of a doughnut—although this is a little difficult to visualize in four-dimensional space. However, we know that we have two independent variables on this set, namely θ 1 and θ 2 , and both are periodic with period 2 π . So this is akin to the two independent circular directions that parametrize the familiar torus in R 3 . Restricted to this torus, the equations now read θ ′ 1 = − ω 1 θ ′ 2 = − ω 2 . It is convenient to think of θ 1 and θ 2 as variables in a square of sidelength 2 π where we glue together the opposite sides θ j = 0 and θ j = 2 π to make the torus. In this square our vector field now has constant slope θ ′ 2 θ ′ 1 = ω 2 ω 1 . Therefore, solutions lie along straight lines with slope ω 2 /ω 1 in this square. When a solution reaches the edge θ 1 = 2 π (say at θ 2 = c ), it instantly reap- pears on the edge θ 1 = 0 with θ 2 coordinate given by c , and then continues onward with slope ω 2 /ω 1 . A similar identification occurs when the solution meets θ 2 = 2 π . So now we have a simplified geometric vision of what happens to these solu- tions on these tori. But what really happens? The answer depends on the ratio ω 2 /ω 1 . If this ratio is a rational number, say n / m , then the solution starting at (θ 1 ( 0 ) , θ 2 ( 0 )) will pass through the torus horizontally exactly m times and vertically n times before returning to its starting point. This is the periodic solution we observed previously. Incidentally, the picture of the straight-line solutions in the θ 1 θ 2 -plane is not at all the same as our depiction of solutions in the x 1 x 2 -plane as shown in Figure 6.6. In the irrational case, something quite different occurs. See Figure 6.7. To understand what is happening here, we return to the notion of a Poincar´ e map discussed in Chapter 1. Consider the circle θ 1 = 0, the left edge of our square

--- PAGE 133 ---
118 Chapter 6 Higher-Dimensional Linear Systems Figure 6.7 A solution with frequency ratio √ 2 projected into the x 1 x 2 -plane, the left curve computed up to time 50 π ; the right, to time 100 π . representation of the torus. Given an initial point on this circle, say θ 2 = x 0 , we follow the solution starting at this point until it next hits θ 1 = 2 π . By our identification, this solution has now returned to the circle θ 1 = 0. The solution may cross the boundary θ 2 = 2 π several times in making this transit, but it does eventually return to θ 1 = 0. So we may define the Poincar´ e map on θ 1 = 0 by assigning to x 0 on this circle the corresponding coordinate of the point of first return. Suppose that this first return occurs at the point θ 2 (τ ) where τ is the time for which θ 1 (τ ) = 2 π . Since θ 1 ( t ) = θ 1 ( 0 ) − ω 1 t , we have τ = 2 π/ω 1 . Thus θ 2 (τ ) = x 0 − ω 2 ( 2 π/ω 1 ) . Therefore, the Poincar´ e map on the circle may be written as f ( x 0 ) = x 0 + 2 π(ω 2 /ω 1 ) mod 2 π , where x 0 = θ 2 ( 0 ) is our initial θ 2 coordinate on the circle. See Figure 6.8. Thus the Poincar´ e map on the circle is just the function that rotates points on the circle by angle 2 π(ω 2 /ω 1 ) . Since ω 2 /ω 1 is irrational, this function is called an irrational rotation of the circle. Definition The set of points x 0 , x 1 = f ( x 0 ) , x 2 = f ( f ( x 0 )) , . . . , x n = f ( x n − 1 ) is called the orbit of x 0 under iteration of f . The orbit of x 0 tracks how our solution successively crosses θ 1 = 2 π as time increases. Proposition. Suppose ω 2 /ω 1 is irrational. Then the orbit of any initial point x 0 on the circle θ 1 = 0 is dense in the circle.

--- PAGE 134 ---
6.2 Harmonic Oscillators 119 x 1 = f ( x 0 ) θ 1 = 0 θ 1 = 2 π x 0 x 2 Figure 6.8 Poincar ´ e map on the circle θ 1 = 0 in the θ 1 θ 2 -torus. Proof: Recall from Chapter 5, Section 6 that a subset of the circle is dense if there are points in this subset that are arbitrarily close to any point whatsoever in the circle. Therefore we must show that, given any point z on the circle and any  > 0, there is a point x n on the orbit of x 0 such that | z − x n | <  where z and x n are measured mod 2 π . To see this, observe first that there must be n , m for which m > n and | x n − x m | <  . Indeed, we know that the orbit of x 0 is not a finite set of points since ω 2 /ω 1 is irrational. Thus there must be at least two of these points where the distance apart is less than  since the circle has finite circumference. These are the points x n and x m (actually, there must be infinitely many such points). Now rotate these points in the reverse direction exactly n times. The points x n and x m are rotated to x 0 and x m − n respectively. We find, after this rotation, that | x 0 − x m − n | <  . Now x m − n is given by rotating the circle through angle ( m − n ) 2 π(ω 2 /ω 1 ) , in which mod 2 π is therefore a rotation of angle less than  . Thus, performing this rotation again, we find | x 2 ( m − n ) − x m − n | <  as well, and, inductively, | x k ( m − n ) − x ( k − 1 )( m − n ) | <  for each k . Thus we have found a sequence of points obtained by repeated rotation through angle ( m − n ) 2 π(ω 2 /ω 1 ) , and each of these points is within  of its predecessor. Thus there must be a point of this form within  of z .  Since the orbit of x 0 is dense in the circle θ 1 = 0, it follows that the straight- line solutions connecting these points in the square are also dense, and so the original solutions are dense in the torus on which they reside. This accounts

--- PAGE 135 ---
120 Chapter 6 Higher-Dimensional Linear Systems for the densely packed solution shown projected into the x 1 x 2 -plane shown in Figure 6.7 when ω 2 /ω 1 = √ 2. Returning to the actual motion of the oscillators, we see that when ω 2 /ω 1 is irrational, the masses do not move in periodic fashion. However, they do come back very close to their initial positions over and over again as time goes on, due to the density of these solutions on the torus. These types of motions are called quasiperiodic motions . In Exercise 7 at the end of this chapter, we investigate a related set of equations, namely a pair of coupled oscillators. 6.3 Repeated Eigenvalues As we saw in the previous chapter, the solution of systems with repeated real eigenvalues reduces to solving systems with matrices that contain blocks of the form         λ 1 λ 1 . . . . . . . . . 1 λ         . Example. Let X ′ =   λ 1 0 0 λ 1 0 0 λ   X . The only eigenvalue for this system is λ , and its only eigenvector is ( 1, 0, 0 ) . We may solve this system as we did in Chapter 3, by first noting that x ′ 3 = λ x 3 , so we must have x 3 ( t ) = c 3 e λ t . Now we must have x ′ 2 = λ x 2 + c 3 e λ t . As in Chapter 3, we guess a solution of the form x 2 ( t ) = c 2 e λ t + α te λ t . Substituting this guess into the differential equation for x ′ 2 , we determine that α = c 3 and find x 2 ( t ) = c 2 e λ t + c 3 te λ t .

--- PAGE 136 ---
6.3 Repeated Eigenvalues 121 Finally, the equation x ′ 1 = λ x 1 + c 2 e λ t + c 3 te λ t suggests the guess x 1 ( t ) = c 1 e λ t + α te λ t + β t 2 e λ t . Solving as before, we find x 1 ( t ) = c 1 e λ t + c 2 te λ t + c 3 t 2 2 e λ t . Altogether, we find X ( t ) = c 1 e λ t   1 0 0   + c 2 e λ t   t 1 0   + c 3 e λ t   t 2 / 2 t 1   , which is the general solution. Despite the presence of the polynomial terms in this solution, when λ < 0, the exponential term dominates and all solu- tions do tend to zero. Some representative solutions when λ < 0 are shown in Figure 6.9. Note that there is only one straight-line solution for this system; this solution lies on the x -axis. Also, the xy -plane is invariant and solutions there behave exactly as in the planar repeated eigenvalue case.  x z Figure 6.9 Phase portrait for repeated real eigenvalues.

--- PAGE 137 ---
122 Chapter 6 Higher-Dimensional Linear Systems Example. Consider the four-dimensional system x ′ 1 = x 1 + x 2 − x 3 x ′ 2 = x 2 + x 4 x ′ 3 = x 3 + x 4 x ′ 4 = x 4 . We may write this system in matrix form as X ′ = AX =     1 1 − 1 0 0 1 0 1 0 0 1 1 0 0 0 1     X . Since A is upper triangular, all of the eigenvalues are equal to 1. Solving ( A − I ) X = 0, we find two independent eigenvectors V 1 = ( 1, 0, 0, 0 ) and W 1 = ( 0, 1, 1, 0 ) . This reduces the possible canonical forms for A to two pos- sibilities. Solving ( A − I ) X = V 1 yields one solution, V 2 = ( 0, 1, 0, 0 ) , and solving ( A − I ) X = W 1 yields another solution, W 2 = ( 0, 0, 0, 1 ) . Thus, we know that the system X ′ = AX may be tranformed into Y ′ = ( T − 1 AT ) Y =     1 1 0 0 0 1 0 0 0 0 1 1 0 0 0 1     Y , where the matrix T is given by T =     1 0 0 0 0 1 1 0 0 0 1 0 0 0 0 1     . Solutions of Y ′ = ( T − 1 AT ) Y therefore are given by y 1 ( t ) = c 1 e t + c 2 te t y 2 ( t ) = c 2 e t y 3 ( t ) = c 3 e t + c 4 te t y 4 ( t ) = c 4 e t .

--- PAGE 138 ---
6.4 The Exponential of a Matrix 123 Applying the change of coordinates T , we find the general solution of the original system: x 1 ( t ) = c 1 e t + c 2 te t x 2 ( t ) = c 2 e t + c 3 e t + c 4 te t x 3 ( t ) = c 3 e t + c 4 te t x 4 ( t ) = c 4 e t .  6.4 The Exponential of a Matrix We turn now to an alternative and elegant approach to solving linear systems using the exponential of a matrix. In a certain sense, this is the more natural way to attack these systems. Recall how we solved the 1 × 1 “system” of linear equations x ′ = ax , where our matrix was now simply ( a ) . We did not go through the process of finding eigenvalues and eigenvectors here (well, actually, we did, but the process was pretty simple). Rather, we just exponentiated the matrix ( a ) to find the general solution x ( t ) = c exp ( at ) . In fact, this process works in the general case where A is n × n . All we need to know is how to exponentiate a matrix. Here’s how: Recall from calculus that the exponential function can be expressed as the infinite series e x = ∞ ∑ k = 0 x k k ! . We know that this series converges for every x ∈ R . Now we can add matrices; we can raise them to the power k ; and we can multiply each entry by 1 / k !. So this suggests that we can use this series to exponentiate them as well. Definition Let A be an n × n matrix. We define the exponential of A to be the matrix given by exp ( A ) = ∞ ∑ k = 0 A k k ! . Of course, we have to worry about what it means for this sum of matrices to converge, but let’s put that off and try to compute a few examples first.

--- PAGE 139 ---
124 Chapter 6 Higher-Dimensional Linear Systems Example. Let A = ( λ 0 0 μ ) . Then we have A k = ( λ k 0 0 μ k ) so that exp ( A ) =       ∞ ∑ k = 0 λ k / k ! 0 0 ∞ ∑ k = 0 μ k / k !       = ( e λ 0 0 e μ ) , as you may have guessed.  Example. For a slightly more complicated example, let A = ( 0 β − β 0 ) . We compute A 0 = I , A 2 = − β 2 I , A 3 = − β 3 ( 0 1 − 1 0 ) , A 4 = β 4 I , A 5 = β 5 ( 0 1 − 1 0 ) , . . . so we find exp ( A ) =           ∞ ∑ k = 0 ( − 1 ) k β 2 k ( 2 k ) ! ∞ ∑ k = 0 ( − 1 ) k β 2 k + 1 ( 2 k + 1 ) ! − ∞ ∑ k = 0 ( − 1 ) k β 2 k + 1 ( 2 k + 1 ) ! ∞ ∑ k = 0 ( − 1 ) k β 2 k ( 2 k ) !           = ( cos β sin β − sin β cos β ) . 

--- PAGE 140 ---
6.4 The Exponential of a Matrix 125 Example. Now let A = ( λ 1 0 λ ) with λ 6 = 0. With an eye toward what comes later, we compute, not exp A , but rather exp ( tA ) . We have ( tA ) k = ( ( t λ) k kt k λ k − 1 0 ( t λ) k ) . Thus we find exp ( tA ) =         ∞ ∑ k = 0 ( t λ) k k ! t ∞ ∑ k = 0 ( t λ) k k ! 0 ∞ ∑ k = 0 ( t λ) k k !         = ( e t λ te t λ 0 e t λ ) .  Note that, in each of these three examples, the matrix exp ( A ) is a matrix with entries that are infinite series. We therefore say that the infinite series of matrices exp ( A ) converges absolutely if each of its individual terms does so. In each of the preceding cases, this convergence was clear. Unfortunately, in the case of a general matrix A , this is not so clear. To prove convergence here, we need to work a little harder. Let a ij ( k ) denote the ij -entry of A k . Let a = max | a ij | . We have | a ij ( 2 ) | = ∣ ∣ ∣ ∣ n ∑ k = 1 a ik a kj ∣ ∣ ∣ ∣ ≤ na 2 | a ij ( 3 ) | = ∣ ∣ ∣ ∣ n ∑ k , ` = 1 a ik a k ` a ` j ∣ ∣ ∣ ∣ ≤ n 2 a 3 . . . | a ij ( k ) | ≤ n k − 1 a k . Thus we have a bound for the ij -entry of the n × n matrix exp ( A ) : ∣ ∣ ∣ ∣ ∞ ∑ k = 0 a ij ( k ) k ! ∣ ∣ ∣ ∣ ≤ ∞ ∑ k = 0 | a ij ( k ) | k ! ≤ ∞ ∑ k = 0 n k − 1 a k k ! ≤ ∞ ∑ k = 0 ( na ) k k ! ≤ exp na , so that this series converges absolutely by the comparison test. Therefore, the matrix exp A makes sense for any A ∈ L ( R n ) .

--- PAGE 141 ---
126 Chapter 6 Higher-Dimensional Linear Systems The following result shows that matrix exponentiation shares many of the familiar properties of the usual exponential function. Proposition. Let A , B, and T be n × n matrices. Then: 1. If B = T − 1 AT, then exp ( B ) = T − 1 exp ( A ) T 2. If AB = BA, then exp ( A + B ) = exp ( A ) exp ( B ) 3. exp ( − A ) = ( exp ( A )) − 1 Proof: The proof of (1) follows from the identities T − 1 ( A + B ) T = T − 1 AT + T − 1 BT and ( T − 1 AT ) k = T − 1 A k T . Therefore, T − 1 ( n ∑ k = 0 A k k ! ) T = n ∑ k = 0 ( T − 1 AT ) k k ! and (1) follows by taking limits. To prove (2), observe that because AB = BA we have by the binomial theorem ( A + B ) n = n ! ∑ j + k = n A j j ! B k k ! . Therefore, we must show that ∞ ∑ n = 0   ∑ j + k = n A j j ! B k k !   =   ∞ ∑ j = 0 A j j !   ( ∞ ∑ k = 0 B k k ! ) . This is not as obvious as it may seem, since we are dealing here with series of matrices, not series of real numbers. So we will prove this in the following lemma, which then proves (2). Putting B = − A in (2) gives (3).  Lemma. For any n × n matrices A and B, we have ∞ ∑ n = 0   ∑ j + k = n A j j ! B k k !   =   ∞ ∑ j = 0 A j j !   ( ∞ ∑ k = 0 B k k ! ) . Proof: We know that each of these infinite series of matrices converges. We just have to check that they converge to each other. To do this, consider the partial sums γ 2 m = 2 m ∑ n = 0   ∑ j + k = n A j j ! B k k !  

--- PAGE 142 ---
6.4 The Exponential of a Matrix 127 and α m =   m ∑ j = 0 A j j !   and β m = ( m ∑ k = 0 B k k ! ) . We need to show that the matrices γ 2 m − α m β m tend to the zero matrix as m → ∞ . Toward that end, for a matrix M = [ m ij ], we let || M || = max | m ij | . We will show that || γ 2 m − α m β m || → 0 as m → ∞ . A computation shows that γ 2 m − α m β m = ∑ ′ A j j ! B k k ! + ∑ ′′ A j j ! B k k ! , where ∑ ′ denotes the sum over terms with indices satisfying j + k ≤ 2 m , 0 ≤ j ≤ m , m + 1 ≤ k ≤ 2 m while ∑ ′′ denotes the sum corresponding to j + k ≤ 2 m , m + 1 ≤ j ≤ 2 m , 0 ≤ k ≤ m . Therefore, || γ 2 m − α m β m || ≤ ∑ ′ ∣ ∣ ∣ ∣ A j j ! ∣ ∣ ∣ ∣ · ∣ ∣ ∣ ∣ B k k ! ∣ ∣ ∣ ∣ + ∑ ′′ ∣ ∣ ∣ ∣ A j j ! ∣ ∣ ∣ ∣ · ∣ ∣ ∣ ∣ B k k ! ∣ ∣ ∣ ∣ . Now ∑ ′ ∣ ∣ ∣ ∣ A j j ! ∣ ∣ ∣ ∣ · ∣ ∣ ∣ ∣ B k k ! ∣ ∣ ∣ ∣ ≤ ( m ∑ j = 0 ∣ ∣ ∣ ∣ A j j ! ∣ ∣ ∣ ∣ )( 2 m ∑ k = m + 1 ∣ ∣ ∣ ∣ B k k ! ∣ ∣ ∣ ∣ ) . This tends to 0 as m → ∞ since, as we saw previously, ∞ ∑ j = 0 ∣ ∣ ∣ ∣ A j j ! ∣ ∣ ∣ ∣ ≤ exp ( n || A || ) < ∞ . Similarly, ∑ ′′ ∣ ∣ ∣ ∣ A j j ! ∣ ∣ ∣ ∣ · ∣ ∣ ∣ ∣ B k k ! ∣ ∣ ∣ ∣ → 0 as m → ∞ . Therefore, lim m →∞ (γ 2 m − α m β m ) = 0, proving the lemma. 

--- PAGE 143 ---
128 Chapter 6 Higher-Dimensional Linear Systems Observe that statement (3) of the proposition implies that exp ( A ) is invert- ible for every matrix A . This is analogous to the fact that e a 6 = 0 for every real number a . There is a very simple relationship between the eigenvectors of A and those of exp ( A ) . Proposition. If V ∈ R n is an eigenvector of A associated with the eigenvalue λ , then V is also an eigenvector of exp ( A ) associated with e λ . Proof: From AV = λ V , we obtain exp ( A ) V = lim n →∞ ( n ∑ k = 0 A k V k ! ) = lim n →∞ ( n ∑ k = 0 λ k k ! V ) = ( ∞ ∑ k = 0 λ k k ! ) V = e λ V .  Now let’s return to the setting of systems of differential equations. Let A be an n × n matrix and consider the system X ′ = AX . Recall that L ( R n ) denotes the set of all n × n matrices. We have a function R → L ( R n ) which assigns the matrix exp ( tA ) to t ∈ R . Since L ( R n ) is identified with R n 2 , it makes sense to speak of the derivative of this function. Proposition. d dt exp ( tA ) = A exp ( tA ) = exp ( tA ) A . In other words, the derivative of the matrix-valued function t → exp ( tA ) is another matrix-valued function A exp ( tA ) . Proof: We have d dt exp ( tA ) = lim h → 0 exp (( t + h ) A ) − exp ( tA ) h = lim h → 0 exp ( tA ) exp ( hA ) − exp ( tA ) h = exp ( tA ) lim h → 0 ( exp ( hA ) − I h ) = exp ( tA ) A .

--- PAGE 144 ---
6.4 The Exponential of a Matrix 129 That the last limit equals A follows from the series definition of exp ( hA ) . Note that A commutes with each term of the series for exp ( tA ) , thus with exp ( tA ) . This proves the proposition.  Now we return to solving systems of differential equations. The following may be considered the fundamental theorem of linear differential equations with constant coefficients. Theorem. Let A be an n × n matrix. Then the solution of the initial value problem X ′ = AX with X ( 0 ) = X 0 is X ( t ) = exp ( tA ) X 0 . Moreover, this is the only such solution. Proof: The preceding proposition shows that d dt ( exp ( tA ) X 0 ) = ( d dt exp ( tA ) ) X 0 = A exp ( tA ) X 0 . Moreover, since exp ( 0 A ) X 0 = X 0 , it follows that this is a solution of the initial value problem. To see that there are no other solutions, let Y ( t ) be another solution satisfying Y ( 0 ) = X 0 and set Z ( t ) = exp ( − tA ) Y ( t ) . Then Z ′ ( t ) = ( d dt exp ( − tA ) ) Y ( t ) + exp ( − tA ) Y ′ ( t ) = − A exp ( − tA ) Y ( t ) + exp ( − tA ) AY ( t ) = exp ( − tA )( − A + A ) Y ( t ) ≡ 0. Therefore, Z ( t ) is a constant. Setting t = 0 shows Z ( t ) = X 0 , so that Y ( t ) = exp ( tA ) X 0 . This completes the proof of the theorem.  Note that this proof is identical to that given in Chapter 1, Section 1.1. Only the meaning of the letter A has changed. Example. Consider the system X ′ = ( λ 1 0 λ ) X .

--- PAGE 145 ---
130 Chapter 6 Higher-Dimensional Linear Systems By the theorem, the general solution is X ( t ) = exp ( tA ) X 0 = exp ( t λ t 0 t λ ) X 0 . But this is precisely the exponential of the matrix we computed earlier. We find that X ( t ) = ( e t λ te t λ 0 e t λ ) X 0 . Note that this agrees with our computations in Chapter 3.  6.5 Nonautonomous Linear Systems Up to this point, virtually all of the linear systems of differential equations that we have encountered have been autonomous. There are, however, certain types of nonautonomous systems that often arise in applications. One such system is of the form X ′ = A ( t ) X , where A ( t ) = [ a ij ( t ) ] is an n × n matrix that depends continuously on time. We will investigate these types of systems further when we encounter the variational equation in subsequent chapters. Here we restrict our attention to a different type of nonautonomous linear system given by X ′ = AX + G ( t ) , where A is a constant n × n matrix and G : R → R n is a forcing term that depends explicitly on t . This is an example of a first-order, linear, nonho- mogeneous system of equations. Example. (The Forced Harmonic Oscillator) If we apply an external force to the harmonic oscillator system, the differential equation governing the motion becomes x ′′ + bx ′ + kx = f ( t ) , where f ( t ) measures the external force. An important special case occurs when this force is a periodic function of time, which corresponds, for example, to moving the table on which the mass-spring apparatus resides back and forth periodically. As a system, the forced harmonic oscillator equation becomes X ′ = ( 0 1 − k − b ) X + G ( t ) , where G ( t ) = ( 0 f ( t ) ) . 

--- PAGE 146 ---
6.5 Nonautonomous Linear Systems 131 For a nonhomogeneous system, the equation that results from dropping the time-dependent term, namely X ′ = AX , is called the homogeneous equation. We know how to find the general solution of this system. Borrowing the nota- tion from the previous section, the solution satisfying the initial condition X ( 0 ) = X 0 is X ( t ) = exp ( tA ) X 0 , so this is the general solution of the homogeneous equation. To find the general solution of the nonhomogeneous equation, suppose that we have one particular solution Z ( t ) of this equation. So Z ′ ( t ) = AZ ( t ) + G ( t ) . If X ( t ) is any solution of the homogeneous equation, then the function Y ( t ) = X ( t ) + Z ( t ) is another solution of the nonhomogeneous equation. This follows since we have Y ′ = X ′ + Z ′ = AX + AZ + G ( t ) = A ( X + Z ) + G ( t ) = AY + G ( t ) . Therefore, since we know all solutions of the homogeneous equation, we can now find the general solution to the nonhomogeneous equation, provided that we can find just one particular solution of this equation. Often one gets such a solution by simply guessing it (in calculus, this method is usually called the method of undetermined coefficients). Unfortunately, guessing a solution does not always work. The following method, called variation of parameters, does work in all cases. However, there is no guarantee that we can actually evaluate the required integrals. Theorem. (Variation of Parameters) Consider the nonhomogeneous equa- tion X ′ = AX + G ( t ) , where A is an n × n matrix and G ( t ) is a continuous function of t. Then X ( t ) = exp ( tA )   X 0 + t ∫ 0 exp ( − sA ) G ( s ) ds   is a solution of this equation satisfying X ( 0 ) = X 0 .

--- PAGE 147 ---
132 Chapter 6 Higher-Dimensional Linear Systems Proof: Differentiating X ( t ) , we obtain X ′ ( t ) = A exp ( tA )   X 0 + t ∫ 0 exp ( − sA ) G ( s ) ds   + exp ( tA ) d dt t ∫ 0 exp ( − sA ) G ( s ) ds = A exp ( tA )   X 0 + t ∫ 0 exp ( − sA ) G ( s ) ds   + G ( t ) = AX ( t ) + G ( t ) .  We now give several applications of this result in the case of the periodically forced harmonic oscillator. Assume first that we have a damped oscillator that is forced by cos t , so the period of the forcing term is 2 π . The system is X ′ = AX + G ( t ) , where G ( t ) = ( 0, cos t ) and A is the matrix A = ( 0 1 − k − b ) with b , k > 0. We claim that there is a unique periodic solution of this system which has period 2 π . To prove this, we must first find a solution X ( t ) satisfy- ing X ( 0 ) = X 0 = X ( 2 π) . By variation of parameters, we need to find X 0 such that X 0 = exp ( 2 π A ) X 0 + exp ( 2 π A ) 2 π ∫ 0 exp ( − sA ) G ( s ) ds . Now the term exp ( 2 π A ) 2 π ∫ 0 exp ( − sA ) G ( s ) ds is a constant vector that we denote by W . Therefore we must solve the equation ( exp ( 2 π A ) − I ) X 0 = − W . There is a unique solution to this equation, since the matrix exp ( 2 π A ) − I is invertible. For if this matrix were not invertible, there would be a nonzero vector V with ( exp ( 2 π A ) − I ) V = 0, or, in other words, the matrix exp ( 2 π A ) would have an eigenvalue 1. But, from the previous section, the eigenvalues of exp ( 2 π A ) are given by exp ( 2 πλ j ) ,

--- PAGE 148 ---
6.5 Nonautonomous Linear Systems 133 where the λ j are the eigenvalues of A . But each λ j has real part less than 0, so the magnitude of exp ( 2 πλ j ) is smaller than 1. Thus the matrix exp ( 2 π A ) − I is indeed invertible, and the unique initial value leading to a 2 π -periodic solution is X 0 = ( exp ( 2 π A ) − I ) − 1 ( − W ) . So let X ( t ) be this periodic solution with X ( 0 ) = X 0 . This solution is called the steady-state solution. If Y 0 is any other initial condition, then we may write Y 0 = ( Y 0 − X 0 ) + X 0 , so the solution through Y 0 is given by Y ( t ) = exp ( tA )( Y 0 − X 0 ) + exp ( tA ) X 0 + exp ( tA ) t ∫ 0 exp ( − sA ) G ( s ) ds = exp ( tA )( Y 0 − X 0 ) + X ( t ) . The first term in this expression tends to 0 as t → ∞ , since it is a solution of the homogeneous equation. Thus every solution of this system tends to the steady state solution as t → ∞ . Physically, this is clear: The motion of the damped (and unforced) oscillator tends to equilibrium, leaving only the motion due to the periodic forcing. We have proved the following theorem. Theorem. Consider the forced, damped harmonic oscillator equation x ′′ + bx ′ + kx = cos t with k , b > 0 . Then all solutions of this equation tend to the steady-state solution, which is periodic with period 2 π .  Now consider a particular example of a forced, undamped harmonic oscillator X ′ = ( 0 1 − 1 0 ) X + ( 0 cos ω t ) , where the period of the forcing is now 2 π/ω with ω 6 = ± 1. Let A = ( 0 1 − 1 0 ) . The solution of the homogeneous equation is X ( t ) = exp ( tA ) X 0 = ( cos t sin t − sin t cos t ) X 0 .

--- PAGE 149 ---
134 Chapter 6 Higher-Dimensional Linear Systems Variation of parameters provides a solution of the nonhomogeneous equation starting at the origin: Y ( t ) = exp ( tA ) t ∫ 0 exp ( − sA ) ( 0 cos ω s ) ds = exp ( tA ) t ∫ 0 ( cos s − sin s sin s cos s ) ( 0 cos ω s ) ds = exp ( tA ) t ∫ 0 ( − sin s cos ω s cos s cos ω s ) ds = 1 2 exp ( tA ) t ∫ 0 ( sin (ω − 1 ) s − sin (ω + 1 ) s cos (ω − 1 ) s + cos (ω + 1 ) s ) ds . Recalling that exp ( tA ) = ( cos t sin t − sin t cos t ) and using the fact that ω 6 = ± 1, evaluation of this integral plus a long computation yields Y ( t ) = 1 2 exp ( tA )    − cos (ω − 1 ) t ω − 1 + cos (ω + 1 ) t ω + 1 sin (ω − 1 ) t ω − 1 + sin (ω + 1 ) t ω + 1    + exp ( tA ) ( (ω 2 − 1 ) − 1 0 ) = 1 ω 2 − 1 ( − cos ω t ω sin ω t ) + exp ( tA ) ( (ω 2 − 1 ) − 1 0 ) . Thus the general solution of this equation is Y ( t ) = exp ( tA ) ( X 0 + ( (ω 2 − 1 ) − 1 0 )) + 1 ω 2 − 1 ( − cos ω t ω sin ω t ) . The first term in this expression is periodic with period 2 π while the second has period 2 π/ω . Unlike the damped case, this solution does not necessarily yield a periodic motion. Indeed, this solution is periodic if and only if ω is a

--- PAGE 150 ---
Exercises 135 rational number. If ω is irrational, the motion is quasiperiodic, just as we saw in Section 6.2. E X E R C I S E S 1. Find the general solution for X ′ = AX where A is given by ( a )   0 0 1 0 1 0 1 0 0   ( b )   1 0 1 0 1 0 1 0 1   ( c )   0 1 0 − 1 0 0 1 1 1   ( d )   0 1 0 1 0 0 1 1 1   ( e )   1 0 1 0 1 0 0 0 1   ( f )   1 1 0 1 1 1 0 1 1   ( g )   1 0 − 1 − 1 1 − 1 0 0 1   ( h )     1 0 0 1 0 1 1 0 0 0 1 0 1 0 0 0     2. Consider the linear system X ′ =   λ 1 0 0 0 λ 2 0 0 0 λ 3   X , where λ 3 < λ 2 < λ 1 < 0. Describe how the solution through an arbitrary initial value tends to the origin. 3. Give an example of a 3 × 3 matrix A for which all nonequilibrium solutions of X ′ = AX are periodic with period 2 π . Sketch the phase portrait. 4. Find the general solution of X ′ =     0 1 1 0 − 1 0 0 1 0 0 0 1 0 0 − 1 0     X . 5. Consider the system X ′ =   0 0 a 0 b 0 a 0 0   X , depending on the two parameters a and b .

--- PAGE 151 ---
136 Chapter 6 Higher-Dimensional Linear Systems (a) Find the general solution of this system. (b) Sketch the region in the ab -plane where this system has different types of phase portraits. 6. Consider the system X ′ =   a 0 b 0 b 0 − b 0 a   X , depending on the two parameters a and b . (a) Find the general solution of this system. (b) Sketch the region in the ab -plane where this system has different types of phase portraits. 7. Coupled Harmonic Oscillators. In this series of exercises you are asked to generalize the material on harmonic oscillators in Section 6.2 to the case where the oscillators are coupled . Suppose there are two masses m 1 and m 2 attached to springs and walls as shown in Figure 6.10. The springs connecting m j to the walls both have spring constants k 1 , while the spring connecting m 1 and m 2 has spring constant k 2 . This coupling means that the motion of either mass affects the behavior of the other. Let x j denote the displacement of each mass from its rest position, and assume that both masses are equal to 1. The differential equations for these coupled oscillators are then given by x ′′ 1 = − ( k 1 + k 2 ) x 1 + k 2 x 2 x ′′ 2 = k 2 x 1 − ( k 1 + k 2 ) x 2 . These equations are derived as follows. If m 1 is moved to the right ( x 1 > 0), the left spring is stretched and exerts a restorative force on m 1 given by − k 1 x 1 . Meanwhile, the central spring is compressed, so it exerts a restorative force on m 1 given by − k 2 x 1 . If the right spring is stretched, then the central spring is compressed and exerts a restorative force on m 1 given by k 2 x 2 (since x 2 < 0). The forces on m 2 are similar. (a) Write these equations as a first-order linear system. (b) Determine the eigenvalues and eigenvectors of the corresponding matrix. k 1 k 2 k 1 m 1 m 2 Figure 6.10 A coupled oscillator.

--- PAGE 152 ---
Exercises 137 (c) Find the general solution. (d) Let ω 1 = √ k 1 and ω 2 = √ k 1 + 2 k 2 . What can be said about the periodicity of solutions relative to the ω j ? Prove this. 8. Suppose X ′ = AX , where A is a 4 × 4 matrix with eigenvalues that are ± i √ 2 and ± i √ 3. Describe this flow. 9. Suppose X ′ = AX , where A is a 4 × 4 matrix with eigenvalues that are ± i and − 1 ± i . Describe this flow. 10. Suppose X ′ = AX , where A is a 4 × 4 matrix with eigenvalues that are ± i and ± 1. Describe this flow. 11. Consider the system X ′ = AX , where X = ( x 1 , . . . , x 6 ) , A =         0 ω 1 − ω 1 0 0 ω 2 − ω 2 0 − 1 1         , and ω 1 /ω 2 is irrational. Describe qualitatively how a solution behaves when, at time 0, each x j is nonzero with the exception that (a) x 6 = 0 (b) x 5 = 0 (c) x 3 = x 4 = x 5 = 0 (d) x 3 = x 4 = x 5 = x 6 = 0 12. Compute the exponentials of the following matrices: ( a ) ( 5 − 6 3 − 4 ) ( b ) ( 2 − 1 1 2 ) ( c ) ( 2 − 1 0 2 ) ( d ) ( 0 1 1 0 ) ( e )   0 1 2 0 0 3 0 0 0   ( f )   2 0 0 0 3 0 0 1 3   ( g )   λ 0 0 1 λ 0 0 1 λ   ( h ) ( i 0 0 − i ) ( i ) ( 1 + i 0 2 1 + i ) ( j )     1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0     13. Find an example of two matrices A , B such that exp ( A + B ) 6 = exp ( A ) exp ( B ) .

--- PAGE 153 ---
138 Chapter 6 Higher-Dimensional Linear Systems 14. Show that if AB = BA , then (a) exp ( A ) exp ( B ) = exp ( B ) exp ( A ) (b) exp ( A ) B = B exp ( A ) 15. Consider the triplet of harmonic oscillators x ′′ 1 = − x 1 x ′′ 2 = − 2 x 2 x ′′ 3 = − ω 2 x 3 , where ω is irrational. What can you say about the qualitative behavior of solutions of this six-dimensional system?

--- PAGE 154 ---
7 Nonlinear Systems In this chapter we begin the study of nonlinear differential equations. In linear (constant coefficient) systems we can always find the explicit solution of any initial value problem; however, this is rarely the case for nonlinear systems. In fact, basic properties such as the existence and uniqueness of solutions, which was so obvious in the linear case, no longer hold for nonlinear systems. As we shall see, some nonlinear systems have no solutions whatsoever to a given initial value problem. On the other hand, there are other systems that have infinitely many differ- ent such solutions. Even if we do find a solution of such a system, this solution need not be defined for all time; for example, the solution may tend to ∞ in finite time. Other questions also arise: For example, what happens if we vary the initial condition of a system ever so slightly? Does the corresponding solu- tion vary continuously? All of this is clear for linear systems, but not at all clear in the nonlinear case. This means that the underlying theory behind nonlinear systems of differential equations is quite a bit more complicated than that for linear systems. In practice, most nonlinear systems that arise are “nice” in the sense that we do have existence and uniqueness of solutions, as well as continuity of solutions when initial conditions are varied and other “natural” properties. Thus we have a choice: Given a nonlinear system, we could simply plunge ahead and either hope that or, if possible, verify that, in each specific case, the system’s solutions behave nicely. Alternatively, we could take a long pause at Differential Equations, Dynamical Systems, and an Introduction to Chaos. DOI: 10.1016/B978-0-12-382010-5.00007-5 c © 2013 Elsevier Inc. All rights reserved. 139

--- PAGE 155 ---
140 Chapter 7 Nonlinear Systems this stage to develop the necessary hypotheses that guarantee that solutions of a given nonlinear system behave nicely. In this book we pursue a compromise route. In this chapter, we spell out in precise detail many of the theoretical results that govern the behavior of solu- tions of differential equations. We present examples of how and when these results fail, but we will not prove these theorems here. Rather, we will postpone all of the technicalities until Chapter 17, primarily because understanding this material demands a firm and extensive background in the principles of real analysis. In subsequent chapters, we will make use of the results stated here, but readers who are primarily interested in applications of differential equations or in understanding how specific nonlinear systems may be ana- lyzed need not get bogged down in these details here. Readers who want the technical details may take a detour to Chapter 17 now. 7.1 Dynamical Systems As mentioned previously, most nonlinear systems of differential equations are impossible to solve analytically. One reason for this is the unfortunate fact that we simply do not have enough functions with specific names that we can use to write down explicit solutions of these systems. Equally problematic is the fact that, as we shall see, higher-dimensional systems may exhibit chaotic behav- ior, a property that makes knowing a particular explicit solution essentially worthless in the larger scheme of understanding the behavior of the system. Thus, to begin to understand these systems we are forced to resort to different means. These are the techniques that arise in the field of dynamical systems. We will use a combination of analytic, geometric, and topological techniques to derive rigorous results about the behavior of solutions of these equations. We begin by collecting together some of the terminology regarding dynam- ical systems that we have introduced at various points in the preceding chapters. A dynamical system is a way of describing the passage in time of all points of a given space S . The space S could be thought of, for example, as the space of states of some physical system. Mathematically, S might be a Euclidean space or an open subset of Euclidean space or some other space such as a surface in R 3 . When we consider dynamical systems that arise in mechanics, the space S will be the set of possible positions and velocities of the system. For the sake of simplicity, we will assume throughout that the space S is Euclidean space R n , although in certain cases the important dynamical behavior will be confined to a particular subset of R n . Given an initial position X ∈ R n , a dynamical system on R n tells us where X is located 1 unit of time later, 2 units of time later, and so on. We denote these

--- PAGE 156 ---
7.1 Dynamical Systems 141 new positions of X by X 1 , X 2 , and so forth. At time zero, X is located at posi- tion X 0 . One unit before time zero, X was at X − 1 . In general the “trajectory” of X is given by X t . If we measure the positions X t using only integer time values, we have an example of a discrete dynamical system, which we shall study in Chapter 15. If time is measured continuously with t ∈ R , we have a continuous dynamical system. If the system depends on time in a continuously differen- tiable manner, we have a smooth dynamical system. These are the three princi- pal types of dynamical systems that arise in the study of systems of differential equations, and they will form the backbone of Chapters 8 through 14. The function that takes t to X t yields either a sequence of points or a curve in R n that represents the life history of X as time runs from −∞ to ∞ . Dif- ferent branches of dynamical systems make different assumptions about how the function X t depends on t . For example, ergodic theory deals with such functions under the assumption that they preserve a measure on R n . Topolog- ical dynamics deals with such functions under the assumption that X t varies only continuously. In the case of differential equations, we will usually assume that the function X t is continuously differentiable. The map φ t : R n → R n that takes X into X t is defined for each t and, from our interpretation of X t as a state moving in time, it is reasonable to expect φ t to have φ − t as its inverse. Also, φ 0 should be the identity function φ 0 ( X ) = X , and φ t (φ s ( X )) = φ t + s ( X ) is also a natural condition. We formalize all of this in the following definition. Definition A smooth dynamical system on R n is a continuously differentiable function φ : R × R n → R n , where φ( t , X ) = φ t ( X ) satisfies 1. φ 0 : R n → R n is the identity function: φ 0 ( X 0 ) = X 0 . 2. The composition φ t ◦ φ s = φ t + s for each t , s ∈ R . Recall that a function is continuously differentiable if all of its partial deriva- tives exist and are continuous throughout its domain. It is traditional to call a continuously differentiable function a C 1 function. If the function is k times continuously differentiable, it is called a C k function. Note that the preceding definition implies that the map φ t : R n → R n is C 1 for each t and has a C 1 inverse φ − t (take s = − t in part 2). Example. For the first-order differential equation x ′ = ax , the function φ t ( x 0 ) = x 0 exp ( at ) gives the solutions of this equation and also defines a smooth dynamical system on R . 

--- PAGE 157 ---
142 Chapter 7 Nonlinear Systems Example. Let A be an n × n matrix. Then the function φ t ( X 0 ) = exp ( tA ) X 0 defines a smooth dynamical system on R n . Clearly, φ 0 = exp ( 0 ) = I and, as we saw in the previous chapter, we have φ t + s = exp (( t + s ) A ) = ( exp ( tA ))( exp ( sA )) = φ t ◦ φ s .  Note that these examples are intimately related to the system of differential equations X ′ = AX . In general, a smooth dynamical system always yields a vector field on R n via this rule: Given φ t , let F ( X ) = d dt ∣ ∣ ∣ ∣ t = 0 φ t ( X ) . Then φ t is just the time t map associated with the flow of X ′ = F ( X ) . Conversely, the differential equation X ′ = F ( X ) generates a smooth dynam- ical system provided the time t map of the flow is well defined and continu- ously differentiable for all time. Unfortunately, this is not always the case. 7.2 The Existence and Uniqueness Theorem We turn now to the fundamental theorem of differential equations, the Exis- tence and Uniqueness Theorem. Consider the system of differential equations X ′ = F ( X ) , where F : R n → R n . Recall that a solution of this system is a function X : J → R n defined on some interval J ⊂ R such that, for all t ∈ J , X ′ ( t ) = F ( X ( t )) . Geometrically, X ( t ) is a curve in R n with a tangent vector X ′ ( t ) that exists for all t ∈ J and equals F ( X ( t )) . As in previous chapters, we think of this vector as being based at X ( t ) , so that the map F : R n → R n defines a vector field on R n . An initial condition or initial value for a solution X : J → R n is a specifi- cation of the form X ( t 0 ) = X 0 , where t 0 ∈ J and X 0 ∈ R n . For simplicity, we usually take t 0 = 0. The main problem in differential equations is to find the solution of any initial value problem —that is, to determine the solution that of the system that satisfies the initial condition X ( 0 ) = X 0 for each X 0 ∈ R n . Unfortunately, nonlinear differential equations may have no solutions that satisfy certain initial conditions.

--- PAGE 158 ---
7.2 The Existence and Uniqueness Theorem 143 Example. Consider the simple first-order differential equation x ′ = { 1 if x < 0 − 1 if x ≥ 0. This vector field on R points to the left when x ≥ 0 and to the right if x < 0. Consequently, there is no solution that satisfies the initial condition x ( 0 ) = 0. Indeed, such a solution must initially decrease since x ′ ( 0 ) = − 1, but for all negative values of x , solutions must increase. This cannot happen. Note fur- ther that solutions are never defined for all time. For example, if x 0 > 0, then the solution through x 0 is given by x ( t ) = x 0 − t , but this solution is only valid for −∞ < t < x 0 for the same reason as before. The problem in this example is that the vector field is not continuous at 0; whenever a vector field is discontinuous we face the possibility that nearby vectors may point in “opposing” directions, thereby causing solutions to halt at these bad points.  Beyond the problem of existence of solutions of nonlinear differential equa- tions, we also must confront the fact that certain equations may have many different solutions to the same initial value problem. Example. Consider the differential equation x ′ = 3 x 2 / 3 . The identically zero function u : R → R given by u ( t ) ≡ 0 is clearly a solution with initial condition u ( 0 ) = 0. But u 0 ( t ) = t 3 is also a solution satisfying this initial condition. Moreover, for any τ > 0, the function given by u τ ( t ) = { 0 if t ≤ τ ( t − τ ) 3 if t > τ is also a solution satisfying the initial condition u τ ( 0 ) = 0. Although the dif- ferential equation in this example is continuous at x 0 = 0, the problems arise because x 2 / 3 is not differentiable at this point.  From these two examples it is clear that, to ensure existence and uniqueness of solutions, certain conditions must be imposed on the function F . In the first example, F was not continuous at the problematic point 0, while in the sec- ond example, F failed to be differentiable at 0. It turns out that the assumption that F is continuously differentiable is sufficient to guarantee both existence and uniqueness of solutions, as we shall see. Fortunately, differential equa- tions that are not continuously differentiable rarely arise in applications, so

--- PAGE 159 ---
144 Chapter 7 Nonlinear Systems the phenomenon of nonexistence or nonuniqueness of solutions with given initial conditions is quite exceptional. The following is the fundamental local theorem of ordinary differential equations. The important proof of this theorem is contained in Chapter 17. The Existence and Uniqueness Theorem. Consider the initial value problem X ′ = F ( X ) , X ( t 0 ) = X 0 , where X 0 ∈ R n . Suppose that F : R n → R n is C 1 . Then, first, there exists a solu- tion of this initial value problem, and second, this is the only such solution. More precisely, there exists an a > 0 and a unique solution, X : ( t 0 − a , t 0 + a ) → R n , of this differential equation satisfying the initial condition X ( t 0 ) = X 0 .  Without dwelling on the details here, the proof of this theorem depends on an important technique known as Picard iteration . Before moving on, we illustrate how the Picard iteration scheme used in the proof of the theorem works in several special examples. The basic idea behind this iterative process is to construct a sequence of functions that converges to the solution of the differential equation. The sequence of functions u k ( t ) is defined inductively by u 0 ( t ) = x 0 , where x 0 is the given initial condition, and then u k + 1 ( t ) = x 0 + t ∫ 0 F ( u k ( s )) ds . Example. Consider the simple differential equation x ′ = x . We will produce the solution of this equation satisfying x ( 0 ) = x 0 . We know, of course, that this solution is given by x ( t ) = x 0 e t . We will construct a sequence of functions u k ( t ) , one for each k , that converges to the actual solution x ( t ) as k → ∞ . We start with u 0 ( t ) = x 0 , the given initial value. Then we set u 1 ( t ) = x 0 + t ∫ 0 F ( u 0 ( s )) ds = x 0 + t ∫ 0 x 0 ds ,

--- PAGE 160 ---
7.2 The Existence and Uniqueness Theorem 145 so that u 1 ( t ) = x 0 + tx 0 . Given u 1 we define u 2 ( t ) = x 0 + t ∫ 0 F ( u 1 ( s )) ds = x 0 + t ∫ 0 ( x 0 + sx 0 ) ds , so that u 2 ( t ) = x 0 + tx 0 + t 2 2 x 0 . You can probably see where this is heading. Inductively, we set u k + 1 ( t ) = x 0 + t ∫ 0 F ( u k ( s )) ds , and so u k + 1 ( t ) = x 0 k + 1 ∑ i = 0 t i i ! . As k → ∞ , u k ( t ) converges to x 0 ∞ ∑ i = 0 t i i ! = x 0 e t = x ( t ) , which is the solution of our original equation.  Example. For an example of Picard iteration applied to a system of differential equations, consider the linear system X ′ = F ( X ) = ( 0 1 − 1 0 ) X with initial condition X ( 0 ) = ( 1, 0 ) . As we have seen, the solution of this initial value problem is X ( t ) = ( cos t − sin t ) . Using Picard iteration, we have U 0 ( t ) = ( 1 0 ) U 1 ( t ) = ( 1 0 ) + t ∫ 0 F ( 1 0 ) ds = ( 1 0 ) + t ∫ 0 ( 0 − 1 ) ds = ( 1 − t )