Home / Math and Physics Files / Math / ODEs / Peter Olver Notes

ne

PDF · 57 pages · 2.2 MB

Textbook chapter by Peter J. Olver (2012 draft), filed among the ODE notes in the archive. It introduces nonlinear discrete dynamical systems u(k+1)=g(u(k)), fixed points and their stability, the cosine iteration example, and the scalar affine case with its convergence criterion |a|<1. The chapter outline also covers bisection, Newton's method, and gradient-descent optimization. Only the opening portion was read.

AI-written summary; may contain errors.

Extracted text (machine-read; may contain errors)

Chapter 19 NonlinearSystems Nonlinearity is ubiquitous in physical phenomena. Fluid an d plasma mechanics, gas dynamics, elasticity, relativity, chemical reactions, co mbustion, ecology, biomechanics, and many, manyother phenomena areallgoverned byinherently no nlinear equations. (Theone notable exception is quantum mechanics, which is a fundamen tally linear theory. Recent attempts at grand uniﬁcation of all fundamental physical th eories, such as string theory and conformal ﬁeld theory, [ 85], do venture into the nonlinear wilderness.) For this reaso n, an ever increasing proportion of modern mathematical resea rch is devoted to the analysis of nonlinear systems. Why, then, have we devoted the overwhelming majority of this text to linear mathe- matics? The facile answer is that nonlinear systems are vast ly more diﬃcult to analyze. In the nonlinear regime, many of the most basic questions remai n unanswered: existence and uniqueness of solutions arenot guaranteed; explicit formu lae arediﬃcult to come by; linear superposition is no longer available; numerical approxima tions are not always suﬃciently accurate; etc., etc. A more intelligent answer is that a thor ough understanding of linear phenomena and linear mathematics isan essential prerequis itefor progress in the nonlinear arena. Therefore, in this introductory text on applied math ematics, we have no choice but to ﬁrst develop the proper linear foundations in suﬃcient de pth before we can realistically confront the untamed nonlinear wilderness. Moreover, many important physical systems are “weakly nonlinear”, in the sense that, while nonlinear e ﬀects do play an essential role, the linear terms tend to dominate the physics, and so, to a ﬁrs t approximation, the system is essentially linear. As a result, such nonlinear phenomen a are best understood as some form of perturbation of their linear approximations. The tr uly nonlinear regime is, even today, only sporadically modeled and even less well underst ood. The advent of powerful computers has fomented a veritable re volution in our under- standing ofnonlinear mathematics. Indeed, manyofthemost importantmodernanalytical techniques drew their inspiration from early computer fora ys into the uncharted nonlinear wilderness. However, despite dramatic advances in both har dware capabilities and so- phisticated mathematical algorithms, many nonlinear syst ems — for instance, Einsteinian gravitation — still remain beyond the capabilities of today ’s computers. Space limitations restrict us to providing just a brief over view of some of the most important ideas, mathematical techniques, and new physica l phenomena that arise when venturing into the nonlinear realm. In this chapter, we star t with iteration of nonlinear functions. Building on our experience with iterative linea r systems, as developed in Chap- ter 10, we will discover that functional iteration, when it c onverges, provides a powerful mechanism for solving equations and for optimization. On th e other hand, even very 12/11/12 1024 c/ci∇cleco√y∇t2012 Peter J. Olver simple nonconvergent nonlinear iterative systems may admi t remarkably complex, chaotic behavior. The second section is devoted to basic solution te chniques for nonlinear systems, and includes bisection, general iteration, and the very pow erful NewtonMethod. The third section is devoted to ﬁnite-dimensional optimization prin ciples, i.e., the minimization or maximization of nonlinear functions. Numerical optimizat ion procedures rely on iterative procedures, and we concentrate on those associated with gra dient descent. 19.1. Iteration of Functions. AsﬁrstnotedinChapter 10, iteration, meaning repeatedapp licationofafunction, can be viewed as a discrete dynamical system in which the continuous time variable has been “quantized” to assume integer values. Even iterating a very simple quadratic scalar func- tion can lead to an amazing variety of dynamical phenomena, i ncluding multiply-periodic solutions and genuine chaos. Nonlinear iterative systems a rise not just in mathematics, but also underlie the growth and decay of biological populat ions, predator-prey interac- tions, spread ofcommunicable diseases such as Aids, and host ofother natural phenomena. Moreover, many numerical solution methods — for systems of a lgebraic equations, ordi- nary diﬀerential equations, partial diﬀerential equation s, and so on — rely on iteration, and so the theory underlies the analysis of convergence and e ﬃciency of such numerical approximation schemes. In general, an iterative system has the form u(k+1)=g(u(k)), (19.1) whereg:Rn→Rnis a real vector-valued function. (One can similarly treat i teration of complex-valued functions g:Cn→Cn, but, for simplicity, we only deal with real systems here.) A solution is a discrete collection of points†u(k)∈Rn, in which the index k= 0,1,2,3,...takes on non-negative integer values. Chapter 10 dealt with the case when g(u) =Auis a linear function, necessarily given by multiplication b y ann×nmatrixA. In this chapter, we enlarge our scope to the nonlinear case. Once we specify the initial iterate, u(0)=c, (19.2) then the resulting solution to the discrete dynamical syste m (19.1) is easily computed: u(1)=g(u(0)) =g(c),u(2)=g(u(1)) =g(g(c)),u(3)=g(u(2)) =g(g(g(c))), ... and so on. Thus, unlike continuous dynamical systems, the ex istence and uniqueness of solutions is not an issue. As long as each successive iterate u(k)lies in the domain of deﬁnition of gone merely repeats the process to produce the solution, u(k)=ktimes/bracehtipdownleft/bracehtipupright/bracehtipupleft/bracehtipdownright g◦···◦g(c), k = 0,1,2,...,(19.3) †The superscripts on u(k)refer to the iteration number, and do notdenote derivatives. 12/11/12 1025 c/ci∇cleco√y∇t2012 Peter J. Olver which is obtained by composing the function gwith itself a total of ktimes. In other words, the solution to a discrete dynamical system correspo nds to repeatedly pushing the gkey on your calculator. For example, entering 0 and then repe atedly hitting the coskey corresponds to solving the iterative system u(k+1)= cosu(k), u(0)= 0. (19.4) The ﬁrst 10 iterates are displayed in the following table: k0 1 2 3 4 5 6 7 8 9 u(k)0 1.540302 .857553 .65429.79348.701369 .76396.722102 .750418 For simplicity, we shall always assume that the vector-valu ed function g:Rn→Rnis deﬁned on all of Rn; otherwise, we must always be careful that the successive it eratesu(k) never leave its domain of deﬁnition, thereby causing the ite ration to break down. To avoid technical complications, we will also assume that gis at least continuous; later results rely on additional smoothness requirements, e.g., continuity o f its ﬁrst and second order partial derivatives. While the solution to a discrete dynamical system is essenti ally trivial, understanding its behavior is deﬁnitely not. Sometimes the solution conve rges to a particular value — the key requirement for numerical solution methods. Someti mes it goes oﬀ to ∞, or, more precisely, the norms†of the iterates are unbounded: /ba∇dblu(k)/ba∇dbl→∞ask→∞. Sometimes the solution repeats itself after a while. And sometimes the iterates behave in a seemingly random, chaotic manner — all depending on the function gand, at times, the initial condition c. Although all of these cases may arise in real-world applica tions, we shall mostly concentrate upon understanding convergence. Deﬁnition 19.1. Aﬁxed point orequilibrium of a discrete dynamical system (19.1) is a vector u⋆∈Rnsuch that g(u⋆) =u⋆. (19.5) We easily see that every ﬁxed point provides a constant solut ion to the discrete dy- namical system, namely u(k)=u⋆for allk. Moreover, it is not hard to prove that any convergent solution necessarily converges to a ﬁxed point. Proposition 19.2. If a solution to a discrete dynamical system converges, lim k→∞u(k)=u⋆, then the limit u⋆is a ﬁxed point. Proof: This is a simple consequence of the continuity of g. We have u⋆= lim k→∞u(k+1)= lim k→∞g(u(k)) =g/parenleftbigg lim k→∞u(k)/parenrightbigg =g(u⋆), the last two equalities following from the continuity of g. Q.E.D. †In view of the equivalence of norms on ﬁnite-dimensional vector spaces , cf. Theorem 3.17, any norm will do here. 12/11/12 1026 c/ci∇cleco√y∇t2012 Peter J. Olver For example, continuing the cosine iteration (19.4), we ﬁnd that the iterates gradually converge to the value u⋆≈.739085, which is the unique solutionto the ﬁxed point equati on cosu=u. Later we will see how to rigorously prove this observed behav ior. Of course, not every solution to a discrete dynamical system will necessarily converge, but Proposition 19.2 says that if it does, then it must conver ge to a ﬁxed point. Thus, a key goal is to understand when a solution converges, and, if s o, to which ﬁxed point — if there is more than one. (In the linear case, only the actual convergence is a signiﬁcant issues since most linear systems admit exactly one ﬁxed poin t, namely u⋆=0.) Fixed points are roughly divided into three classes: •asymptotically stable , with the property that all nearby solutions converge to it, •stable, with the property that all nearby solutions stay nearby, an d •unstable, almost all of whose nearby solutions diverge away from the ﬁ xed point. Thus, from a practical standpoint, convergence of the itera tes of a discrete dynamical system requires asymptoticstability of the ﬁxed point. Exa mples will appear in abundance in the following sections. Scalar Functions As always, the ﬁrst step is to thoroughly understand the scal ar case, and so we begin with a discrete dynamical system u(k+1)=g(u(k)), u(0)=c, (19.6) in which g:R→Ris a continuous, scalar-valued function. As noted above, we will assume, for simplicity, that gis deﬁned everywhere, and so we do not need to worry about whet her the iterates u(0),u(1),u(2),...are all well-deﬁned. The linear case g(u) =auwas treated in Section 10.1 — following equation (10.2). The simplest “nonlinear” case is that of an aﬃne function g(u) =au+b, (19.7) leading to an aﬃne discrete dynamical system u(k+1)=au(k)+b. (19.8) The only ﬁxed point is the solution to u⋆=g(u⋆) =au⋆+b, namely, u⋆=b 1−a. (19.9) The formula for u⋆requires that a/ne}ationslash= 1, and, indeed, the case a= 1 has no ﬁxed point, as the reader can easily conﬁrm; see Exercise . Since we already know the value of u⋆, we can readily analyze the diﬀerences e(k)=u(k)−u⋆, (19.10) 12/11/12 1027 c/ci∇cleco√y∇t2012 Peter J. Olver u⋆ 1u⋆ 2u⋆ 3 Figure 19.1. Fixed Points. between successive iterates and the ﬁxed point. Observe tha t, the smaller e(k)is, the closer u(k)is to the desired ﬁxed point. In many applications, the itera teu(k)is viewed as an approximation to the ﬁxed point u⋆, and so e(k)is interpreted as the errorin thekth iterate. Subtracting the ﬁxed point equation (19.9) from th e iteration equation (19.8), we ﬁnd u(k+1)−u⋆=a(u(k)−u⋆). Therefore the errors e(k)are related by a linear iteration e(k+1)=ae(k),and hence e(k)=ake(0). (19.11) Therefore, as we already demonstrated in Section 10.1, the s olutions to this scalar linear iteration converge: e(k)−→0 and hence u(k)−→u⋆,if and only if |a|<1. This is the criterion for asymptotic stability of the ﬁxed point, or, equivalently, convergence oftheaﬃneiterativesystem(19.8). Themagnitudeof adeterminestherateofconvergence, and the closer it is to 0, the faster the iterates approach the ﬁxed point. Example 19.3. The aﬃne function g(u) =1 4u+2 leads to the iterative scheme u(k+1)=1 4u(k)+2. Starting with the initial condition u(0)= 0, the ensuing values are k1 2 3 4 5 6 7 8 u(k)2.0 2.5 2.625 2.6562 2 .6641 2 .6660 2 .6665 2 .6666 12/11/12 1028 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 19.2. Tangent Line Approximation. Thus, after 8 iterations, the iterates have produced the ﬁxe d pointu⋆=8 3to 4 decimal places. The rate of convergence is1 4, and indeed |e(k)|=|u(k)−u⋆|=/parenleftbig1 4/parenrightbigk|u(0)−u⋆|=8 3/parenleftbig1 4/parenrightbigk−→0 as k−→ ∞. Let us now turn to the fully nonlinear case. First note that th e ﬁxed points of g(u) correspond to the intersections of its graph with the graph o f the function i(u) =u. For instanceFigure19.1showsthegraphofafunctionthathas3ﬁ xedpoints, labeled u⋆ 1,u⋆ 2,u⋆ 3. In general, near any point in its domain, a (smooth) nonlinea r function can be well approximated by its tangent line, which repre4sents the gra ph of an aﬃne function; see Figure 19.2. Therefore, if we are close to a ﬁxed point u⋆, then we might expect the iterative system based on the nonlinear function g(u) to behave very much like that of its aﬃne tangent line approximation. And, indeed, this intuiti on turns out to be essentially correct. This result forms our ﬁrst concrete example of linearization , in which the analysis of a nonlinear system is based on its linear (or, more precise ly, aﬃne) approximation. The explicit formula for the tangent line to g(u) near the ﬁxed point u=u⋆=g(u⋆) is g(u)≈g(u⋆)+g′(u⋆)(u−u⋆)≡au+b, (19.12) where a=g′(u⋆), b =g(u⋆)−g′(u⋆)u⋆=/parenleftbig 1−g′(u⋆)/parenrightbig u⋆. Note that u⋆=b/(1−a) remains a ﬁxed point for the aﬃne approximation: au⋆+b= u⋆. According to the preceding discussion, the convergence of the iterates for the aﬃne approximation is governed by the size of the coeﬃcient a=g′(u⋆). This observation inspires the basic stability criterion for ﬁxed points of sc alar iterative systems. Theorem 19.4. Letg(u)be a continuously diﬀerentiable scalar function. Suppose u⋆=g(u⋆)is a ﬁxed point. If |g′(u⋆)|<1, thenu⋆is an asymptotically stable ﬁxed point, and hence any sequence of iterates u(k)which starts out suﬃciently close to u⋆will converge to u⋆. On the other hand, if |g′(u⋆)|>1, thenu⋆is an unstable ﬁxed point, and 12/11/12 1029 c/ci∇cleco√y∇t2012 Peter J. Olver the only iterates which converge to it are those that land exa ctly on it, i.e., u(k)=u⋆for somek≥0. Proof: The goal is to prove that the errors e(k)=u(k)−u⋆between the iterates and the ﬁxed point tend to 0 as k→∞. To this end, we try to estimate e(k+1)in terms of e(k). According to (19.6) and the Mean Value Theorem C.2 from calc ulus, e(k+1)=u(k+1)−u⋆=g(u(k))−g(u⋆) =g′(v)(u(k)−u⋆) =g′(v)e(k),(19.13) for some vlying between u(k)andu⋆. By continuity, if |g′(u⋆)|<1 at the ﬁxed point, then we can choose δ >0 and|g′(u⋆)|< σ <1 such that the estimate |g′(v)|≤σ <1 whenever |v−u⋆|< δ (19.14) holds in a (perhaps small) interval surrounding the ﬁxed poi nt. Suppose |e(k)|=|u(k)−u⋆|< δ. Then the point vin (19.13), which is closer to u⋆thanu(k), satisﬁes (19.14). Therefore, |u(k+1)−u⋆| ≤σ|u(k)−u⋆|,and hence|e(k+1)| ≤σ|e(k)|.(19.15) In particular, since σ <1, we have|u(k+1)−u⋆|< δ, and hence the subsequent iterate u(k+1)also lies in the interval where (19.14) holds. Repeating the argument, we conclude that, provided the initial iterate satisﬁes |e(0)|=|u(0)−u⋆|< δ, the subsequent errors are bounded by e(k)≤σke(0),and hence e(k)=|u(k)−u⋆| −→0 ask→∞, which completes the proof of the theorem in the stable case. The proof in unstable case is left as Exercise for the reader. Q.E.D. Remark: The constant σgoverns the rate of convergence of the iterates to the ﬁxed point. The closer the iterates are to the ﬁxed point, the smal ler we can choose δin (19.14), and hence the closer we can choose σto|g′(u⋆)|. Thus, roughly speaking, |g′(u⋆)|governs the speed of convergence, once the iterates get close to the ﬁ xed point. This observation will be developed more fully in the following subsection. Remark: The cases when g′(u⋆) =±1 arenotcovered by the theorem. For a linear system, such ﬁxed points are stable, but not asymptotically stable. For nonlinear systems, more detailed knowledge of the nonlinear terms is required i n order to resolve the status — stable or unstable — of the ﬁxed point. Despite their importa nce in certain applications, we will not try to analyze such borderline cases any further h ere. 12/11/12 1030 c/ci∇cleco√y∇t2012 Peter J. Olver m u Figure 19.3. Planetary Orbit. Example 19.5. Given constants ǫ,m, the trigonometric equation u=m+ǫsinu (19.16) isknownas Kepler’s equation . Itarisesinthestudyofplanetarymotion,inwhich0 < ǫ <1 represents the eccentricity of an elliptical planetary orbit, uis theeccentric anomaly , deﬁned as the angle formed at the center of the ellipse by the p lanet and the major axis, andm= 2πt/Tis itsmean anomaly , which is the time, measured in units of T/(2π) whereTis the period of the orbit, i.e., the length of the planet’s ye ar, since perihelion or point of closest approach to the sun; see Figure 19.3. The solutions to Kepler’s equation are the ﬁxed points of the discrete dynamical system based on the function g(u) =m+ǫsinu. Note that |g′(u)|=|ǫcosu|=|ǫ|<1, (19.17) which automaticallyimplies that the as yet unknown ﬁxed poi nt is stable. Indeed, Exercise implies that condition (19.17) is enough to prove the existe nce of a unique stable ﬁxed point. In the particular case m=ǫ=1 2, the result of iterating u(k+1)=1 2+1 2sinu(k) starting with u(0)= 0 is k1 2 3 4 5 6 7 8 9 u(k).5.7397.8370.8713.8826.8862.8873.8877.8878 After 13 iterations, we have converged suﬃciently close to t he solution (ﬁxed point) u⋆= .887862 to have computed its value to 6 decimal places. Remark: In Exercise C.3.57, we 12/11/12 1031 c/ci∇cleco√y∇t2012 Peter J. Olver ug(u) u⋆L+(u) L−(u) Figure 19.4. Graph of a Contraction. outline Bessel’s explicit Fourier series solution to Keple r’s equation that led him to the original deﬁnition of Bessel functions. Inspection of the proof of Theorem 19.4 reveals that we never really used the diﬀer- entiability of g, except to verify the inequality |g(u)−g(u⋆)|≤σ|u−u⋆|for some ﬁxed σ <1. (19 .18) A function that satisﬁes (19.18) for all nearby uis called a contraction at the point u⋆. Anyfunction g(u) whose graph lies between the two lines L±(u) =g(u⋆)±σ(u−u⋆) for some σ <1, for allusuﬃciently close to u⋆, i.e., such that|u−u⋆|< δfor some δ >0, deﬁnes a contraction, and hence ﬁxed point iteration starting with |u(0)−u⋆|< δwill converge to u⋆; see Figure 19.4. In particular, Exercise asks you to prove that any function that is diﬀerentiable at u⋆with|g′(u⋆)|<1 deﬁnes a contraction at u⋆. Example 19.6. The simplest truly nonlinear example is a quadratic polynom ial. The most important case is the so-called logistic map g(u) =λu(1−u), (19.19) whereλ/ne}ationslash= 0 is a ﬁxed non-zero parameter. (The case λ= 0 is completely trivial. Why?) In fact, an elementary change of variables can make any quadr atic iterative system into one involving a logistic map; see Exercise . The ﬁxed points of the logistic map are the solutions to the qu adratic equation u=λu(1−u),or λu2−λu+1 = 0. 12/11/12 1032 c/ci∇cleco√y∇t2012 Peter J. Olver Using the quadratic formula, we conclude that g(u) has two ﬁxed points: u⋆ 1= 0, u⋆ 2= 1−1 λ. Let us apply Theorem 19.4 to determine their stability. The d erivative is g′(u) =λ−2λu, and so g′(u⋆ 1) =λ, g′(u⋆ 2) = 2−λ. Therefore, if|λ|<1, the ﬁrst ﬁxed point is stable, while if 1 < λ <3, the second ﬁxed point is stable. For λ <−1 orλ >3 neither ﬁxed point is stable, and we expect the iterates to not converge at all. Numerical experiments with this example show that it is the s ource of an amazingly diverse range of behavior, depending upon the value of the pa rameterλ. In the accompa- nying Figure 19.5, we display the results of iteration start ing with initial point u(0)=.5 for several diﬀerent values of λ; in each plot, the horizontal axis indicates the iterate num - berkand the vertical axis the iterate valoue u(k)fork= 0,...,100. As expected from Theorem 19.4, the iterates converge to one of the ﬁxed points in the range−1< λ <3, except when λ= 1. For λa little bit larger than λ1= 3, the iterates do not converge to a ﬁxed point. But it does not take long for them to settle down, switching back and forth between two particular values. This behavior indicates the exitence of a (stable) period2 orbitfor the discrete dynamical system, in accordance with the fo llowing deﬁnition. Deﬁnition 19.7. Aperiodkorbitof a discrete dynamical system is a solution that satisﬁesu(n+k)=u(n)for alln= 0,1,2,.... The (minimal)periodis the smallest positive value ofkfor which this condition holds. Thus, a ﬁxed point u(0)=u(1)=u(2)=··· is a period 1 orbit. A period 2 orbit satisﬁes u(0)=u(2)=u(4)=···andu(1)=u(3)=u(5)=···, butu(0)/ne}ationslash=u(1), as otherwise the minimal period would be 1. Similarly, a per iod 3 orbit has u(0)=u(3)=u(6)=···, u(1)=u(4)=u(7)=···, u(2)=u(5)=u(8)=···, withu(0),u(1),u(2)distinct. Stability of a period korbit implies that nearby iterates converge to this periodic solution. For the logistic map, the period 2 orbit persists until λ=λ2≈3.4495, after which the iterates alternate between four values — a period 4 orbit . This again changes at λ=λ3≈3.5441, after which the iterates end up alternating between ei ght values. In fact, there is an increasing sequence of values 3 =λ1< λ2< λ3< λ4<···, where, for any λn< λ≤λn+1, the iterates eventually follow a period 2norbit. Thus, as λ passes through each value λnthe period of the orbit goes from 2nto 2·2n= 2n+1, and the 12/11/12 1033 c/ci∇cleco√y∇t2012 Peter J. Olver 20 40 60 80 1000.20.40.60.81 λ= 1.020 40 60 80 1000.20.40.60.81 λ= 2.020 40 60 80 1000.20.40.60.81 λ= 3.0 20 40 60 80 1000.20.40.60.81 λ= 3.420 40 60 80 1000.20.40.60.81 λ= 3.520 40 60 80 1000.20.40.60.81 λ= 3.55 20 40 60 80 1000.20.40.60.81 λ= 3.620 40 60 80 1000.20.40.60.81 λ= 3.720 40 60 80 1000.20.40.60.81 λ= 3.8 Figure 19.5. Logistic Iterates. discrete dynamical system experiences a bifurcation . The bifurcation values λnare packed closer and closer together as nincreases, piling up on an eventual limiting value λ⋆= lim n→∞λn≈3.5699, at which point the orbit’s period has, so to speak, become inﬁ nitely large. The entire phenomena is known as a period doubling cascade . Interestingly, the ratios of the distances between success ive bifurcation points ap- proaches a well-deﬁned limit, λn+2−λn+1 λn+1−λn−→4.6692... , (19.20) known as Feigenbaum’s constant . In the 1970’s, the American physicist Mitchell Feigen- baum, [65], discovered that similar period doubling cascades appear in a broad range of discrete dynamical systems. Even more remarkably, in almos t all cases, the corresponding ratios of distances between bifurcation points has the samelimiting value. Feigenbaum’s experimental observations were rigorously proved by Oscar Lanford in 1982, [ 125]. 12/11/12 1034 c/ci∇cleco√y∇t2012 Peter J. Olver 2.5 3 3.5 40.20.40.60.81 Figure 19.6. The Logistic Map. Afterλpasses the limiting value λ⋆, all hell breaks loose. The iterates become com- pletely chaotic†, moving at random over the interval [0 ,1]. But this is not the end of the story. Embedded within this chaotic regime are certain s mall ranges of λwhere the system settles down to a stable orbit, whose period is no long er necessarily a power of 2. In fact, there exist values of λfor which the iterates settle down to a stable orbit of periodkforanypositive integer k. For instance, as λincreases past λ3,⋆≈3.83, a period 3 orbit appears over a small range of values, after which, as λincreses slightly further, there is a period doubling cascade where period 6, 12, 24 ,...orbits successively appear, each persisting on a shorter and shorter range of parameter v alues, until λpasses yet an- other critical value where chaos breaks out yet again. There is a well-prescribed order in which the periodic orbits make their successive appearance , and each odd period korbit is followed by a very closely spaced sequence of period doubl ing bifurcations, of periods 2nkforn= 1,2,3,..., after which the iterates revert to completely chaotic beha vior until the next periodic case emerges. The ratios of distances betw een bifurcation points always have the same Feigenbaum limit (19.20). Finally, these peri odic and chaotic windows all pile up on the ultimate parameter value λ⋆ ⋆= 4. And then, when λ >4, all the iterates go oﬀ to∞, and the system ceases to be interesting. The reader is encouraged to write a simple computer program a nd perform some numerical experiments. In particular, Figure 19.6 shows th e asymptotic behavior of the iterates for values of the parameter in the interesting rang e 2< λ <4. The horizontal axis isλ, and the marked points show the ultimate fate of the iteratio n for the given value of λ. For instance, each point the single curve lying above the sm aller values of λrepresents a stable ﬁxed point; this bifurcates into a pair o f curves representing stable †The term “chaotic” does have a precise mathematical deﬁnition, [ 56], but the reader can take it more ﬁguratively for the purposes of this elementary exposition . 12/11/12 1035 c/ci∇cleco√y∇t2012 Peter J. Olver period 2 orbits, which then bifurcates into 4 curves represe nting period 4 orbits, and so on. Chaotic behavior is indicated by a somewhat random patte rn of points lying above the value ofλ. To plot this ﬁgure, we ran the logistic iteration u(n)for 0≤n≤100, discarded the ﬁrst 50 points, and then plotted the next 50 iterates u(51),...,u(100). Investigation of the ﬁne detailed structure of the logistic map requires yet m ore iterations with increased numerical accuracy. In addition one should discard more of t he initial iterates so as to give the system enough time to settle down to a stable periodic orb it or, alternatively, continue in a chaotic manner. Remark: So far, wehaveonlylookedatrealscalariterativesystems . Complexdiscrete dynamical systems display yet more remarkable and fascinat ing behavior. The complex version of the logistic iteration equation leads to the just ly famous Julia and Mandelbrot sets, [128], with their stunning, psychedelic fractal structure, [ 152]. The rich range of phenomena in evidence, even in such extreme ly simple nonlinear iterative systems, is astounding. While intimations ﬁrst a ppeared in the late nineteenth century research of the inﬂuential French mathematician He nri Poincar´ e, serious investi- gations were delayed until the advent of the computer era, wh ich precipitated an explosion of research activity in the area of dynamical systems. Simil ar period doubling cascades and chaos are found in a broad range of nonlinear systems, [ 7], and are often encountered in physical applications, [ 133]. A modern explanation of ﬂuid turbulence is that it is a (very complicated) form of chaos, [ 7]. Quadratic Convergence Let us now return to the more mundane case when the iterates co nverge to a stable ﬁxedpointofthediscretedynamicalsystem. Inapplication s,weusetheiteratestocompute a precise†numerical value for the ﬁxed point, and hence the eﬃciency of the algorithm depends on the speed of convergence of the iterates. According to the remark following the proof Theorem 19.4, th e convergence rate of an iterative system is essentially governed by the magnitud e of the derivative |g′(u⋆)|at the ﬁxed point. The basic inequality (19.15) for the errors e(k)=u(k)−u⋆, namely |e(k+1)|≤σ|e(k)|, is known as a linear convergence estimate . It means that, once the iterates are close to the ﬁxed point, the error decreases by a factor of (at least) σ≈|g′(u⋆)|at each step. If thekthiterateu(k)approximates the ﬁxed point u⋆correctly to mdecimal places, so its error is bounded by |e(k)|< .5×10−m, then the ( k+1)stiterate satisﬁes the error bound |e(k+1)|≤σ|e(k)|< .5×10−mσ=.5×10−m+log10σ. †The degree of precision is to be speciﬁed by the user and the applicat ion. 12/11/12 1036 c/ci∇cleco√y∇t2012 Peter J. Olver More generally, for any j >0, |e(k+j)|≤σj|e(k)|< .5×10−mσj=.5×10−m+jlog10σ, which means that the ( k+j)thiterateu(k+j)has at least‡ m−jlog10σ=m+jlog10σ−1 correct decimal places. For instance, if σ=.1 then each new iterate produces one new decimal place of accuracy (atleast), whileif σ=.9 then it typicallytakes22 ≈−1/log10.9 iterates to produce just one additional accurate digit! Thismeansthatthereisahugeadvantage—particularlyinth eapplicationofiterative methods to the numerical solution of equations — to arrange t hat|g′(u⋆)|be as small as possible. The fastest convergence rate of all will occur whe ng′(u⋆) = 0. In fact, in such a happy situation, the rate of convergence is not just slightl y, but dramatically faster than linear. Theorem 19.8. Suppose that g∈C2, andu⋆=g(u⋆)is a ﬁxed point such that g′(u⋆) = 0. Then, for all iterates u(k)suﬃciently close to u⋆, the errors e(k)=u(k)−u⋆ satisfy the quadratic convergence estimate |e(k+1)| ≤τ|e(k)|2(19.21) for some constant τ >0. Proof: Just as that of the linear convergence estimate (19.15), th e proof relies on approximating g(u) by a simpler function near the ﬁxed point. For linear conver gence, an aﬃne approximation suﬃced, but here we require a higher orde r approximation. Thus, we replace the mean value formula (19.13) by the ﬁrst order Tayl or expansion g(u) =g(u⋆)+g′(u⋆)(u−u⋆)+1 2g′′(w)(u−u⋆)2, (19.22) where the ﬁnal error term depends on an (unknown) point wthat lies between uandu⋆; see (C.10) for details. At a ﬁxed point, the constant term is g(u⋆) =u⋆. Furthermore, under our hypothesis g′(u⋆) = 0, and so (19.22) reduces to g(u)−u⋆=1 2g′′(w)(u−u⋆)2. Therefore, |g(u)−u⋆|≤τ|u−u⋆|2, (19.23) whereτis chosen so that 1 2|g′′(w)|≤τ (19.24) for allwsuﬃciently close to u⋆. Therefore, the magnitude of τis governed by the size of thesecond derivative of the iterative function g(u) near the ﬁxed point. We use the inequality (19.23) to estimate the error |e(k+1)|=|u(k+1)−u⋆|=|g(u(k))−g(u⋆)| ≤τ|u(k)−u⋆|2=τ|e(k)|2, which establishes the quadratic convergence estimate (19. 21). Q.E.D. ‡Note that since σ <1, the logarithm log10σ−1=−log10σ >0 is positive. 12/11/12 1037 c/ci∇cleco√y∇t2012 Peter J. Olver Let us see how the quadratic estimate (19.21) speeds up the co nvergence rate. Fol- lowing our earlier argument, suppose u(k)is correct to mdecimal places, so |e(k)|< .5×10−m. Then (19.21) implies that |e(k+1)|< .5×(10−m)2τ=.5×10−2m+log10τ, and sou(k+1)has 2m−log10τaccurate decimal places. If τ≈|g′′(u⋆)|is of moderate size, we have essentially doubledthe number of accurate decimal places in just a single iterate! A second iteration will double the number of accura te digits yet again. Thus, the convergence of a quadratic iteration scheme is extremely rapid, and, barring round-oﬀ errors, one can produce any desired number of digits of accur acy in a very short time. For example, if we start with an initial guess that is accurate in the ﬁrst decimal digit, then a linear iteration with σ=.1 will require 49 iterations to obtain 50 decimal place accur acy, whereas a quadratic iteration (with τ= 1) will only require 6 iterations to obtain 26= 64 decimal places of accuracy! Example 19.9. Consider the function g(u) =2u3+3 3u2+3. There is a unique (real) ﬁxed point u⋆=g(u⋆), which is the real solution to the cubic equation 1 3u3+u−1 = 0. Note that g′(u) =2u4+6u2−6u 3(u2+1)2=6u/parenleftbig1 3u3+u−1/parenrightbig 3(u2+1)2, and hence g′(u⋆) = 0 vanishes at the ﬁxed point. Theorem 19.8 implies that the iterations will exhibit quadratic convergence to the root. Indeed, we ﬁ nd, starting with u(0)= 0, the following values: k 1 2 3 u(k)1.00000000000000 .833333333333333 .817850637522769 4 5 6 .817731680821982 .817731673886824 .817731673886824 The convergence rate is dramatic: after only 5 iterations, w e have produced the ﬁrst 15 decimal places of the ﬁxed point. In contrast, the linearly c onvergent scheme based on /tildewideg(u) = 1−1 3u3takes 29 iterations just to produce the ﬁrst 5 decimal places of the same solution. 12/11/12 1038 c/ci∇cleco√y∇t2012 Peter J. Olver In practice, the appearance of a quadratically convergent ﬁ xed point is a matter of luck. The construction of quadratically convergent iterat ive methods for solving equations will be the focus of the following Section 19.2. Vector–Valued Iteration Extending the preceding analysis to vector-valued iterati ve systems is not especially diﬃcult. We will build on our experience with linear iterati ve systems, and so the reader is advised to review the basic concepts and results from Chap ter 10 before proceeding to the nonlinear systems presented here. We begin by ﬁxing a norm /ba∇dbl·/ba∇dblonRn. Since we will also be computing the asso- ciated matrix norm /ba∇dblA/ba∇dbl, as deﬁned in Theorem 10.20, it may be more convenient for computations to adopt either the 1 or the ∞norms rather than the standard Euclidean norm. We begin by deﬁning the vector-valued counterpart of the bas ic linear convergence condition (19.18). Deﬁnition 19.10. A function g:Rn→Rnis acontraction at a point u⋆∈Rnif there exists a constant 0 ≤σ <1 such that /ba∇dblg(u)−g(u⋆)/ba∇dbl ≤σ/ba∇dblu−u⋆/ba∇dbl (19.25) for allusuﬃciently close to u⋆, i.e.,/ba∇dblu−u⋆/ba∇dbl< δfor some ﬁxed δ >0. Remark: The notion of a contraction depends on the underlying choic e of matrix norm. Indeed, the linear function g(u) =Auif and only if/ba∇dblA/ba∇dbl<1, which implies thatAis a convergent matrix. While every convergent matrix satis ﬁes/ba∇dblA/ba∇dbl<1 insome matrix norm, and hence deﬁnes a contraction relative to that norm, it may very well have /ba∇dblA/ba∇dbl>1inaparticularnorm, violatingthecontactioncondition; see(10.42)foranexplicit example. Theorem 19.11. Ifu⋆=g(u⋆)is a ﬁxed point for the discrete dynamical system (19.1)andgis a contraction at u⋆, thenu⋆is an asymptotically stable ﬁxed point. Proof: The proof is a copy of the last part of the proof of Theorem 19. 4. We write /ba∇dblu(k+1)−u⋆/ba∇dbl=/ba∇dblg(u(k))−g(u⋆)/ba∇dbl≤σ/ba∇dblu(k)−u⋆/ba∇dbl, using the assumed estimate (19.25). Iterating this basic in equality immediately demon- strates that /ba∇dblu(k)−u⋆/ba∇dbl ≤σk/ba∇dblu(0)−u⋆/ba∇dblfor k= 0,1,2,3,... . (19.26) Sinceσ <1, the right hand side tends to 0 as k→∞, and hence u(k)→u⋆.Q.E.D. In most interesting situations, the function gis diﬀerentiable, and so can be approxi- mated by its ﬁrst order Taylor polynomial (C.14): g(u)≈g(u⋆)+g′(u⋆)(u−u⋆) =u⋆+g′(u⋆)(u−u⋆). (19.27) 12/11/12 1039 c/ci∇cleco√y∇t2012 Peter J. Olver Here g′(u) = ∂g1 ∂u1∂g1 ∂u2...∂g1 ∂un ∂g2 ∂u1∂g2 ∂u2...∂g2 ∂un ............ ∂gn ∂u1∂gn ∂u2...∂gn ∂un , (19.28) denotes the n×nJacobian matrix of the vector-valued function g, whose entries are the partial derivatives of its individual components. Since u⋆is ﬁxed, the the right hand side of (19.27) is an aﬃne function of u. Moreover, u⋆remains a ﬁxed point of the aﬃne approximation. Proposition 10.44 tells us that iteration o f the aﬃne function will converge tothe ﬁxedpoint ifandonlyifitscoeﬃcient matrix, namely g′(u⋆), isaconvergent matrix, meaning that its spectral radius ρ(g′(u⋆))<1. This observation motivates the following theorem and corollary. Theorem 19.12. Letu⋆be a ﬁxed point for the discrete dynamical system u(k+1)= g(u(k)). If the Jacobian matrix norm /ba∇dblg′(u⋆)/ba∇dbl<1, thengis a contraction at u⋆, and hence the ﬁxed point u⋆is asymptotically stable. Proof: The ﬁrst order Taylor expansion (C.14) of g(u) at the ﬁxed point u⋆takes the form g(u) =g(u⋆)+g′(u⋆)(u−u⋆)+R(u−u⋆), (19.29) where the remainder term satisﬁes lim u→u⋆R(u−u⋆) /ba∇dblu−u⋆/ba∇dbl= 0. Letε >0 be such that σ=/ba∇dblg′(u⋆)/ba∇dbl+ε <1. Choose 0 < δ <1 such that/ba∇dblR(u−u⋆)/ba∇dbl≤ε/ba∇dblu−u⋆/ba∇dblwhenever/ba∇dblu−u⋆/ba∇dbl≤δ. For suchu, we have, by the Triangle Inequality, /ba∇dblg(u)−g(u⋆)/ba∇dbl ≤ /ba∇dblg′(u⋆)(u−u⋆)/ba∇dbl+/ba∇dblR(u−u⋆)/ba∇dbl ≤/parenleftbig /ba∇dblg′(u⋆)/ba∇dbl+ε/parenrightbig /ba∇dblu−u⋆/ba∇dbl=σ/ba∇dblu−u⋆/ba∇dbl, which establishes the contraction inequality (19.25). Q.E.D. Corollary 19.13. If the Jacobian matrix g′(u⋆)is a convergent matrix, meaning that its spectral radius satisﬁes ρ/parenleftbig g′(u⋆)/parenrightbig <1, thenu⋆is an asymptotically stable ﬁxed point. Proof: Corollary 10.32 assures us that /ba∇dblg′(u⋆)/ba∇dbl<1 in some matrix norm. Using this norm, the result immediately follows from the theorem. Q.E.D. 12/11/12 1040 c/ci∇cleco√y∇t2012 Peter J. Olver Theorem 19.12 tells us that initial values u(0)that are suﬃciently near a stable ﬁxed pointu⋆are guaranteed to converge to it. In the linear case, closene ss of the initial data to the ﬁxed point was not, in fact, an issue; all stable ﬁxed po ints are, in fact, globally stable. For nonlinear iteration, it is of critical importan ce, and one does not typically expect iteration starting with far away initial data to conv erge to the desired ﬁxed point. An interesting (and diﬃcult) problem is to determine the so- calledbasin of attraction of a stable ﬁxed point, deﬁned as the set of all initial data that e nds up converging to it. As in the elementary logistic map (19.19), initial values that li e outside a basin of attraction can lead to divergent iterates, periodic orbits, or even exhibi t chaotic behavior. The full range of possible phenomena is a topic of contemporary research in dynamical systems theory, [152], and in numerical analysis, [ 7]. Example 19.14. Consider the function g(u,v) =/parenleftbigg −1 4u3+9 8u+1 4v3 3 4v−1 2uv/parenrightbigg . There are four (real) ﬁxed points; stability is determined b y the size of the eigenvalues of the Jacobian matrix g′(u,v) =/parenleftigg 9 8−3 4u2−1 2v 3 4v23 4−1 2u/parenrightigg at each of the ﬁxed points. The results are summarized in the f ollowing table: ﬁxed point u⋆ 1=/parenleftbigg 0 0/parenrightbigg u⋆ 2=/parenleftigg 1√ 2 0/parenrightigg u⋆ 3=/parenleftigg −1√ 2 0/parenrightigg u⋆ 4=/parenleftigg −1 2 1 2/parenrightigg Jacobian matrix/parenleftigg 9 80 03 4/parenrightigg /parenleftigg 3 40 03 4−1 2√ 2/parenrightigg /parenleftigg 3 40 03 4+1 2√ 2/parenrightigg /parenleftigg 15 16−1 4 3 161/parenrightigg eigenvalues 1.125, .75.75, .396447 1 .10355, .75.96875±.214239i spectral radius 1.125 .75 1 .10355 .992157 Thus,u⋆ 2andu⋆ 4are stable ﬁxed points, whereas u⋆ 1andu⋆ 3are both unstable. Indeed, starting with u(0)= (.5,.5)T, it takes 24 iterates to converge to u⋆ 2with 4 signiﬁcant decimal digits, whereas starting with u(0)= (−.7,.7)T, it takes 1049 iterates to converge to within 4 digits of u⋆ 4; the slower convergence rate is predicted by the larger Jaco bian spectral radius. The two basins of attraction are plotted in Figure 19.7. The stable ﬁxed points are indicated by black dots. The light gray region con tainsu⋆ 2and indicates all the points that converge to it; the darker gray indicates points converging, more slowly, to u⋆ 4. All other initial points, except u⋆ 1andu⋆ 3, have rapidly unbounded iterates: /ba∇dblu(k)/ba∇dbl→∞. The smaller the spectral radius or matrix norm of the Jacobia n matrix at the ﬁxed point, the faster the nearby iterates will converge to it. As in the scalar case, quadratic 12/11/12 1041 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 19.7. Basins of Attraction. convergence will occur when the Jacobian matrix g′(u⋆) = O is the zero matrix†, i.e., allﬁrst order partial derivatives of the components of gvanish at the ﬁxed point. The quadratic convergence estimate /ba∇dblu(k+1)−u⋆/ba∇dbl ≤τ/ba∇dblu(k)−u⋆/ba∇dbl2(19.30) is a consequence of the second order Taylor expansion at the ﬁ xed point. Details of the proof are left as an exercise. Of course, in practice we don’t know the norm or spectral radi us of the Jacobian matrixg′(u⋆) because we don’t know where the ﬁxed point is. This apparent diﬃculty can be easily circumvented by requiring that /ba∇dblg′(u)/ba∇dbl<1 for all u— or, at least, for all uin a domain Ω containing the ﬁxed point. In fact, this hypothe sis can be used to prove the exitence and uniqueness of asymptotically stable ﬁxed p oints. Rather than work with the Jacobian matrix, let us return to the contraction condit ion (19.25), but now imposed uniformly on an entire domain. Deﬁnition 19.15. A function g:Rn→Rnis called a contraction mapping on a domain Ω⊂Rnif (a) it maps Ω to itself, so g(u)∈Ω whenever u∈Ω, and (b) there exists a constant 0≤σ <1 such that /ba∇dblg(u)−g(v)/ba∇dbl ≤σ/ba∇dblu−v/ba∇dblfor all u,v∈Ω. (19.31) In other words, applying a contraction mapping reduces the m utual distance between points. In view of Theorem 19.12, we can detect contraction m appings by lloking at the norm of their Jacobian matrix. Lemma 19.16. Ifg:Ω→Ωand/ba∇dblg′(u)/ba∇dbl<1for allu∈Ω, thengis a contraction mapping. †Having zero spectral radius is not suﬃcient for quadratic convergenc e; see Exercise . 12/11/12 1042 c/ci∇cleco√y∇t2012 Peter J. Olver Ωg Ω Figure 19.8. A Contraction Mapping. So, as its name indicates, a contraction mapping eﬀectively shrinks the size of its domain; see Figure 19.8. As the iterations proceed, the succ essive image domains become smaller and smaller. If the original domain is closed and bou nded, then it is forced to shrink down to a single point, which is the unique ﬁxed point o f the iterative system, leading to the Contraction Mapping Theorem . Theorem 19.17. Ifgisacontractionmappingonaclosedbounded domain Ω⊂Rn, thengadmits a unique ﬁxed point u⋆∈Ω. Moreover, starting with any initial point u(0)∈Ω, the iterates u(k+1)=g(u(k))necessarily converge to the ﬁxed point: u(k)→u⋆. The proof is outlined in Exercise . 19.2. Solution of Equations and Systems. Solvingnonlinearequationsandsystemsofequationsis, of course, aproblemofutmost importanceinmathematicsanditsmanifoldapplications. I tcanalsobeextremelydiﬃcult. Indeed, ﬁnding a complete set of (numerical) solutions to a c omplicated nonlinear system can be an almost insurmountable challlenge. In its most gene ral version, we are given a collection of mfunctions f1,...,fmdepending upon nvariables u1,...,un, and are asked to determine all possible solutions u= (u1,u2,...,un)Tto the system f1(u1,...,un) = 0, ... fm(u1,...,un) = 0. (19.32) In many applications, the number of equations equals the num ber of unknowns, m=n, in which case one expects both existence and uniqueness of so lutions. This point will be discussed in further detail below. Here are some prototypical examples: (a) Find the roots of the quintic polynomial equation u5+u+1 = 0. (19.33) Graphing the left hand side of the equation, as in Figure 19.9 , convinces us that there is just one real root, lying somewhere between −1 and−.5. While there are explicit algebraic formulas for the roots of quadratic, cubic, and quartic poly nomials, a famous theorem†due †A modern proof of this fact relies on Galois theory, [ 75]. 12/11/12 1043 c/ci∇cleco√y∇t2012 Peter J. Olver -1 -0.5 0.5 1 -4-3-2-1123 Figure 19.9. Graph of u5+u+1. to the Norwegian mathematician Nils Henrik Abel in the early 1800’s states that there is nosuch formula for generic ﬁfth order polynomial equations. (b) Any ﬁxed point equation u=g(u) has the form (19.35) where f(u) =u−g(u). For example, the trigonometric Kepler equation u−ǫsinu=m arises in the study of planetary motion, cf. Example 19.5. He reǫ,mare ﬁxed constants, and we seek a corresponding solution u. (c) Suppose we are given chemical compounds A,B,Cthat react to produce a fourth compound Daccording to 2A+B←→D, A +3C←→D. Leta,b,cbe the initial concentrations of the reagents A,B,Cinjected into the reaction chamber. If udenotes the concentration of Dproduced by the ﬁrst reaction, and vthat by the second reaction, then the ﬁnal equilibrium concentra tions a⋆=a−2u−v, b⋆=b−u, c⋆=c−3v, d⋆=u+v, of the reagents will be determined by solving the nonlinear s ystem (a−2u−v)2(b−u) =α(u+v),(a−2u−v)(c−3v)3=β(u+v),(19.34) whereα,βare the known equilibrium constants of the two reactions. Our immediate goal is to develop numerical algorithms for so lving such nonlinear equations. Unfortunately, there is no direct universal sol ution method for nonlinear sys- tems comparable to Gaussian elimination. As a result, numer ical solution techniques rely almost exclusively on iterative algorithms. This section p resents the principal methods for numerically approximating the solution(s) to a nonlinear s ystem. We shall only discuss general purpose algorithms; specialized methods for solvi ng particular classes of equations, e.g., polynomial equations, can be found in numerical analy sis texts, e.g., [ 27,35,153]. Of course, the most important specialized methods — those de signed for solving linear systems — will continue to play a critical role, even in the no nlinear regime. 12/11/12 1044 c/ci∇cleco√y∇t2012 Peter J. Olver a b u⋆f(u) Figure 19.10. Intermediate Value Theorem. The Bisection Method We begin, as always, with the scalar case. Thus, we are given a real-valued function f:R→R, and seek its roots, i.e., the real†solution(s) to the scalar equation f(u) = 0. (19.35) Our immediate goal is to develop numerical algorithms for so lving such nonlinear scalar equations. The most primitive algorithm, and the only one that is guaranteed to work in all cases, is the Bisection Method. While it has an iterative ﬂavor, it cannot be properly classed as a method governed by functional iteration as deﬁn ed in the preceding section, and so must be studied directly in its own right. The starting point is the Intermediate Value Theorem, which we state in simpliﬁed form. See Figure 19.10 for an illustration, and [ 9] for a proof. Lemma 19.18. Letf(u)be a continuous scalar function. Suppose we can ﬁnd two pointsa < bwhere the values of f(a)andf(b)takeopposite signs, so either f(a)<0and f(b)>0, orf(a)>0andf(b)<0. Then there exists at least one point a < u⋆< bwhere f(u⋆) = 0. Note that if f(a) = 0 or f(b) = 0, then ﬁnding a root is trivial. If f(a) andf(b) have the same sign, then there may or may not be a root in between. Fi gure 19.11 plots the functions u2+ 1,u2andu2−1, on the interval −2≤u≤2. The ﬁrst has two simple roots; the second has a single double root, while the third ha s no root. We also note that continuity of the function on the entire interval [ a,b] is an essential hypothesis. For example, the function f(u) = 1/usatisﬁesf(−1) =−1 andf(1) = 1, but there is no root to the equation 1 /u= 0. Note carefully that the Lemma 19.18 does notsay there is a unique root between a andb. There may be many roots, or even, in pathological examples, inﬁnitely many. All the theorem guarantees is that, under the stated hypotheses , there is at least one root. †Complex roots to complex equations will be discussed later. 12/11/12 1045 c/ci∇cleco√y∇t2012 Peter J. Olver -2 -1 1 2 -2-112345 -2 -1 1 2 -2-112345 -2 -1 1 2 -2-112345 Figure 19.11. Roots of Quadratic Functions. Once we are assured that a root exists, bisection relies on a “ divide and conquer” strategy. The goal is to locate a root a < u⋆< bbetween the endpoints. Lacking any additionalevidence, onetacticwouldbetotrythemidpoint c=1 2(a+b) asaﬁrst guess for the root. If, by some miracle, f(c) = 0, then we are done, since we have found a solution! Otherwise (and typically) we look at the sign of f(c). There are two possibilities. If f(a) andf(c) are of opposite signs, then the Intermediate Value Theorem tells us that there is a rootu⋆lying between a < u⋆< c. Otherwise, f(c) andf(b) must have opposite signs, and so there is a root c < u⋆< b. In either event, we apply the same method to the interval in which we are assured a root lies, and repeat the procedure. Each iteration halves the length of the interval, and chooses the half in which a root is sure to lie. (There may, of course, be a root in the other half interval, but as we cannot b e sure, we discard it from further consideration.) The root we home in on lies trapped i n intervals of smaller and smaller width, and so convergence of the method is guarantee d. Figure bisect illustrates the steps in a particular example. Example 19.19. The roots of the quadratic equation f(u) =u2+u−3 = 0 can be computed exactly by the quadratic formula: u⋆ 1=−1+√ 13 2≈1.302775... , u⋆ 2=−1−√ 13 2≈−2.302775... . Let usseehow onemight approximatethem byapplyingtheBise ctionAlgorithm. Westart the procedure by choosing the points a=u(0)= 1, b=v(0)= 2, noting that f(1) =−1 andf(2) = 3 have opposite signs and hence we are guaranteed that th ere is at least one root between 1 and 2. In the ﬁrst step we look at the midpoint of the interval [1 ,2], which is 1 .5, and evaluate f(1.5) =.75. Since f(1) =−1 andf(1.5) =.75 have opposite signs, we know that there is a root lying between 1 and 1 .5. Thus, we take u(1)= 1 and v(1)= 1.5 as the endpoints of the next interval, and continue. The nex t midpoint is at 1.25, where f(1.25) =−.1875 has the opposite sign to f(1.5) =.75, and so a root lies between u(2)= 1.25 andv(2)= 1.5. The process is then iterated as long as desired — or, more practically, as long as your computer’s precision does not become an issue. The table displays the result of the algorithm, rounded oﬀ to four decimal places. After 14 iterations, the Bisection Method has correctly com puted the ﬁrst four decimal 12/11/12 1046 c/ci∇cleco√y∇t2012 Peter J. Olver k u(k)v(k)w(k)=1 2(u(k)+v(k))f(w(k)) 0 1 2 1 .5 .75 1 1 1 .5 1 .25−.1875 2 1.25 1 .5 1 .375 .2656 3 1.25 1 .375 1 .3125 .0352 4 1.25 1 .3125 1 .2813−.0771 5 1.2813 1 .3125 1 .2969−.0212 6 1.2969 1 .3125 1 .3047 .0069 7 1.2969 1 .3047 1 .3008−.0072 8 1.3008 1 .3047 1 .3027−.0002 9 1.3027 1 .3047 1 .3037 .0034 10 1.3027 1 .3037 1 .3032 .0016 11 1.3027 1 .3032 1 .3030 .0007 12 1.3027 1 .3030 1 .3029 .0003 13 1.3027 1 .3029 1 .3028 .0001 14 1.3027 1 .3028 1 .3028−.0000 digits of the positive root u⋆ 1. A similar bisection starting with the interval from u(1)=−3 tov(1)=−2 will produce the negative root. A formal implementation of the Bisection Algorithm appears in the accompanying pseudocode program. The endpoints of the kthinterval are denoted by u(k)andv(k). The midpoint is w(k)=1 2/parenleftbig u(k)+v(k)/parenrightbig , and the key decision is whether w(k)should be the right or left hand endpoint of the next interval. The integer n, governing the number of iterations, is to be prescribed in accordance with how accur ately we wish to approximate the root u⋆. The algorithm produces two sequences of approximations u(k)andv(k)that both converge monotonically to u⋆, one from below and the other from above: a=u(0)≤u(1)≤u(2)≤ ··· ≤ u(k)−→u⋆←−v(k)≤ ··· ≤ v(2)≤v(1)≤v(0)=b. andu⋆is trapped between the two. Thus, the root is trapped inside a sequence of intervals [u(k),v(k)] of progressively shorter and shorter length. Indeed, the l ength of each interval is exactly half that of its predecessor: v(k)−u(k)=1 2(v(k−1)−u(k−1)). Iterating this formula, we conclude that v(n)−u(n)=/parenleftbig1 2/parenrightbign(v(0)−u(0)) =/parenleftbig1 2/parenrightbign(b−a)−→0 as n−→ ∞. The midpoint w(n)=1 2(u(n)+v(n)) 12/11/12 1047 c/ci∇cleco√y∇t2012 Peter J. Olver The Bisection Method start iff(a)f(b)<0setu(0)=a,v(0)=b else print “Bisection Method not applicable” fork= 0ton−1 setw(k)=1 2(u(k)+v(k)) iff(w(k)) = 0, stop; print u⋆=w(k) iff(u(k))f(w(k))<0, setu(k+1)=u(k),v(k+1)=w(k) else set u(k+1)=w(k),v(k+1)=v(k) nextk printu⋆=w(n)=1 2(u(n)+v(n)) end lies within a distance |w(n)−u⋆|≤1 2(v(n)−u(n)) =/parenleftbig1 2/parenrightbign+1(b−a) oftheroot. Consequently,ifwedesiretoapproximatethero otwithinaprescribedtolerance ε, we should choose the number of iterations nso that /parenleftbig1 2/parenrightbign+1(b−a)< ε, orn >log2b−a ε−1. (19.36) Summarizing: Theorem 19.20. Iff(u)is a continuous function, with f(a)f(b)<0, then the Bisection Method starting with u(0)=a, v(0)=b, will converge to a solution u⋆to the equation f(u) = 0lyingbetween aandb. Afternsteps, themidpoint w(n)=1 2(u(n)+v(n)) will be within a distance of ε= 2−n−1(b−a)from the solution. For example, in the case of the quadratic equation in Example 19.19, after 14 itera- tions, we have approximated the positive root to within ε=/parenleftbig1 2/parenrightbig15(2−1)≈3.052×10−5, reconﬁrmingourobservationthatwehaveaccuratelycomput editsﬁrstfourdecimalplaces. If we are in need of 10 decimal places, we set our tolerance to ε=.5×10−10, and so, according to (19.36), must perform n= 34>33.22≈log22×1010−1 successive bisections†. †This assumes we have suﬃcient precision on the computer to avoid rou nd-oﬀ errors. 12/11/12 1048 c/ci∇cleco√y∇t2012 Peter J. Olver Example 19.21. As noted at the beginning of this section, the quintic equati on f(u) =u5+u+1 = 0 hasonerealroot,whosevaluecanbereadilycomputedbybise ction. Westartthealgorithm withtheinitialpoints u(0)=−1, v(0)= 0,notingthat f(−1) =−1<0whilef(0) = 1>0 are of opposite signs. In order to compute the root to 6 decima l places, we set ε=.5×10−6 in (19.36), and so need to perform n= 20>19.93≈log22×106−1 bisections. Indeed, the algorithm produces the approximation u⋆≈−.754878 to the root, and the displayed digits are guaranteed to be accurate. Fixed Point Methods The Bisection Method has an ironclad guarantee to converge t o a root of the function — provided it can be properly started by locating two points w here the function takes opposite signs. This may be tricky if the function has two ver y closely spaced roots and is, say, negative only for a very small interval between them , and may be impossible for multiple roots, e.g., the root u⋆= 0 of the quadratic function f(u) =u2. When applicable, its convergence rate is completely predictabl e, but not especially fast. Worse, it has no immediately apparent extension to systems of equat ions, since there is noobvious counterpart to the Intermediate Value Theorem for vector-v alued functions. Most other numerical schemes for solving equations rely on s ome form of ﬁxed point iteration. Thus, we seek to replace the system of equations f(u) =0with a ﬁxed point systemu=g(u), that leads to the iterative solution scheme u(k+1)=g(u(k)). For this to work, there are two key requirements: (a) The solution u⋆to the equation f(u) =0is also a ﬁxed point for g(u), and (b)u⋆is, in fact a stable ﬁxed point, meaning that the Jacobian g′(u⋆) is a convergent matrix, or, slightly more restrictively, /ba∇dblg′(u⋆)/ba∇dbl<1 for a prescribed matrix norm. If both conditions hold, then, provided we choose the initial iterate u(0)=csuﬃciently close to u⋆, the iterates u(k)→u⋆will converge to the desired solution. Thus, the key to the practical use of functional iteration for solving equ ations is the proper design of an iterative system — coupled with a reasonably good initial gu ess for the solution. Before implementing general procedures, let us discuss a na¨ ıve ex ample. Example 19.22. To solve the cubic equation f(u) =u3−u−1 = 0 (19 .37) we note that f(1) =−1 whilef(2) = 5, and so there is a root between 1 and 2. Indeed, the Bisection Method leads to the approximate value u⋆≈1.3247 after 17 iterations. Let us try to ﬁnd the same root by ﬁxed point iteration. As a ﬁrs t, na¨ ıve, guess, we rewrite the cubic equation in ﬁxed point form u=u3−1 =/tildewideg(u). Starting with the initial guess u(0)= 1.5, successive approximations to the solution are found by iterating u(k+1)=/tildewideg(u(k)) = (u(k))3−1, k = 0,1,2,... . 12/11/12 1049 c/ci∇cleco√y∇t2012 Peter J. Olver However, their values u(0)= 1.5, u(1)= 2.375, u(2)= 12.396, u(3)= 1904, u(4)= 6.9024×109, u(5)= 3.2886×1029, ... rapidlybecomeunbounded, andsofailtoconverge. Thiscoul d, infact, havebeenpredicted by the convergence criterion in Theorem 19.4. Indeed, /tildewideg′(u) =−3u2and so|/tildewideg′(u)|>3 for allu≥1, including the root u⋆. This means that u⋆is an unstable ﬁxed point, and the iterates cannot converge to it. On the other hand, we can rewrite the equation (19.37) in the a lternative iterative form u=3√ 1+u=g(u). In this case 0≤g′(u) =1 3(1+u)2/3≤1 3foru >0. Thus, the stability condition (19.14) is satisﬁed, and we an ticipate convergence at a rate of at least1 3. (The Bisection Method converges more slowly, at rate1 2.) Indeed, the ﬁrst few iterates u(k+1)=3√ 1+u(k)are 1.5,1.35721,1.33086,1.32588,1.32494,1.32476,1.32473, and we have converged to the root, correct to four decimal pla ces, in only 6 iterations. Newton’s Method Our immediate goal is to design an eﬃcient iterative scheme u(k+1)=g(u(k)) whose iterates converge rapidly to the solution of the given scala r equation f(u) = 0. As we learned in Section 19.1, the convergence of the iteration is governed by the magnitude of its derivative at the ﬁxed point. At the very least, we shou ld impose the stability criterion|g′(u⋆)|<1, and the smaller this quantity can be made, the faster the it erative scheme converges. if we are able to arrange that g′(u⋆) = 0, then the iterates will converge quadraticallyfast, leading, asnotedinthediscussionfol lowingTheorem19.8, toadramatic improvement in speed and eﬃciency. Now, the ﬁrst condition requires that g(u) =uwhenever f(u) = 0. A little thought will convince you that the iterative function should take th e form g(u) =u−h(u)f(u), (19.38) whereh(u) is a reasonably nice function. If f(u⋆) = 0, then clearly u⋆=g(u⋆), and so u⋆ is a ﬁxed point. The converse holds provided h(u)/ne}ationslash= 0 is never zero. For quadratic convergence, the key requirement is that the d erivative of g(u) be zero at the ﬁxed point solutions. We compute g′(u) = 1−h′(u)f(u)−h(u)f′(u). Thus,g′(u⋆) = 0 at a solution to f(u⋆) = 0 if and only if 0 = 1−h′(u⋆)f(u⋆)−h(u⋆)f′(u⋆) = 1−h(u⋆)f′(u⋆). 12/11/12 1050 c/ci∇cleco√y∇t2012 Peter J. Olver Consequently, we should require that h(u⋆) =1 f′(u⋆)(19.39) to ensure a quadratically convergent iterative scheme. Thi s assumes that f′(u⋆)/ne}ationslash= 0, which means that u⋆is asimple root off. For here on, we leave aside multiple roots, which require a diﬀerent approach, to be outlined in Exercis e. Of course, there are many functions h(u) that satisfy (19.39), since we only need to specify its value at a single point. The problem is that we do n ot know u⋆— after all this is what we are trying to compute — and so cannot compute the val ue of the derivative offthere. However, we can circumvent this apparent diﬃculty by a simple device: we impose equation (19.39) at all points, setting h(u) =1 f′(u), (19.40) which certainly guarantees that it holds at the solution u⋆. The result is the function g(u) =u−f(u) f′(u), (19.41) and the resulting iteration scheme is known as Newton’s Method , which, as the name suggests, dates back to the founder of the calculus. To this d ay, Newton’s Method remains themost important general purpose algorithm for solving equat ions. It starts with an initial guess u(0)to be supplied by the user, and then successively computes u(k+1)=u(k)−f(u(k)) f′(u(k)). (19.42) As long as the initial guess is suﬃciently close, the iterate su(k)are guaranteed to converge, quadratically fast, to the (simple) root u⋆of the equation f(u) = 0. Theorem 19.23. Suppose f(u)∈C2is twice continuously diﬀerentiable. Let u⋆be a solution to the equation f(u⋆) = 0such that f′(u⋆)/ne}ationslash= 0. Given an initial guess u(0) suﬃciently close to u⋆, the Newton iteration scheme (19.42)converges at a quadratic rate to the solution u⋆. Proof: Bycontinuity, if f′(u⋆)/ne}ationslash= 0, then f′(u)/ne}ationslash= 0 forall usuﬃciently closeto u⋆, and hence the Newton iterative function (19.41) is well deﬁned a nd continuously diﬀerentiable nearu⋆. Sinceg′(u) =f(u)f′′(u)/f′(u)2, we have g′(u⋆) = 0 when f(u⋆) = 0, as promised by our construction. The quadratic convergence of the resul ting iterative scheme is an immediate consequence of Theorem 19.8. Q.E.D. Example 19.24. Consider the cubic equation f(u) =u3−u−1 = 0, that we already solved in Example 19.22. The function used in the Newton iteration is g(u) =u−f(u) f′(u)=u−u3−u−1 3u2−1, 12/11/12 1051 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 1 -0.25-0.2-0.15-0.1-0.050.050.1 Figure 19.12. The function f(u) =u3−3 2u2+5 9u−1 27. which is well-deﬁned as long as u/ne}ationslash=±1√ 3. We will try to avoid these singular points. The iterative procedure u(k+1)=g(u(k)) =u(k)−(u(k))3−u(k)−1 3(u(k))2−1 with initial guess u(0)= 1.5 produces the following values: 1.5,1.34783,1.32520,1.32472, and we have computed the root to 5 decimal places after only th ree iterations. The quadratic convergence of Newton’s Method implies that, rou ghly, each new iterate doubles the number of correct decimal places. Thus, to compute the ro ot accurately to 40 decimal placeswouldonlyrequire3further iterations†. Thisunderscores thetremendous advantage that the Newton algorithm oﬀers over competing methods. Example 19.25. Consider the cubic polynomial equation f(u) =u3−3 2u2+5 9u−1 27= 0. Since f(0) =−1 27, f/parenleftbig1 3/parenrightbig =1 54, f/parenleftbig2 3/parenrightbig =−1 27, f (1) =1 54, the Intermediate Value Lemma 19.18 guarantees that there ar e three roots on the interval [0,1]: one between 0 and1 3, the second between1 3and2 3, and the third between2 3and 1. The graph in Figure 19.12 reconﬁrms this observation. Since we are dealing with a cubic polynomial, there are no other roots. (Why?) †Thisassumes weareworkinginasuﬃcientlyhighprecisionarithmet icsoastoavoid round-oﬀ errors. 12/11/12 1052 c/ci∇cleco√y∇t2012 Peter J. Olver It takes sixteen iterations of the Bisection Method startin g with the three subintervals/bracketleftbig 0,1 3/bracketrightbig ,/bracketleftbig1 3,2 3/bracketrightbig and/bracketleftbig2 3,1/bracketrightbig , to produce the roots to six decimal places: u⋆ 1≈.085119, u⋆ 2≈.451805, u⋆ 3≈.963076. Incidentally, if we start with the interval [0 ,1] and apply bisection, we converge (perhaps surprisingly) to the largest root u⋆ 3in 17 iterations. Fixed point iteration based on the formulation u=g(u) =−u3+3 2u2+4 9u+1 27 can be used to ﬁnd the ﬁrst and third roots, but not the second r oot. For instance, starting withu(0)= 0 produces u⋆ 1to 5 decimal places after 23 iterations, whereas starting wi th u(0)= 1 produces u⋆ 3to 5 decimal places after 14 iterations. The reason we cannot produce u⋆ 2is due to the magnitude of the derivative g′(u) =−3u2+3u+4 9 at the roots, which is g′(u⋆ 1)≈0.678065, g′(u⋆ 2)≈1.18748, g′(u⋆ 3)≈0.551126. Thus,u⋆ 1andu⋆ 3are stable ﬁxed points, but u⋆ 2is unstable. However, because g′(u⋆ 1) and g′(u⋆ 3) are both bigger than .5, this iterative algorithm actually converges slowerthan ordinary bisection! Finally, Newton’s Method is based upon iteration of the rati onal function g(u) =u−f(u) f′(u)=u−u3−3 2u2+5 9u−1 27 3u2−3u+5 9. Starting with an initial guess of u(0)= 0, the method computes u⋆ 1to 6 decimal places after only 4 iterations; starting with u(0)=.5, it produces u⋆ 2to similar accuracy after 2 iterations; while starting with u(0)= 1 produces u⋆ 3after 3 iterations — a dramatic speed up over the other two methods. Newton’s Method has a very pretty graphical interpretation , that helps us understand what is going on and why it converges so fast. Given the equati onf(u) = 0, suppose we know an approximate value u=u(k)for a solution. Nearby u(k), we can approximate the nonlinear function f(u) by its tangent line y=f(u(k))+f′(u(k))(u−u(k)). (19.43) As long as the tangent line is not horizontal — which requires f′(u(k))/ne}ationslash= 0 — it crosses the axis at u(k+1)=u(k)−f(u(k)) f′(u(k)), 12/11/12 1053 c/ci∇cleco√y∇t2012 Peter J. Olver u(k)u(k+1)f(u) Figure 19.13. Newton’s Method. which represents a new, and, presumably more accurate, appr oximation to the desired root. The procedure is illustrated pictorially in Figure 19 .13. Note that the passage from u(k)tou(k+1)is exactly the Newton iteration step (19.42). Thus, Newtoni an iteration is the same as the approximation of function’s root by those of i ts successive tangent lines. Given a suﬃciently accurate initial guess, Newton’s Method will rapidly produce highly accurate values for the simple roots to the equation i n question. In practice, barring some kind of special exploitable structure, Newton’s Metho d is the root-ﬁnding algorithm of choice. The one caveat is that we need to start the process r easonably close to the root we are seeking. Otherwise, there is no guarantee that a p articular set of iterates will converge, although if they do, the limiting value is necessa rily a root of our equation. The behavior of Newton’s Method as we change parameters and vary the initial guess is very similar to the simpler logistic map that we studied in Sectio n 19.1, including period dou- bling bifurcations and chaotic behavior. The reader is invi ted to experiment with simple examples, some of which are provided in Exercise ; further details can be found in [ 152]. Example 19.26. For ﬁxed values of the eccentricity ǫ, Kepler’s equation u−ǫsinu=m (19.44) can be viewed as a implicit equation deﬁning the eccentric an omalyuas a function of the mean anomaly m. To solve Kepler’s equation by Newton’s Method, we introduc e the iterative function g(u) =u−u−ǫsinu−m 1−ǫcosu. Notice that when |ǫ|<1, the denominator never vanishes and so the iteration remai ns well-deﬁned everywhere. Starting with a suﬃciently close i nitial guess u(0), we are assured that the method will quickly converge to the solution. Fixing the eccentricity ǫ, we can employ the method of continuation to determine how the solution u⋆=h(m) depends upon the mean anomaly m. Namely, we start at 12/11/12 1054 c/ci∇cleco√y∇t2012 Peter J. Olver 0.2 0.4 0.6 0.8 10.20.40.60.811.21.4 Figure 19.14. The Solution to the Kepler Equation for Eccentricity ǫ=.5. m=m0= 0 with the obvious solution u⋆=h(0) = 0. Then, to compute the solution at successive closely spaced values 0 < m1< m2< m3<···, we use the previously computed value as an initial guess u(0)=h(mk) for the value of the solution at the next mesh point mk+1, and run the Newton scheme until it converges to a suﬃciently accurate approximation to the value u⋆=h(mk+1). As long as mk+1is reasonably close to mk, Newton’s Method will converge to the solution quite quickly . The continuation method will quickly produce the values of uat the sample points. Intermediate values can either be determined by an interpol ation scheme, e.g., a cubic spline ﬁt of the data, or by running the Newton scheme using th e closest known value as an initial condition. A plot for 0 ≤m≤1 using the value ǫ=.5 appears in Figure 19.14. Systems of Equations Let us now turn our attention to nonlinear systems of equatio ns. We shall only consider the case when there are the same number of equations as unknowns: f1(u1,...,un) = 0, ... fn(u1,...,un) = 0. (19.45) We shall rewrite the system in vector form f(u) =0, (19.46) wheref:Rn→Rnis a vector-valued function of nvariables. In practice, we do not necessarilyrequirethat fbedeﬁnedonallof Rn, althoughthisdoessimplifytheexposition. We shall only consider solutions that are separated from any others. More formally: Deﬁnition 19.27. A solution u⋆to a system f(u) =0is called isolated if there existsδ >0 such that f(u)/ne}ationslash=0for allusatisfying 0 </ba∇dblu−u⋆/ba∇dbl< δ. Example 19.28. Consider the planar equation x2+y2= (x2+y2)2. 12/11/12 1055 c/ci∇cleco√y∇t2012 Peter J. Olver Rewriting the equation in polar coordinates as r=r2or r(r−1) = 0, we immediately see that the solutions consist of the origin x=y= 0 and all points on the unit circle r2=x2+y2= 1. Only the origin is an isolated solution, since every solu tion lying on the circle has plenty of other points on the circle th at lie arbitrarily close to it. Typically, solutions to a system of nequations in nunknowns are isolated, although this is not always true. For example, if Ais a singular n×nmatrix, then the solutions to the homogeneous linear system Au=0form a nontrivial subspace, and so are not isolated. Nonlinear systems with non-isolated solutions can similar ly be viewed as exhibiting some form of degeneracy. In general, the numerical computation o f non-isolated solutions, e.g., solving the implicit equations for a curve or surface, is a mu ch more diﬃcult problem, and we will not attempt to discuss these issues inthis introduct ory presentation. (However, our continuation approach to the Kepler equation in Example 19. 26 indicates how one might proceed in such situations.) In the case of a single scalar equation, the simple roots, mea ning those for which f′(u⋆)/ne}ationslash= 0, are the easiest to compute. In higher dimensions, the rol e of the derivative of the function is played by the Jacobian matrix (19.28), and this motivates the following deﬁnition. Deﬁnition 19.29. A solution u⋆to a system f(u) =0is called nonsingular if the associated Jacobian matrix is nonsingular there: det f′(u⋆)/ne}ationslash= 0. Note that the Jacobian matrix is square if and only if the syst em has the same number of equations as unknowns, which is thus one of the requiremen ts for a solution to be nonsingular in our sense. Moreover, the Inverse Function Th eorem from multivariable calculus, [ 9,129], implies that a nonsingular solution is necessarily isola ted. Theorem 19.30. Every nonsingular solution u⋆to a system f(u) =0is isolated. Being the multivariate counterparts of simple roots also me ans that nonsingular solu- tions of systems are the most amenable to practical computat ion. Computing non-isolated solutions, as well as isolated solutions with a singular Jac obian matrix, is a considerable challenge, and practical algorithms remain much less well d eveloped. For this reason, we focus exclusively on numerical solution techniques for non singular solutions. Now, let us turn to numerical solution techniques. The ﬁrst r emark is that, unlike the scalar case, proving existence of a solution to a system of eq uations is often a challenging issue. There is no counterpart to the Intermediate Value Lem ma 19.18 for vector-valued functions. It is not hard to ﬁnd vector-valued functions who se entries take on both positive and negative values, but admit no solutions; a simple exampl e can be found in Exercise . This precludes any simple analog of the Bisection Method for nonlinear systems in more than one unknown. On the other hand, Newton’s Method can be straightforwardly adapted to compute nonsingular solutions to systems of equations, and is themost widely used method for this 12/11/12 1056 c/ci∇cleco√y∇t2012 Peter J. Olver purpose. The derivation proceeds in very similar manner to t he scalar case. First, we replace the system (19.46) by a ﬁxed point system u=g(u) (19 .47) having the same solutions. By direct analogy with (19.38), a ny (reasonable) ﬁxed point method will take the form g(u) =u−L(u)f(u), (19.48) whereL(u) is ann×nmatrix-valued function. Clearly, if f(u) =0theng(u) =u; conversely, if g(u) =u, thenL(u)f(u) =0. If we further require that the matrix L(u) be nonsingular, i.e., det L(u)/ne}ationslash= 0, then every ﬁxed point of the iterator (19.48) will be a solution to the system (19.46) and vice versa. According to Theorem 19.12, the speed of convergence (if any ) of the iterative method u(k+1)=g(u(k)) (19 .49) is governed by the matrix norm (or, more precisely, the spect ral radius) of the Jacobian matrixg′(u⋆) at the ﬁxed point. In particular, if g′(u⋆) = O (19 .50) is the zero matrix, then the method converges quadratically fast. Let’s ﬁgure out how this can be arranged. Computing the derivative using the matrix v ersion of the Leibniz rule for the derivative of a matrix product, (9.40), we ﬁnd g′(u⋆) = I−L(u⋆)f′(u⋆), (19.51) where I is the n×nidentity matrix. (Fortunately, all the terms that involve d erivatives of the entries of L(u) go away since f(u⋆) =0by assumption; details are relegated to Exercise .) Therefore, the quadratic convergence criterion (19.50) holds if and only if L(u⋆)f′(u⋆) = I,and hence L(u⋆) =f′(u⋆)−1(19.52) should be the inverse of the Jacobian matrix of fat the solution, which, fortuitously, was already assumed to be nonsingular. Asinthescalarcase, wedon’tknow thesolution u⋆, but wecanarrangethatcondition (19.52) holds by setting L(u) =f′(u)−1 everywhere — or at least everywhere that fhas a nonsingular Jacobian matrix. The resulting ﬁxed point system u=g(u) =u−f′(u)−1f(u), (19.53) leads to the quadratically convergent Newton iteration scheme u(k+1)=u(k)−f′(u(k))−1f(u(k)). (19.54) All it requires is that we guess an initial value u(0)that is suﬃciently close to the de- sired solution u⋆. We are then guaranteed, by Exercise , that the iterates u(k)converge quadratically fast to u⋆. 12/11/12 1057 c/ci∇cleco√y∇t2012 Peter J. Olver Theorem 19.31. Letu⋆be a nonsingular solution to the system f(u) =0. Then, provided u(0)is suﬃciently close to u⋆, the Newton iteration scheme (19.54)converges at a quadratic rate to the solution: u(k)→u⋆. Example 19.32. Consider the pair of simultaneous cubic equations f1(u,v) =u3−3uv2−1 = 0, f2(u,v) = 3u2v−v3= 0. (19.55) It is not diﬃcult to prove that there are precisely three solu tions: u⋆ 1=/parenleftbigg 1 0/parenrightbigg , u⋆ 2=/parenleftbigg −.5 .866025.../parenrightbigg , u⋆ 3=/parenleftbigg −.5 −.866025.../parenrightbigg .(19.56) The Newton scheme relies on the Jacobian matrix f′(u) =/parenleftbigg 3u2−3v2−6uv 6uv3u2−3v2/parenrightbigg . Since det f′(u) = 9(u2+v2) is non-zero except at the origin, all three solutions are no n- singular, and hence, for a suﬃciently close initial value, N ewton’s Method will converge to the nearby solution. We explicitly compute the inverse Jaco bian matrix: f′(u)−1=1 9(u2+v2)/parenleftbigg 3u2−3v26uv −6uv3u2−3v2/parenrightbigg . Hence, in this particular example, the Newton iterator (19. 53) is g(u) =/parenleftbigg u v/parenrightbigg −1 9(u2+v2)/parenleftbigg 3u2−3v26uv −6uv3u2−3v2/parenrightbigg/parenleftbigg u3−3uv2−1 3u2v−v3/parenrightbigg . A complete diagram of the three basins of attraction, consis ting of points whose New- ton iterates converge to each of the three roots, has a remark ably complicated, fractal-like structure, as illustrated in Figure 19.15. In this plot, the xandycoordinates run from −1.5 to 1.5. The points in the black region all converge to u⋆ 1; those in the light gray region all converge to u⋆ 2; while those in the dark gray region all converge to u⋆ 3. The closer one is to the root, the sooner the iterates converge. On the inter faces between the basins of attraction are points for which the Newton iterates fail to c onverge, but exhibit a random, chaotic behavior. However, round-oﬀ errors will cause such iterates to fall into one of the basins, making it extremely diﬃcult to observe such behavio over the long run. Remark: The alert reader may notice that in this example, we are in fa ct merely computing the cube roots of unity, i.e., equations (19.55) a re the real and imaginary parts of the complex equation z3= 1 when z=u+ iv. Example 19.33. A robot arm consists of two rigid rods that are joined end-to- end to a ﬁxed point in the plane, which we take as the origin 0. The arms are free to rotate, and the problem is to conﬁgure them so that the robot’s hand en ds up at the prescribed position a= (a,b)T. The ﬁrst rod has length ℓand makes an angle αwith the horizontal, soitsendisatposition v1= (ℓcosα,ℓsinα)T. Thesecondrodhaslength mandmakesan 12/11/12 1058 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 19.15. Computing the Cube Roots of Unity by Newton’s Method. αβ ℓm v1v1+v2 Figure 19.16. Robot Arm. angleβwith the horizontal, and so is represented by the vector v2= (mcosβ,msinβ)T. The hand at the end of the second arm is at position v1+v2, and the problem is to ﬁnd values for the angles α,βso thatv1+v2=a; see Figure 19.16. To this end, we need to solve the system of equations ℓcosα+mcosβ=a, ℓ sinα+msinβ=b, (19.57) for the angles α,β. To ﬁnd the solution, weshall apply Newton’sMethod. First, w ecompute the Jacobian 12/11/12 1059 c/ci∇cleco√y∇t2012 Peter J. Olver matrix of the system with respect to α,β, which is f′(α,β) =/parenleftbigg −ℓsinα−msinβ ℓcosα mcosβ/parenrightbigg , with inverse f′(α,β)−1=1 ℓmsin(β−α)/parenleftbigg −ℓsinα msinβ −ℓcosα mcosβ/parenrightbigg . As a result, the Newton iteration equation (19.54) has the ex plicit form /parenleftbigg α(k+1) β(k+1)/parenrightbigg =/parenleftbigg α(k) β(k)/parenrightbigg − −1 ℓmsin(β(k)−α(k))/parenleftbigg −ℓcosα(k)msinβ(k) −ℓcosα(k)msinβ(k)/parenrightbigg/parenleftbigg ℓcosα(k)+mcosβ(k)−a ℓsinα(k)+msinβ(k)−b/parenrightbigg . when running the iteration, one must be careful to avoid poin ts at which α(k)−β(k)= 0 orπ, i.e., where the robot arm has straightened out. As an example, let us assume that the rods have lengths ℓ= 2,m= 1, and the desired location of the hand is at a= (1,1)T. We start with an initial guess of α(0)= 0, β(0)=1 2π, so the ﬁrst rod lies along the x–axis and the second is perpendicular. The ﬁrst few Newton iterates are given in the accompanying table. The ﬁrst column is the iterate numberk; the second and third columns indicate the angles α(k),β(k)of the rods. The fourth and ﬁfth give the position ( x(k),y(k))Tof the joint or elbow, while the ﬁnal two indicate the position ( z(k),w(k))Tof the robot’s hand. kα(k)β(k)x(k)y(k)z(k)w(k) 0.0000 1 .5708 2 .0000 .0000 2 .0000 1 .0000 1.0000 2 .5708 2 .0000 .0000 1 .1585 .5403 2.3533 2 .8642 1 .8765 .6920 .9147 .9658 3.2917 2 .7084 1 .9155 .5751 1 .0079 .9948 4.2987 2 .7176 1 .9114 .5886 1 .0000 1 .0000 5.2987 2 .7176 1 .9114 .5886 1 .0000 1 .0000 Observe thattherobot hasrapidlyconverged to oneofthetwo possible conﬁgurations. (Can you ﬁgure out what the second equilibrium is?) In genera l, convergence depends on the choice of initial conﬁguration, and the Newton iterates do not always settle down to a ﬁxed point. For instance, if /ba∇dbla/ba∇dbl> ℓ+m, there is no possible solution, since the arms are too short for the hand to reach to desired location; thus, no choice of initial conditions will lead to a convergent scheme and the robot arm ﬂaps around in a chaotic manner. Now that we have gained a little experience with Newton’s Met hod for systems of equations, some supplementary remarks are in order. As we le arned back in Chapter 1, except perhaps in very low-dimensional situations, one sho uld not directly invert a matrix, 12/11/12 1060 c/ci∇cleco√y∇t2012 Peter J. Olver but rather use Gaussian elimination, or, in favorable situa tions, a linear iterative scheme, e.g., Jacobi, Gauss–Seidel or even SOR. So a better strategy is to leave the Newton system (19.54) in unsolved, implicit form f′(u(k))v(k)=−f(u(k)),u(k+1)=u(k)+v(k). (19.58) Given the iterate u(k), we compute the Jacobian matrix f′(u(k)) and the right hand side −f(u(k)), and then use our preferred linear systems solver to ﬁnd v(k). Adding u(k)to the result immediately yields the updated approximation u(k+1)to the solution. The main bottleneck in the implementation of the Newton sche me, particularly for large systems, is solving the linear system in (19.58). The c oeﬃcient matrix f′(u(k)) must berecomputedateachstepoftheiteration,andhenceknowin gthesolutiontothe kthlinear system does not appear to help us solve the subsequent system . Pereforming a complete Gaussian elimination at every step will tend to slow down the algorithm, particularly in high dimensional situations involving many equations in ma ny unknowns. One simple dodge for speeding up the computation is to note th at, once we start converging, u(k)will be very close to u(k−1)and so we will probably not go far wrong by usingf′(u(k−1)) in place of the updated Jacobian matrix f′(u(k)). Since we have already solved the linear system with coeﬃcient matrix f′(u(k−1)), we know its LUfactorization, and hence can use Forward and Back Substitution to quickly so lve the modiﬁed system f′(u(k−1))v(k+1)=−f(u(k)),u(k+1)=u(k)+v(k). (19.59) Ifu(k+1)is still close to u(k−1), we can continue to use f′(u(k−1)) as the coeﬃcient matrix when proceeding on to the next iterate u(k+2). We proceed in this manner until there has been a notable change in the iterates, at which stage we ca n revert to solving the correct, unmodiﬁed linear system (19.58) by Gaussian Elimi nation. This strategy may dramatically reduce the total amount of computation requir ed to approximate the solution to a prescribed accuracy. The down side is that this quasi-Newton scheme is usually only linearly convergent, and so does not home in on the root a s fast as the unmodiﬁed implementation. The user needs to balance the trade-oﬀ betw een speed of convergence versus amount of time needed to solve the linear system at eac h step in the process. See [153] for further discussion. 19.3. Optimization. We have already noted the importance of quadratic minimizat ion principles for char- acterizing the equilibrium solutions of linear systems of p hysical signiﬁcance. In nonlinear systems, optimization — either maximization or minimizati on — retains its centrality, and the wealth of practical applications has spawned an enti re subdiscipline of applied mathematics. Physical systems naturally seek to minimize t he potential energy function, and so determination of the possible equilibrium conﬁgurat ions requires solving a non- linear minimization principle. Engineering design is guid ed by a variety of optimization constraints, such as performance, longevity, safety, and c ost. Non-quadratic minimization principles also arise in the ﬁtting of data by schemes that go beyond the simple linear least squares approximation method discussed in Section 4.3. Add itional applications naturally 12/11/12 1061 c/ci∇cleco√y∇t2012 Peter J. Olver appear in economics and ﬁnancial mathematics — one often wis hes to minimize expenses or maximize proﬁts, in biological and ecological systems, i n pattern recognition and signal processing, in statistics, and so on. In this section, we wil l describe the basic mathematics underlying simple nonlinear optimization problems along w ith basic numerical techniques. The Objective Function Throughout this section, the real-valued function F(u) =F(u1,...,un) to be op- timized — the energy, cost, entropy, performance, etc. — wil l be called the objective function. As such, it depends upon one or more variables u= (u1,u2,...,un)Tthat belong to a prescribed subset Ω ⊂Rn. Deﬁnition 19.34. A pointu⋆∈Ω is aglobal minimum of the objective function F(u⋆)≤F(u) for all u∈Ω. (19.60) The minimum is called strictif F(u⋆)< F(u) for u⋆/ne}ationslash=u∈Ω. (19.61) The point u⋆is called a ( strict)local minimum if the relevant inequality holds just for pointsu∈Ω nearby u⋆, i.e., satisfying/ba∇dblu−u⋆/ba∇dbl< δfor some δ >0. In particular, strict local minima are isolated. The deﬁnition of a maximum — local or global — is the same, but with the reversed inequality: F(u⋆)≥F(u) or, in the strict case, F(u⋆)> F(u). Alternatively, a maximum ofF(u) is the same as a minimum of the negative −F(u). Therefore, every result that appliestominimizationofafunctioncaneasilybetranslat edintoaresult onmaximization, which allowsus to concentrate exclusively on the minimizat ionproblem without any loss of generality. Wewilluse extremum asashorthandtermfor eitheraminimum oramaximum. Remark: As we already noted in Section 4.1, anysystem of equations can be readily converted into a minimization principle. Given a system f(u) =0, we introduce the objective function F(u) =/ba∇dblf(u)/ba∇dbl2, (19.62) where/ba∇dbl·/ba∇dblis any convenient norm on Rn. By the basic properties of the norm, the minimum value is F(u) = 0, and this is achieved if and only if f(u) =0, i.e., at a solution to the system. More generally, if there is no solution to the s ystem, the minimizer(s) of F(u) play the role of a least squares solution, at least for an inn er product-based norm, along with the extensions to more general norms. In contrast to the rather diﬃcult question of existence of so lutions to systems of equations, there is a general theorem that guarantees the ex istence of minima (and, hence, maxima) for a broad class of optimization problems. Theorem 19.35. IfF:Ω→Ris continuous, and Ω⊂Rnis acompact, meaning closed and bounded, subset, then Fhas at least one global minimum u⋆∈Ω. 12/11/12 1062 c/ci∇cleco√y∇t2012 Peter J. Olver -1 -0.5 0.5 1246 Figure 19.17. The function 8 u3+5u2−6u. Proof: Letm⋆= min{F(u)|u∈Ω}, which may, a priori, be−∞. Choose points u(1),u(2),...∈Ω, such that F(u(k))→mask→∞. By the basic properties of compact sets, [158], there is a convergent subsequence u(ki)→u⋆∈Ω. By continuity, F(u⋆) =F/parenleftbigg lim ki→∞u(ki)/parenrightbigg = lim ki→∞F/parenleftig u(ki)/parenrightig =m⋆, and hence u⋆is a minimizer. Q.E.D. Although Theorem 19.35 assures us of the existence of a globa l minimum of any continuous function on a bounded domain, it does not guarant ee uniqueness, nor does it indicate how to go about ﬁnding it. Just as with the solutio n of nonlinear systems of equations, it is quite rare that one can extract explicit f ormulae for the minima of non-quadratic functions. Our goal, then, is to formulate pr actical algorithms that can accurately compute the minima of general nonlinear functio ns. The most na¨ ıve algorithm, but one that is often successful i n small scale problems, [153,opt], is to select a reasonably dense set of sample points u(k)in the domain and choose the one that provides the smallest value for F(u(k)). If the points are suﬃciently densely distributed and the function is not too wild, this wi ll give a reasonable approxima- tionto the minimum. The algorithm can bespeeded up by appeal ing to moresophisticated means of selecting the sample points. In the rest of this section, we will discuss optimization str ategies that exploit the diﬀerential calculus. Let us ﬁrst review the basic procedur e for optimizing functions that you learned in ﬁrst and second year calculus. As you no doubt r emember, there are two diﬀerent possible types of minima. An interior minimum occurs at an interior point of the domain of deﬁnition of the function, whereas a boundary minimum occurs on its boundary ∂Ω. Interior local minima are easier to ﬁnd, and, to keep the pr esentation simple, we shall focus our eﬀorts on them. Let us begin with a simple scalar exa mple. Example 19.36. Let us optimize the scalar objective function F(u) = 8u3+5u2−6u on the domain−1≤u≤1. To locate the minimum, the ﬁrst step is to look at the critical 12/11/12 1063 c/ci∇cleco√y∇t2012 Peter J. Olver pointswhere the derivative vanishes: F′(u) = 24u2+10u−6 = 0,and hence u=1 3,−3 4. To ascertain the local nature of the two critical points, we a pply the second derivative test. SinceF′′(u) = 48u+10, we have F′′/parenleftbig1 3/parenrightbig = 26>0,whereas F′′/parenleftbig −3 4/parenrightbig =−26<0. We conclude that1 3is a local minimum, while3 4is a local maximum. To ﬁnd the global minimum and maximum on the interval [ −1,1], we must also take into account the boundary points ±1. Comparing the function values at the four points, F(−1) = 3, F/parenleftbig1 3/parenrightbig =−31 27≈−1.148, F/parenleftbig −3 4/parenrightbig =63 16= 3.9375, F(1) = 7, we see that1 3is the global minimum, whereas 1 is the global maximum — which occurs on the boundary of the interval. This is borne out by the graph of the function, as displayed in Figure 19.17. The Gradient As you ﬁrst learn in multi-variable calculus, [ 9,129], the interior extrema — minima and maxima — of a smooth function F(u) =F(u1,...,un) are necessarily critical points , meaning places where the gradient of Fvanishes. The standard gradient is the vector ﬁeld whose entries are its ﬁrst order partial derivatives: ∇F(u) =/parenleftbigg∂F ∂u1, ... ,∂F ∂un/parenrightbiggT . (19.63) Letus, inpreparationforthemoregeneralminimizationpro blemsoverinﬁnite-dimensional function spaces to be treated in Chapter 21, reformulate the deﬁnition of the gradient in a more intrinsic manner. An important but subtle point is tha t the gradient operator, in fact, relies upon the introduction of an inner product on the underlying vector space. The version (19.63) is, in fact, based upon on the Euclidean dot p roduct on Rn. Altering the inner product will change the formula for the gradient! Deﬁnition 19.37. LetVbe an inner product space. The gradient of a function F:V→Rat a point u∈Vis the vector∇F(u)∈Vthat satisﬁes /an}b∇acketle{t∇F(u);v/an}b∇acket∇i}ht=d dtF(u+tv)/vextendsingle/vextendsingle/vextendsingle/vextendsingle t=0for all v∈V. (19.64) Remark: The function Fdoes not have to be deﬁned on all of the space Vin order for this deﬁnition to make sense. The quantity displayed in the preceding formula is known as t hedirectional derivative ofFwith respect to v∈V, and typically denoted by ∂F/∂v. Thus, by deﬁnition, the directional derivative equals the inner product with the gr adient vector. The directional 12/11/12 1064 c/ci∇cleco√y∇t2012 Peter J. Olver derivative measures the rate of change of Fin the direction of the vector v, scaled in proportion to its length. In the Euclidean case, when F(u) =F(u1,...,un) is a function of nvariables, deﬁned foru= (u1,u2,...,un)T∈Rn, we can use the chain rule to compute d dtF(u+tv) =d dtF(u1+tv1,...,un+tvn) =∂F ∂u1(u+tv)v1+···+∂F ∂un(u+tv)vn.(19.65) Settingt= 0, the right hand side of (19.64) reduces to d dtF(u+tv)/vextendsingle/vextendsingle/vextendsingle/vextendsingle t=0=∂F ∂u1(u)v1+···+∂F ∂un(u)vn=∇F(u)·v. Therefore, the directional derivative equals the Euclidea n dot product between the usual gradient of the function (19.63) and the direction vector v, justifying (19.64) in the Eu- clidean case. Remark: In this chapter, we will only deal with the standard Euclide an dot product, which results in the usual gradient formula (19.63). If we in troduce an alternative inner product on Rn, then the notion of gradient, as deﬁned in (19.64) will chang e. Details are outlined in Exercise . A function F(u) iscontinuously diﬀerentiable if and only if its gradient ∇F(u) is a continuously varying vector-valued function of u. This is equivalent to the requirement that its ﬁrst order partial derivatives ∂F/∂uiare all continuous. As usual, we use C1(Ω) to denote the vector space of all continuously diﬀerentiabl e scalar-valued functions deﬁned on a domain Ω⊂Rn. From now on, all objective functions are assumed to be conti nuously diﬀerentiable on their domain of deﬁnition. Ifu(t) represents a parametrized curve contained within the doma in of deﬁnition of F(u), then the instantaneous rate of change in the scalar quanti tyFas we move along the curve is given by d dtF(u(t)) =/angbracketleftbigg ∇F(u);du dt/angbracketrightbigg , (19.66) which is the directional derivative of Fwith respect to the velocity or tangent vector v=/squaresmallsoliduto the curve. For instance, suppose F(u1,u2) represents the elevation of a mountain range at position u= (u1,u2)T. If we travel through the mountains along the path u(t) = (u1(t),u2(t))T, then our instantaneous rate of ascent or descent (19.66) is equal to the dot product of our velocity vector/squaresmallsolidu(t) with the gradient of the elevation function. This observation leads to an important interpretation of th e gradient vector. Theorem 19.38. The gradient∇F(u)of a scalar function F(u)points in the direc- tion of its steepest increase at the point u. The negative gradient, −∇F(u), which points in the opposite direction, indicates the direction of steep est decrease. 12/11/12 1065 c/ci∇cleco√y∇t2012 Peter J. Olver Thus, when Frepresents elevation, ∇Ftells us the direction that is steepest uphill, while−∇Fpoints directly downhill — the direction water will ﬂow. Sim ilarly, if F represents the temperature of a solid body, then ∇Ftells us the direction in which it is heating up the quickest. Heat energy (like water) will ﬂow in the opposite, coldest direction, namely that of the negative gradient vector −∇F. But you need to be careful in how you interpret Theorem 19.38. Clearly, the faster you move along a curve, the faster the function F(u) will vary, and one needs to take this into account when comparing the rates of change along di ﬀerent curves. The easiest way to eﬀect the comparison is to assume that the tangent vect ora=/squaresmallsoliduhas unit norm, so/ba∇dbla/ba∇dbl= 1, which means that we are passing through the point u(t) with unit speed. Once this is done, Theorem 19.38 is an immediate consequence of the Cauchy–Schwarz inequality (3.18). Indeed, /vextendsingle/vextendsingle/vextendsingle/vextendsingle∂F ∂a/vextendsingle/vextendsingle/vextendsingle/vextendsingle=|a·∇F|≤/ba∇dbla/ba∇dbl/ba∇dbl∇F/ba∇dbl=/ba∇dbl∇F/ba∇dbl,when/ba∇dbla/ba∇dbl= 1, with equality if and only if apoints in the same direction as the gradient. Therefore, the maximum rate of change is when a=∇F//ba∇dbl∇F/ba∇dblis the unit vector in the gradient direction, while the minimum is achieved when a=−∇F//ba∇dbl∇F/ba∇dblpoints in the opposite direction. Critical Points Thus, the only points at which the gradient fails to indicate directions of increase/de- crease of the objective fulnction are where it vanishes. Suc h points play a critical role in the analysis, whence the following deﬁnition. Deﬁnition 19.39. A pointu⋆is called a critical point of the objective function F(u) if ∇F(u⋆) =0. (19.67) Let usprovethatalllocalminimaareindeed criticalpoints . Themostimportantthing about this proof is that it only relies on the intrinsic deﬁni tion of gradient, and therefore applies to any function on any inner product space. Moreover , even though the gradient will change if we alter the underlying inner product, the req uirement that it vanish at a local minimum does not. Theorem 19.40. Every local (interior)minimum u⋆of a continuously diﬀerentiable function F(u)is a critical point :∇F(u⋆) =0. Proof: Let0/ne}ationslash=v∈Rnbe any vector. Consider the scalar function g(t) =F(u⋆+tv) =F(u⋆ 1+tv1,...,u⋆ n+tvn), (19.68) wheret∈Ris suﬃciently small to ensure that u⋆+tvremains strictly inside the domain ofF. Note that gmeasures the values of Falong a straight line passing through u⋆in the direction prescribed by v. Sinceu⋆is a local minimum, F(u⋆)≤F(u⋆+tv),and hence g(0)≤g(t) 12/11/12 1066 c/ci∇cleco√y∇t2012 Peter J. Olver for alltsuﬃciently close to zero. In other words, g(t), as a function of the single variable t, has a local minimum at t= 0. By the basic calculus result on minima of functions of one variable, the derivative of g(t) must vanish at t= 0. Therefore, by the deﬁnition (19.64) of gradient, 0 =g′(0) =d dtF(u⋆+tv)/vextendsingle/vextendsingle/vextendsingle/vextendsingle t=0=/an}b∇acketle{t∇F(u⋆);v/an}b∇acket∇i}ht. We conclude that the gradient vector ∇F(u⋆) at the critical point must be orthogonal to everyvectorv∈Rn, which is only possible if ∇F(u⋆) =0. Q.E.D. Thus, provided the objective function is continuously diﬀe rentiable, every interior minimum, both local and global, is necessarily a critical po int. The converse is not true; critical points can be maxima; they can also be saddle points or of some degenerate form. The basic analytical method†for determining the (interior) minima of a given function is to ﬁrst ﬁnd all its critical points by solving the system of eq uations (19.67). Each critical point then needs to be examined more closely — as it could be ei ther a minimum, or a maximum, or neither. Example 19.41. Consider the function F(u,v) =u4−2u2+v2, whichisdeﬁnedandcontinuouslydiﬀerentiableonallof R2. Since∇F=/parenleftbig 4u3−4u,2v/parenrightbigT, its critical points are obtained by solving the pair of equat ions 4u3−4u= 0,2v= 0. The solutions to the ﬁrst equation are u= 0,±1, while the second equation requires v= 0. Therefore, Fhas three critical points: u⋆ 1=/parenleftbigg 0 0/parenrightbigg , u⋆ 2=/parenleftbigg 1 0/parenrightbigg , u⋆ 3=/parenleftbigg −1 0/parenrightbigg . (19.69) Inspecting its graph in Figure 19.18, we suspect that the ﬁrs t critical point u⋆ 1is a saddle point, whereas the other two appear to be local minima, havin g the same value F(u⋆ 2) = F(u⋆ 3) =−1. This will be conﬁrmed once we learn how to rigorously disti nguish critical points. The student should also pay attention to the distinction bet ween local minima and global minima. In the absence of theoretical justiﬁcation, the only practical way to deter- mine whether or not a minimum is global is to ﬁnd all the diﬀere nt local minima, including those on the boundary, and see which one gives the smallest va lue. If the domain is un- bounded, one must also worry about the asymptotic behavior o f the objective function for largeu. †Numerical methods are discussed below. 12/11/12 1067 c/ci∇cleco√y∇t2012 Peter J. Olver Figure 19.18. The Function u4−2u2+v2. The Second Derivative Test The statusofcriticalpoint —minimum, maximum, orneither — canoftenberesolved by analyzing the second derivative of the objective functio n at the critical point. Let us ﬁrst review the one variable second derivative test you lear ned in ﬁrst year calculus. Proposition 19.42. Letg(t)∈C2beascalarfunction, andsupposethat t⋆acritical point:g′(t⋆) = 0. Ift⋆is a local minimum, then g′′(t⋆)≥0. Conversely, if g′′(t⋆)>0, thent⋆is a strict local minimum. Similarly, g′′(t⋆)≤0is required at a local maximum, whileg′′(t⋆)<0implies that t⋆is a strict local maximum. The proof of this result relies on the fact that we can approxi mate the function by its quadratic Taylor polynomial near the critical point: g(t)≈g(t⋆)+1 2(t−t⋆)2g′′(t⋆), sinceg′(t⋆) = 0, and so the linear terms in the Taylor polynomial vanish. Ifg′′(t⋆)/ne}ationslash= 0, then the quadratic Taylor polynomial has a minimum or maximu m att⋆according to the sign of the second derivative, and this provides the key to th e proof. In the borderline case, when g′′(t⋆) = 0, the second derivative test is inconclusive, and the poi nt could be either maximum or minimum or neither. One must analyze the hi gher order terms in the Taylor expansion to resolve the staus of the critical point; see Exercise . In multi-variate calculus, the “second derivative” of a fun ctionF(u) =F(u1,...,un) is represented by the n×nHessian†matrix, whose entries are its second order partial †Named after the early eighteenth century German mathematician Ludwi g Otto Hesse. 12/11/12 1068 c/ci∇cleco√y∇t2012 Peter J. Olver derivatives: ∇2F(u) = ∂2F ∂u2 1∂2F ∂u1∂u2...∂2F ∂u1∂un ∂2F ∂u2∂u1∂2F ∂u2 2...∂2F ∂u2∂un ............ ∂2F ∂un∂u1∂2F ∂un∂u2...∂2F ∂u2 n , (19.70) We will always assume that F(u)∈C2has continuous second order partial derivatives. In this case, its mixed partial derivatives are equal: ∂2F/∂ui∂uj=∂2F/∂uj∂ui, cf. [9,129]. As a result, the Hessian is a symmetric matrix: ∇2F(u) =∇2F(u)T. The second derivative test for a local minimum of scalar func tion relies on the positiv- ity of its second derivative. For a function of several varia bles, the corresponding condition is that the Hessian matrix be positive deﬁnite, as in Deﬁniti on 3.20. More speciﬁcally: Theorem 19.43. LetF(u) =F(u1,...,un)∈C2(Ω)be a real-valued, twice con- tinuously diﬀerentiable function deﬁned on an open domain Ω⊂Rn. Ifu⋆∈Ωis a (local, interior )minimum for F, then it is necessarily a critical point, so ∇F(u⋆) =0. Moreover, the Hessian matrix (19.70)must be positive semi-deﬁnite at the minimum, so ∇2F(u⋆)≥0. Conversely, if u⋆is a critical point with positive deﬁnite Hessian matrix ∇2F(u⋆)>0, thenu⋆is a strict local minimum of F. A maximum requires a negative semi-deﬁnite Hessian matrix. If, moreover, the Hes- sian at the critical point is negative deﬁnite, then the crit ical point is a strict local maxi- mum. If the Hessian matrix is indeﬁnite, then the critical po int is a saddle point — neither minimum nor maximum. In the borderline case, when the Hessia n is only positive or nega- tive semi-deﬁnite at the critical point, the second derivat ive test is inconclusive. Resolving the nature of the critical point requires more detailed know ledge of the objective function, e.g., its higher order derivatives. We defer the proof of Theorem 19.43 until the end of this secti on. Example 19.44. As a ﬁrst, elementary example, consider the quadratic funct ion F(u,v) =u2−2uv+3v2. To minimize F, we begin by computing its gradient ∇F(u,v) =/parenleftbigg 2u−2v −2u+6v/parenrightbigg . Solving the pair of equations ∇F=0, namely 2u−2v= 0,−2u+6v= 0, 12/11/12 1069 c/ci∇cleco√y∇t2012 Peter J. Olver we see that the only critical point is the origin u=v= 0. To test whether the origin is a maximum or minimum, we further compute the Hessian matrix H=∇2F(u,v) =/parenleftbigg FuuFuv FuvFvv/parenrightbigg =/parenleftbigg 2−2 −2 6/parenrightbigg . Using the methods of Section 3.5, we easily prove that the Hes sian matrix is positive deﬁnite. Therefore, by Theorem 19.43, u⋆=0is a strict local minimum of F. Indeed, we recognize F(u,v) to be, in fact, a homogeneous positive deﬁnite quadratic form, which can be written in the form F(u,v) =uTKu,where K=/parenleftbigg 1−1 −1 3/parenrightbigg =1 2H,u=/parenleftbigg u v/parenrightbigg . Positive deﬁniteness of the coeﬃcient matrix Kimplies that F(u,v)>0 for all u= (u,v)T/ne}ationslash=0, and hence 0is, in fact, a global minimum. In general, any quadratic function Q(u) =Q(u1,...,un) can be written in the form Q(u) =uTKu−2bTu+c=m/summationdisplay i,j=1kijuiuj−2n/summationdisplay i=1biui+c, (19.71) whereK=KTis a symmetric n×nmatrix,b∈Rnis a ﬁxed vector, and c∈Ris a scalar. A straightforward computation produces the formul a for its gradient and Hessian matrix: ∇Q(u) = 2Ku−2b,∇2Q(u) = 2K. (19.72) As a result, the critical points of the quadratic function ar e the solutions to the linear systemKu=b. IfKis nonsingular, there is a unique critical point u⋆, which is a strict local minimum if and only if K >0 is positive deﬁnite. In fact, Theorem 4.1 tells us that, in the positive deﬁnite case, u⋆is a strict globalminimum for Q(u). Thus, the algebraic approach of Chapter 4 provides additional, global informat ion that cannot be gleaned directly from the local, multivariable calculus Theorem 19 .43. But algebra is only able to handle quadratic minimization problems with ease. The anal ytical classiﬁcation of minima and maxima of more complicated objective functions necessa rily relies the gradient and Hessian criteria of Theorem 19.43. Example 19.45. The function F(u,v) =u2+v2−v3has gradient ∇F(u,v) =/parenleftbigg 2u 2v−3v2/parenrightbigg . The critical point equation ∇F=0has two solutions: u⋆ 1=/parenleftbigg 0 0/parenrightbigg andu⋆ 2=/parenleftbigg 0 2 3/parenrightbigg . The Hessian matrix of the objective function is ∇2F(u,v) =/parenleftbigg 2 0 0 2−6v/parenrightbigg . 12/11/12 1070 c/ci∇cleco√y∇t2012 Peter J. Olver u2+v2−v3u2+v4u2+v3 Figure 19.19. Critical Points. At the ﬁrst critical point, the Hessian ∇2F(0,0) =/parenleftbigg 2 0 0 2/parenrightbigg is positive deﬁnite. Therefore, the origin is a strict local minimum. On the other hand, ∇2F/parenleftbig 0,2 3/parenrightbig =/parenleftbigg 2 0 0−2/parenrightbigg is indeﬁnite, and hence u⋆ 2=/parenleftbigg 0 2 3/parenrightbigg a saddle point. The function is graphed in Figure 19.19, with the critical points indicated by the small solid balls. The origin is, in fact, only a local minimum, since F(0,0) = 0, whereas F(0,v)<0 for all v >1. Thus, this particular function has no global minimum or maximum on R2. Next, consider the function F(u,v) =u2+v4,with gradient ∇F(u,v) =/parenleftbigg 2u 4v3/parenrightbigg . The only critical point is the origin u=v= 0. The origin is a strict global minimum becauseF(u,v)>0 =F(0,0) for all ( u,v)/ne}ationslash= (0,0)T. However, its Hessian matrix ∇2F(u,v) =/parenleftbigg 2 0 0 12v2/parenrightbigg is only positive semi-deﬁnite at the origin, ∇2F(0,0) =/parenleftbigg 2 0 0 0/parenrightbigg . On the other hand, the origin u=v= 0 is also the only critical point for the function F(u,v) =u2+v3with∇F(u,v) =/parenleftbigg 2u 3v2/parenrightbigg . The Hessian matrix is ∇2F(u,v) =/parenleftbigg 2 0 0 6v/parenrightbigg ,and so∇2F(0,0) =/parenleftbigg 2 0 0 0/parenrightbigg is the same positive semi-deﬁnite matrix at the critical poi nt. However, in this case (0 ,0) is not a local minimum; indeed F(0,v)<0 =F(0,0) whenever v <0, and so there exist points arbitrarily close to the origin whe reFtakes on smaller values. As illustrated in Figure 19.19, the origin is, in fact, a dege nerate saddle point. 12/11/12 1071 c/ci∇cleco√y∇t2012 Peter J. Olver Finally, the function F(u,v) =u2−2uv+v2has gradient∇F(u,v) =/parenleftbigg 2u−2v −2u+2v/parenrightbigg , and so every point u=vis a critical point. The Hessian matrix ∇2F(u,v) =/parenleftbigg FuuFuv FuvFvv/parenrightbigg =/parenleftbigg 2−2 −2 2/parenrightbigg is positive semi-deﬁnite everywhere. Since F(u,u) = 0, while F(u,v) = (u−v)2>0 whenu/ne}ationslash=v, each of these critical points is a non-isolated (and hence n on-strict) local minimum. Thus, comparing the three preceding examples, we s ee that a semi-deﬁnite Hessian is unable to distinguish between diﬀerent types of d egenerate critical points. Finally, the reader should always keep in mind that ﬁrst and s econd derivative tests only determine the local behavior of the function near the cr itical point. They cannot be used to determine whether or not we are at a global minimum. This requires some additional analysis, and, often, a fair amount of ingenuity . Proof of Theorem 19.43 : We return to the proof of Theorem 19.40. Given a local minimum u⋆, the scalar function g(t) =F(u⋆+tv) in (19.68) has a local minimum at t= 0. As noted above, basic calculus tells us that its derivati ves att= 0 must satisfy g′(0) = 0, g′′(0)≥0. (19.73) The ﬁrst condition leads to the critical point equation ∇F(u⋆) =0. A straightforward chain rule calculation produces the formula g′′(0) =n/summationdisplay i,j=1∂2F ∂ui∂uj(u⋆)vivj=vT∇2F(u⋆)v. As a result, the second condition in (19.73) requires that vT∇2F(u⋆)v≥0. Sincethisconditionisrequiredforeverydirection v∈Rn,theHessianmatrix ∇2F(u⋆)≥0 satisﬁes the criterion for positive semi-deﬁniteness, pro ving the ﬁrst part of the theorem. The proof of the converse relies†on the second order Taylor expansion of the function: F(u) =F(u⋆)+∇F(u⋆)·v+1 2vT∇2F(u⋆)v+S(v,u⋆) =F(u⋆)+1 2vT∇2F(u⋆)v+S(v,u⋆),where v=u−u⋆, (19.74) †Actually, it is not hard to prove the ﬁrst part using the ﬁrst order T aylor expansion without resorting to the scalar function g; see Exercise . On the other hand, when we look at inﬁnite- dimensional minimization problems arising in the calculus of variation s, we will no longer have the luxury of appealilng to the ﬁnite-dimensional Taylor expansion, whe reas the previous argument continues to apply in general contextx. 12/11/12 1072 c/ci∇cleco√y∇t2012 Peter J. Olver at the critical point, whence ∇F(u⋆) =0. The remainder term in the Taylor formula goes to 0 asu→u⋆at a rate faster than quadratic: S(v,u⋆) /ba∇dblv/ba∇dbl2−→0 as v−→0. (19.75) Assuming∇2F(u⋆), there is a constant C >0 such that vT∇2F(u⋆)v≥C/ba∇dblv/ba∇dbl2for all v∈Rn. This is a consequence of the Theorem 3.17 on the equivalence o f norms, coupled with the fact that every positive deﬁnite matrix deﬁnes a norm. By (19 .75), we can ﬁnd δ >0 such that |S(v,u⋆)|<1 2C/ba∇dblv/ba∇dbl2whenever 0 </ba∇dblv/ba∇dbl=/ba∇dblu−u⋆/ba∇dbl< δ. But then the Taylor formula (19.74) implies that, for all usatisfying the preceding in- equality, 0<1 2vT∇2F(u⋆)v+S(v,u⋆) =F(u)−F(u⋆), which implies u⋆is a strict local minimum of F(u). Q.E.D. Constrained Optimization and Lagrange Multipliers In many applications, the function to be minimized is subjec t to constraints. For instance, ﬁnding boundary minima requires constraining th e minima to the boundary of the domain. Another example would be to ﬁnd the minimum tempe rature on the surface of the earth. Assuming the earth is a perfect sphere of radius R, the temperature function T(u,v,w) is then to be minimized subject to the constraint u2+v2+w2=R2. Let us focus on ﬁnding the minimum value of an objective funct ionF(u) =F(u,v,w) when its arguments ( u,v,w) are constrained to lie on a regular surface S⊂R3. Suppose u⋆= (u⋆,v⋆,w⋆)T∈Sis a (local) minimum for the constrained objective function . Let u(t) = (u(t),v(t),w(t))T⊂Sbe any curve contained within the surface that passes through the minimum, with u(0) =u⋆. Then the scalar function g(t) =F(u(t)) must have a local minimum at t= 0, and hence, in view of (19.66), 0 =g′(0) =d dtF(u(t))/vextendsingle/vextendsingle/vextendsingle/vextendsingle t=0=∇F(u(0))·/squaresmallsolidu(0) =∇F(u⋆)·/squaresmallsolidu(0). (19.76) Thus, the gradient of the objective function at the surface m inimum must be orthogonal to the tangent vector to the curve. Since the curve was constr ained to lies entirely in S, its tangent vector/squaresmallsolidu(0) is tangent to the surface at the point u⋆. Since every tangent vector to the surface is tangent to some curve contained in the surface ,∇F(u⋆) must be orthogonal to every tangent vector, and hence point in the normal direct ion to the surface. Thus, a contrained critical point u⋆∈Sof a function on a surface is deﬁned so that ∇F(u⋆) =λn, (19.77) 12/11/12 1073 c/ci∇cleco√y∇t2012 Peter J. Olver wherendenotes the normal to the surface at the point u⋆. The scalar factor λis known as theLagrange multiplier in honor or Lagrange, one of the pioneers of constrained opti miza- tion. The value of the Lagrange multiplier is not ﬁxed a prior i., but must be determined by solving the critical point system (19.77). The same reaso ning applies to local maxima, which are also constrained critical points. The nature of a c onstrained critical point — local minimum, local maximum, local saddle point, etc. — is ﬁ xed by a constrained second derivative test; details are discussed in Exercise . Example 19.46. Our problem isto ﬁnd theminimum valueof theobjectivefunct ion F(u,v,w) =u2−2w3whenu,v,warerestrictedtotheunitsphere S= (u2+v2+w2= 1). Theradial vector n= (u,v,w)Tisnormalto thesphere, andso thecriticalpoint condition (19.77) is ∇F= 2u 0 −6w2 =λ u v w . Thus, we must solve the system of equations 2u=λu,0 =λv,−6w2=λw, subject to u2+v2+w2= 1, for the unknowns u,v,wandλ. This needs to be done carefully to avoid missing any cases. First, ifu/ne}ationslash= 0, then λ= 2,v= 0, and either w= 0 whence u=±1, orw=−1 3and so u=±√ 1−w2=±2√ 2 3. On the other hand, if u= 0, then either λ= 0,w= 0 and so v=±1, orv= 0,w=±1, andλ=∓6. Collecting these together, we discover that there are a total of 8 critical points on the unit sphere: u⋆ 1= 1 0 0 ,u⋆ 2= −1 0 0 ,u⋆ 3= 2√ 2 3 0 −1 3 ,u⋆ 4= −2√ 2 3 0 −1 3 , u⋆ 5= 0 1 0 ,u⋆ 6= 0 −1 0 ,u⋆ 7= 0 0 1 ,u⋆ 8= 0 0 −1 . Since the unit sphere is closed and bounded, we are assured th atFhas a global maximum and a global minimum when restricted to S, which are both to be found among our candidate critical points. Thus, we merely compute the valu e of the objective function at each critical point, F(u⋆ 1) = 1, F(u⋆ 2) = 1, F(u⋆ 3) =22 27, F(u⋆ 4) =26 27, F(u⋆ 5) = 0, F(u⋆ 6) = 0, F(u⋆ 7) =−2, F(u⋆ 8) = 2. Therefore, u⋆ 7must be the global minimum and u⋆ 8the global maximum of the objective function restricted to the unit sphere. The status of the oth er six critical points — con- strainned local maximum, minimum, or neither — is less evide nt, and a full classiﬁcation requires the second derivative test outlined in Exercise . 12/11/12 1074 c/ci∇cleco√y∇t2012 Peter J. Olver If the surface is given as the level set of a function G(u,v,w) =c, (19.78) then at any point u⋆∈S, the gradient vector ∇G(u⋆) points in the normal direction to the surface, and hence, provided n=∇G(u⋆)/ne}ationslash=0, the surface critical point condition can be rewritten as ∇F(u⋆) =λ∇G(u⋆), (19.79) or, in full detail, the critical point ( u⋆,v⋆,w⋆)Tmust satisfy ∂F ∂u(u,v,w) =λ∂G ∂u(u,v,w), ∂F ∂v(u,v,w) =λ∂G ∂v(u,v,w), ∂F ∂w(u,v,w) =λ∂G ∂w(u,v,w).(19.80) Thus, to ﬁnd the contrained critical points, one needs to sol ve the combined system (19.78, 80) of 4 equations for the four unknowns u,v,wand the Lagrange multiplier λ. Formally, one can reformulate the problem as an unconstrain ed optimization problem by introducing the augmented objective function E(u,v,w,λ )=F(u,v,w)−λ/parenleftbig G(u,v,w)−c/parenrightbig . (19.81) Thecriticalpointsoftheaugmentedfunction arewhereitsg radient, withrespect toallfour arguments, vanishes. Setting the partial derivatives with respect to u,v,wto 0 reproduces the system (19.80), while it partial derivative with respec t toλreporduces the constraint (19.78). IfF(u) is deﬁned on a closed subdomain Ω ⊂Rn, then its minima may also occur at boundary points u∈∂Ω. When the boundary is smooth, there is an analogous critica l point condition for local boundary extrema. Theorem 19.47. LetΩ⊂Rnbe a domain with smooth boundary ∂Ω. Suppose F(u)is co;ntinuously diﬀerentiable at all points in Ω = Ω∪∂Ω. If the boundary point u0∈∂Ωis a(local)minimum for Fwhen restricted to the closed domain Ω, then the gradient vector∇F(u0)is either 0or points inside the domain in the normal direction to ∂Ω. See Figure bmin for a sketch of the geometrical conﬁguration. Proof: Letu(t)⊂∂Ω be any curve that is entirely contained in the boundary, wit h u(0) =u0. Then the scalar function g(t) =F(u(t)) must have a local minimum at t= 0, and hence, in view of (19.66), 0 =g′(0) =d dtF(u(t))/vextendsingle/vextendsingle/vextendsingle/vextendsingle t=0=/angbracketleftbig ∇F(u(0));/squaresmallsolidu(0)/angbracketrightbig . Since the curve lies entirely in ∂Ω, its tangent vector/squaresmallsolidu(0) is tangent to the boundary at the point u0; moreover, we can realize any such tangent vector by a suitab le choice of curve. We conclude that ∇F(u0) is orthogonal to every tangent vector to ∂Ω at the point 12/11/12 1075 c/ci∇cleco√y∇t2012 Peter J. Olver u0, and hence must point in the normal direction to the boundary . Moreover, if non-zero, it cannot point outside Ω since then −∇F(u0), which is the direction of decrease in the objective function, would point inside the domain, which wo uld preclude u0from being a local minimum. Q.E.D. The same ideas can be applied to optimization problems invol ving functions of several varaibles subject to one or more constraints. Suppose the ob jective function F(u) = F(u1,...,un) is subject to the constraints G1(u) =c1,···Gk(u) =ck. (19.82) A pointusatisfying the constraints is termed regularif the corresponding gradient vectors ∇G1(u),...,∇Gk(u) are linearly independent. (Irregular points are more tric ky, and must be handled separately.) A regular constrained critical poi nt necessarily satisﬁes the vector equation ∇F(u) =λ1∇G1(u)+···+λ1∇G1(u), (19.83) where the unspeciﬁed scalars λ1,...,λkare called the Lagrange multipliers for the con- strained optimization problem. The critical points are thu s found by solving the combined system (19.82–83) for the n+kvaraibles u1,...,unandλ1,...,λk. As in (19.81) we can reformulate this as an uncontrained optimization probl em for the augmented objective function E(u,λ) =F(u)−k/summationdisplay i=1λi/parenleftbig Gi(u)−ci/parenrightbig . (19.84) The gradient with respect to ureproduces the critical point system (19.83), while its gradient with respect to to λ= (λ1,...,λk) recovers the constraints (19.82). Theorem 19.48. Every regular constrained local minimum and maximum is a con - strained critical point. Example 19.49. The problem is to ﬁnd the point or points on the intersection o f the elliptical cylinders u2+4v2= 1, u2+9w2= 4, (19.85) that is the closest to the origin. Thus, we seek to minimize th e squared†distance function F(u,v,w) =u2+v2+w2 subject to the constraints G(u,v,w) =u2+4v2= 1, H (u,v,w) =u2+9w2= 4, The augmented objective function (19.81) is E(u,v,w,λ,µ )=u2+v2+w2−λ/parenleftbig u2+4v2−1/parenrightbig +µ/parenleftbig u2+9w2−4/parenrightbig . †Any distance minimizer also minimizes the squared distance; we w ork with the latter in order to avoid square roots in the computation. 12/11/12 1076 c/ci∇cleco√y∇t2012 Peter J. Olver To ﬁnd its critical points, we set all its partial derivative s to zero: ∂E ∂u= 2u+2λu+2µu= 0,∂E ∂v= 2v+8λv= 0,∂E ∂w= 2w+18λw= 0, while the partial derivatives with respect to the Lagrange m ultipliers λ,µreproduce the two constraints (19.85). Thus, eitheru= 0 orλ+µ=−1,eitherv= 0 orλ=−1 4,and either w= 0 orµ=−1 18. Thus, at least one of u,v,wmust be zero. If u= 0, then v=±1 2, w=±2 3; ifv= 0, thenu=±1, w=±1√ 3; while there are no real solutions to the constraints when w= 0. The ﬁrst four critical points,/parenleftbig 0,±1 2,±2 3/parenrightbigT, all lie a distance5 6≈.8333 from the origin, while the second four,/parenleftbig ±1,0,±2 3/parenrightbigT, are further away, at distance√ 2 3≈1.1547. Thus, the closest points on intersection of the cylinders are the ﬁ rst four, while the furthest points from the origin are the last four. (The latter comment relies on the fact that the intersection is a bounded subset of R3.) Remark: A second derivative test for constrained minima and maxima can be found in [129]. Numerical Minimization of Scalar Functions In practical optimization, one typically bypasses the prel iminary characterization of minima as critical points, and instead implements a direct i terative procedure that con- structs a sequence of successively better approximations t o the desired minimum. As the computation progresses, the approximations are adjusted s o that the objective function is made smaller and smaller, which, we hope, will ensure that we are converging to some form of minimum. As always, to understand the issues involved, it is essentia l to consider the simplest scalar situation. Thus, we are given the problem of minimizi ng a scalar function F(u) on a bounded interval a≤u≤b. The minimum value can either be at an endpoint or an interior minimum. Let us ﬁrst state a result that plays a simi lar role to the Intermediate Value Lemma 19.18 that formed the basis of the Bisection Meth od for locating roots. Lemma 19.50. Suppose that F(u)is deﬁned and continuous for all a≤u≤b. Suppose that we can ﬁnd a point a < c < b such that F(c)< F(a)andF(c)< F(b). Then F(u)has a minimum at some point a < u⋆< b. The proof is an easy consequence of Theorem 19.35. Therefore , if we ﬁnd three points a < c < b satisfying the conditions of the lemma, we are assured of the existence of a local minimum for the function between the two endpoints. Once thi s is done, we can design an algorithm to home in on the minimum u⋆. We choose another point, say dbetween aandcand evaluate F(d). IfF(d)< F(c), thenF(d)< F(a) also, and so the points a < d < c satisfy the hypotheses of Lemma 19.50. Otherwise, if F(d)> F(c) then the pointsd < c < b satisfy the hypotheses of the lemma. In either case, a local m inimum has been narrowed down to a smaller interval, either [ a,c] or [d,b]. In the unlikely even 12/11/12 1077 c/ci∇cleco√y∇t2012 Peter J. Olver thatF(d) =F(c), one can try another point instead — unless the objective fu nction is constant, one will eventually ﬁnd a suitable value of d. Iterating the method will produce a sequence of progressively smaller and smaller intervals i n which the minimum is trapped, and, just like the Bisection Method, the endpoints of the int ervals get closer and closer to u⋆. The one question is how to choose the point d. We described the algorithm when it was selected to lie between aandc, but one could equally well try a point between candb. To speed up the algorithm, it makes sense to place din the larger of the two subintervals [ a,c] and [c,b]. One could try placing din the midpoint of the interval, but a more inspired choice is to place it at a fraction θ=5√ 2−1 2≈.61803 of the way along the interval, i.e., at θa+ (1−θ)cif[a,c] is the longer interval. (One could equally well take the point (1−θ)a+θc.) The result is the Golden Section Method . At each stage, the length of the interval has been reduced by a factor of θ, so the convergence rate is linear, although a bit slower than bisection. Another strategy is to use an interpolating polynomial pass ing through the three points on the graph of F(u) and use its minimum value as the next approximation to the minimum. According to Exercise 4.4.24, the minimizing valu e occurs at d=ms−nt s−t, where s=F(c)−F(a) c−a, t =F(b)−F(c) b−c, m =a+c 2, n =c+b 2. As long as a < c < b satisfy the hypothesis of Lemma 19.50, we are assured that th e quadratic interpolant has a minimum (and not a maximum!), an d that the minimum remains between the endpoints of the interval: a < d < b . If the length of the interval is small, the minimum value should be a good approximation to th e minimizer u⋆ofF(u) itself. Once dis determined, the algorithm proceeds as before. In this cas e, convergence is not quite guaranteed, or, in unfavorable situations, could be much slower than the Golden Section Method. One can even try using the method when the fun ction values do not satisfy the hypothesis of Lemma 19.50, although now the new p ointdwill not necessarily lie between aandb. Worse, the quadratic interpolant may have a maximum at d, and one ends up going in the wrong direction, which can even happen in the minimizing case due to the discrepancy between it and the objective function F(u). Thus, this method must be handled with care. A ﬁnal idea is to focus not on the objective function F(u) but rather its derivative f(u) =F′(u). The critical points of Fare the roots of f(u) = 0, and so one can use one of the solution methods, e.g., Bisection or Newton’s Method , to ﬁnd the critical points. Of course, one must then take care that the critical point u⋆is indeed a minimum, as it could equally well be a maximum of the original objective fun ction. (It will probably not be an inﬂection point, as these do not correspond to simple ro ots off(u).) The status of the critical point can be checked by looking at the sign of F′′(u⋆) =f′(u⋆); indeed, if we use Newton’s Method we will be computing the derivative at ea ch stage of the algorithm, and can stop looking if the derivative turns out to be of the wr ong sign. 12/11/12 1078 c/ci∇cleco√y∇t2012 Peter J. Olver Gradient Descent Now, let us turn our attention to multi-dimensional optimiz ation problems. We are seeking to minimize a (smooth) scalar objective function F(u) =F(u1,...,un). According to Theorem 19.38, at any given point uin the domain of deﬁnition of F, the negative gradient vector−∇F(u), if nonzero, points in the direction of the steepest decrea se in F. Thus, to minimize F, an evident strategy is to “walk downhill”, and, to be eﬃcien t, walk downhill as fast as possible, namely in the direction −∇F(u). After walking in this direction for a little while, we recompute the gradient, and this tells us the new direction to head downhill. With luck, we will eventually end up at the b ottom of the valley, i.e., at a (local) minimum value of the objective function. This simple idea forms the basis of the Gradient Descent Method for minimizing the objective function F(u). In a numerical implementation, we start the iterative pro cedure with an initial guess u(0), and let u(k)denote the kthapproximation to the minimum u⋆. To compute the next approximation, we set out from u(k)in the direction of the negative gradient, and set u(k+1)=u(k)−tk∇F(u(k)). (19.86) for some positive scalar tk>0. We are free to adjust tkso as to optimize our descent path, and this is the key to the success of the method. If∇F(u(k))/ne}ationslash=0, then, at least when tk>0 is suﬃciently small, F(u(k+1))< F(u(k)), (19.87) and sou(k+1)is, presumably, a better approximation to the desired minim um. Clearly, we cannot choose tktoo large or we run the risk of overshooting the minimum and re versing the inequality (19.87). Think of walking downhill in the Swi ss Alps. If you walk too far in a straight line, which is what happens as tkincreases, then you might very well miss the valley and end up higher than you began — not a good strateg y for descending to the bottom! On the other hand, if we choose tktoo small, taking very tiny steps, then the method may end up converging to the minimum much too slowly to be of practical use. How should we choose an optimal value for the factor tk? Keep in mind that the goal is to minimize F(u). Thus, a good strategy would be to set tkequal to the value of t >0 that minimizes the scalar objective function g(t) =F/parenleftbig u(k)−t∇F(u(k))/parenrightbig (19.88) obtained by restricting F(u) to the ray emanating from u(k)that lies in the negative gradient direction. Physically, this corresponds to setti ng oﬀ in a straight line in the direction of steepest decrease, and continuing on until we c annot go down any further. Barring luck, we will not have reached the actual bottom of th e valley, but must then readjust our direction and continue on down the hill in a seri es of straight line paths. Inpractice, onecanrarelycomputetheminimizingvalue t⋆of(19.88)exactly. Instead, we employ one of the scalar minimization algorithms present ed in the previous subsection. Note that we only need to look for a minimum among positive val ues oft >0, since our choice of the negative gradient direction assures us that, a t least for tsuﬃciently small and positive, g(t)< g(0). 12/11/12 1079 c/ci∇cleco√y∇t2012 Peter J. Olver A more sophisticated approach is to employ the second order T aylor polynomial to approximate the function and then use its minimum (assuming such exists) as the next approximation to the desired minimizer. Speciﬁcally, if u(k)is the current approximation to the minimum, then we approximate F(u)≈c(k)+(u−u(k))Tg(k)+1 2(u−u(k))TH(k)(u−u(k))(19.89) nearu(k), where c(k)=F(u(k)),g(k)=∇F(u(k)), H(k)=∇2F(u(k)), (19.90) are, respectively, the function value, the gradient and the Hessian at the current iterate. If u⋆is a stric local minimum, then ∇2F(u⋆) is positive deﬁnite, and hence, assuming u(k)is close, so is H(k)=∇2F(u(k)). Thus, by Theorem 4.1, the quadratic Taylor approximation has a unique minimum value u(k+1)which satisﬁes the linear system H(k)(u(k+1)−u(k)) =−g(k). (19.91) The solution serves to deﬁne the next approximation u(k+1). While not guaranteed to converge, the method does perform well in all reasonable sit uations. 12/11/12 1080 c/ci∇cleco√y∇t2012 Peter J. Olver