odq

NonlinearOrdinaryDifferentialEquations by Peter J. Olver University of Minnesota 1. Introduction. These notes are concerned with initial value problems for sy stems of ordinary dif- ferential equations. Here our emphasis will be on nonlinear phenomena and properties, particularly those with physical relevance. Finding a solu tion to a differential equation may not be so important if that solution never appears in the p hysical model represented by the system, or is only realized in exceptional circumstan ces. Thus, equilibrium solu- tions, which correspond to configurations in which the physi cal system does not move, only occur in everyday situations if they are stable. An unst able equilibrium will not ap- pear in practice, since slight perturbations in the system o r its physical surroundings will immediately dislodge the system far away from equilibrium. Of course, very few nonlinear systems can be solved explicit ly, and so one must typ- ically rely on a numerical scheme to accurately approximate the solution. Basic methods for initial value problems, beginning with the simple Euler scheme, and working up to the extremely popular Runge–Kutta fourth order method, wil l be the subject of the final section of the chapter. However, numerical schemes do not al ways give accurate results, and we breifly discuss the class of stiff differential equation s, which present a more serious challenge to numerical analysts. Without some basic theoretical understanding of the nature of solutions, equilibrium points, and stability properties, one would not be able to un derstand when numerical so- lutions (even those provided by standard well-used package s) are to be trusted. Moreover, when testing a numerical scheme, it helps to have already ass embled a repertoire of nonlin- ear problems in which one already knows one or more explicit a nalytic solutions. Further tests and theoretical results can be based on first integrals (also known as conservation laws) or, more generally, Lyapunov functions. Although we h ave only space to touch on these topics briefly, but, we hope, this will whet the reader’ s appetite for delving into this subject in more depth. The references [ 2,9,13,15,17] can be profitably consulted. 2. First Order Systems of Ordinary Differential Equations. Let us begin by introducing the basic object of study in discr ete dynamics: the initial value problem for a first order system of ordinary differentia l equations. Many physical applications lead to higher order systems of ordinary differ ential equations, but there is a simple reformulation that will convert them into equivalen t first order systems. Thus, we do not lose any generality by restricting our attention to th e first order case throughout. Moreover, numerical solution schemes for higher order init ial value problems are entirely based on their reformulation as first order systems. 11/17/13 1 c/ci∇cleco√y∇t2013 Peter J. Olver Scalar Ordinary Differential Equations As always, when confronted with a new problem, it is essentia l to fully understand the simplest case first. Thus, we begin with a single scalar, fi rst order ordinary differential equation du dt=F(t,u). (2.1) In many applications, the independent variable trepresents time, and the unknown func- tionu(t) is some dynamical physical quantity. Throughout this chap ter, all quantities are assumed to be real. (Results on complex ordinary differen tial equations can be found in [14].) Under appropriate conditions on the right hand side (to b e formalized in the following section), the solution u(t) is uniquely specified by its value at a single time, u(t0) =u0. (2.2) The combination (2.1–2) is referred to as an initial value problem , and our goal is to devise both analytical and numerical solution strategies. A differential equation is called autonomous if the right hand side does not explicitly depend upon the time variable: du dt=F(u). (2.3) All autonomous scalar equations can be solved by direct inte gration. We divide both sides byF(u), whereby 1 F(u)du dt= 1, and then integrate with respect to t; the result is /integraldisplay1 F(u)du dtdt=/integraldisplay dt=t+k, wherekis the constant of integration. The left hand integral can be evaluated by the change of variables that replaces tbyu, whereby du= (du/dt)dt, and so /integraldisplay1 F(u)du dtdt=/integraldisplaydu F(u)=G(u), whereG(u) indicates a convenient anti-derivative†of the function 1 /F(u). Thus, the solution can be written in implicit form G(u) =t+k. (2.4) If we are able to solve the implicit equation (2.4), we may the reby obtain the explicit solution u(t) =H(t+k) (2 .5) †Technically, a second constant of integration should appear here, but th is can be absorbed into the previous constant k, and so proves to be unnecessary. 11/17/13 2 c/ci∇cleco√y∇t2013 Peter J. Olver 0.5 1 1.5 2 -1-0.50.511.52 Figure 1. Solutions to/squaresmallsolidu=u2. in terms of the inverse function H=G−1. Finally, to satisfy the initial condition (2.2), we sett=t0in the implicit solution formula (2.4), whereby G(u0) =t0+k. Therefore, the solution to our initial value problem is G(u)−G(u0) =t−t0,or, explicitly, u(t) =H/parenleftbig t−t0+G(u0)/parenrightbig .(2.6) Remark: A more direct version of this solution technique is to rewri te the differential equation (2.3) in the “separated form” du F(u)=dt, in which all terms involving u, including its differential du, are collected on the left hand side of the equation, while all terms involving tand its differential are placed on the right, and then formally integrate both sides, leading to the same i mplicit solution formula: G(u) =/integraldisplaydu F(u)=/integraldisplay dt=t+k. (2.7) Before completing our analysis of this solution method, let us run through a couple of elementary examples. Example 2.1. Consider the autonomous initial value problem du dt=u2, u (t0) =u0. (2.8) To solve the differential equation, we rewrite it in the separ ated form du u2=dt,and then integrate both sides: −1 u=/integraldisplaydu u2=t+k. 11/17/13 3 c/ci∇cleco√y∇t2013 Peter J. Olver Solving the resulting algebraic equation for u, we deduce the solution formula u=−1 t+k. (2.9) To specify the integration constant k, we evaluate uat the initial time t0; this implies u0=−1 t0+k,so that k=−1 u0−t0. Therefore, the solution to the initial value problem is u=u0 1−u0(t−t0). (2.10) Figure 1 shows the graphs of some typical solutions. Astapproaches the critical value t⋆=t0+1/u0from below, the solution “blows up”, meaning u(t)→ ∞ast→t⋆. The blow-up time t⋆depends upon the initial data — the largeru0>0 is, the sooner the solution goes off to infinity. If the initia l data is negative, u0<0, the solution is well-defined for all t > t0, but has a singularity in the past, at t⋆=t0+1/u0< t0. The only solution that exists for all positive and negative time is the constant solution u(t)≡0, corresponding to the initial condition u0= 0. In general, the constant equilibrium solutions to an autonomous ordinary differential equation, also known as its fixed points , play a distinguished role. If u(t)≡u⋆is a constant solution, then du/dt≡0, and hence the differential equation (2.3) implies that F(u⋆) = 0. Therefore, the equilibrium solutions coincide with the rootsof the function F(u). In point of fact, since we divided by F(u), the derivation of our formula for the solution (2.7) assumed that we were notat an equilibrium point. In the preceding example, our final solutionformula (2.10)happens toincludetheequilibrium solutionu(t)≡0, corresponding tou0= 0, but this is a lucky accident. Indeed, the equilibrium sol ution does notappear in the “general” solution formula (2.9). One must typically ta ke extra care that equilibrium solutions do not elude us when utilizing this basic integrat ion method. Example 2.2. Although a population of people, animals, or bacteria consi sts of individuals, the aggregatebehavior can often be effectivel y modeled by a dynamical system that involves continuously varying variables. As first prop osed by the English economist Thomas Malthus in 1798, the population of a species grows, ro ughly, in proportion to its size. Thus, the number of individuals N(t) at time tsatisfies a first order differential equation of the form dN dt=ρN, (2.11) where the proportionality factor ρ=β−δmeasures the rate of growth, namely the difference between the birth rate β≥0 and the death rate δ≥0. Thus, if births exceed deaths,ρ >0, and the population increases, whereas if ρ <0, more individuals are dying and the population shrinks. In the very simplest model, the growth rate ρis assumed to be independent of the population size, and (2.11) reduces to a simple linear ordin ary differential equation whose 11/17/13 4 c/ci∇cleco√y∇t2013 Peter J. Olver solutions satisfy the Malthusian exponential growth law N(t) =N0eρt, whereN0=N(0) is the initial population size. Thus, if ρ >0, the population grows without limit, while if ρ <0, the population dies out, so N(t)→0 ast→ ∞, at an exponentially fast rate. The Malthusian population model provides a reasonably accurat e description of the behavior of an isolated population in an environment with unlimited r esources. Inamorerealisticscenario, thegrowthratewilldependupo nthesizeofthepopulation aswellasexternal environmental factors. For example, int hepresence oflimitedresources, relatively small populations will increase, whereas an exc essively large population will have insufficient resources to survive, and so its growth rate will be negative. In other words, the growth rate ρ(N)>0 whenN < N⋆, whileρ(N)<0 whenN > N⋆, where the carrying capacity N⋆>0 depends upon the resource availability. The simplest clas s of functions that satifies these two inequalities are of the for mρ(N) =µ(N⋆−N), where µ >0 is a positive constant. This leads us to the nonlinear popul ation model dN dt=µN(N⋆−N). (2.12) In deriving this model, we assumed that the environment is no t changing over time; a dynamical environment would require a more complicated non -autonomous differential equation. Before analyzing the solutions to the nonlinear population model, let us make a pre- liminary change of variables, and set u(t) =N(t)/N⋆, so that urepresents the size of the population in proportion to the carrying capacity N⋆. A straightforward computation shows that u(t) satisfies the so-called logistic differential equation du dt=λu(1−u), u (0) =u0, (2.13) whereλ=N⋆µ, and, for simplicity, we assign the initial time to be t0= 0. The logistic differentialequationcanbeviewedasthecontinuouscounte rpartofthelogisticmapstudied in my Notes on Nonlinear Systems. However, unlike its discre te namesake, the logistic differential equation is quite sedate, and its solutions eas ily understood. First, there are two equilibrium solutions: u(t)≡0 andu(t)≡1, obtained by setting the right hand side of the equation equal to zero. The first rep resents a nonexistent populationwithnoindividualsandhencenoreproduction. T hesecondequilibriumsolution corresponds to a static population N(t)≡N⋆that is at the ideal size for the environment, so deaths exactly balance births. In all other situations, t he population size will vary over time. To integrate the logistic differential equation, we proceed as above, first writing it in the separated form du u(1−u)=λdt. Integrating both sides, and using partial fractions, λt+k=/integraldisplaydu u(1−u)=/integraldisplay/bracketleftbigg1 u+1 1−u/bracketrightbigg du= log/vextendsingle/vextendsingle/vextendsingle/vextendsingleu 1−u/vextendsingle/vextendsingle/vextendsingle/vextendsingle, 11/17/13 5 c/ci∇cleco√y∇t2013 Peter J. Olver 2 4 6 8 10 -1-0.50.511.52 Figure 2. Solutions to u′=u(1−u). wherekis a constant of integration. Therefore u 1−u=ceλt,where c=±ek. Solving for u, we deduce the solution u(t) =ceλt 1+ceλt. (2.14) The constant of integration is fixed by the initial condition . Solving the algebraic equation u0=u(0) =c 1+cyields c=u0 1−u0. Substituting the result back into the solution formula (2.1 4) and simplifying, we find u(t) =u0eλt 1−u0+u0eλt. (2.15) The resulting solutions are illustrated in Figure 2. Intere stingly, while the equilibrium solutions are not covered by the integration method, they re appear in the final solution formula, corresponding to initial data u0= 0 and u0= 1 respectively. However, this is a lucky accident, and cannot be anticipated in more complicat ed situations. When using the logistic equation to model population dynami cs, the initial data is assumedtobepositive, u0>0. Astime t→ ∞,thesolution(2.15)tendstotheequilibrium valueu(t)→1 — which corresponds to N(t)→N⋆approaching the carrying capacity in the original population model. For small initial values u0≪1 the solution initially grows at an exponential rate λ, corresponding to a population with unlimited resources. However, as the population increases, the gradual lack of re sources tends to slow down 11/17/13 6 c/ci∇cleco√y∇t2013 Peter J. Olver the growth rate, and eventually the population saturates at the equilibrium value. On the other hand, if u0>1, the population is too large to be sustained by the availabl e resources, and so dies off until it reaches the same saturatio n value. If u0= 0, then the solution remains at equilibrium u(t)≡0. Finally, when u0<0, the solution only exists for a finite amount of time, with u(t)−→ −∞ as t−→t⋆=1 λlog/parenleftbigg 1−1 u0/parenrightbigg . Of course, thisfinal case does appear in the physical world, s ince we cannot have a negative population! The separation of variables method used to solve autonomous equations can be stra- ightforwardly extended to a special class of non-autonomou s equations. A separable ordi- nary differential equation has the form du dt=a(t)F(u), (2.16) in which the right hand side is the product of a function of tand a function of u. To solve the equation, we rewrite it in the separated form du F(u)=a(t)dt. Integrating both sides leads to the solution in implicit for m G(u) =/integraldisplaydu F(u)=/integraldisplay a(t)dt=A(t)+k. (2.17) The integration constant kis then fixed by the initial condition. And, as before, one mus t properly account for any equilibrium solutions, when F(u) = 0. Example 2.3. Let us solve the particular initial value problem du dt= (1−2t)u, u (0) = 1. (2.18) We begin by writing the differential equation in separated fo rm du u= (1−2t)dt. Integrating both sides leads to logu=/integraldisplaydu u=/integraldisplay (1−2t)dt=t−t2+k, wherekis the constant of integration. We can readily solve for u(t) =cet−t2, wherec=±ek. The latter formula constitutes the general solution to the differential equation, and happens to include the equilibrium solution u(t)≡0 whenc= 0. The given initial condition requires that c= 1, and hence u(t) =et−t2is the unique solution to the initial value problem. The solution is graphed in Figure 3. 11/17/13 7 c/ci∇cleco√y∇t2013 Peter J. Olver 0.5 11.5 22.5 30.250.50.7511.251.5 Figure 3. Solution to the Initial Value Problem/squaresmallsolidu= (1−2t)u,u(0) = 1. First Order Systems Afirst order system of ordinary differential equations has the general form du1 dt=F1(t,u1,...,un),···dun dt=Fn(t,u1,...,un).(2.19) The unknowns u1(t),...,un(t) are scalar functions of the real variable t, which usually represents time. We shall write the system more compactly in vector form du dt=F(t,u), (2.20) whereu(t) = (u1(t),...,un(t))T, andF(t,u) = (F1(t,u1,...,un),...,Fn(t,u1,...,un))T is a vector-valued function of n+ 1 variables. By a solution to the differential equation, we mean a vector-valued function u(t) that is defined and continuously differentiable on an interval a < t < b , and, moreover, satisfies the differential equation on its in terval of definition. Each solution u(t) serves to parametrize a curve C⊂Rn, also known as a trajectory ororbitof the system. In this chapter, we shall concentrate on initial value probl ems for such first order systems. The general initial conditions are u1(t0) =a1, u2(t0) =a2,··· un(t0) =an, (2.21) or, in vectorial form, u(t0) =a (2.22) Heret0is a prescribed initial time, while the vector a= (a1,a2,...,an)Tfixes the initial position of the desired solution. In favorable situations, as described below, the initial conditions serve to uniquely specify a solution to the differ ential equations — at least for nearby times. The general issues of existence and uniquenss of solutions will be addressed in the following section. 11/17/13 8 c/ci∇cleco√y∇t2013 Peter J. Olver A system of differential equations is called autonomous if the right hand side does not explicitly depend upon the time t, and so takes the form du dt=F(u). (2.23) One important class of autonomous first order systems are the steady state fluid flows. HereF(u) =vrepresents the fluid velocity vector field at the position u. The solution u(t) to the initial value problem (2.23,22) describes the motio n of a fluid particle that starts at position aat timet0. The differential equation tells us that the fluid velocity at each point on the particle’s trajectory matches the prescri bed vector field. Anequilibrium solution is constant: u(t)≡u⋆for allt. Thus, its derivative must vanish,du/dt≡0, and hence, every equilibrium solution arises as a solution to the system of algebraic equations F(u⋆) =0 (2.24) prescribed by the vanishing of the right hand side of the syst em (2.23). Example 2.4. Apredator-prey system isa simplified ecological model oftwo species: the predators which feed on the prey. For example, the predat ors might be lions roaming the Serengeti and the prey zebra. We let u(t) represent the number of prey, and v(t) the number of predators at time t. Both species obey a population growth model of the form (2.11), and so the dynamical equations can be written as du dt=ρu,dv dt=σv, (2.25) where the growth rates ρ,σmay depend upon the other species. The more prey, i.e., the largeruis, the faster the predators reproduce, while a lack of prey w ill cause them to die off. On the other hand, the more predators, the faster the prey are consumed and the slower their net rate of growth. If we assume that the environment has unlimited resources fo r the prey, which, bar- ring drought, is probably valid in the case of the zebras, the n the simplest model that incorporates these assumptions is the Lotka–Volterra system du dt=αu−δuv,dv dt=−βv+γuv, (2.26) corresponding to growth rates ρ=α−δv,σ=−β+γu. The parameters α,β,γ,δ > 0 are all positive, and their precise values will depend upon the species involved and how they interact, as indicated by field data, combined with, perhaps, educated guesses. In particular, αrepresents the unrestrained growth rate of the prey in the ab sence of predators, while −βrepresents the rate that the predators die off in the absence o f their prey. The nonlinear terms model the interaction of the two sp ecies: the rate of increase in the predators is proportional to the number of available p rey, while the rate of decrese in the prey is proportional to the number of predators. The in itial conditions u(t0) =u0, v(t0) =v0represent the initial populations of the two species. 11/17/13 9 c/ci∇cleco√y∇t2013 Peter J. Olver We will discuss the integration of the Lotka–Volterra syste m (2.26) in Section 4. Here, let us content ourselves with determining the possible equi libria. Setting the right hand sides of the system to zero leads to the nonlinear algebraic s ystem 0 =αu−δuv=u(α−δv),0 =−βv+γuv=v(−β+γu). Thus, there are two distinct equilibria, namely u⋆ 1=v⋆ 1= 0, u⋆ 2=β/γ, v⋆ 2=α/δ. Thefirstistheuninteresting(or,rathercatastropic)situ ationwheretherearenoanimals— no predators and no prey. The second is a nontrivial solution in which both populations maintain a steady value, for which the birth rate of the prey i s precisely sufficient to continuously feed the predators. Is this a feasible solutio n? Or, to state the question more mathematically, is this a stable equilibrium? We shall deve lop the tools to answer this question below. Higher Order Systems A wide variety of physical systems are modeled by nonlinear s ystems of differential equations depending upon second and, occasionally, even hi gher order derivatives of the unknowns. But there is an easy device that will reduce any hig her order ordinary differ- ential equation or system to an equivalent first order system . “Equivalent” means that each solution to the first order system uniquely corresponds to a solution to the higher order equation and vice versa. The upshot is that, for all pra ctical purposes, one only needs to analyze first order systems. Moreover, the vast majo rity of numerical solution algorithms are designed for first order systems, and so to num erically integrate a higher order equation, one must place it into an equivalent first ord er form. We have already encountered the main idea in our discussion o f the phase plane approach to second order scalar equations d2u dt2=F/parenleftbigg t,u,du dt/parenrightbigg . (2.27) Weintroduceanew dependent variable v=du dt. Sincedv dt=d2u dt2, thefunctions u,vsatisfy the equivalent first order system du dt=v,dv dt=F(t,u,v). (2.28) Conversely, it is easy to check that if u(t) = (u(t),v(t))Tis any solution to the first order system, then its first component u(t) defines a solution to the scalar equation, which establishes their equivalence. The basic initial conditio nsu(t0) =u0, v(t0) =v0, for the firstordersystemtranslateintoapairofinitialcondition su(t0) =u0,/squaresmallsolidu(t0) =v0,specifying the value of the solution and its first order derivative for th e second order equation. Similarly, given a third order equation d3u dt3=F/parenleftbigg t,u,du dt,d2u dt2/parenrightbigg , 11/17/13 10 c/ci∇cleco√y∇t2013 Peter J. Olver we set v=du dt, w =dv dt=d2u dt2. The variables u,v,wsatisfy the equivalent first order system du dt=v,dv dt=w,dw dt=F(t,u,v,w). The general technique should now be clear. Example 2.5. The forced van der Pol equation d2u dt2+(u2−1)du dt+u=f(t) (2 .29) arises in the modeling of an electrical circuit with a triode whose resistance changes with the current. It also arises in certain chemical reactions an d wind-induced motions of structures. To convert the van der Pol equation into an equiv alent first order system, we setv=du/dt, whence du dt=v,dv dt=f(t)−(u2−1)v−u, (2.30) is the equivalent phase plane system. Example 2.6. The Newtonian equations for a mass mmoving in a potential force field are a second order system of the form md2u dt2=−∇F(u) in which u(t) = (u(t),v(t),w(t))Trepresents the position of the mass and F(u) = F(u,v,w) the potential function. In components, md2u dt2=−∂F ∂u, md2v dt2=−∂F ∂v, md2w dt2=−∂F ∂w.(2.31) Forexample, a planetmovinginthesun’s gravitationalfield satisfies theNewtoniansystem for the gravitational potential F(u) =−α /ba∇dblu/ba∇dbl=−α√ u2+v2+w2, (2.32) whereαdepends on the masses and the universal gravitational const ant. (This simplified model ignores any additional interplanetary forces.) Thus , the mass’ motion in such a gravitational force field follows the solution to the second order Newtonian system md2u dt2=−∇F(u) =−αu /ba∇dblu/ba∇dbl3=α (u2+v2+w2)3/2 u v w . 11/17/13 11 c/ci∇cleco√y∇t2013 Peter J. Olver The same system of ordinary differential equations describe s the motion of a charged particle in a Coulomb electric force field, where the sign of αis positive for attracting opposite charges, and negative for repelling like charges. To convert the second order Newton equations into a first orde r system, we set v=/squaresmallsolidu to be the mass’ velocity vector, with components p=du dt, q=dv dt, r=dw dt, and so du dt=p,dv dt=q,dw dt=r, (2.33) dp dt=−1 m∂F ∂u(u,v,w),dq dt=−1 m∂F ∂v(u,v,w),dr dt=−1 m∂F ∂w(u,v,w). One of Newton’s greatest acheivements was to solve this syst em in the case of the central gravitationalpotential(2.32),andtherebyconfirm theval idityofKepler’slawsofplanetary motion. Finally, we note that there is a simple device that will conve rt any non-autonomous system into an equivalent autonomous system involving one a dditional variable. Namely, one introduces an extra coordinate u0=tto represent the time, which satisfies the el- ementary differential equation du0/dt= 1 with initial condition u0(t0) =t0. Thus, the original system (2.19) can be written in the autonomous form du0 dt= 1,du1 dt=F1(u0,u1,...,un),···dun dt=Fn(u0,u1,...,un).(2.34) For example, the autonomous form of the forced van der Pol sys tem (2.30) is du0 dt= 1,du1 dt=u2,du2 dt=f(u0)−(u2 1−1)u2−u1,(2.35) in which u0represents the time variable. 3. Existence, Uniqueness, and Continuous Dependence. It goes without saying that there is no general analytical me thod that will solve all differential equations. Indeed, even relatively simple firs t order, scalar, non-autonomous ordinary differential equations cannot be solved in closed f orm. For example, the solution to the particular Riccati equation du dt=u2+t (3.1) cannot be written in terms of elementary functions, althoug h it can be solved in terms of Airy functions, [ 25]. TheAbel equation du dt=u3+t (3.2) 11/17/13 12 c/ci∇cleco√y∇t2013 Peter J. Olver fares even worse, since its general solution cannot be writt en in terms of even standard special functions — although power series solutions can be t ediously ground out term by term. Understanding when a given differential equation ca n be solved in terms of elementary functions or knownspecial functions isanactiv earea ofcontemporary research, [3]. In this vein, we cannot resist mentioning that the most imp ortant class of exact solutiontechniquesfordifferentialequationsarethoseba sedonsymmetry. Anintroduction can be found in the author’s graduate level monograph [ 26]; see also [ 5,16]. Existence Before worrying about how to solve a differential equation, e ither analytically, qual- itatively, or numerically, it behooves us to try to resolve t he core mathematical issues of existence and uniqueness. First, does a solution exist? If, not, it makes no sense trying to find one. Second, is the solution uniquely determined? Other wise, the differential equation probably has scant relevance for physical applications sin ce we cannot use it as a predictive tool. Since differential equations inevitably have lots of s olutions, the only way in which we can deduce uniqueness is by imposing suitable initial (or boundary) conditions. Unlike partial differential equations, which must be treate d on a case-by-case basis, there are complete general answers to both the existence and uniqueness questions for initial value problems for systems of ordinary differential equations. (Boundary value problems are more subtle.) While obviously important, we wi ll not take the time to present the proofs of these fundamental results, which can b e found in most advanced textbooks on the subject, including [ 2,13,15,17]. Let us begin by stating the Fundamental Existence Theorem fo r initial value problems associated with first order systems of ordinary differential equations. Theorem 3.1. LetF(t,u)be a continuous function. Then the initial value problem† du dt=F(t,u), u(t0) =a, (3.3) admits a solution u=f(t)that is, at least, defined for nearby times, i.e., when |t−t0|< δ for some δ >0. Theorem 3.1 guarantees that the solution to the initial valu e problem exists — at least for times sufficiently close to the initial instant t0. This may be the most that can be said, although in many cases the maximal interval α < t < β of existence of the solution might be much larger — possibly infinite, −∞< t <∞, resulting in a global solution . The interval of existence of a solution typically depends up on both the equation and the particular initial data. For instance, even though its righ t hand side is defined everywhere, the solutions to the scalar initial value problem (2.8) only exist up until time 1 /u0, and so, the larger the initial data, the shorter the time of exist ence. In this example, the only global solution is the equilibrium solution u(t)≡0. It is worth noting that this short-term †IfF(t,u) is only defined on a subdomain Ω ⊂Rn+1, then we must assume that the point (t0,a)∈Ω specifying the initial conditions belongs to its domain of definition . 11/17/13 13 c/ci∇cleco√y∇t2013 Peter J. Olver existence phenomenon does not appear in the linear regime, w here, barring singularities in the equation itself, solutions to a linear ordinary differ ential equation are guaranteed to exist for all time. In practice, one always extends a solutions to its maximal in terval of existence. The Existence Theorem 3.1 implies that there are only two possib le ways in whcih a solution cannot be extended beyond a time t⋆: Either (i) the solution becomes unbounded: /ba∇dblu(t)/ba∇dbl → ∞ast→t⋆, or (ii) if the right hand side F(t,u) is only defined on a subset Ω ⊂Rn+1, then the solution u(t) reaches the boundary ∂Ω ast→t⋆. If neither occurs in finite time, then the solution is necessa rily global. In other words, a solution to an ordinary differential equation cannot sudden ly vanish into thin air. Remark: The existence theorem can be readily adapted to any higher o rder system of ordinary differential equations through the method of con verting it into an equivalent first order system by introducing additional variables. The appropriate initial conditions guaranteeing existence are induced from those of the corres ponding first order system, as in the second order example (2.27) discussed above. Uniqueness and Smoothness Asimportant asexistenceisthequestionof uniqueness. Doe stheinitialvalueproblem have more than one solution? If so, then we cannot use the diffe rential equation to predict the future behavior of the system from its current state. Whi le continuity of the right hand side of the differential equation will guarantee that a s olution exists, it is not quite sufficient to ensure uniqueness of the solution to the initial value problem. The difficulty can be appreciated by looking at an elementary example. Example 3.2. Consider the nonlinear initial value problem du dt=5 3u2/5, u (0) = 0. (3.4) Since the right hand side is a continuous function, Theorem 3 .1 assures us of the existence of a solution — at least for tclose to 0. This autonomous scalar equation can be easily solved by the usual method: /integraldisplay3 5du u2/5=u3/5=t+c, and so u= (t+c)5/3. Substituting into the initialconditionimpliesthat c= 0, andhence u(t) =t5/3isa solution to the initial value problem. On the other hand, since the right hand side of the differentia l equation vanishes at u= 0, the constant function u(t)≡0 is anequilibrium solutionto the differential equation. (Hereisanexamplewhere theintegrationmethodfailstorec over theequilibriumsolution.) Moreover, the equilibrium solution has the same initial val ueu(0) = 0. Therefore, we have constructed two different solutions to the initial value pro blem (3.4). Uniqueness is not 11/17/13 14 c/ci∇cleco√y∇t2013 Peter J. Olver 0.5 1 1.5 20.511.52 Figure 4. Solutions to the Differential Equation/squaresmallsolidu=5 3u2/5. valid! Worse yet, there are, in fact, an infinitenumber of solutions to the initial value problem. For anya >0, the function u(t) =/braceleftbigg0, 0≤t≤a, (t−a)5/3, t≥a,(3.5) is differentiable everywhere, even at t=a. (Why?) Moreover, it satisfies both the differ- ential equation and the initial condition, and hence defines a solution to the initial value problem. Several of these solutions are plotted in Figure 4. Thus, to ensure uniqueness of solutions, we need to impose a m ore stringent condition, beyond mere continuity. The proof of the following basic uni queness theorem can be found in the above references. Theorem 3.3. IfF(t,u)∈C1is continuously differentiable, then there exists one and only one solution†to the initial value problem (3.3). Thus, the difficulty with the differential equation (3.4) is th at the function F(u) = 5 3u2/5, although continuous everywhere, is not differentiable at u= 0, and hence the Uniqueness Theorem 3.3 does not apply. On the other hand, F(u) is continuously differ- entiable away from u= 0, and so any nonzero initial condition u(t0) =u0/ne}ationslash= 0 will produce a unique solution — for as long as it remains away from the prob lematic value u= 0. Blanket Hypothesis : From now on, all differential equations must satisfy the uni que- ness criterion that their right hand side is continuously di fferentiable. While continuous differentiability is sufficient to guarante e uniqueness of solutions, the smoother the right hand side of the system, the smoother t he solutions. Specifically: Theorem 3.4. IfF∈Cnforn≥1, then any solution to the system/squaresmallsolidu=F(t,u)is of classu∈Cn+1. IfF(t,u)is an analytic function, then all solutions u(t)are analytic. †As noted earlier, we extend all solutions to their maximal interval of e xistence. 11/17/13 15 c/ci∇cleco√y∇t2013 Peter J. Olver The basic outline of the proof of the first result is clear: Con tinuity of u(t) (which is a basic prerequiste of any solution) implies continuity o fF(t,u(t)), which means/squaresmallsoliduis continuous and hence u∈C1. This in turn implies F(t,u(t)) =/squaresmallsoliduis a continuously differentiable of t, and so u∈C2. And so on, up to order n. The proof of analyticity follows from a detailed analysis of the power series solutio ns, [14]. Indeed, the analytic result underlies the method of power series solutions of ord inary differential equations, [2,13]. Uniqueness has a number of particularly important conseque nces for the solutions to autonomous systems, i.e., those whoe right hand side does no t explicitly depend upon t. Throughout the remainder of this section, we will deal with a n autonomous system of ordinary differential equations du dt=F(u),where F∈C1, (3.6) whose right hand side is defined and continuously differentia ble for all uin a domain Ω⊂Rn. As a consequence, each solution u(t) is, on its interval of existence, uniquely determined by its initial data. Autonomy of the differential equation is an essential hy- pothesis for the validity of the following properties. The first result tells us that the solution trajectories of an autonomous system do not vary over time. Proposition 3.5. Ifu(t)is the solution to the autonomous system (3.6)with initial condition u(t0) =u0, then the solution to the initial value problem /tildewideu(t1) =u0is/tildewideu(t) = u(t−t1+t0). Proof: Let/tildewideu(t) =u(t−t1+t0), where u(t) is the original solution. In view of the chain rule and the fact that t1andt0are fixed, d dt/tildewideu(t) =du dt(t−t1+t0) =F/parenleftbig u(t−t1+t0)/parenrightbig =F/parenleftbig /tildewideu(t)/parenrightbig , and hence /tildewideu(t) is also a solution to the system (3.6). Moreover, /tildewideu(t1) =u(t0) =u0 has the indicated initial conditions, and hence, by uniquen ess, must be the one and only solution to the latter initial value problem. Q.E.D. Note that the two solutions u(t) and/tildewideu(t) parametrize the samecurve inRn, differing only by an overall “phase shift”, t1−t0, in their parametrizations. Thus, all solutions passing through the point u0follow the same trajectory, irrespective of the time they arrive there. Indeed, not only is the trajectory the same, bu t the solutions have identical speeds at each point along the trajectory curve. For instanc e, if the right hand side of (3.6) represents the velocity vector field of steady state fluid flow , Proposition 3.5 implies that the stream lines — the paths followed by the individual fluid p articles — do not change in time, even though the fluid itself is in motion. This, indee d, is the meaning of the term “steady state” in fluid mechanics. 11/17/13 16 c/ci∇cleco√y∇t2013 Peter J. Olver One particularly important consequence of uniqueness is th at a solution u(t) to an autonomous system is either stuck at an equilibrium for all t ime, or is always in motion. In other words, either/squaresmallsolidu≡0, in the case of equilibrium, or, otherwise,/squaresmallsolidu/ne}ationslash=0wherever defined. Proposition 3.6. Letu⋆be an equilibrium for the autonomous system (3.6), so F(u⋆) =0. Ifu(t)is any solution such that u(t⋆) =u⋆at some time t⋆, thenu(t)≡u⋆is the equilibrium solution. Proof: We regard u(t⋆) =u⋆as initial data for the given solution u(t) at the initial timet⋆. SinceF(u⋆) =0, the constant function u⋆(t)≡u⋆is a solution of the differential equation that satisfies the same initial conditions. Theref ore, by uniqueness, it coincides with the solution in question. Q.E.D. In other words, it is mathematically impossible for a soluti on to reach an equilibrium position in a finite amount of time — although it may well appro ach equilibrium in an asymptotic fashion as t→ ∞; see Proposition 3.9 below for details. Physically, this ob ser- vation has the interesting and physically counterintuitiv e consequence that a mathematical system never actually attains an equilibrium position! Eve n at very large times, there is always some very slight residual motion. In practice, thoug h, once the solution gets suffi- ciently close to equilibrium, we are unable to detect the mot ion, and the physical system has, in all but name, reached its stationary equilibrium con figuration. And, of course, the inherent motion of the atoms and molecules not included in su ch a simplified model would hide any infinitesimal residual effects of the mathematical s olution. Without uniqueness, the result is false. For example, the function u(t) = (t−t⋆)5/3is a solution to the scalar ordinary differential equation (3.4) that reaches the equil ibrium point u⋆= 0 in a finite timet=t⋆. Continuous Dependence Inareal-worldapplications, initialconditionsarealmos tnever knownexactly. Rather, experimental and physical errors will only allow us to say th at their values are approxi- mately equal to those in our mathematical model. Thus, to ret ain physical relevance, we need to be sure that small errors in our initial measurements do not induce a large change in the solution. A similar argument can be made for any physic al parameters, e.g., masses, charges, spring stiffnesses, frictional coefficients, etc., that appear in the differential equa- tion itself. A slight change in the parameters should not hav e a dramatic effect on the solution. Mathematically, what we are after is a criterion of continuous dependence of solutions upon both initial data and parameters. Fortunately, the des ired result holds without any additional assumptions, beyond requiring that the paramet ers appear continuously in the differential equation. We state both results in a single theo rem. Theorem 3.7. Consider an initial value problem problem du dt=F(t,u,µ), u(t0) =a(µ), (3.7) 11/17/13 17 c/ci∇cleco√y∇t2013 Peter J. Olver in which the differential equation and/or the initial condit ions depend continuously on one or more parameters µ= (µ1,...,µk). Then the unique†solution u(t,µ)depends continuously upon the parameters. Example 3.8. Let us look at a perturbed version du dt=αu2, u (0) =u0+ε, of the initial value problem that we considered in Example 2. 1. We regard εas a small perturbation of our original initial data u0, andαas a variable parameter in the equation. The solution is u(t,ε) =u0+ε 1−α(u0+ε)t. (3.8) Note that, where defined, this is a continuous function of bot h parameters α,ε. Thus, a small change in the initial data, or in the equation, produce s a small change in the solution — at least for times near the initial time. Continuous dependence does not preclude nearby solutions from eventually becoming far apart. Indeed, the blow-up time t⋆= 1//bracketleftbig α(u0+ε)/bracketrightbig for the solution (3.8) depends upon both the initial data and the parameter in the equation. Thus, as we approach the singularity, solutions that started out very close to each o ther will get arbitrarily far apart; see Figure 1 for an illustration. An even simpler example is the linear model of exponential gr owth/squaresmallsolidu=αuwhen α >0. A very tiny change in the initial conditions has a negligib le short term effect upon the solution, but over longer time intervals, the difference s between the two solutions will be dramatic. Thus, the “sensitive dependence” of solutions on initial conditions already appears in very simple linear equations. For similar reason s, sontinuous dependence does notprevent solutions from exhibiting chaotic behavior. Furth er development of these ideas can be found in [ 1,8] and elsewhere. As an application, let us show that if a solution to an autonom ous system converges to a single limit point, then that point is necessarily an equ ilibrium solution. Keep in mind that, owing to uniqueness of solutions, the limiting equili brium cannot be mathematically achieved in finite time, but only as a limit as time goes to infin ity. Proposition 3.9. Letu(t)be a solution to the autonomous system/squaresmallsolidu=F(u), with F∈C1, such that lim t→∞u(t) =u⋆. Thenu⋆is an equilibrium solution, and so F(u⋆) =0. Proof: Letv(t,a)denotethesolutiontotheinitialvalueproblem/squaresmallsolidv=F(v),v(0) =a. (We use a different letter to avoid confusion with the given so lutionu(t).) Theorem 3.7 implies that v(t,a) is a continuous function of the initial position a. and hence v(t,u(s)) is a continuous function of s∈R. Since lim s→∞u(s) =u⋆, we have lim s→∞v(t,u(s)) =v(t,u⋆). †We continue to impose our blanket uniqueness hypothesis. 11/17/13 18 c/ci∇cleco√y∇t2013 Peter J. Olver Ontheotherhand, sincethesystemisautonomous,Propositi on3.5impliesthat v(t,u(s)) = u(t+s), and hence lim s→∞v(t,u(s)) = lim s→∞u(t+s) =u⋆. Equating the preceding two limit equations, we conclude tha tv(t,u⋆) =u⋆for allt, and hence the solution with initial value v(0) =u⋆is an equilibrium solution. Q.E.D. The same conclusion holds if we run time backwards: if lim t→−∞u(t) =u⋆, thenu⋆is also an equilibrium point. When they exist, solutions that s tart and end at equilibrium pointsplayaparticularlyroleinthedynamics, andareknow nasheteroclinic ,or,ifthestart and end equilibria are the same, homoclinic orbits . Of course, limiting equilibrium points arebutoneofthepossiblelongtermbehaviorsofsolutionst ononlinearordinarydifferential equations, which can also become unbounded infinite or infini te time, or approach periodic orbits, known as limit cycles , or become completely chaotic, depending upon the nature of the system and the initial conditions. Resolving the long te rm behavior os solutions is one of the many challenges awaiting the detailed analysis of any nonlinear ordinary differential equation. 4. Stability. Once a solution to a system of ordinary differential equation s has settled down, its limiting value is an equilibrium solution; this is the conte nt of Proposition 3.9. However, not all equilibria appear in this fashion. The only steady st ate solutions that one directly observes in a physical system are the stable equilibria. Uns table equilibria are hard to sustain, and will disappear when subjected to even the tinie st perturbation, e.g., a breath of air, or outside traffic jarring the experimental apparatus . Thus, finding the equilibrium solutions to a system of ordinary differential equations is o nly half the battle; one must thenunderstand theirstabilitypropertiesinordertochar acterizethosethatcanberealized in normal physical circumstances. We will focus our attention on autonomous systems /squaresmallsolidu=F(u) whose right hand sides are at least continuously differentia ble, so as to ensure the unique- ness of solutions to the initial value problem. If everysolution that starts out near a given equilibrium solution tends to it, the equilibrium is called asymptotically stable . If the solutions that start out nearby stay nearby, then the equili brium is stable. More formally: Definition 4.1. An equilibrium solution u⋆to an autonomous system of first order ordinary differential equations is called •stableif for every (small) ε >0, there exists a δ >0 such that every solution u(t) having initial conditions within distance δ >/ba∇dblu(t0)−u⋆/ba∇dblof the equilibrium rmains within distance ε >/ba∇dblu(t)−u⋆/ba∇dblfor allt≥t0. •asymptotically stable if it is stable and, in addition, there exists δ0>0 such that whenever δ0>/ba∇dblu(t0)−u⋆/ba∇dbl, thenu(t)→u⋆ast→ ∞. 11/17/13 19 c/ci∇cleco√y∇t2013 Peter J. Olver δ εu⋆ u(t0) Stabilityδ0u⋆ u(t0) Asymptotic Stability Figure 5. Stability of Equilibria. Thus, although solutions nearby a stable equilibrium may dr ift slightly farther away, they must remain relatively close. In the case of asymptotic stability, they will eventually return to equilibrium. This is illustrated in Figure 5 Example 4.2. As we saw, the logistic differential equation du dt=λu(1−u) has two equilibrium solutions, corresponding to the two roo ts of the quadratic equation λu(1−u) = 0. The solution graphs in Figure 1 illustrate the behavior of the solutions. Observe that the first equilibrium solution u⋆ 1= 0 is unstable, since all nearby solutions go away from it at an exponentially fast rate. On the other han d, the other equilibrium solution u⋆ 2= 1 is asymptotically stable, since any solution with initia l condition 0 < u0 tends to it, again at an exponentially fast rate. Example 4.3. Consider anautonomous(meaningconstantcoefficient) homog eneous linear planar system du dt=au+bv,dv dt=cu+dv, with coefficient matrix A=/parenleftbigg a b c d/parenrightbigg . The origin u⋆=v⋆= 0 is an evident equilibrium, solution, and, moreover, is the only equilibrium provided Ais nonsingular. According to the results in [ 27; Section9.3], the stability of the origin depends upon the e igenvalues ofA: It is (globally) asymptotically stable if and only if both e igenvalues are real and negative, and is stable, but not asymptotically stable if an d only if both eigenvalues are purely imaginary, or if 0 is a double eigenvalue and so A= O. In all other cases, the origin is an unstable equilibrium. Later, we will see how this simpl e linear analysis has a direct bearing on the stability question for nonlinear planar syst ems. 11/17/13 20 c/ci∇cleco√y∇t2013 Peter J. Olver Stability of Scalar Differential Equations Before looking at any further examples, we need to develop so me basic mathematical tools for investigating the stability of equilibria. We beg in at the beginning. The stability analysis for first order scalar ordinary differential equati ons du dt=F(u) (4 .1) is particularly easy. The first observation is that all non-e quilibrium solutions u(t) are strictly monotone functions, meaning they are either always increasing or alw ays decreas- ing. Indeed, when F(u)>0, then (4.1) implies that the derivative/squaresmallsolidu>0, and hence u(t) is increasing at such a point. Vice versa, solutions are dec reasing at any point where F(u)<0. SinceF(u(t))depends continuously on t, anynon-monotonesolutionwouldhave topass throughanequilibrium valuewhere F(u⋆) = 0, inviolationofProposition3.6. This proves the claim. As a consequence of monotonicity, there are only three possi ble behaviors for a non- equilibrium solution: (a) it becomes unbounded at some finite time: |u(t)| → ∞ast→t⋆; or (b) it exists for all t≥t0, but becomes unbounded as t→ ∞; or (c) it exists for all t≥t0and has a limiting value, u(t)→u⋆ast→ ∞, which, by Proposition 3.9 must be an equilibrium point. Let us look more carefully at the last eventuality. Suppose u⋆is an equilibrium point, so F(u⋆) = 0. Suppose that F(u)>0 for all ulying slightly below u⋆, i.e., on an interval of the form u⋆−δ < u < u⋆. Any solution u(t) that starts out on this interval, u⋆−δ < u(t0)< u⋆must be increasing. Moreover, u(t)< u⋆for alltsince, according to Proposition 3.6, the solution cannot pass through the equil ibrium point. Therefore, u(t) is a solution of type ( c). It must have limiting value u⋆, since by assumption, this is the only equilibrium solution it can increase to. Therefore, in this situation, the equilibrium pointu⋆isasymptotically stable from below : solutions that start out slightly below return to it in the limit. On the other hand, if F(u)<0 for all uslightly below u⋆, then any solution that starts out in this regime will be monotonicall y decreasing, and so will move downwards, away from the equilibrium point, which is thus unstable from below . By the same reasoning, if F(u)<0 foruslightly above u⋆, then solutions starting out there will be monotonically decreasing, bounded from be low byu⋆, and hence have no choice but to tend to u⋆in the limit. Under this condition, the equilibrium point isasymptotically stable from above . The reverse inequality, F(u)>0, corresponds to solutions that increase away from u⋆, which is hence unstable from above . Combining the two stable cases produces the basic asymptotic stabilit y criterion for scalar ordinary differential equations. Theorem 4.4. A equilibrium point u⋆of an autonomous scalar differential equation is asymptotically stable if and only if F(u)>0foru⋆−δ < u < u⋆andF(u)<0for u⋆< u < u⋆+δ, for some δ >0. 11/17/13 21 c/ci∇cleco√y∇t2013 Peter J. Olver uF(u) u⋆uF(u) u⋆ u tu⋆ Stable Equilibriumu tu⋆ Unstable Equilibrium Figure 6. Equilibria of Scalar Ordinary Differential Equations. In other words, if F(u) switches sign from positive to negative as uincreases through the equilibrium point, then the equilibrium is asymptotica lly stable. If the inequalities are reversed, and F(u) goes from negative to positive, then the equilibrium point is unstable. The two cases are illustrated in Figure 6. An equilibrium poi nt where F(u) is of one sign on both sides, e.g., the point u⋆= 0 forF(u) =u2, is stable from one side, and unstable from the other. Example 4.5. Consider the differential equation du dt=u−u3. (4.2) Solving the algebraic equation F(u) =u−u3= 0, we find that the equation has three equilibria: u⋆ 1=−1,u⋆ 2= 0,u⋆ 3= +1, As uincreases, the graph of the function F(u) = u−u3switches from positive to negative at the first equilibrium p ointu⋆ 1=−1, which proves its stability. Similarly, the graph goes back from po sitive to negative at u⋆ 2= 0, 11/17/13 22 c/ci∇cleco√y∇t2013 Peter J. Olver -1 -0.5 0.5 1 -0.75-0.5-0.250.250.50.75 0.5 1 1.5 2 -2-1.5-1-0.50.511.52 Figure 7. Stability of/squaresmallsolidu=u−u3. establishing the instability of the second equilibrium. Th e final equilibrium u⋆ 3= +1 is stable because F(u) again changes from negative to positive there. With this information in hand, we are able to completely char acterize the behavior of all solutions to the system. Any solution with negative init ial condition, u0<0, will end up, asymptotically, at the first equilibrium, u(t)→ −1 ast→ ∞. Indeed, if u0<−1, then u(t) is monotonically increasing to −1, while if −1< u0<0, the solution is decreasing towards −1. On the other hand, if u0>0, the corresponding solution ends up at the other stable equilibrium, u(t)→+1; those with 0 < u0<1 are monotonically increasing, while those with u0>1 are decreasing. The only solution that does not end up at eit her −1 or +1 as t→ ∞is the unstable equilibrium solution u(t)≡0. Any perturbation of it, no matter how tiny, will force the solutions to choose o ne of the stable equilibria. Representative solutions are plotted in Figure 7. Note that all the curves, with the sole exception of the horizontal axis, converge to one of the stab le solutions ±1, and diverge from the unstable solution 0 as t→ ∞. Thus, the sign of the function F(u) nearby an equilibrium determines its stability. In most instances, this can be checked by looking at the deriv ative of the function at the equilibrium. If F′(u⋆)<0, then we are in the stable situation, where F(u) goes from positive to negative with increasing u, whereas if F′(u⋆)>0, then the equilibrium u⋆ unstable on both sides. Theorem 4.6. Letu⋆be a equilibrium point for a scalar ordinary differential equ a- tion/squaresmallsolidu=F(u). IfF′(u⋆)<0, thenu⋆is asymptotically stable. If F′(u⋆)>0, thenu⋆is unstable. For instance, in the preceding example, F′(u) = 1−3u2, and its value at the equilibria are F′(−1) =−2<0, F′(0) = 1>0, F′(1) =−2<0. 11/17/13 23 c/ci∇cleco√y∇t2013 Peter J. Olver The signs reconfirm our conclusion that ±1 are stable equilibria, while 0 is unstable. In the borderline case when F′(u⋆) = 0, the derivative test is inconclusive, and further analysisisneeded toresolvethestatusoftheequilibriump oint. Forexample, theequations /squaresmallsolidu=u3and/squaresmallsolidu=−u3bothsatisfy F′(0) = 0attheequilibriumpoint u⋆= 0. But, according to the criterion of Theorem 4.4, the former has an unstable eq uilibrium, while the latter’s is stable. Thus, Theorem 4.6 is not as powerful as the direct a lgebraic test in Theorem 4.4. But it does have the advantage of being a bit easier to use. Mor e significantly, unlike the algebraic test, it can be directly generalized to systems of ordinary differential equations. Linearization and Stability In higher dimensional situations, we can no longer rely on si mple monotonicity prop- erties, and a more sophisticated approach to stability issu es is required. The key idea is already contained in the second characterization of stable equilibria in Theorem 4.6. The derivative F′(u⋆) determines the slope of the tangent line, which is a linear a pproximation to the function F(u) near the equilibrium point. In a similar fashion, a vector- vallued function F(u) is replaced by its linear approximation near an equilibriu m point. The basic stability criteria for the resulting linearized differenti al equation were established in [ 27; Section9.2]. and, in most situations, the linearized stabi lity or instability carries over to the nonlinear regime. Let us first revisit the scalar case du dt=F(u) (4 .3) from this point of view. Linearization of a scalar function at a point means to replace it by its tangent line approximation F(u)≈F(u⋆)+F′(u⋆)(u−u⋆) (4 .4) Ifu⋆is an equilibrium point, then F(u⋆) = 0, and so the first term disappears. Therefore, we anticipate that, near the equilibrium point, the solutio ns to the nonlinear ordinary differential equation (4.3) will be well approximated by its linearization du dt=F′(u⋆)(u−u⋆). Let us rewrite the linearized equation in terms of the deviat ionv(t) =u(t)−u⋆of the solution from equilibrium. Since u⋆is fixed,dv/dt=du/dt, and so the linearized equation takes the elementary form dv dt=av, where a=F′(u⋆) (4 .5) is the value of the derivative at the equilibrium point. Note that the original equilibrium pointu⋆corresponds to the zero equilibrium point v⋆= 0 of the linearized equation (4.5). We already know that the linear differential equation (4.5) h as an asymptotically stable equilibrium at v⋆= 0 if and only if a=F′(u⋆)<0, while for a=F′(u⋆)>0 the origin is unstable. In this manner, the linearized stability criteri on reproduces that established in Theorem 4.6. 11/17/13 24 c/ci∇cleco√y∇t2013 Peter J. Olver The same linearization technique can be applied to analyze t he stability of an equi- librium solution u⋆to a first order autonomous system /squaresmallsolidu=F(u). (4.6) We approximate the function F(u) near an equilibrium point, where F(u⋆) =0, by its first order Taylor polynomial: F(u)≈F(u⋆)+F′(u⋆)(u−u⋆) =F′(u⋆)(u−u⋆). (4.7) Here,F′(u⋆) denotes its n×nJacobian matrix at the equilibrium point. Thus, for nearby solutions, we expect that the deviationfrom equilibrium, v(t) =u(t)−u⋆, will be governed by the linearized system dv dt=Av,where A=F′(u⋆). (4.8) Now, we already know the complete stability criteria for lin ear systems, [ 27; Section 9.2]. The zero equilibrium solution to (4.8) is asymptotica lly stable if and only if all the eigenvalues of the coefficient matrix A=F′(u⋆) have negative real part. In contrast, if one or more of the eigenvalues has positive real part, then th e zero solution is unstable. Indeed, it can be proved, [ 13,15], that these linearized stability criteria are also valid i n the nonlinear case. Theorem 4.7. Letu⋆be an equilibrium point for the first order ordinary different ial equation/squaresmallsolidu=F(u). If all of the eigenvalues of the Jacobian matrix F′(u⋆)have negative real part, Reλ <0, thenu⋆is asymptotically stable. If, on the other hand, F′(u⋆)has one or more eigenvalues with positive real part, Reλ >0, thenu⋆is an unstable equilibrium. Intuitively, the additional nonlinear terms in the full sys tem should only slightly per- turb the eigenvalues, and hence, at least for those with nonz ero real part, not alter their effect on the stability of solutions. The borderline case occ urs when one or more of the eigenvalues of F′(u⋆) is either 0 or purely imaginary, i.e., Re λ= 0, while all other eigenval- ues have negative real part. In such situations, the lineari zed stability test is inconclusive, and we need more detailed information (which may not be easy t o come by) to resolve the status of the equilibrium. Example 4.8. The second order ordinary differential equation md2θ dt2+µdθ dt+κsinθ= 0 (4 .9) describes the damped oscillations of a rigid pendulum that r otates on a pivot subject to a uniform gravitational force in the vertical direction. The unknown function θ(t) measures the angle of the pendulum from the vertical, as illustratedi n Figure 8. The constant m >0 is the mass of the pendulum bob, µ >0 is the coefficient of friction, assumed here to be strictly positive, and κ >0 represents the gravitational force. 11/17/13 25 c/ci∇cleco√y∇t2013 Peter J. Olver θ Figure 8. The Pendulum. In order to study the equilibrium solutions and their stabil ity, we must first convert the equation into a first order system. Setting u(t) =θ(t), v(t) =dθ dt, we find du dt=v,dv dt=−αsinu−βv,where α=κ m, β=µ m,(4.10) are both positive constants. The equilibria occur where the right hand sides of the first order system (4.10) simultaneously vanish, that is, v= 0,−αsinu−βv= 0,and hence u= 0,±π,±2π, ... . Thus, the system has infinitely many equilibrium points: u⋆ k= (kπ,0) where k= 0,±1,±2,...is any integer. (4 .11) The equilibrium point u⋆ 0= (0,0) corresponds to u=θ= 0,v=/squaresmallsolid θ= 0, which means that the pendulum is at rest at the bottom of its arc. Our physi cal intuition leads us to expect this to describe a stable configuration, as the fricti onal effects will eventually damp out small nearby motions. The next equilibrium u⋆ 1= (π,0) corresponds to u=θ=π, v=/squaresmallsolid θ= 0, which means that the pendulum is sitting motionless at th e top of its arc. This is a theoretically possible equilibrium configuration, but highly unlikely to be observed in practice, and is thus expected to be unstable. Now, since u=θis an angular variable, equilibria whose uvalues differ by an integer multiple of 2 πdefine the same physical configuration, and hence should have identical stability pr operties. Therefore, all the remaining equilibria u⋆ kphysically correspond to one or the other of these two possib ilities: whenk= 2jis even, the pendulum is at the bottom, while when k= 2j+1 is odd, the pendulum is at the top. Let us now confirm our intuition by applying the linearizatio n stability criterion of Theorem 4.7. The right hand side of the system (4.10), namely F(u,v) =/parenleftbigg v −αsinu−βv/parenrightbigg ,has Jacobian matrix F′(u,v) =/parenleftbigg 0 1 −αcosu−β/parenrightbigg . 11/17/13 26 c/ci∇cleco√y∇t2013 Peter J. Olver Figure 9. The Underdamped Pendulum. At the bottom equilibrium u⋆ 0= (0,0), the Jacobian matrix F′(0,0) =/parenleftbigg 0 1 −α−β/parenrightbigg has eigenvalues λ=−β±/radicalbig β2−4α 2. Under our assumption that α,β >0, both eigenvalues have negative real part, and hence the origin is a stable equilibrium. If β2<4α— theunderdamped case — the eigenvalues are complex, and hence, in the terminology of Section 2dode, the origin is a stable focus . In the phase plane, the solutions spiral in to the focus, whic h corresponds to a pendulum with damped oscillations of decreasing magnitude. On the ot her hand, if β2>4α, then the system is overdamped . Both eigenvalues are negative, and the origin is a stable node . In this case, the solutions decay exponentially fast to 0. Physically, this would be like a pendulum moving in a vat of molasses. The exact same analysi s applies at all even equilibria u⋆ 2j= (2jπ,0) — which really represent the same bottom equilibrium poin t. On the other hand, at the top equilibrium u⋆ 1= (π,0), the Jacobian matrix F′(0,0) =/parenleftbigg 0 1 α−β/parenrightbigg has eigenvalues λ=−β±/radicalbig β2+4α 2. In this case, one of the eigenvalues is real and positive whil e the other is negative. The linearized system has an unstable saddle point, and hence th e nonlinear system is also unstable at this equilibrium point. Any tiny perturbation o f an upright pendulum will dislodge it, causing it to swing down, and eventually settle into a damped oscillatory motion converging on one of the stable bottom equilibria. The complete phase portrait of an underdamped pendulum appe ars in Figure 9. Note that,asadvertised, almostallsolutionsendupspiralingi ntothestableequilibria. Solutions with a large initial velocity will spin several times around the center, but eventually the cumulative effect of frictional forces wins out and the pendu lum ends up in a damped oscillatory mode. Each of the the unstable equilibria has th e same saddle form as its linearizations, with two very special solutions, correspo nding to the stable eigenline of the linearization, in which the pendulum spins around a few time s, and, in the t→ ∞limit, ends up standing upright at the unstable equilibrium positi on. However, like unstable 11/17/13 27 c/ci∇cleco√y∇t2013 Peter J. Olver Figure 10. Phase Portrait of the van der Pol System. equilibria, such solutions are practically impossible to a chieve in a physical environment as any tiny perturbation will cause the pendulum to sightly d eviate and then end up eventually decaying into the usual damped oscillatory moti on at the bottom. A deeper analysis demonstrates the local structural stability of any nonlinear equi- librium whose linearization is structurally stable, and he nce has no eigenvalues on the imaginary axis: Re λ/ne}ationslash= 0. Structural stability means that, not only are the stabil ity properties of the equilibrium dictated by the linearized ap proximation, but, nearby the equilibrium point, all solutions to the nonlinear system ar e slight perturbations of solu- tionsto the corresponding linearized system, and so, close to theequilibrium point, the two phase portraits have the same qualitative features. Thus, s table foci of the linearization remain stable foci of the nonlinear system; unstable saddle points remain saddle points, although the eigenlines become slightly curved as they depa rt from the equilibrium. Thus, the structural stability of linear systems, as discussed at the end of [ 27; Section9.3] also carries over to the nonlinear regime near an equilibrium. A m ore in depth discussion of these issues can be found, for instance, in [ 13,15]. , Example 4.9. Consider the unforced van der Pol system du dt=v,dv dt=−(u2−1)v−u. (4.12) that we derived in Example 2.5. The only equilibrium point is at the origin u=v= 0. Computing the Jacobian matrix of the right hand side, F′(u,v) =/parenleftbigg 0 1 2uv−1 1/parenrightbigg ,and hence F′(0,0) =/parenleftbigg 0 1 −1 1/parenrightbigg . 11/17/13 28 c/ci∇cleco√y∇t2013 Peter J. Olver Figure 11. Phase Portrait for/squaresmallsolidu=u(v−1),/squaresmallsolidv= 4−u2−v2.. The eigenvalues of F′(0,0) are1 2±i√ 3 2, and correspond to an unstable focus of the lin- earized system near the equilibrium point. Therefore, the o rigin is an unstable equilibrium for nonlinear van der Pol system, and all non-equilibrium so lutions starting out near 0 eventually spiral away. On the other hand, it can be shown that solutions that are suffic iently far away from the origin spiral in towards the center. So what happens to th e solutions? As illustrated in the phase plane portrait sketched in Figure 10, all non-equi librium solutions spiral towards a stable periodic orbit, known as a limit cycle for the system. Any non-zero initialdata will eventually end up closely following the limit cycle orbit as it periodically circles around the origin. A rigorous proof of the existence of a limit cycle rel ies on the more sophisticated Poincar´ e–Bendixson Theory for planar autonomous systems, discussed in detail in [ 13]. Example 4.10. The nonlinear system du dt=u(v−1),dv dt= 4−u2−v2, has four equilibria: (0 ,±2) and (±√ 3,1). Its Jacobian matrix is F′(u,v) =/parenleftbigg v−1u −2u−2v/parenrightbigg . A table of the eigenvalues at the equilibrium points and thei r stability follows: These results are reconfirmed by the phase portrait drawn in Figure 11 11/17/13 29 c/ci∇cleco√y∇t2013 Peter J. Olver Equilibrium Point Jacobian matrix Eigenvalues Stability (0,2)/parenleftbigg 1 0 0−4/parenrightbigg 1,−4unstable saddle (0,−2)/parenleftbigg −3 0 0 6/parenrightbigg −3,6unstable saddle (√ 3,1)/parenleftbigg 0−√ 3 2√ 3−2/parenrightbigg −1±i√ 5stable focus (−√ 3,1)/parenleftbigg 0−√ 3 2√ 3−2/parenrightbigg −1±i√ 5stable focus Conservative Systems When modeling a physical system that includes some form of da mping — due to fric- tion, viscosity, or dissipation — linearization will usual ly suffice to resolve the stability or instability of equilibria. However, when dealing with cons ervative systems, when damp- ing is absent and energy is preserved, the linearization tes t is often inconclusive, and one must rely on more sophisticated stability criteria. In such situations, one can often exploit conservation of energy, appealing to our general philosoph y that minimizers of an energy function should be stable (but not necessarily asymptotica lly stable) equilibria. By saying that energy is conserved , we mean that it remains constant as the solution evolves. Conserved quantities are also known as first integrals for the system of ordinary differential equations. Additional well-known examples in clude the laws of conservation of mass, and conservation of linear and angular momentum. Let u s mathematically formulate the general definition. Definition 4.11. Afirst integral of an autonomous system/squaresmallsolidu=F(u) is a real-valued function I(u) which is constant on solutions. In other words, for each solution u(t) to the differential equation, I(u(t)) =cfor all t, (4.13) wherecis a fixed constant, which will depend upon which solution is b eing monitored. The value of cis fixed by the initial data since, in particular, c=I(u(t0)) =I(u0). Or, to rephrase this condition in another, equivalent manner, e very solution to the dynamical system is constrained to move along a single level set {I(u) =c}of the first integral, namely the level set that contains the initial data u0. Note first that any constant function, I(u)≡c0, is trivially a first integral, but this provides no useful information whatsoever about the soluti ons, and so is uninteresting. We 11/17/13 30 c/ci∇cleco√y∇t2013 Peter J. Olver willcallanyautonomoussystemthatpossessesanontrivial firstintegral I(u)aconservative system. How do we find first integrals? In applications, one often appe als to the underlying physical principles such as conservation of energy, moment um, or mass. Mathematically, the most convenient way to check whether a function is consta nt is to verify that its derivative is identically zero. Thus, differentiating (4.1 3) with respect to tand invoking the chain rule leads to the basic condition 0 =d dtI(u(t)) =∇I(u(t))·du dt=∇I(u(t))·F(u(t)). (4.14) The final expression can be identified as the directional deri vative of I(u) with respect to the vector field v=F(u) that specifies the differential equation. Writing out (4.14 ) in detail, we find that a first integral I(u1,...,un) must satisfy a first order linear partial differential equation: F1(u1,...,un)∂I ∂u1+···+Fn(u1,...,un)∂I ∂un= 0. (4.15) As such, it looks harder to solve than the original ordinary d ifferential equation! Often, one falls back on either physical intuition, intelligent gu esswork, or, as a last resort, a lucky guess. A deeper fact, due to the pioneering twentieth c entury mathematician Emmy Noether,cf.[ 24,26], isthatfirstintegralsandconservationlawsaretheresul tofunderlying symmetry properties of the differential equation. Like many nonlinear methods, it remains the subject of contemporary research. Let us specialize to planar autonomous systems du dt=F(u,v),dv dt=G(u,v). (4.16) According to (4.15), any first integral I(u,v) must satisfy the linear partial differential equation F(u,v)∂I ∂u+G(u,v)∂I ∂v= 0. (4.17) This nonlinear first order partial differential equation can be solved by the method of characteristics. Consider the auxiliary first order scalar ordinary differential equation† dv du=G(u,v) F(u,v)(4.18) forv=h(u). Note that (4.18) can be formally obtained by dividing the s econd equation in the original system (4.16) by the first, and then canceling th e timedifferentials dt. Suppose we can write the general solution to the scalar equation (4.1 8) in the implicit form I(u,v) =c, (4.19) †We assume that F(u,v)/negationslash≡0. Otherwise, I(u) =uis itself a first integral, and the system reduces to a scalar equation for v. 11/17/13 31 c/ci∇cleco√y∇t2013 Peter J. Olver wherecis a constant of integration. We claim that the function I(u,v) is a first integral of the original system (4.16). Indeed, differentiating (4.1 9) with respect to u, and using the chain rule, we find 0 =d duI(u,v) =∂I ∂u+dv du∂I ∂v=∂I ∂u+G(u,v) F(u,v)∂I ∂v. Clearing the denominator, we conclude that I(u,v) solves the partial differential equation (4.17), which justifies our claim. Example 4.12. As an elementary example, consider the linear system du dt=−v,dv dt=u. (4.20) To construct a first integral, we form the auxiliary equation (4.18), which is dv du=−u v. This first order ordinary differential equation can be solved by separating variables: vdv=−udu, and hence1 2u2+1 2v2=c, wherecis the constant of integration. Therefore, by the preceding result, I(u,v) =1 2u2+1 2v2 is a first integral. The level sets of I(u,v) are the concentric circles centered at the origin, and we recover the fact that the solutions of (4.20) go around the circles. The origin is a stable equilibrium — a center. This simple example hints at the importance of first integral s in stability theory. The following key result confirms our general philosophy that en ergy minimizers, or, more generally, minimizers of first integrals, are stable equili bria. Theorem 4.13. LetI(u)be a first integral for the autonomous system/squaresmallsolidu=F(u). Ifu⋆is a strict local extremum — minimum or mximum — of I, thenu⋆is a stable equilibrium point for the system. Remark: At first sight, the fact that strict maxima are also stable eq uilibria appears to contradict our intuition. However, energy functions typ ically do not have local maxima. Indeed, physical energy is the sum of kinetic and potential c ontributions. While potential energy can admit maxima, e.g., the pendulum at the top of its a rc, these are only unstable saddle points for the full energy function, since the kineti c energy can always be increased by moving a bit faster. Proof: We first prove that u⋆is an equilibrium. Indeed, the solution u(t) with initial condition u(t0) =u⋆must maintain the value of I(u(t)) =I(u⋆). But, by definition of a strict local minimum, I(u)> I(u⋆) for allunearu⋆, and hence, by continuity, the solution has no choice but to remain at the point u⋆. 11/17/13 32 c/ci∇cleco√y∇t2013 Peter J. Olver To prove stability, set M(r) = max{I(u)| /ba∇dblu−u⋆/ba∇dbl ≤r}, m (r) = min{I(u)| /ba∇dblu−u⋆/ba∇dbl=r}. Thus,M(r) is the maximum value of the first integral over a ball†of radius rcentered at the minimum, while m(r) is the minimum over its boundary sphere of radius r. SinceIis continuous, so are mandM. Sinceu⋆is a strict local minimum, M(r)≥m(r)> m(0) = M(0) =I(u⋆) for any 0 < r < εsufficiently small. For each ε >0, we can choose a δ >0 such that M(δ)< m(ε). Then, whenever /ba∇dblu(t0)−u⋆/ba∇dbl ≤δ, thenI(u(t)) =I(u(t0))≤M(δ)< m(ε). Since m(ε) is the minimum possible value for I(u) when/ba∇dblu(t)−u⋆/ba∇dbl=ε, the solution u(t) cannot cross the sphere of radius εat anyt, and so /ba∇dblu(t)−u⋆/ba∇dbl< εfor allt≥t0. Hence, we have fulfilled the stability criteria of Definition 4.1. Q.E.D. Example 4.14. Consider the specific predator-prey system du dt= 2u−uv,dv dt=−9v+3uv, (4.21) modeling populations of, say, lions and zebra, and a special case of (2.26). According to Example 2.4, there are two possible equilibria: u⋆ 1=v⋆ 1= 0, u⋆ 2= 3, v⋆ 2= 2. Let us first try to determine their stability by linearizatio n. The Jacobian matrix for the system is F′(u,v) =/parenleftbigg 2−v−u 3v3u−9/parenrightbigg . At the first, trivial equilibrium, F′(0,0) =/parenleftbigg 2 0 0−9/parenrightbigg ,with eigenvalues 2 and −9. Since there is one positive and one negative eigenvalue, the origin is an unstable saddle point. On the other hand, at the nonzero equilibrium, the Jac obian matrix F′(3,2) =/parenleftbigg 0−3 6 0/parenrightbigg ,has purely imaginary eigenvalues ±3√ 2 i. So the linearized system has a stable center. However, as pur ely imaginary eigenvalues is a borderline situation, Theorem 4.7 cannot be applied. Thus , the linearization stability test isinconclusive . It turns out that the predator-prey model is a conservative s ystem. To find a first integral, we need to solve the auxiliary equation (4.18), wh ich is dv du=−9v+3uv 2u−uv=−9/u+3 2/v−1. †We write as if the norm is the Euclidean norm, but any other norm will wor k equally well for this proof. 11/17/13 33 c/ci∇cleco√y∇t2013 Peter J. Olver 2 4 6 8123456 Figure 12. Phase Portrait and Solution of the Predator-Prey System. Fortunately, this is a separable first order ordinary differe ntial equation. Integrating, 2 logv−v=/integraldisplay/parenleftbigg2 v−1/parenrightbigg dv=/integraldisplay/parenleftbigg −9 u+3/parenrightbigg du=−9logu+3u+c, wherecis the constant of integration. Writing the solution in the f orm (4.19), we conclude that I(u,v) = 9 log u−3u+2 logv−v=c, is a first integral of the system. The solutions to (4.21) must stay on the level sets of I(u,v). Note that ∇I(u,v) =/parenleftbigg 9/u−3 2/v−1/parenrightbigg ,and hence ∇I(3,2) =0, which shows that the second equilibrium is a critical point. (On the other hand, I(u,v) is not defined at the unstable zero equilibrium.) Moreover, the Hessian matrix at the critical point, ∇2I(3,2) =/parenleftbigg −3 0 0−1/parenrightbigg , is negative definite, and hence u⋆ 2= (3,2)Tis a strict local maximum of the first integral I(u,v). Thus, Theorem 4.13 proves that the equilibrium point is a s table center. The first integral serves to completely characterize the qua litative behavior of the sys- tem. In the physically relevant region, i.e., the upper righ t quadrant Q={u >0, v >0} where both populations are positive, all of the level sets of the first integral are closed curves encircling the equilibrium point u⋆ 2= (3,2)T. The solutions move in a counter- clockwise direction around the closed level curves, and hen ce all non-equilibrium solutions in the positive quadrant are periodic. The phase portrait is illustrated in Figure 12, along with a typicl periodic solution. Thus, in such an idealized e cological model, for any initial conditions starting with some zebra and lions, i.e., where u(t0),v(t0)>0, the populations will maintain a balance over the long term, each varying peri odically between maximum and minimum values. Observe also that the maximum and minimu m values of the two 11/17/13 34 c/ci∇cleco√y∇t2013 Peter J. Olver populations are not achieved simultaneously. Starting wit h a small number of predators, the number of prey will initially increase. The predators th en have more food availbable, and so also start to increase in numbers. At a certain critica l point, the predators are sufficiently numerous as to kill prey faster than they can repr oduce. At this point, the prey population has reached its maximum, and begins to decli ne. But it takes a while for the predator population to feel the effect, and so it conti nues to increase. Eventually the increasingly rapid decline in the number of prey begins t o affect the predators. After the predators reach their maximum number, both populations are in decline. Eventually, enough predators have died off so as to relieve the pressure on the prey, whose population bottoms out, and then slowly begins to rebound. Later, the nu mber of predators also reaches a minimum, at which point the entire growth and decay cycle starts over again. In contrast to a linear system, the period of the population c ycle is not fixed, but depends upon how far away from the stable equilibrium the sol ution orbit lies. Near equilibrium, the solutions are close to those of the lineari zed system which, in view of its eigenvalues ±3i√ 2, are periodic of frequency 3√ 2 and period√ 2π/3. However, solutions that are far away from equilibrium have much longer periods, and so greater imbalances between predator and prey populations leads to longer perio ds, and more radically varying numbers. Understanding the mechanisms behind these popula tion cycles is of increasing important in the ecological management of natural resource s. Example 4.15. In our next example, we look at the undamped oscillations of a pendulum. When we set the friction coefficient µ= 0, the nonlinear second order ordinary differential equation (4.9) reduces to md2θ dt2+κsinθ= 0. (4.22) As before, we convert this into a first order system du dt=v,dv dt=−αsinu, (4.23) where u(t) =θ(t), v (t) =dθ dt, α =κ m. The equilibria, u⋆ k= (kπ,0) for k= 0,±1,±2,... , are the same as in the damped case, i.e., the pendulum is eithe r at the top ( keven) or the bottom ( kodd) of the circle. Let us see what the linearization stability test tells us. In this case, the Jacobian matrix of (4.23) is F′(u,v) =/parenleftbigg 0 1 −αcosu0/parenrightbigg . At the top equilibria F′(u⋆ 2j+1) =F′/parenleftbig (2j+1)π,0/parenrightbig =/parenleftbigg 0 1 α0/parenrightbigg has real eigenvalues ±√α, 11/17/13 35 c/ci∇cleco√y∇t2013 Peter J. Olver Figure 13. The Undamped Pendulum. and hence these equilibria are unstable saddle points, just as in the damped version. On the other hand, at the bottom equilibria F′(u⋆ 2j) =F′(2jπ,0) =/parenleftbigg 0 1 −α0/parenrightbigg ,has purely imaginary eigenvalues ±i√α. Without the benefit of damping, the linearization test is inc onclusive, and the stability of the bottom equilibria remains in doubt. Since we are dealing with a conservative system, the total en ergy of the pendulum, namely E(u,v) =1 2mv2+κ(1−cosu) =m 2/parenleftbiggdθ dt/parenrightbigg2 +κ(1−cosθ) (4 .24) should provide us with a first integral. Note that Eis a sum of two terms, which represent, respectively, the kinetic energy due to the pendulum’s moti on, and the potential energy† due to the height of the pendulum bob. To verify that E(u,v) is indeed a first integral, we compute dE dt=∂E ∂udu dt+∂E ∂vdv dt= (κsinu)v+(mv)(−αsinu) = 0,since α=κ m. Therefore, Eis indeed constant on solutions, reconfirming the physical b asis of the model. The phase plane solutions to the pendulum equation will move along the level sets of the energy function E(u,v), which are plotted in Figure 13. Its critical points are the equilibria, where ∇E(u) =/parenleftbigg κsinu mv/parenrightbigg =0,and hence u=u⋆ k= (kπ,0), k= 0,±1,±2,... . To characterize the critical points, we appeal to the second derivative test, and so evaluate the Hessian ∇2E(u,v) =/parenleftbigg κcosu0 0m/parenrightbigg . †In a physical system, the potential energy is only defined up to an addi tive constant. Here we have fixed the zero energy level to be at the bottom of the pendulum ’s arc. 11/17/13 36 c/ci∇cleco√y∇t2013 Peter J. Olver At the bottom equilibria, ∇2E(u⋆ 2j) =∇2E(2jπ,0) =/parenleftbigg κ0 0m/parenrightbigg is positive definite, since κandmare positive constants. Therefore, the bottom equilibria are strict local minima of the energy, and so Theorem 4.13 gua rantees their stability. Each stable equilibrium is surrounded by a family of closed o val-shaped level curves, and hence forms a center. Each oval corresponds to a periodic solution‡of the system, in which the pendulum oscillates back and forth symmetrically around the bottom of its arc. Near the equilibrium, the period is close to that of the linea rized system, namely 2 π/√α as prescribed by the eigenvalues of the Jacobian matrix. Thi s fact underlies the use of pendulum-based clocks for keeping time, first recognized by Galileo. A grandfather clock is accurate because the amplitude of its pendulum’s oscilla tions is kept relatively small. However, moving further away from the equilibrium point in t he phase plane, we find that the periodic solutions with very large amplitude oscil lations, in which the pendulum becomes nearly vertical, have much longer periods, and so wo uld lead to inaccurate time- keeping. The large amplitude limit of the periodic solutions is of par ticular interest. The pair of trajectories connecting two distinct unstable equilibr ia are known as the homoclinic orbits. Physically, a homoclinic orbit corresponds to a pendulum t hat starts out just shy of vertical, goes through exactly one full rotation, and eve ntually (as t→ ∞) ends up vertical again. The homoclinic orbits play an essential rol e in the analysis of the chaotic behavior of a periodically forced pendulum, [ 1,8,23]. Finally, the level sets lying above and below the “cat’s-eye s” formed by the homoclinic and periodic orbits are known as the running orbits . Sinceu=θis a 2πperiodic angular variable, a running orbit solution ( u(t),v(t))T= (θ(t),/squaresmallsolid θ(t))T, in fact, also corresponds to a periodic physical motion, in which the pendulum spins ar ound and around its pivot. The larger the total energy E(u,v), the farther away from the u–axis the running orbit lies, and the faster the pendulum spins. In summary, the qualitativebehavior of a solution to the pen dulum equation is almost entirely characterized by its energy: •E= 0, stable equilibria, •0< E <2κ, periodic oscillating orbits, •E= 2κ, unstable equilibria and homoclinic orbits, •E >2κ, running orbits. Example 4.16. The system governing the dynamical rotations of a rigid soli d body around its center of mass are known as the Euler equations of rigid body mechanics, in honor of the prolific eighteenth century Swiss mathematicia n Leonhard Euler, cf. [ 11]. ‡More precisely, a family of periodic solutions indexed by their in itial condition on the oval, and differing only by a phase shift: u(t−δ). 11/17/13 37 c/ci∇cleco√y∇t2013 Peter J. Olver Figure 14. The Rigid Body Phase Portrait. The eigenvectors of the positive definite inertia tensor of t he body prescribe its three mutually orthogonal principal axes . The corresponding eigenvalues I1,I2,I3>0 are called theprincipal moments of inertia . Letu1(t),u2(t),u3(t) denote the angular momenta of the body around its three principal axes. In the absence of ex ternal forces, the dynamical system governing the body’s rotations around its center of m ass takes the symmetric form du1 dt=I2−I3 I2I3u2u3,du2 dt=I3−I1 I1I3u1u3,du3 dt=I1−I2 I1I2u1u2.(4.25) Th Euler equations model, for example, the dynamics of a sate llite spinning in its orbit above the earth. The solution will prescribe the angular mot ions of the satellite around its center of mass, but not the overall motion of the center of mass as the satellite orbits the earth. Letusassumethatthemomentsofinertiaarealldifferent, wh ichweplaceinincreasing order 0< I1< I2< I3. The equilibria of the Euler system (4.25) are where the righ t hand sides simultaneously vanish, which requires that either u2=u3= 0, oru1=u3= 0, or u1=u2= 0. In other words, every point on the three coordinate axes i s an equilibrium configuration. Since the variables represent angular momen ta, these equilibria correspond to the body spinning around one of its principal axes at a fixed angular velocity. Let us analyze the stability of these equilibrium configurat ions. The linearization test fails completely — as it must whenever dealing with a non-iso lated equilibrium. But the Euler equations turn out to admit two independent first integ rals: E(u) =1 2/parenleftbiggu2 1 I1+u2 2 I2+u2 3 I3/parenrightbigg , A (u) =1 2/parenleftbig u2 1+u2 2+u2 3/parenrightbig .(4.26) 11/17/13 38 c/ci∇cleco√y∇t2013 Peter J. Olver The first is the total kinetic energy, while the second is the t otal angular momentum. The proof that dE/dt= 0 =dA/dtfor any solution u(t) to (4.25) is left as an exercise for the reader. Since both EandAare constant, the solutions to the system are constrained to move along a common level set C={E=e, A=a}. Thus, the solution trajectories are the curves obtained by intersecting the sphere Sa={A(u) =a}of radius√ 2awith the ellipsoid Le={E(u) =e}. In Figure 14, we have graphed the solution trajectories on a fixed sphere. (To see the figure, make sure the left hand period ic orbits are perceived on the back side of the sphere.) The six equilibria on the sphere are at its intersections with the coordinate axes. Those on the xandzaxes are surrounded by closed periodic orbits, and hence are stable equilibria; indeed, they are, respecti vely, local minima and maxima of the energy when restricted to the sphere. On the other hand, t he two equilibria on the y axis have the form of unstable saddle points, and are connect ed by four distinct homoclinic orbits. We conclude that a body that spins around either of it s principal axes with the smallest or the largest moment of inertia is stable, whereas one that spins around the axis corresponding to the intermediate moment of inertia is unstable. This mathematical deduction can be demonstrated physically by flipping a solid rectangular object, e.g., this book, up into the air. It is easy to arrange it to spin around it s long axis or its short axis in a stable manner, but it will balk at attempts to make it rota te around its middle axis! Lyapunov’s Method Systems that incorporate damping, viscosity and/or fricti onal effects do not typically possessnon-constantfirstintegrals. Fromaphysicalstand point, thedampingwillcausethe total energy of the system to be a decreasing function of time . Asymptotically, the system returns to a (stable) equilibrium, and the extra energy has b een dissipated away. However, this physical principle captures important mathematical i mplications for the behavior of solutions. It leads to a useful alternative method for estab lishing stability of equilibria, even in cases when the linearization stability test is incon clusive. The nineteenth century Russian mathematician Alexander Lyapunov was the first to pi npoint the importance of such functions in dynamics. Definition 4.17. ALyapunov function for the first order autonomous system/squaresmallsolidu= F(u) is a continuous real-valued function L(u) that is non-increasing on all solutions u(t), meaning that L(u(t))≤L(u(t0)) for all t > t0. (4.27) Astrict Lyapunov function satisfies the strict inequality L(u(t))< L(u(t0)) for all t > t0, (4.28) whenever u(t) is anon-equilibrium solutionto the system. (Clearly, the Lyapunov function must be constant on an equilibrium solution.) We can characterize continuously differentiable Lyapunov f unctions by using the ele- mentary calculus results that a scalar function is non-incr easing if and only if its derivative is non-negative, and is strictly decreasing if its derivati ve is strictly less than 0. We can 11/17/13 39 c/ci∇cleco√y∇t2013 Peter J. Olver compute the derivative of L(u(t)) by applying the same chain rule computation used to establish (4.14). As a result, we establish the basic criter ia for Lyapunov functions. Proposition 4.18. A continuously differentiable function L(u)is a Lyapunov func- tion for the system/squaresmallsolidu=F(u)if and only if it satisfies the Lyapunov inequality d dtL(u(t)) =∇L(u)·F(u)≤0 for all solutions u(t). (4.29) Furthermore, L(u)is a strict Lyapunov function if and only if d dtL(u(t)) =∇L(u)·F(u)<0 whenever F(u)/ne}ationslash=0. (4.30) The main result on stability and instability of equilibria o f a system that possesses a Lyapunov function follows. Theorem 4.19. LetL(u)be a Lyapunov function for the autonomous system/squaresmallsolidu= F(u). Ifu⋆is a strict local minimum of L, thenu⋆is a stable equilibrium point. If L(u)is a strict Lyapunov function, then u⋆is an asymptotically stable equilibrium. On the other hand, any critical point of a strict Lyapunov function L(u)which is not a local minimum is an unstable equilibrium point. In particular, local maxima of strict Lyapunov functions ar enotstable equilibria. In outline, the proof relies on the fact that the Lyapunov fun ction is non-increasing on solutions, and hence any solution that starts out sufficientl y near a minimum value can- not go far away, demonstrating stability. Figuratively, if you start near the bottom of a valley and never walk uphill you stay near the bottom. In the s trict case, the Lyapunov function must be strictly decreasing on the non-equilibriu m solutions, which thus must go to the minimum value of Last→ ∞. Again, starting near the bottom of a valley and always walking downhill takes you eventually to the bott om. On the other hand, if a non-equilibrium solution starts near an equilibrium poin t that is not a local minimum, then the fact that the Lyapunov function is steadily decreas ing implies that the solution must move further and further away from the equilibrium, whi ch suffices to demonstrate instability. Unlike first integrals, which can, at least in principle, be s ystematically constructed by solving a first order partial differential equation, findin g Lyapunov functions is much more of an art form, usually requiring some combination of ph ysical intuition and inspired guesswork. Example 4.20. Return to the planar system du dt=v,dv dt=−αsinu−βv, describing the damped oscillations of a pendulum, as in (4.1 0). Physically, we expect that the damping will cause a continual decrease in the total ener gy in the system, which, by (4.24), is E(u,v) =1 2mv2+κ(1−cosu). 11/17/13 40 c/ci∇cleco√y∇t2013 Peter J. Olver Let us prove that Eis, indeed, a Lyapunov function. We compute its time derivat ive, whenu(t),v(t) is a solution to the damped system. Recalling that α=κ/m,β=µ/m, we find dE dt=∂E ∂udu dt+∂E ∂vdv dt= (κsinu)v+(mv)(−αsinu−βv) =−µv2≤0, since we are assuming that the frictional coefficient µ >0. Therefore, the energy satisfies the Lyapunov stability criterion (4.29). Consequently, Th eorem 4.19 re-establishes the stability of the energy minima u⋆ 2k= 2kπ,v⋆ 2k= 0, where the damped pendulum is at the bottom of the arc. In fact, since dE/dt < 0 except when v= 0, with a little more work, the Lyapunov criterion can be used to establish their asympt otic stability. Hamiltonian and Poisson Systems One particularly important class of conservative systems w ere first introduced in the work of William Rowan Hamilton and Sim´ eon–Dennis Poisson i n the early nineteenth century. Their research was concerned with methods for solv ing the conservative systems arising in classical mechanics. Remarkably, Hamiltonian s ystems turned out to be the mathematical vehicle that unlocked the modern world of suba tomic quantum mechanics. Definition 4.21. APoisson system is a first order system of ordinary differential equations of the form /squaresmallsolidu=J∇H(u),where JT=−J (4.31) is a constant†skew-symmetric matrix. The real-valued function H(u) is known as the Hamiltonian function for the system. The Hamiltonian function is automatically a first integral f or the Poisson system (4.31). Indeed, to verify (4.14), we compute dH dt=∇H·du dt=∇H·F=∇HTF=∇HTJ∇H= 0. The final equality follows from the fact that Jis skew-symmetric, and hence vTJv= 0 for any vector v. In applications, the Hamiltonian function often represen ts the total energy of the system, and so a Poisson system necessarily conserves energy. Thus, friction and damping are not typically included in the Poisson framework . Conservation of the Hamiltonian function means that the sol utions of the Poisson system (4.31) move along its level sets {H(u) =c}. In particular, if a level set is a single point u0, then the corresponding solution cannot move and hence is an equilibrium solution: u(t)≡u0. Ifu0is an isolated maximum or minimum of H, then it is necessarily stable, since the nearby level sets remain close to u0and hence nearby solutions can never †Nonconstant Poisson matrices have additional, more subtle requiremen ts, [26; Chapter 6]. 11/17/13 41 c/ci∇cleco√y∇t2013 Peter J. Olver Figure 15. A Stable Planar Hamiltonian System. get far away from the equilibrium solution. In Figure 15 we pl ot the elliptical level sets of a typical stable quadratic Hamiltonian in the plane — the s olutions move periodically around the ellipses, all with the same period. The simplest example is when J=/parenleftbigg 0−1 1 0/parenrightbigg is a 2×2 matrix. In this case, we setu(t) = (p(t),q(t))T, and the Poisson system (4.31) corresponding to a Hamiltoni an function H(p,q) takes the form dp dt=−∂H ∂q,∂q ∂t=∂H ∂p. (4.32) For example, consider the undamped pendulum equation (4.23 ). Let us introduce the position variable q=θrepresenting the angular coordinate of the pendulum. Let p=m/squaresmallsolid θ be its angular momentum. The pendulum energy function (4.24 ) is rewritten in terms of the position and momentum variables: H(p,q) =p2 2m+κ(1−cosq). The resulting Hamiltonian system is (4.32) dp dt=κsinq,∂q ∂t=p m. (4.33) If we convert back to u=qand velocity v=/squaresmallsolid θ=p/mwe immediately recover the phase plane form (4.23) of the pendulum equation. Thus, we recover the constancy of the energy first integral directly from the general Hamiltonian framew ork. A Poisson system is called Hamiltonian if detJ/ne}ationslash= 0 is nonsingular. Since Jis a skew-symmetric matrix, this can only happen if its size, whi ch is also the dimension of 11/17/13 42 c/ci∇cleco√y∇t2013 Peter J. Olver the underlying space is even: n= 2k. The most important case, generalizing the two- dimensional version (4.32), is when J=/parenleftbigg O−I I O/parenrightbigg where O denotes the k×kzero matrix and I denotes the k×kidentity matrix. In this case, the variables are split into two vector components: u= (p,q)T= (p1,...,pk,q1,...,qk)T. We find ∇H(p,q) =/parenleftbigg∂H ∂p1, ...∂H ∂pk,∂H ∂q1, ...∂H ∂qk/parenrightbiggT and so (4.31) takes on the explicit form dp1 dt=−∂H ∂q1, ...dpk dt=−∂H ∂qk,dq1 dt=∂H ∂p1, ...dqk dt=∂H ∂pk.(4.34) In physical applications, the pvariables typically represent momenta, while the qvariables represent positions of the objects in the system; this was al ready noted in the pendulum example. Conservativemechanicalsystemscanalwaysberepresented inHamiltonianform, [ 28]. The prototypical example is to take the Hamiltonian functio n of the particular form H(p,q) =p2 1 2m1+···+p2 k 2mk+V(q1,...,qk), (4.35) where each mi>0 is a positive constant. The canonical system (4.34) has the form dp1 dt=−∂V ∂q1, ...dpk dt=−∂V ∂qk,dq1 dt=p1 m1, ...dqk dt=pk mk, which we can rewrite in vectorial form Mdq dt=p,dp dt=−∇V(q), (4.36) whereM= diag(m1,...,mn). Eliminating the pvariables, we find that (4.36) reduces to a second order system of ordinary differential equations in t he Newtonian form Md2q dt2=−∇V(q). The vector qrepresents the position vector for the mechanical system, w hileV(q) is the potential energy. The constants miare masses, and each pi=mi/squaresmallsolidqirepresents the momentum variable associated with the position variable qi. The Hamiltonian function is the total energy, the first term being the kinetic energy sinc e the term p2 i 2mi=mi 2/parenleftbiggdqi dt/parenrightbigg2 is precisely one half mass times velocity squared. 11/17/13 43 c/ci∇cleco√y∇t2013 Peter J. Olver Example 4.22. Consideraplanetarysystemconsistingof nbodiesinthree-dimensional space. We represent the planets as point masses, where the ithplanet is concentrated at position qi(t) = (xi(t),yi(t),zi(t))T. The planet’s linear momentum vector is pi(t) = mi/squaresmallsolidqi(t) =/parenleftbig/squaresmallsolidxi(t),/squaresmallsolidyi(t),/squaresmallsolidzi(t)/parenrightbigT, and its kinetic energy is /ba∇dblpi/ba∇dbl2 2mi=mi 2/ba∇dbl/squaresmallsolidqi/ba∇dbl2=mi 2/parenleftbig/squaresmallsolidx2 i+/squaresmallsolidy2 i+/squaresmallsolidz2 i/parenrightbig . The Newtonian gravitational potential between two point ma sses is proportional to the inverse square of their distance; more specifically V(qi,qj) =Gmimj /ba∇dblqi−qj/ba∇dbl2(4.37) represents the gravitational potential between planets iandj, where mi,mjare their masses, and Gthe universal gravitational constant. The Hamiltonian of t he system is the total kinetic plus potential energy, which is obtained by su mming the individual contribu- tions from (pairs of) planets: H(p,q) =n/summationdisplay i=1/ba∇dblpi/ba∇dbl2 2mi+/summationdisplay 1≤i<j≤nGmimj /ba∇dblqi−qj/ba∇dbl2. The resulting Hamiltonian system is known as the n–body problem , and it has been the subject of intense research since the time of Newton. For a tw o body system, e.g., a planet and a sun, the motion is well understood; the smaller body mov es around either an ellipse with the sum at a focus, or a parabola, or a hyperbola — the latt er two cases apply to comets. Even today, there are fundamental unanswered quest ions about the nature of solutions to this basic physical system for 3 or more bodies. Recent developments have been the discovery of systems of planets in which one or more g oes off to infinity in finite time, [20], and systems of planets that move around a figure 8 and even mo re complicated curves, [6,22]. Remark: For systems in Hamiltonian form, Theorem 4.13 tells us that the strict lo- cal minima (and maxima) of the Hamiltonian function are nece ssarily stable equilibrium points! This reconfirms our general observation. Saddle poi nts are typically unstable. The Hamiltonian form of a classical mechanical system plays a profound role in the constructionofitsquantummechanicsversion. Thereisaca nonicalquantizationprocedure that replaces the nonlinear Hamiltonian system of ordinary differential equations with a linear, quantum mechanics partialdifferential equationkn ownastheSchr¨ odinger equation. For example, the quantum mechanical system describing an el ectron orbiting a proton in a hydrogen atom is the quantized version of the classical 2–b ody problem describing the motion of a single planet around the sun. For details, we refe r the reader to any basic text in quantum mechanics, e.g., [ 10,19,21]. 11/17/13 44 c/ci∇cleco√y∇t2013 Peter J. Olver 5. Numerical Methods. Since we have no hope of solving the vast majority of different ial equations in explicit, analytic form, the design of suitable numerical algorithms for accurately approximating solutions is essential. The ubiquity of differential equati ons throughout mathematics and its applications has driven the tremendous research effort d evoted to numerical solution schemes, some dating back to the beginnings of the calculus. Nowadays, one has the luxury of choosing from a wide range of excellent software pa ckages that provide reliable and accurate results for a broad range of systems, at least fo r solutions over moderately long time periods. However, all of these packages, and the un derlying methods, have their limitations, and it is essential that one be able to to recogn ize when the software is working as advertised, and when it produces spurious results! Here i s where the theory, particularly the classification of equilibria and their stability proper ties, as well as first integrals and Lyapunov functions, can play an essential role. Explicit so lutions, when known, can also be used as test cases for tracking the reliability and accura cy of a chosen numerical scheme. In this section, we survey the most basic numerical methods f or solving initial value problems. For brevity, we shall only consider so-called sin gle step schemes, culminating in the very popular and versatile fourth order Runge–Kutta Met hod. This should only serve as a extremely basic introduction to the subject, and many ot her important and useful methods can be found in more specialized texts, [ 12,18]. It goes without saying that some equations are more difficult to accurately approximate than o thers, and a variety of more specialized techniques are employed when confronted with a recalcitrant system. But all of the more advanced developments build on the basic schemes and ideas laid out in this section. Euler’s Method The key issues confronting the numerical analyst of ordinar y differential equations already appear in the simplest first order ordinary different ial equation. Our goal is to calculate a decent approxiomation to the (unique) solution to the initial value problem du dt=F(t,u), u (t0) =u0. (5.1) Tokeepmatterssimple, wewillfocusour attentiononthesca larcase; however, allformulas and results written in a manner that can be readily adapted to first order systems — just replace the scalar functions u(t) andF(t,u) by vector-valued functions uandF(t,u) throughout. (The time t, of course, remains a scalar.) Higher order ordinary differe ntial equations are inevitably handled by first converting them in to an equivalent first order system, as discussed in Section 2, and then applying the nume rical scheme. The very simplest numerical solution method is named after L eonhard Euler — al- though Newton and his contemporaries were well aware of such a simple technique. Euler’s Method is rarely used in practice because much more efficient a nd accurate techniques can be implemented with minimal additional work. Nevertheless , the method lies at the core of the entire subject, and must be thoroughly understood bef ore progressing on to the more sophisticated algorithms that arise in real-world com putations. 11/17/13 45 c/ci∇cleco√y∇t2013 Peter J. Olver Starting at the initial time t0, we introduce successive mesh points (or sample times) t0< t1< t2< t3<···, continuing on until we reach a desired final time tn=t⋆. The mesh points should be fairly closely spaced. To keep the analysis as simple as possible, w e will always use a uniform step size, and so h=tk+1−tk>0, (5.2) does not depend on kand is assumed to be relatively small. This assumption serve s to simplify the analysis, and does not significantly affect the u nderlying ideas. For a uniform step size, the kthmesh point is at tk=t0+kh. More sophisticated adaptive methods, in which the step size is adjusted in order to maintain accuracy of the numerical solution, can be found in more specialized texts, e.g., [ 12,18]. Our numerical algorithm will recursively compute approximations uk≈u(tk), fork= 0,1,2,3,..., to the sampled values of the solution u(t) at the chosen mesh points. Our goal is to make the errorEk=uk−u(tk) in the approximation at each time tkas small as possible. If required, the values of the solution u(t) between mesh points may be computed by a subsequent interpo lation procedure, e.g., cubic splines. As you learned in first year calculus, the simplest approxima tion to a (continuously differentiable) function u(t) is provided by its tangent line or first order Taylor polynom ial. Thus, near the mesh point tk u(t)≈u(tk)+(t−tk)du dt(tk) =u(tk)+(t−tk)F(tk,u(tk)), in which we replace the derivative du/dtof the solution by the right hand side of the governing differential equation (5.1). In particular, the a pproximate value of the solution at the subsequent mesh point is u(tk+1)≈u(tk)+(tk+1−tk)F(tk,u(tk)). (5.3) This simple idea forms the basis of Euler’s Method. Since in practice we only know the approximation ukto the value of u(tk) at the current mesh point, we are forced to replace u(tk) by its approximation ukin the preceding formula. We thereby convert (5.3) into the iterative scheme uk+1=uk+(tk+1−tk)F(tk,uk). (5.4) In particular, when based on a uniform step size (5.2), Euler’s Method takes the simple form uk+1=uk+hF(tk,uk). (5.5) As sketched in Figure 16, the method starts off approximating the solution reasonably well, but gradually loses accuracy as the errors accumulate . Euler’s Method is the simplest example of a one-step numerical scheme for integrating anordinarydifferentialequation. Thisreferstothefactth atthesucceeding approximation, uk+1≈u(tk+1), depends only upon the current value, uk≈u(tk), which is one mesh point or “step” behind. 11/17/13 46 c/ci∇cleco√y∇t2013 Peter J. Olver t0t1t2t3u0u1u2u3u(t) Figure 16. Euler’s Method. To begin to understand how Euler’s Method works in practice, let us test it on a problem we know how to solve, since this will allow us to preci sely monitor the resulting errors in our numerical approximation to the solution. Example 5.1. The simplest “nontrivial” initial value problem is du dt=u, u (0) = 1, whose solution is, of course, the exponential function u(t) =et. SinceF(t,u) =u, Euler’s Method (5.5) with a fixed step size h >0 takes the form uk+1=uk+huk= (1+h)uk. This is a linear iterative equation, and hence easy to solve: uk= (1+h)ku0= (1+h)k is our proposed approximation to the solution u(tk) =etkat the mesh point tk=kh. Therefore, the Euler scheme to solve the differential equati on, we are effectively approxi- mating the exponential by a power function: etk=ekh≈(1+h)k When we use simply tto indicate the mesh time tk=kh, we recover, in the limit, a well-known calculus formula: et= lim h→0(1+h)t/h= lim k→∞/parenleftbigg 1+t k/parenrightbiggk . (5.6) A reader familiar with the computation of compound interest will recognize this particular approximation. As the time interval of compounding, h, gets smaller and smaller, the amount in the savings account approaches an exponential. How good is the resulting approximation? The error E(tk) =Ek=uk−etk 11/17/13 47 c/ci∇cleco√y∇t2013 Peter J. Olver measures the difference between the true solution and its num erical approximation at time t=tk=kh. Let us tabulate the error at the particular times t= 1,2 and 3 for various values of the step size h. The actual solution values are e1=e= 2.718281828 ... , e2= 7.389056096 ... , e3= 20.085536912 ... . In this case, the numerical approximation always underesti mates the true solution. h E(1) E(2) E(3) .1 −.125 −.662 −2.636 .01 −.0134 −.0730 −.297 .001 −.00135 −.00738 −.0301 .0001 −.000136 −.000739 −.00301 .00001 −.0000136 −.0000739 −.000301 Some key observations: •For a fixed step size h, the further we go from the initial point t0= 0, the larger the magnitude of the error. •On the other hand, the smaller the step size, the smaller the e rror at a fixed value of t. The trade-off is that more steps, and hence more computationa l effort†is required to produce the numerical approximation. For instance, we ne edk= 10 steps of sizeh=.1, butk= 1000 steps of size h=.001 to compute an approximation to u(t) at time t= 1. •The error is more or less in proportion to the step size. Decre asing the step size by a factor of1 10decreases the error by a similar amount, but simultaneously increases the amount of computation by a factor of 10. The final observation indicates that the Euler Method is of first order , which means that the error depends linearlyon the step size h. More specifically, at a fixed time t, the error is bounded by |E(t)|=|uk−u(t)| ≤C(t)h, when t=tk=kh, (5.7) for some positive C(t)>0 that depends upon the time, and the initial condition, but n ot on the step size. Example 5.2. The solution to the initial value problem du dt=/parenleftbig 1−4 3t/parenrightbig u, u (0) = 1, (5.8) †In this case, there happens to be an explicit formula for the numeri cal solution which can be used to bypass the iterations. However, in almost any other situation, one cannot compute the approximation ukwithout having first determined the intermediate values u0,...,uk−1. 11/17/13 48 c/ci∇cleco√y∇t2013 Peter J. Olver 0.5 11.5 22.5 30.20.40.60.811.21.4 h=.10.5 11.5 22.5 30.20.40.60.811.21.4 h=.01 Figure 17. Euler’s Method for/squaresmallsolidu=/parenleftbig 1−4 3t/parenrightbig u. was found in Example 2.3 by the method of separation of variab les: u(t) = exp/parenleftbig t−2 3t2/parenrightbig . (5.9) Euler’s Method leads to the iterative numerical scheme uk+1=uk+h/parenleftbig 1−4 3tk/parenrightbig uk, u0= 1. In Figure 17 we compare the graphs of the actual and numerical solutions for step sizes h=.1 and.01. In the former plot, we expliticly show the mesh points, bu t not in the latter, since they are too dense; moreover, the graphs of the numerical and true solutions are almost indistinguishable at this resolution. Thefollowingtableliststhenumerical errors E(tk) =uk−u(tk)betweenthecomputed and actual solution values u(1) = 1.395612425 ... , u (2) =.513417119 ... , u (3) =.049787068 ... , for several different step sizes: h E(1) E(2) E(3) .1000 .07461761 .03357536 −.00845267 .0100 .00749258 .00324416 −.00075619 .0010 .00074947 .00032338 −.00007477 .0001 .00007495 .00003233 −.00000747 As in the previous example, each decrease in step size by a fac tor of 10 leads to one additional decimal digit of accuracy in the computed soluti on. 11/17/13 49 c/ci∇cleco√y∇t2013 Peter J. Olver Taylor Methods In general, the order of a numerical solution method governs both the accuracy of its approximations and the speed at which they converge to the tr ue solution as the step size is decreased. Although the Euler Method is simple and easy to implement, it is only a first order scheme, and therefore of limited utility in serio us computations. So, the goal is to devise simple numerical methods that enjoy a much higher o rder of accuracy. Our derivation of the Euler Method was based on a first order Ta ylor approximation to the solution. So, an evident way to design a higher order me thod is to employ a higher order Taylor approximation. The Taylor series expansion fo r the solution u(t) at the succeeding mesh point tk+1=tk+hhas the form u(tk+1) =u(tk+h) =u(tk)+hdu dt(tk)+h2 2d2u dt2(tk)+h3 6d3u dt3(tk)+···.(5.10) As we just saw, we can evaluate the first derivative term throu gh use of the underlying differential equation: du dt=F(t,u). (5.11) The second derivative term can be found by differentiating th e equation with respect to t. Invoking the chain rule†, d2u dt2=d dtdu dt=d dtF(t,u(t)) =∂F ∂t(t,u)+∂F ∂u(t,u)du dt =∂F ∂t(t,u)+∂F ∂u(t,u)F(t,u)≡F(2)(t,u).(5.12) This operation is known as the total derivative , indicating that that we must treat the second variable uas a function of twhen differentiating. Substituting (5.11–12) into (5.10) and truncating at order h2leads to the Second Order Taylor Method uk+1=uk+hF(tk,uk)+h2 2F(2)(tk,uk) =uk+hF(tk,uk)+h2 2/parenleftbigg∂F ∂t(tk,uk)+∂F ∂u(tk,uk)F(tk,uk)/parenrightbigg ,(5.13) in which, as before, we replace the solution value u(tk) by its computed approximation uk. The resulting method is of second order, meaning that the er ror function satisfies the quadratic error estimate |E(t)|=|uk−u(t)| ≤C(t)h2when t=tk=kh. (5.14) †We assume throughout that Fhas as many continuous derivatives as needed. 11/17/13 50 c/ci∇cleco√y∇t2013 Peter J. Olver Example 5.3. Let us explicitly formulate the second order Taylor Method f or the initial value problem (5.8). Here du dt=F(t,u) =/parenleftbig 1−4 3t/parenrightbig u, d2u dt2=d dtF(t,u) =−4 3u+/parenleftbig 1−4 3t/parenrightbigdu dt=−4 3u+/parenleftbig 1−4 3t/parenrightbig2u, and so (5.13) becomes uk+1=uk+h/parenleftbig 1−4 3tk/parenrightbig uk+1 2h2/bracketleftbig −4 3uk+/parenleftbig 1−4 3tk/parenrightbig2uk/bracketrightbig , u0= 1. Thefollowingtableliststheerrorsbetween thevaluescomp uted by thesecond order Taylor scheme and the actual solution values, as given in Example 5. 2. h E(1) E(2) E(3) .100 .00276995 −.00133328 .00027753 .010 .00002680 −.00001216 .00000252 .001 .00000027 −.00000012 .00000002 Observe that, in accordance withthe quadratic error estima te (5.14), a decrease in the step size by a factor of1 10leads in an increase in accuracy of the solution by a factor1 100, i.e., an increase in 2 significant decimal places in the numerical a pproximation of the solution. Higher order Taylor methods are obtained by including furth er terms in the expansion (5.10). For example, to derive a third order Taylor method, w e include the third order term (h3/6)d3u/dt3in the Taylor expansion, where we evaluate the third derivat ive by differentiating (5.12), and so d3u dt3=d dtd2u dt2=d dtF(2)(t,u) =∂F(2) ∂t+∂F(2) ∂udu dt=∂F(2) ∂t+F∂F(2) ∂u =∂2F ∂t2+2F∂2F ∂t∂u+F2∂2F ∂u2+∂F ∂t∂F ∂u+F/parenleftbigg∂F ∂u/parenrightbigg2 ≡F(3)(t,u),(5.15) where we continue to make use of the fact that du/dt=F(t,u) is provided by the right hand side of the differential equation. The resulting third o rder Taylor method is uk+1=uk+hF(tk,uk)+h2 2F(2)(tk,uk)+h3 6F(3)(tk,uk), (5.16) where the last two summand are given by (5.12), (5.15), respe ctively. The higher order expressions are even worse, and a good symbolic manipulatio n system is almost essential for accurate computation. Whereas higher order Taylor methods are easy to motivate, th ey are rarely used in practice. There are two principal difficulties: 11/17/13 51 c/ci∇cleco√y∇t2013 Peter J. Olver •Owing to their dependence upon the partial derivatives of F(t,u), the right hand side of the differential equation needs to be rather smooth. •Even worse, the explicit formulae become exceedingly compl icated, even for relatively simplefunctions F(t,u). EfficientevaluationofthemultiplicityoftermsintheTay lor approximation and avoidance of round off errors become signi ficant concerns. As a result, mathematicians soon abandoned the Taylor serie s approach, and began to look elsewhere for high order, efficient integration methods. Error Analysis Before pressing on, we need to engage in a more serious discus sion of the error in a numerical scheme. A general one-step numerical method can be written in the form uk+1=G(h,tk,uk), (5.17) whereGis a prescribed function of the current approximate solutio n valueuk≈u(tk), the timetk, and the step size h=tk+1−tk, which, for illustrative purposes, we assume to be fixed. (We leave the discussion of multi-step methods , in which Gcould also depend upon the earlier values uk−1,uk−2,..., to more advanced texts, e.g., [ 12,18].) In any numerical integration scheme there are, in general, t hree sources of error. •The first is the local error committed in the current step of the algorithm. Even if we have managed to compute a completely accurate value of the so lutionuk=u(tk) at timetk, the numerical approximation scheme (5.17) is almost certa inly not exact, and will therefore introduce an error into the next computed valueuk+1≈u(tk+1). Round-off errors, resulting from the finite precision of the c omputer arithmetic, will also contribute to the local error. •The second is due to the error that is already present in the cu rrent approximation uk≈ u(tk). The local errors tend to accumulate as we continue to run th e iteration, and the net result is the global error , which is what we actually observe when compaing the numerical apporximation with the exact solution. •Finally, if the initial condition u0≈u(t0) is not computed accurately, this initial error will also make a contribution. For example, if u(t0) =π, then we introduce some initial error by using a decimal approximation, say π≈3.14159. The third error source will, for simplicity, be ignored in ou r discussion, i.e., we will assume u0=u(t0) is exact. Further, for simplicity we will assume that round -off errors do not play any significant role — although one must always keep them in mind when analyzing the computation. Since the global error is entirely due to th e accumulation of successive local errors, we must first understand the local error in deta il. Tomeasurethelocalerroringoingfrom tktotk+1,wecomparetheexactsolutionvalue u(tk+1) with its numerical approximation (5.17) under the assumpt ion that the current computed value is correct: uk=u(tk). Of course, in practice this is never the case, and so the local error is an artificial quantity. Be that as it may, in most circumstances the local error is ( a) easy to estimate, and, ( b) provides a reliable guide to the global accuracy of the numerical scheme. To estimate the local error, we assume that the step size his small 11/17/13 52 c/ci∇cleco√y∇t2013 Peter J. Olver and approximate the solution u(t) by its Taylor expansion† u(tk+1) =u(tk)+hdu dt(tk)+h2 2d2u dt2(tk)+··· =uk+hF(tk,uk)+h2 2F(2)(tk,uk)+h3 6F(3)(tk,uk)+···.(5.18) In the second expression, we have employed (5.12,15) and the ir higher order analogs to evaluate the derivative terms, and then invoked our local ac curacy assumption to replace u(tk) byuk. On the other hand, a direct Taylor expansion, in h, of the numerical scheme produces uk+1=G(h,tk,uk) =G(0,tk,uk)+h∂G ∂h(0,tk,uk)+h2 2∂2G ∂h2(0,tk,uk)+h3 6∂3G ∂h3(0,tk,uk)+···. (5.19) The local error is obtained by comparing these two Taylor exp ansions. Definition 5.4. A numerical integration method is of ordernif the Taylor expan- sions (5.18,19) of the exact and numerical solutions agree u p to order hn. For example, the Euler Method uk+1=G(h,tk,uk) =uk+hF(tk,uk), is already in the form of a Taylor expansion — that has no terms involving h2,h3,.... Comparing with the exact expansion (5.18), we see that the co nstant and order hterms are the same, but the order h2terms differ (unless F(2)≡0). Thus, according to the definition, the Euler Method is a first order method. Similarl y, the Taylor Method (5.13) is a second order method, because it was explicitly designed to match the constant, hand h2terms in the Taylor expansion of the solution. The general Ta ylor Method of order n setsG(h,tk,uk) to be exactly the order nTaylor polynomial, differing from the full Taylor expansion at order hn+1. Under fairly general hypotheses, it can be proved that, if th e numerical scheme has ordernas measured by the local error, then the global error is bounded by a multiple of hn. In other words, assuming no round-off or initial error, the c omputed value ukand the solution at time tkcan be bounded by |uk−u(tk)| ≤M hn, (5.20) where theconstant M >0 may depend onthe time tkand theparticular solution u(t). The error bound (5.20) serves to justify our numerical observat ions. For a method of order n, decreasing the step size by a factor of1 10will decrease the overall error by a factor of about 10−n, and so, roughly speaking, we anticipate gaining an additio nalndigits of accuracy — †In ouranalysis, weassume thatthedifferentialequation, and hencethe solution, has sufficient smoothness to justify the relevant Taylor approximation. 11/17/13 53 c/ci∇cleco√y∇t2013 Peter J. Olver at least up until the point that round-off errors begin to play a role. Readers interested in a more complete error analysis of numerical integration sche mes should consult a specialized text, e.g., [ 12,18]. The bottom line is the higher its order, the more accurate the numerical scheme, and hence thelargerthestepsizethatcanbeusedtoproducethes olutiontoadesiredaccuracy. However, this must be balanced with the fact that higher orde r methods inevitably require more computational effort at each step. If the total amount of computation has also decreased, then the high order method is to be preferred over a simpler, lower order method. Our goal now is to find another route to the design of hi gher order methods that avoids the complications inherent in a direct Taylor expans ion. An Equivalent Integral Equation The secret to the design of higher order numerical algorithm s is to replace the dif- ferential equation by an equivalent integral equation. By w ay of motivation, recall that, in general, differentiation is a badly behaved process; a rea sonable function can have an unreasonable derivative. On the other hand, integration am eliorates; even quite nasty functions have relatively well-behaved integrals. For the same reason, accurate numerical integration is relatively painless, whereas numerical diff erentiation should be avoided un- less necessary. While we have not dealt directly with integr al equations in this text, the subject has been extensively developed by mathematicians, [7], and has many important physical applications. Conversion of an initial value problem (5.1) to an integral e quation is straightforward. We integrate both sides of the differential equation from the initial point t0to a variable timet. The Fundamental Theorem of Calculus is used to explicitly e valuate the left hand integral: u(t)−u(t0) =/integraldisplayt t0/squaresmallsolidu(s)ds=/integraldisplayt t0F(s,u(s))ds. Rearranging terms, we arrive at the key result. Lemma 5.5. The solution u(t)to the the integral equation u(t) =u(t0)+/integraldisplayt t0F(s,u(s))ds (5.21) coincides with the solution to the initial value problemdu dt=F(t,u), u(t0) =u0. Proof: Our derivation already showed that the solution to the init ial value problem satisfies the integral equation (5.21). Conversely, suppos e thatu(t) solves the integral equation. Since u(t0) =u0is constant, the Fundamental Theorem of Calculus tells us th at the derivative of the right hand side of (5.21) is equal to the integrand, sodu dt=F(t,u(t)). Moreover, at t=t0, the upper and lower limits of the integral coincide, and so i t vanishes, whenceu(t) =u(t0) =u0has the correct initial conditions. Q.E.D. 11/17/13 54 c/ci∇cleco√y∇t2013 Peter J. Olver Left Endpoint Rule Trapezoid Rule Midpoint Rule Figure 18. Numerical Integration Methods. Observe that, unlike the differential equation, the integra l equation (5.21) requires no additional initial condition — it is automatically built in to the equation. The proofs of the fundamental existence and uniqueness Theorems 3.1 and 3 .3 for ordinary differential equations are, in fact, based on the integral equation refor mulation of the initial value problem; see [ 13,15] for details. The integral equation reformulation is equally valid for sy stems of first order ordinary differential equations. As noted above, u(t) andF(t,u(t)) become vector-valued func- tions. Integrating a vector-valued function is accomplish ed by integrating its individual components. Complete details are left to the exercises. Implicit and Predictor–Corrector Methods From this point onwards, we shall abandon the original initi al value problem, and turn our attention to numerically solving the equivalent in tegral equation (5.21). Let us rewrite the integral equation, starting at the mesh point tkinstead of t0, and integrating until time t=tk+1. The result is the basic integral formula u(tk+1) =u(tk)+/integraldisplaytk+1 tkF(s,u(s))ds (5.22) that (implicitly) computes the value of the solution at the s ubsequent mesh point. Com- paring this formula with the Euler Method uk+1=uk+hF(tk,uk),where h=tk+1−tk, and assuming for the moment that uk=u(tk) is exact, we discover that we are merely approximating the integral by /integraldisplaytk+1 tkF(s,u(s))ds≈hF(tk,u(tk)). (5.23) This is the Left Endpoint Rule for numerical integration — th at approximates the area under the curve g(t) =F(t,u(t)) between tk≤t≤tk+1by the area of a rectangle whose heightg(tk) =F(tk,u(tk))≈F(tk,uk) is prescribed by the left-hand endpoint of the graph. As indicated in Figure 18, this is a reasonable, but no t especially accurate method of numerical integration. 11/17/13 55 c/ci∇cleco√y∇t2013 Peter J. Olver Infirstyearcalculus,younodoubtencounteredmuchbetterm ethodsofapproximating the integral of a function. One of these is the Trapezoid Rule , which approximates the integralofthefunction g(t)bytheareaofatrapezoidobtainedbyconnectingthetwopoi nts (tk,g(tk))and(tk+1,g(tk+1))onthegraphof gbyastraightline, asinthesecondFigure18. Let us therefore try replacing (5.23) by the more accurate tr apezoidal approximation /integraldisplaytk+1 tkF(s,u(s))ds≈1 2h/bracketleftbig F(tk,u(tk))+F(tk+1,u(tk+1))/bracketrightbig . (5.24) Substituting this approximationinto the integral formula (5.22), and replacing the solution valuesu(tk),u(tk+1) by their numerical approximations, leads to the (hopefull y) more accurate numerical scheme uk+1=uk+1 2h/bracketleftbig F(tk,uk)+F(tk+1,uk+1)/bracketrightbig , (5.25) known as the Trapezoid Method . It is an implicit scheme , since the updated value uk+1 appears on both sides of the equation, and hence is only define d implicitly. Example 5.6. Consider the differential equation/squaresmallsolidu=/parenleftbig 1−4 3t/parenrightbig ustudied in Exam- ples 5.2 and 5.3. The Trapezoid Method with a fixed step size htakes the form uk+1=uk+1 2h/bracketleftbig/parenleftbig 1−4 3tk/parenrightbig uk+/parenleftbig 1−4 3tk+1/parenrightbig uk+1/bracketrightbig . In this case, we can explicit solve for the updated solution v alue, leading to the recursive formula uk+1=1+1 2h/parenleftbig 1−4 3tk/parenrightbig 1−1 2h/parenleftbig 1−4 3tk+1/parenrightbiguk=1+1 2h−2 3htk 1−1 2h+2 3h(tk+h)uk. (5.26) Implementing this scheme for three different step sizes give s the following errors between the computed and true solutions at times t= 1,2,3. h E(1) E(2) E(3) .100 −.00133315 .00060372 −.00012486 .010 −.00001335 .00000602 −.00000124 .001 −.00000013 .00000006 −.00000001 The numerical data indicates that the Trapezoid Method is of second order. For each reduction in step size by1 10, the accuracy in the solution increases by, roughly, a facto r of 1 100=1 102; that is, the numerical solution acquires two additional ac curate decimal digits. The main difficulty with the Trapezoid Method (and any other im plicit scheme) is immediately apparent. The updated approximate value for th e solution uk+1appears on both sides of the equation (5.25). Only for very simple funct ionsF(t,u) can one expect to solve (5.25) explicitly for uk+1in terms of the known quantities tk,ukandtk+1=tk+h. The alternative is to employ a numerical equation solver, su ch as the bisection algorithm 11/17/13 56 c/ci∇cleco√y∇t2013 Peter J. Olver or Newton’s Method, to compute uk+1. In the case of Newton’s Method, one would use the current approximation ukas a first guess for the new approximation uk+1. The resulting scheme requires some work to program, but can be effective in c ertain situations. An alternative, less involved strategy is based on the follo wing far-reaching idea. We already know a half-way decent approximation to the solutio n valueuk+1— namely that provided by the more primitive Euler scheme /tildewideuk+1=uk+hF(tk,uk). (5.27) Let’s use this estimated value in place of uk+1on the right hand side of the implicit equation (5.25). The result uk+1=uk+1 2h/bracketleftbig F(tk,uk)+F(tk+h,/tildewideuk+1)/bracketrightbig =uk+1 2h/bracketleftbig F(tk,uk)+F/parenleftbig tk+h,uk+hF(tk,uk)/parenrightbig/bracketrightbig .(5.28) is known as the Improved Euler Method . It is a completely explicit scheme since there is no need to solve any equation to find the updated value uk+1. Example 5.7. For our favorite equation/squaresmallsolidu=/parenleftbig 1−4 3t/parenrightbig u, the Improved Euler Method begins with the Euler approximation /tildewideuk+1=uk+h/parenleftbig 1−4 3tk/parenrightbig uk, and then replaces it by the improved value uk+1=uk+1 2h/bracketleftbig/parenleftbig 1−4 3tk/parenrightbig uk+/parenleftbig 1−4 3tk+1/parenrightbig /tildewideuk+1/bracketrightbig =uk+1 2h/bracketleftbig/parenleftbig 1−4 3tk/parenrightbig uk+/parenleftbig 1−4 3(tk+h)/parenrightbig/parenleftbig uk+h/parenleftbig 1−4 3tk/parenrightbig uk/parenrightbig/bracketrightbig =/bracketleftBig/parenleftbig 1−2 3h2/parenrightbig/bracketleftbig 1+h/parenleftbig 1−4 3tk/parenrightbig/bracketrightbig +1 2h2/parenleftbig 1−4 3tk/parenrightbig2/bracketrightBig uk. Implementing this scheme leads to the following errors in th e numerical solution at the indicated times. The Improved Euler Method performs compar ably to the fully implicit scheme (5.26), and significantly better than the original Eu ler Method. h E(1) E(2) E(3) .100 −.00070230 .00097842 .00147748 .010 −.00000459 .00001068 .00001264 .001 −.00000004 .00000011 .00000012 The Improved Euler Method is the simplest of a large family of so-called predictor– corrector algorithms . In general, one begins a relatively crude method — in this ca se the Euler Method — to predicta first approximation /tildewideuk+1to the desired solution value uk+1. One then employs a more sophisticated, typically implicit, method to correctthe original prediction, by replacing the required update uk+1on the right hand side of the implicit 11/17/13 57 c/ci∇cleco√y∇t2013 Peter J. Olver scheme by the less accurate prediction /tildewideuk+1. The resulting explicit, corrected value uk+1 will, provided the method has been designed with due care, be an improved approximation to the true solution. The numerical data in Example 5.7 indicates that the Improve d Euler Method is of second order since each reduction in step size by1 10improves the solution accuracy by, roughly, a factor of1 100. To verify this prediction, we expand the right hand side of ( 5.28) ina Taylorseries in h, and then compare, term by term, withthesolutionexpansion (5.18). First†, F/parenleftbig tk+h,uk+hF(tk,uk)/parenrightbig =F+h/parenleftbig Ft+F Fu/parenrightbig +1 2h2/parenleftbig Ftt+2F Ftu+F2Fuu/parenrightbig +···, where all the terms involving Fand its partial derivatives on the right hand side are evaluated at tk,uk. Substituting into (5.28), we find uk+1=uk+hF+1 2h2/parenleftbig Ft+F Fu/parenrightbig +1 4h3/parenleftbig Ftt+2F Ftu+F2Fuu/parenrightbig +···.(5.29) The two Taylor expansions (5.18) and (5.29) agree in their or der 1,handh2terms, but differ at order h3. This confirms our experimental observation that the Improv ed Euler Method is of second order. We can design a range of numerical solution schemes by implem enting alternative numerical approximations to the basic integral equation (5 .22). For example, the Midpoint Ruleapproximatestheintegralofthefunction g(t)bytheareaoftherectanglewhoseheight is the value of the function at the midpoint: /integraldisplaytk+1 tkg(s)ds≈hg/parenleftbig tk+1 2h/parenrightbig ,where h=tk+1−tk. (5.30) See Figure 18 for an illustration. The Midpoint Rule is known to have the same order of accuracy as the Trapezoid Rule, [ 4]. Substituting into (5.22) leads to the approximation uk+1=uk+/integraldisplaytk+1 tkF(s,u(s))ds≈uk+hF/parenleftbig tk+1 2h,u/parenleftbig tk+1 2h/parenrightbig/parenrightbig . Of course, we don’t know the value of the solution u/parenleftbig tk+1 2h/parenrightbig at the midpoint, but can predict it through a straightforward adaptation of the basi c Euler approximation: u/parenleftbig tk+1 2h/parenrightbig ≈uk+1 2hF(tk,uk). The result is the Midpoint Method uk+1=uk+hF/parenleftbig tk+1 2h,uk+1 2hF(tk,uk)/parenrightbig . (5.31) A comparison of the terms in the Taylor expansions of (5.18) a nd (5.31) reveals that the Midpoint Method is also of second order. †We use subscripts to indicate partial derivatives to save space. 11/17/13 58 c/ci∇cleco√y∇t2013 Peter J. Olver Runge–Kutta Methods The Improved Euler and Midpoint Methods are the most element ary incarnations of a general class of numerical schemes for ordinary differenti al equations that were first sys- tematically studied by the German mathematicians Carle Run ge and Martin Kutta in the late nineteenth century. Runge–Kutta Methods are by far the most popular and powerful general-purpose numerical methods for integrating ordina ry differential equations. While not appropriate in all possible situations, Runge–Kutta sc hemes are surprisingly robust, performing efficiently and accurately in a wide variety of pro blems. Barring significant complications, they are the method of choice in most basic ap plications. They comprise the engine that powers most computer software for solving ge neral initial value problems for systems of ordinary differential equations. The most general Runge–Kutta Method takes the form uk+1=uk+hm/summationdisplay i=1ciF(ti,k,ui,k), (5.32) wheremcounts the number of termsin the method. Each ti,kdenotes a point in the kth mesh interval, so tk≤ti,k≤tk+1. The second argument ui,k≈u(ti,k) can be viewed as an approximation to the solution at the point ti,k, and so is computed by a similar, but simpler formula of the same type. There is a lot of flexibil ity in the design of the method, through choosing the coefficients ci, the times ti,k, as well as the scheme (and all parameters therein) used to compute each of the intermed iate approximations ui,k. As always, the orderof the method is fixed by the power of hto which the Taylor expansions of the numerical method (5.32) and the actual solution (5.18 ) agree. Clearly, the more terms we include in the Runge–Kutta formula (5.32), the more free parameters available to match terms in the solution’s Taylor series, and so the hig her the potential order of the method. Thus, the goalis toarrange theparameters so that th emethod hasa highorder of accuracy, while, simultaneously, avoiding unduly complic ated, and hence computationally costly, formulae. Both the Improved Euler and Midpoint Methods are instances o f a family of two term Runge–Kutta Methods uk+1=uk+h/bracketleftbig aF(tk,uk)+bF/parenleftbig tk,2,uk,2/parenrightbig/bracketrightbig =uk+h/bracketleftbig aF(tk,uk)+bF/parenleftbig tk+λh,uk+λhF(tk,uk)/parenrightbig/bracketrightbig ,(5.33) based on the current mesh point, so tk,1=tk, and one intermediate point tk,2=tk+λh with 0≤λ≤1. The basic Euler Method is used to approximate the solution value uk,2=uk+λhF(tk,uk) attk,2. The Improved Euler Method sets a=b=1 2andλ= 1, while the Midpoint Method corresponds to a= 0, b= 1, λ=1 2. The range of possible values for a,bandλis found by matching the Taylor expansion uk+1=uk+h/bracketleftbig aF(tk,uk)+bF/parenleftbig tk+λh,uk+λhF(tk,uk)/parenrightbig/bracketrightbig =uk+h(a+b)F(tk,uk)+h2bλ/bracketleftbigg∂F ∂t(tk,uk)+F(tk,uk)∂F ∂u(tk,uk)/bracketrightbigg +···. 11/17/13 59 c/ci∇cleco√y∇t2013 Peter J. Olver (in powers of h) of the right hand side of (5.33) with the Taylor expansion (5 .18) of the solution, namely u(tk+1) =uk+hF(tk,uk)+h2 2[Ft(tk,uk)+F(tk,uk)Fu(tk,uk)]+···, to as high an order as possible. First, the constant terms, uk, are the same. For the order hand order h2terms to agree, we must have, respectively, a+b= 1, bλ =1 2. Therefore, setting a= 1−µ, b =µ,andλ=1 2µ,whereµis arbitrary†, leads to the following family of two term, second order Runge –Kutta Methods: uk+1=uk+h/bracketleftbigg (1−µ)F(tk,uk)+µF/parenleftbigg tk+h 2µ,uk+h 2µF(tk,uk)/parenrightbigg/bracketrightbigg .(5.34) The case µ=1 2corresponds to the Improved Euler Method (5.28), while µ= 1 yields the Midpoint Method (5.31). Unfortunately, none of these me thods are able to match all of the third order terms in the Taylor expansion for the so lution, and so we are left with a one-parameter family of two step Runge–Kutta Methods , all of second order, that include the Improved Euler and Midpoint schemes as particul ar instances. The methods with1 2≤µ≤1 all perform more or less comparably, and there is no special reason to prefer one over the other. To construct a third order Runge–Kutta Method, we need to tak e at least m≥3 terms in (5.32). A rather intricate computation (best done w ith the aid of computer algebra) will produce a range of valid schemes; the results c an be found in [ 12,18]. The algebraic manipulations are rather tedious, and we leave a c omplete discussion of the available options to a more advanced treatment. In practica l applications, a particularly simple fourth order, four term formula has become the most us ed. The method, often abbreviated as RK4, takes the form uk+1=uk+h 6/bracketleftbig F(tk,uk)+2F(t2,k,u2,k)+2F(t3,k,u3,k)+F(t4,k,u4,k)/bracketrightbig ,(5.35) where the times and function values are successively comput ed according to the following procedure: t2,k=tk+1 2h, u2,k=uk+1 2hF(tk,uk), t3,k=tk+1 2h, u3,k=uk+1 2hF(t2,k,u2,k), t4,k=tk+h, u4,k=uk+hF(t3,k,u3,k).(5.36) The four term RK4 scheme (5.35–36) is, in fact, a fourth order method. This is confirmed by demonstrating that the Taylor expansion of the right hand side of (5.35) in powers of †Although we should restrict µ≥1 2in order that 0 ≤λ≤1. 11/17/13 60 c/ci∇cleco√y∇t2013 Peter J. Olver hmatches all of the terms in the Taylor series for the solution (5.18) up to and including those of order h4, and hence the local error is of order h5. This is not a computation for the faint-hearted — bring lots of paper and erasers, or, better y et, a good computer algebra package! The RK4 scheme is but one instance of a large family o f fourth order, four term Runge–Kutta Methods, and by far the most popular owing to its relative simplicity. Example 5.8. Application of the RK4 Method (5.35–36) to our favorite init ial value problem (5.8) leads to the following errors at the indicated times: h E(1) E(2) E(3) .100 −1.944×10−71.086×10−64.592×10−6 .010 −1.508×10−111.093×10−103.851×10−10 .001 −1.332×10−15−4.741×10−141.932×10−14 The accuracy is phenomenally good — much better than any of ou r earlier numerical schemes. Each decrease in the step size by a factor of1 10results in about 4 additional decimal digits of accuracy in the computed solution, in comp lete accordance with its status as a fourth order method. Actually, it is not entirely fair to compare the accuracy of t he methods using the same step size. Each iteration of the RK4 Method requires fou r evaluations of the func- tionF(t,u), and hence takes the same computational effort as four Euler iterations, or, equivalently, two Improved Euler iterations. Thus, the mor e revealing comparison would be between RK4 at step size h, Euler at step size1 4h, and Improved Euler at step size1 2h, as these involve roughly the same amount of computational eff ort. The resulting errors E(1) at time t= 1 are listed in the following table. Thus, even taking computational effort into account, the Run ge–Kutta Method con- tinues to outperform its rivals. At a step size of .1, it is almost as accurate as the Im- proved Euler Method with step size .0005, and hence 200 times as much computation, while the Euler Method would require a step size of approxima tely.24×10−6, and would be 4,000,000 times as slow as Runge–Kutta! With a step size of .001, RK4 computes a solution value that is near the limits imposed by machine acc uracy (in single precision arithmetic). The superb performance level and accuracy of t he RK4 Method immediately explains its popularity for a broad range of applications. h Euler Improved Euler Runge–Kutta 4 .11.872×10−2−1.424×10−4−1.944×10−7 .011.874×10−3−1.112×10−6−1.508×10−11 .001 1.870×10−4−1.080×10−8−1.332×10−15 11/17/13 61 c/ci∇cleco√y∇t2013 Peter J. Olver 1234567812345678 Euler Method, h=.011234567812345678 Euler Method, h=.001 1234567812345678 Improved Euler Method, h=.011234567812345678 RK4 Method, h=.01 Figure 19. Numerical Solutions of Predator–Prey Model. Example 5.9. Asnotedearlier, by writingthe function valuesasvectors uk≈u(tk), one can immediately use all of the preceding methods to integ rate initial value problems for first order systems of ordinary differential equations/squaresmallsolidu=F(u). Consider, by way of example, the Lotka–Volterra system du dt= 2u−uv,dv dt=−9v+3uv, (5.37) analyzed in Example 4.14. To find a numerical solution, we wri teu= (u,v)Tfor the solutionvector, while F(u) = (2u−uv,−9v+3uv)Tistherighthandsideofthesystem. The Euler Method with step size his given by u(k+1)=u(k)+hF(u(k)), or, explicitly, as a first order nonlinear iterative system u(k+1)=u(k)+h(2u(k)−u(k)v(k)), v(k+1)=v(k)+h(−9v(k)+3u(k)v(k)). 11/17/13 62 c/ci∇cleco√y∇t2013 Peter J. Olver 5 10 15 20 25 -0.8-0.6-0.4-0.2 Euler Method h=.00110 20 30 40 50 -0.004-0.003-0.002-0.001 Improved Euler Method h=.0110 20 30 40 50 -0.00001-7.5·10-6-5·10-6-2.5·10-62.5·10-65·10-6 RK4 Method h=.01 Figure 20. Numerical Evaluation of Lotka–Volterra First Integral. The Improved Euler and Runge–Kutta schemes are implemented in a similar fashion. Phase portraits of the three numerical algorithms starting with initial conditions u(0)= v(0)= 1.5, and up to time t= 25 in the case of the Euler Method, and t= 50 for the other two, appear in Figure 19. Recall that the solution is su pposed to travel periodically around a closed curve, which is the level set I(u,v) = 9 log u−3u+2 logv−v=I(1.5,1.5) =−1.53988 of the first integral. The Euler Method spirals away from the e xact periodic solution, whereas the Improved Euler and RK4 Methods perform rather we ll. Since we do not have an analyticformula for the solution, we arenot ableto measu re the error exactly. However, the known first integral is supposed to remain constant on the solution trajectories, and so one means of monitoring the accuracy of the solution is to t rack the variation in the numerical values of I(u(k),v(k)). These are graphed in Figure 20; the Improved Euler keeps the value within .0005, while in the RK4 solution, the first integral only exper iences change in its the fifth decimal place over the indicated time p eriod. Of course, the ononger one continues to integrate, the more error will gradually cr eep into the numerical solution. Still, for most practical purposes, the RK4 solution is indi stinguishable from the exact solution. In practical implementations, it is important to monitor th e accuracy of the numer- ical solution, so to gauge when to abandon an insufficiently pr ecise computation. Since accuracy is dependent upon the step size h, one may try adjusting hso as stay within a preassigned error. Adaptive methods , allow one to change the step size during the course of the computation, in response to some estimation of the ove rall error. Insufficiently ac- curate numerical solutions would necessitate a suitable re duction in step size (or increase in the order of the scheme). On the other hand, if the solution is more accurate than the application requires, one could increase the step size so as to reduce the total amount of computational effort. How might one decide when a method is giving inaccurate resul ts, since one presum- ably does not know the true solution and so has nothing to dire ctly test the numerical approximation against? One useful idea is to integrate the d ifferential equation using two different numerical schemes, usually of different orders of a ccuracy, and then compare the results. If the two solution values are reasonably close, th en one is usually safe in as- suming that the methods are both giving accurate results, wh ile in the event that they 11/17/13 63 c/ci∇cleco√y∇t2013 Peter J. Olver differ beyond some preassigned tolerance, then one needs to r e-evaluate the step size. The required adjustment to the step size relies on a more detaile d analysis of the error terms. Several well-studied methods are employed in practical sit uations; the most popular is the Runge–Kutta–Fehlberg Method, which combines a fourth and a fifth order Runge–Kutta scheme for error control. Details can be found in more advanc ed treatments of the subject, e.g., [12,18]. Stiff Differential Equations While the fourth order Runge–Kutta Method with a sufficiently small step size will successfully integrateabroadrangeofdifferentialequati ons—atleastovernotundulylong time intervals — it does occasionally experience unexpecte d difficulties. While we have not developed sufficiently sophisticated analytical tools t o conduct a thorough analysis, it will be instructive to look at why a breakdown might occur in a simpler context. Example 5.10. The elementary linear initial value problem du dt=−250u, u (0) = 1, (5.38) is an instructive and sobering example. The explicit soluti on is easily found; it is a very rapidly decreasing exponential: u(t) =e−250t. u(t) =e−250twith u(1)≈2.69×10−109. Thefollowingtablegivestheresultofapproximatingtheso lutionvalue u(1)≈2.69×10−109 at timet= 1 using three of our numerical integration schemes for vari ous step sizes: h Euler Improved Euler RK4 .16.34×10133.99×10242.81×1041 .014.07×10171.22×10211.53×10−19 .001 1.15×10−1256.17×10−1082.69×10−109 The results are not misprints! When the step size is .1, the computed solution values are perplexingly large, and appear to represent an exponent ially growing solution — the complete opposite of the rapidly decaying true solution. Re ducing the step size beyond a critical threshold suddenly transforms the numerical solu tion to an exponentially decaying function. Only the fourth order RK4 Method with step size h=.001 — and hence a total of 1,000 steps — does a reasonable job at approximating the correc t value of the solution att= 1. You may well ask, what on earth is going on? The solution could n’t be simpler — why is it so difficult to compute it? To understand the basic issue, let us analyze how the Euler Method handles such simple differential equations. Conside r the initial value problem du dt=λu, u (0) = 1, (5.39) 11/17/13 64 c/ci∇cleco√y∇t2013 Peter J. Olver which has an exponential solution u(t) =eλt. (5.40) As in Example 5.1, the Euler Method with step size hrelies on the iterative scheme uk+1= (1+λh)uk, u0= 1, with solution uk= (1+λh)k. (5.41) Ifλ >0, theexact solution(5.40)isexponentiallygrowing. Sinc e1+λh >1, thenumerical iteratesarealsogrowing, albeitat a somewhat slower rate. Inthiscase, there isno inherent surprise with the numerical approximation procedure — in th e short run it gives fairly accurate results, but eventually trails behind the exponen tially growing solution. On the other hand, if λ <0, then the exact solution is exponentially decaying and positive. But now, if λh <−2, then 1 + λh <−1, and the iterates (5.41) grow expo- nentially fast in magnitude, with alternating signs. In thi s case, the numerical solution is nowhere close to the true solution; this explains the prev iously observed pathological behavior. If −1<1+λh <0, the numerical solutions decay in magnitude, but continue to alternate between positive and negative values. Thus, to correctly model the qualitative features of the solution and obtain a numerically respectab le approximation, we need to choose the step size hso as to ensure that 0 <1+λh, and hence h <−1/λwhenλ <0. For the value λ=−250 in the example, then, we must choose h <1 250=.004 in order that the Euler Method give a reasonable numerical answer. A simil ar, but more complicated analysis applies to any of the Runge–Kutta schemes. Thus, the numerical methods for ordinary differential equat ions exhibit a form of conditional stability. Paradoxically, the larger negativ eλis — and hence the faster the solutiontendstoatrivialzeroequilibrium—the moredifficultandexpensivethenumerical integration. The system (5.38) is the simplest example of what is known as a stiff differential equation. In general, an equation or system is stiff if it has one or more very rapidly decayingsolutions. Inthecaseofautonomous(constantcoe fficient)linearsystems/squaresmallsolidu=Au, stiffness occurs whenever the coefficient matrix Ahas an eigenvalue with a large negative real part: Re λ≪0, resulting in a very rapidly decaying eigensolution. It on ly takes one such eigensolution to render the equation stiff, and ruin the numerical computation of even the well behaved solutions! Curiously, the component of the actual solution corresponding to such large negative eigenvalues is almost irrelevant, as it becomes almost instanteously tiny. However,thepresenceofsuchaneigenvaluecontinues torenderthenumericalsolution to the system very difficult, even to the point of exhausting an y available computing resources. Stiff equations require more sophisticated nume rical procedures to integrate, and we refer the reader to [ 12,18] for details. Most of the other methods derived above also suffer from insta bility due to stiffness of the ordinary differential equation for sufficiently large neg ativeλ. Interestingly, stability for solvingthetrivialtest scalarordinarydifferentialeq uation(5.39)suffices tocharacterize acceptable step sizes h, depending on the size of λ, which, in the case of systems, becomes the eigenvalue. The analysis is not so difficult, owing to the i nnate simplicity of the test 11/17/13 65 c/ci∇cleco√y∇t2013 Peter J. Olver ordinary differential equation (5.39). A significant except ion, which also illustrates the test for behavior under rapidly decaying solutions, is the T rapezoid Method (5.25). Let us analyzethebehavio of theresulting numerical solutionto ( 5.39). Substituting f(t,u) =λu into the Trapezoid iterative equation (5.25), we find uk+1=uk+1 2h/bracketleftbig λuk+λuk+1/bracketrightbig , which we solve for uk+1=1+1 2hλ 1−1 2hλuk≡µuk. Thus, the behavior of the solution is entirely determined by the size of the coefficient µ. Ifλ >0, thenµ >1 and the numerical solution is exponentially growing, as lo ng as the denominator is positive, which requires h <2/λto be sufficiently small. In other words, rapidly growing exponential solutions require reasonably small step sizes to accurately compute, which is not surprising. On the other hand, if λ <0, then|µ|<1,no matter how large negative λgets! (But we should also restrict h <−2/λto be sufficiently small, as otherwise µ <0 and the numerical solution would have oscillating signs, e ven though it is decaying, and hence vanishing small. If this were part o f a larger system, such minor oscillations would not worry us because they would be unnoti ceable in the long run.) Thus, the Trapezoid Method is notaffected by very large negative exponents, and hence not subject to the effects of stiffness. The TrapezoidMethod is the simplest example of an Astable method. More precisely, a numerical solution method is called Astableif the zero solution is asymptotically stable for the iterative equation resulting from the numerical sol ution to the ordinary differential equation/squaresmallsolidu=λufor allλ <0. The big advantage of Astable methods is that they are not affected by stiffness. Unfortunately, Astable methods are few and far between. In fact, they are all implicit one-step methods! No explicit Runge–Kutta Method is Astable; see [18] for a proof of this disappointing result. Moreover, multis tep methods also suffer from the lack of Astability and so are all prone to the effects of stiffness. Stil l, when confronted with a seriously stiff equation, one should discard the sophi sticated methods and revert to a low order, but Astable scheme like the Trapezoid Method. 11/17/13 66 c/ci∇cleco√y∇t2013 Peter J. Olver References [1] Alligood, K.T., Sauer, T.D., and Yorke, J.A., Chaos. An Introduction to Dynamical Systems, Springer-Verlag, New York, 1997. [2] Birkhoff, G., and Rota, G.–C., Ordinary Differential Equations , Blaisdell Publ. Co., Waltham, Mass., 1962. [3] Bronstein, M., Symbolic integration I :Transcendental Functions , Springer–Verlag, New York, 1997. [4] Burden, R.L., and Faires, J.D., Numerical Analysis , Seventh Edition, Brooks/Cole, Pacific Grove, CA, 2001. [5] Cantwell, B.J., Introduction to Symmetry Analysis , Cambridge University Press, Cambridge, 2003. [6] Chenciner, A., and Montgomery, R., A remarkable periodic s olution of the three-body problem in the case of equal masses, Ann. Math. 152(2000), 881–901. [7] Courant, R., and Hilbert, D., Methods of Mathematical Physics , vol. I, Interscience Publ., New York, 1953. [8] Devaney, R.L., An Introduction to Chaotic Dynamical Systems , Addison–Wesley, Redwood City, Calif., 1989. [9] Diacu, F., An Introduction to Differential Equations , W.H. Freeman and Co., New York, 2000. [10] Dirac, P.A.M., The Principles of Quantum Mechanics , 3rd ed., Clarendon Press, Oxford, 1947. [11] Goldstein, H., Classical Mechanics , Second Edition, Addison–Wesley, Reading, Mass., 1980. [12] Hairer, E., Nørsett, S.P., and Wanner, G., Solving Ordinary Differential Equations , 2nd ed., Springer–Verlag, New York, 1993–1996. [13] Hale, J.K., Ordinary Differential Equations , Second Edition, R.E. Krieger Pub. Co., Huntington, N.Y., 1980. [14] Hille, E., Ordinary Differential Equations in the Complex Domain , John Wiley & Sons, New York, 1976. [15] Hirsch, M.W., and Smale, S., Differential Equations, Dynamical Systems, and Linear Algebra , Academic Press, New York, 1974. [16] Hydon, P.E., Symmetry Methods for Differential Equations , Cambridge Texts in Appl. Math., Cambridge University Press, Cambridge, 2000. [17] Ince, E.L., Ordinary Differential Equations , Dover Publ., New York, 1956. [18] Iserles, A., A First Course in the Numerical Analysis of Differential Equa tions, Cambridge University Press, Cambridge, 1996. [19] Landau, L.D., and Lifshitz, E.M., Quantum Mechanics (Non-relativistic Theory ), Course of Theoretical Physics, vol. 3, Pergamon Press, New Y ork, 1977. 11/17/13 67 c/ci∇cleco√y∇t2013 Peter J. Olver [20] Mather, J.N., and McGehee, R., Solutions of the collinear four body problem which become unbounded in finite time , Dynamical Systems, Theory and Applications; J. Moser, ed., Springer, Berlin, 1975, pp. 573–597.. [21] Messiah, A., Quantum Mechanics , John Wiley & Sons, New York, 1976. [22] Montaldi, J., and Steckles, K., Classification of symmetry groups for planar n-body choreographies, preprint, 2013, arXiv:1305.0470v2 . [23] Moon, F.C., Chaotic Vibrations , John Wiley & Sons, New York, 1987. [24] Noether, E., Invariante Variationsprobleme, Nachr. K¨ onig. Gesell. Wissen. G¨ ottingen, Math.–Phys. Kl. (1918), 235–257. (See Kosmann-Schwarzbach, Y.,The Noether Theorems. Invariance and Conservation Laws in t he Twentieth Century, Springer, New York, 2011, for an English translation.) [25] Olver, F.W.J., Lozier, D.W., Boisvert, R.F., and Clark, C. W., eds., NIST Handbook of Mathematical Functions , Cambridge University Press, Cambridge, 2010. [26] Olver, P.J., Applications of Lie Groups to Differential Equations , 2nd ed., Graduate Texts in Mathematics, vol. 107, Springer–Verlag, New York, 1993. [27] Olver, P.J., and Shakiban, C., Applied Linear Algebra , Prentice–Hall, Inc., Upper Saddle River, N.J., 2006. [28] Whittaker, E.T., A Treatise on the Analytical Dynamics of Particles and Rigid Bodies, Cambridge University Press, Cambridge, 1937. 11/17/13 68 c/ci∇cleco√y∇t2013 Peter J. Olver