The Method of Lagrange Multipliers

1 The Method of Lagrange Multipliers Phil Lucht Rimrock Digital Technology, Salt Lake City, Utah 84103 last update: Oct 25, 2016 Maple code is available upon request. Co mments and errata are welcome. The material in this document is copyrighted by the author. The graphics look ratty in Windows Adobe PDF viewers wh en not scaled up, but look just fine in this excellent freeware viewer: https://www.tracker-software.co m /product/pdf-xchange-editor . The table of contents has live links. Most PDF view ers provide these links as bookmarks on the left. Overview an d Summary........................................................................................................... .............. 2 1. Finding the stationary points of a function subject to constraints ................................................. 4 2. The Method of Lagrange Multipliers.......................................................................................... ....11 3. Three simple examples...................................................................................................................... 15 3.1 Example 1: 3-dimensional hemisp here with one simple constraint............................................. 15 3.2 Example 2: 4-dimensional hemisp here with two simple constraints ........................................... 20 3.3 Example 3: N-dimensional hemisphere of radius R in EN with C simple constraints.................. 23 3.4 Example 1 Revisited: What is the half-width of the peak?.......................................................... 24 4. Example 4: Extremal distances between an ellipse and a circle ................................................... 27 5. Example 5: The Boltzmann factor in Statistical Mechanics ......................................................... 34 5.1 Statement of The Boltzmann Extremum Problem ........................................................................ 34 5.2 Solution of The Boltzmann Extremum Problem........................................................................... 36 5.3 More details of the solution ............................................................................................... ........... 40 5.4 Example: N particles in 3 states : a numerical example................................................................ 46 5.5 Example: The Maxwell-Boltzmann Distribution.......................................................................... 55 Section 6. Matrix proof of the Method of Lagrange Multipliers ...................................................... 68 6.1. The R matrix and its Rank ................................................................................................. .......... 68 6.2. Proof of Theorem 3 ( ⇒) .............................................................................................................. 71 (a) Proof for N = 6 and S = 4 .................................................................................................. ........ 71 (b) Proof for general S ≤ N ............................................................................................................. 75 Appendix A: Matrix rank equals the number of independent columns and rows.......................... 79 Appendix B: Determinants....................................................................................................... ............ 83 B.1 Definition of the determinant .............................................................................................. ......... 83 B.2 The permutation group, the permutation tensor ε, and expansions for det(M) ............................ 85 B.3 Minors, Cofactors, and the Cofactor Expansions of det(M) ........................................................ 91 B.4 Expressions for the inverse matrix M-1 and Cramer's Rule ......................................................... 93 Appendix C: The Intersection of Constraint Surfaces ...................................................................... 96 Appendix D: Lagrange Multipli ers in Classical Mechanics............................................................ 101 References............................................................................................................................................ 116 2 Overview and Summary The Method of Lagrange Multipliers is used to determ ine the stationary points (including extrema) of a real function f( r) subject to some number of (holonomic) constraints. The main purpose of this document is to provide a solid derivation of the method and thus to show why the method works. A geometric derivation is presen ted in Section 1 while an alternate matrix-based derivation appears in Section 6. A secondary purpose is to provide some interesting and illustrative examples . Examples 1, 2 and 3 are geometric examples in 3, 4 and N dimensions. Example 4 concerns extremal distances between a circle and an ellipse. Extended Example 5 involves th e Boltzmann distribution of statistical mechanics where the variables N i are state population counts. Finally, A ppendix D gives several examples from classical Lagrangian dynamics (with all supporting material) in which a functional F( φ) is rendered stationary rather than a function like f( r). In these examples the variables are functions φi known as generalized coordinates. Here is a concise summary of the document: Section 1 defines a stationary point r of a function f( r) subject to constraints. It then proves Theorem 1 which says that iff r is such a stationary point, then constants λ i must exist so that a certain gradient statement is true. The λi are the "Lagrange Multipliers". Section 2 presents the Method of Lagrange Multipliers as Theorem 2. This theorem is really just a recasting of Theorem 1, so Theorem 1 ⇔ Theorem 2 and the derivation of Theorem 2 is trivial. Section 3 gives three detailed examples involving hemisphe res in 3,4 and N dimensions, with 1, 2 and S constraints respectively. Section 4 does a numerical example to find the extremal distances between a circle and an ellipse which are coplanar. Section 5 presents a more sophisticated physics exampl e associated with the name Boltzmann wherein one determines the particle populations of energy le vels in the context of statistical mechanics and thermal equilibrium. This self-contained section first tr eats the case of discrete energy levels with a numerical example involving three levels. It then examines a case with continuous energy levels and derives the well-known Maxwell-Boltzmann distribution. Numbers are obtained for a small box of helium atoms. It is shown for example that such atoms at room temperature have a mean speed of 2,780 mph. Section 6 derives Theorem 1 in a manner totally different from the derivation presented in Section 1. The function f( r) and the constraint functions are treated as variables of a certain transformation which has a "differential matrix" called R( r). Theorem 3 is then proven, claiming that iff r is a stationary point, then the rank of the matrix R( r) falls below full rank and this in turn implies Theorem 1. 3 Appendix A proves for a general n x m matrix that column rank = row rank = rank. This fact is used in the proof of Theorem 3 in Section 6. Appendix B derives a batch of facts about determinants some of which are used in Appendix A. This is a standalone section that might be a useful deriva tion resource for workers who deal with matrices. Appendix C gives examples of intersecting constraint surfaces and discusses the notion of eliminating variables using constraint equations, as is done in the proof Section 6. It concludes with an explanation of why a constraint equation a( r) = 0 with r = (x 1,x2.....xN) can be interpreted as a surface in EN. Appendix D presents two applications of Lagrange Multip liers in classical Lagrangian mechanics. References (not many) are provided on the last page. The general tone is more that of an engineer or ph ysicist, not that of an abstract mathematician. Maple code is used where it seems useful to implement calculations or display graphs. When an earlier equation is quoted, its equation number is put in italics. We often use the abbreviations ∂ if ≡ ∂f/∂xi ≡ ∂f ∂xi . Sometimes we also use f i ≡ ∂f ∂xi . 4 1. Finding the stationary points of a function subject to constraints Here we develop a m ethod for finding the stationary points of a function f( r) subject to a set of constraint conditions a i(r) = 0. This method does involve "Lagrange mu ltipliers", but the Method of Lagrange Multipliers is usually stated in the equivalent form shown in Section 2. Nevertheless, the proof of the method is presented below in this section. In the following r = (x 1,x2.....xN) is a point in Euclidian space EN. Preliminary Facts about Surfaces Fact 1: The equation F(x 1,x2.....xN) = K defines an N-1 dimensional surface in EN . (1.1) For example, F(x,y,z) = x 2+y2+z2 with F = R2 describes a spherical 2D surface of radius R in E3. F(x,y) = x2+y2 with F= R2 describes a 1D surface (curve) in E2, a circle . See the discussion starting at (C.16) in Appendix C for more details about Fact 1. Fact 2: The vector ∇ F(r) is always normal to the surface F( r) = K at point r . (1.2) Proof : Start at point r on the surface F( r) = K and move a small distance d r in an arbitrary direction along the surface. Since one stays on the surface, F( r+dr) = K. Then dF = F( r+dr) - F(r ) = 0. But one knows that dF = ∇F • dr so, in order that dF = 0 for an arbitrary dr displacement on the surface, ∇ F(r) must be locally normal to the surface at point r . Definition: As K takes a set of different values K 1, K2.....KM, the equation F(x 1,x2.....xn) = K describes a family of so-called level surfaces (level sets) . If M is large and the range of the K i is small, these surfaces will be closely spaced. Adjace nt surfaces then have nearly the same shape, although in general all the surfaces in the family of surfaces do not have the same shape. (1.3) Example : F(x,y) = x2+2y4 and F= K for K = 1 to 10. Here "level surfaces" means "level curves". 5 (1.4) Fact 3: The gradient ∇F points in the direction in which F changes most rapidly. (1.5) Proof : For tiny displacement d r, dF = ∇F • dr is largest positive when d r points in the direction of ∇ F ("uphill"). And dF = ∇ F • dr is largest negative when d r points opposite the direction of ∇F ("downhill"). For the above example ∇F = 2x x^ + 8y3y^ we can superpose this gradient field on the above level curves using a Maple fieldplot routine, (1.6) Wherever a gradient arrow base lies on a level curve, the arrow is normal to the curve and points in the direction in which the function x2+2y increases most rapidly. Thinking of z = f(x,y) = x2+2y4, the blue lines above are the level curves of an oblong bowl in to which the viewer is looking, and the gradient arrows all point uphill on the bowl surface. The arro ws are longer where the level curves lines (contour lines, topo lines) are closer together. 6 The Intersection Constraint Surface A(r ) = 0 Assume that a point r = (x1, x2 ... xN) in EN is required to satisfy C constraint equations, a i(r) = 0 i = 1,2...C . // constraints (1.7) According to Fact 1 above, each constraint equation defines a surface of dimension N-1 in E N. Our point r must lie simultaneously on all C constraint surfaces. The intersection of all these constraint surfaces is in fact an intersection constraint surface of dimension N-C in E N whose constraint function we shall call A(r). In symbolic form we write ( A ( r) = 0) = (a 1(r) = 0) ∩ (a2(r) = 0) ∩ .....∩ (aC(r) = 0) . (1.8) dimension : N-C N-1 N-1 N-1 Barring degenerate cases (which we shall ignore in th e following general discussion) each extra constraint lowers the dimensionality of the intersecting constraint surface by 1 degree of freedom. We assume that the intersection constraint surface is not null, otherwi se our problem of interest would have no solutions. The intersection surface could in general have multiple disjoint pieces. If there are no constraints, all points in E N are "allowed". If there is one constraint a 1(r) = 0, then A( r) = a1(r). For two constraint surfaces a 1(r) = 0 and a 2(r) = 0, one must compute the intersecting constraint surface (A( r) = 0) which we symbolically write as (a 1(r) = 0) ∩ (a2(r) = 0). Example: N = 3 and C = 2 constraints. Suppose in E3 there are two constraints as follows : a 1(x1,x2,x3) = x12 + x22 + x32 - 22 = 0 a2(x1,x2,x3) = (x1-2)2 + x22 + x32 - 22 = 0 . ( 1 . 9 ) Each of these constraints describes a surface of dimension 2 in E 3 (a sphere of radius 2). These two spheres intersect in a circle x 22 + x32 = 3, so the intersection surface has dimension 1 in E3. This agrees with N-C = 1. Here is a crude drawing, where the inte rsection constraint surface (edge-on circle) is red, (1.10) 7 In this example one has A( r) = x22 + x32 - 3 and then A(r ) = 0 is the equation of the intersection constraint surface. If the origin of the second sphere were moved to x 1= 4 (slide right sphere to the right 2 units), one has a degenerate case where the intersec tion constraint surface is just a point and so has dimension 0. Further separation of the spheres results in a null intersection surface. Setting x 1 = 0 aligns the two spheres for another degenerate case, resulti ng in the intersection surface having dimension 2. Definitions: The tangent space and the perp space at r Consider a point r located on (A( r)=0). As noted above, this intersection constraint surface (A( r)=0) is a surface of dimension N-C embedded in E N. Therefore at r one can construct a se t of N-C independent vectors d r such that r+dr also lies on the intersection constraint surface, so both A( r) = 0 and A(r+dr ) = 0. These N-C independent vectors span a vector space (dimension N-C) at point r on the surface known as the tangent space at point r. These vectors d r are all "tangent to" the surface. (1.11) There is another vector space within E N at point r of dimension N - (N-C) = C which is orthogonal to the tangent space, and we shall call it the perp space at r. If d r is any vector in the tangent space at r, and N is a vector in the perp space at r, then N • dr = 0. So one can find C linearly independent vectors N such that N • dr = 0 and all these N vectors are normal to surface (A( r)=0) at r . (1.12) As the point r moves on the surface (A( r)=0), the axes of both the tangent space and the perp space also move, perhaps rotating a bit. The union of the set of N-C basis vectors of the tangent space at r and the set of C basis vectors of the perp space at r forms a basis for E N. Summarizing the above: Fact 4 : Let M ≡ (A(r)=0) embedded in E N represent the intersection constraint surface of C constraint surfaces (a i(r) = 0). The dimension of surface M, and the dimension of the tangent space at any r on M, is N-C, while the dimension of the perp space at any r on M is C. (1.13) Fact 5: If A • dr = 0 for all d r in the tangent space at r , then A lies in the perp space at r. (1.14) Usually in such discussions one regards the surf aces involved as having smoothness properties which make them be manifolds . If the intersecting constraint surface is a manifold M, then the tangent space at r on M is often denoted as T rM and the perp space as (T rM)⊥. For example , in Fig 1.10 for E3 we started with two constraint s pheres and the intersection constraint surface is the edge-on circle shown in red. Since this is a surface of dimension 1, the tangent space has dimension 1 and the perp space has dimension 2. At any point r on a circle embedded in E3, there is only one ±d r direction which keeps one on the circle, but one could construct two lin early independent normal vectors at r . On the other hand, if either sphere were the only constraint in the problem, the tangent space would have dimension 2 and the perp space dimensi on 1. At a point on the sphere there is only one normal vector direction in E3. 8 Fact 6 : If d r lies in the tangent space of a point r on the intersection constraint surface, then ∇ai(r) • dr = 0 and so the vectors ∇ai(r) lie in the perp space at r. ( 1 . 1 5 ) Proof : Let r be a point on (A( r)=0) and let d r be a vector in the tangent space at r on (A( r)=0). This means that both r and r + dr lie on the intersection constraint surface (A( r)=0), so dA = A( r+dr) = A( r) = 0. Since this surface is the intersection of all the (a i(r)=0) component surfaces, it is also true that da i= ai(r+dr) = ai(r) = 0 for each surface i. But da i = ∇ai(r) • dr so we have 0 = ∇ai(r) • dr. Therefore, if dr is any vector in the tangent space, then ∇ai(r) • dr = 0. From Fact 5 one concludes that ∇ai(r) must lie in the perp space at r on (A( r)=0). Comment: Avoiding gradient confusion Let r be a point in EN. Consider the equation u = f( r) which maps EN to E1. Construct a space EN+1 having coordinates ( r,u). Let g( r,u) = f( r) - u be a function on this new EN+1. The equation u = f( r) is the same as the equation g( r,u) = 0. The equation g( r,u) = 0 is, according to Fact 1, a surface of dimension N in EN+1. From Fact 2, a normal to this surface is given by ∇(N+1)g(r,u) which is a vector in EN+1. For any reasonable smooth surface, this normal vector cannot be null. On the other hand, the vector ∇(N)f(r) is a vector in EN and there is no reason why this vect or cannot be null. We write this as ∇f(r) below. If this paragraph seems confusing, reread it thinking of N = 2 so u = f( r) = f(x,y) is a regular surface in E3. The Stationary Point Problem STATIONARY POINT (NO CONSTRAINTS) We know that if we move from r to r+dr, where d r is any dr, the differential change in f( r) is given by df(d r) = f( r+dr) - f(r ) = ∇ f(r) • dr, dr = any differential vector in E N . (1.16) If for some r it happens that ∇ f(r) = 0, then the above says df = 0 and so r is a stationary point (also known as a critical point ) of the function f(r ). That is to say, the phrase " r is a stationary point of f( r)" means that df(d r) = 0 for any d r. In this problem there are no constraints specifying legal values for r. All points r in EN are legal. For r in E2 the reader is no doubt familiar with the fact that, for u = f( r) plotted in E3, such a stationary point could occur at a hill top, or a valley bottom, or a saddle point (which is neither a maximum or a minimum). For fun, go look up "monkey saddle". STATIONARY POINT (WITH CONSTRAINTS) Let r lie on the intersection constraint surface (A( r)=0). Now the only "allowed" d r are vectors in the tangent space for this r, because any other d r will violate at least one constraint since it will take us off the intersecting constraint surface. We modify the previous equation to say df(d r) = f( r+dr) - f(r ) = ∇ f(r) • dr , dr = any diff. vector lying in the tangent space at r (1.17) 9 If one could find a value of r such that ∇ f(r) = 0 and A(r) = 0, one would have df = 0 and such an r would be a stationary point of "f( r) subject to A(r) = 0", but this situation is generally very unlikely. Instead, we define a stationary point as follows: Definition : If r lies on A(r ) = 0, then r is a " stationary point of f( r) subject to A(r ) = 0" if df(d r) = 0 for any d r in the tangent space of (A( r)=0) at point r. In simple terms, r is a stationary point if df(dr ) = 0 for all "constraint-allowed" dr displacements. (1.18) Since df(d r) = ∇ f(r) • dr, this requirement will certainly be me t in the special case just noted that ∇ f(r) = 0 and A(r) = 0, but the above definition is more general and there will be stationary points for which ∇f(r) ≠ 0. Fact 7: If r is a stationary point of f( r), then ∇f(r) must lie in the perp space of (A( r)=0) at r. (1.19) Proof : If r is a stationary point, then by definition df = ∇f(r) • dr = 0 for any d r in the tangent space at r. Then by Fact 5 ∇f(r) must lie in the perp space at r. In the special case ∇ f(r) = 0 and A (r) = 0 we have noted that r is a stationary point, but we now look for stationary points for which ∇f(r) ≠ 0. Consider the set of C+1 vectors { ∇a i(r), ∇f(r)}. From Facts 6 and 7 we know that all C+1 vectors must lie in the perp space at r. According to Fact 4, the dimensionality of the perp space at r is C. Therefore, the vectors { ∇ai(r), ∇f(r)} for i = 1 to C must be linearly dependent (see (6.1.6) below) -- one cannot have C+1 independent vectors in a space of dimension C. This means that one can find a set of constants λ'i (not all 0) for i = 1 to C such that λ'0∇f(r) + Σi=1C λ'i ∇ai(r) = 0 . ( 1 . 2 0 ) If λ'0 = 0, then Σi=1C λ'i ∇ai(r) = 0 which means the ∇ai(r) are linearly dependent at point r . We explicitly rule out this case by requiring that the a i(r) are such that the C vectors ∇ai(r) are linearly independent near any point r where we might have a stationary point. So assuming then that λ'0 ≠ 0, define λi ≡ λ'i/λ'0 to get this new version of (1.20), ∇f(r) + Σi=1C λi ∇ai(r) = 0 . ( 1 . 2 1 ) We have just proven the following theorem: If a point r is a stationary point for f( r) subject to constraints a i(r) = 0, then the following are true : 1. ai(r) = 0 for i = 1,2...C (point r must satisfy all the constraints ) 2. There must exist C constants λi such that ∇f(r) + Σi=1C λi ∇ai(r) = 0 . (1.22) 10 Conversely, suppose we can find a set of λ i such that the above two conditions are met for some r. Let dr be a vector in the tangent space at r. We then compute df = ∇f • dr = [- Σi=1C λi ∇ai(r)] • dr = - Σi=1C λi [ ∇ai(r) • dr ] . (1.23) But we know from Fact 6 that ∇ai(r) • dr = 0 for d r in the tangent space, and therefore df = 0 for any such d r, and therefore r is a stationary point. Thus we improve the above Theorem to get Theorem 1 : A point r is a "stationary point for f( r) subject to constraints a i(r) = 0" ⇔ (1.24) (a) a i(r) = 0 for i = 1,2...C (point r must satisfy all the constraints ) (b) There exist C constants λi such that ∇f(r) + Σi=1C λi ∇ai(r) = 0 . Corollary : If r satisfies all the constraints [ r lies on (A( r)=0)] but r is not a stationary point of f( r), then there exists no equation of the form ∇f(r) + Σi=1C λi ∇ai(r) = 0. That is to say, there is no set of real numbers λi which makes this equation be true. In this case the vectors { ∇ai(r), ∇f(r)} are linearly independent. (1.25) The constants λ i are referred to as Lagrange Multipliers . 11 2. The Method of Lagrange Multipliers The object is to solve this problem: Find the stationary points r of function f( r) subject to C constraints a i(r) = 0 . (2.1) Define the following function of N+C variables r = (x 1, x2....xN) and λ = (λ1,λ2....λC), H( r,λ) ≡ f(r) + Σi=1Cλiai(r) = f(r) + λ • a(r) . (2.2) We then have the following theorem concerning this function H : Theorem 2 : Method of Lagrange Multipliers (2.3) The stationary points r in EN of the function f(r ) subject to C constraints a i(r) = 0 are the same as the r values obtained from finding stationary points ( r,λ) in EN+C of the unconstrained function H(r ,λ) . Proof: Compute the unconstrained stationary points ( r,λ) of H by setting all N+C partial derivatives to 0, 0 = ∂H/∂λ i = ai(r) for coordinates λi i = 1..C 0 = ∂ H/∂xk = ∂f(r)/∂xk+ Σi=1C λi∂ai(r)/∂xk for coordinates x k k = 1..N . (2.4) Rewrite as 1. ai(r) = 0 f o r i = 1 , 2 . . . C 2. ∇f(r) + Σi=1C λi ∇ai(r) = 0 . // and therefore this equations exists (2.5) But we have already proven in Theorem 1 (1.24) th at the conditions (2.5) are necessary and sufficient for r to be a stationary point of function f( r) subject to C constraints a i(r) = 0, so that then concludes the proof of Theorem 2. We now write out these equations in more detail: a 1(r) = 0 a2(r) = 0 .... a C(r) = 0 / / C e q u a t i o n s ( 2 . 6 ) ∂ 1f(r) + λ1∂1a1(r) + λ2∂1a2(r) + ..... + λC ∂1aC(r) = 0 ∂2f(r) + λ1∂2a1(r) + λ2∂2a2(r) + ..... + λC ∂2aC(r) = 0 .... ∂Nf(r) + λ1∂Na1(r) + λ2∂Na2(r) + ..... + λC ∂NaC(r) = 0 . // N equations (2.7) 12 Since all the functions in these equations are known (since f( r) and all the constraints are known), it is in principle possible to solve the equations for r and the λ i. Then those solution stationary points r need to be studied more to see if they really solve the problem of interest (max, min, saddle, etc). In general the various functions of r are non-linear and not just polynomials, and there are likely to be multiple candidate solutions. A numeric solution may be requ ired. The Method of Lagrange Multipliers provides no silver bullet for solving these equations. Most examples one sees in texts and on the web have very simple functions allowing a straightforward analytic solution. Since the λi appear linearly in (2.7), it is possible to obtain expressions for the λi as follows. First, write out the equations (2.7) in matrix form, - ⎝⎜⎛ ⎠⎟⎞∂1f ∂2f ... ∂Nf = ⎝⎜⎛ ⎠⎟⎞∂1a1 ∂1a2 ... ∂1aC ∂2a1 ∂2a2 ... ∂2aC ... ... ... ... ∂Na1 ∂Na2 ... ∂NaC ⎝⎜⎛ ⎠⎟⎞λ1 λ2 ... λC ( 2 . 8 ) where the matrix has C columns (one column for each constraint) and N rows. In general one must have N ≥ C+1 so the matrix is "tall" -- it has more rows th an columns. Recall that the intersection constraint surface (A( r)=0) has dimension N-C. If one were to a llow N < C, this surface would have negative dimension so the problem is overconstrained. If N = C, then (A(r )=0) has dimension 0 so the problem is constrained to a single point r in EN. This point either is or is not a stationary point according to Theorem 1 so there is no real stationarity problem to solve. Taking only the first C rows of the above matrix equation, one finds - ⎝⎜⎛ ⎠⎟⎞∂1f ∂2f ... ∂Cf = ⎝⎜⎛ ⎠⎟⎞∂1a1 ∂1a2 ... ∂1aC ∂2a1 ∂2a2 ... ∂2aC ... ... ... ... ∂Ca1 ∂Ca2 ... ∂CaC ⎝⎜⎛ ⎠⎟⎞λ1 λ2 ... λC ( 2 . 9 ) where the matrix is now square, CxC. Assuming this matrix has non-zero determinant, one can invert (2.9) to obtain expressions for the Lagrange multipliers λi(r) as functions of r. These λi functions can be inserted into equations (2.7) and then one can proceed to solve equations (2.6) and (2.7) for r. By making different versions of (2.9) taking di fferent sets of C rows from (2.8), one may produce several different sets of {λ i(r)} corresponding to different solutions r. Each stationary point will have its own set of Lagrange Multipliers { λi}, a fact which follows from Theorem 2. This will be clearly illustrated in Example 4 below. Stationary Points vs. Extremal Points The Method of Lagrange Multipliers identifi es stationary points of f some of which may be extremal points. We have not mentioned the domain of f : E N → reals, but in most problems there is some restricted domain of interest for r in f(r ) and an extremum might occur on the boundary of that domain, in which case it won't be picked up by the Method of Lagra nge Multipliers. The following 1D function f(x) plotted in red illustrates this point and shows various cases which can arise with stationary points : 13 ( 2 . 1 0 ) Here the domain of interest is the portion of the x axis lying between the two thin vertical bars. Point a is the true minimum and occurs at a boundary. Point b is the true global maximum, while point e is only a local maximum. Point d wants to be the global mi nimum but is trumped by point a. Point c has ∂ xf = 0 but is neither a local minimum nor a local maximum -- it is the 1D analog of a saddle point. Had points b and e the same height, there would be two equal global maxima. All these situations can occur for general f(r). Comments 1. For the Lagrange Multiplier method to be valid, the functions f(r ) and ai(r) must be continuous and differentiable in all arguments x i and the derivatives must also be continuous. The functions are therefore C0 and C1. This ensures that the gradients in (2 .4) are smooth continuous functions. The constraints are also assumed to be independent in that none can be written as a linear combination of the others, such as a 3(r) = 2a1(r) + a2(r). In this case if a 1(r) = 0 and a 2(r) = 0, then a 3(r) = 0 and a 3 adds nothing new. 2. The λ i terms in (2.5) often appear with minus signs in place of plus signs. This removes the minus signs in (2.8) and (2.9). These signs are just a c onvention and have no effect on the solution points r. 3. The function H in (2.2) is sometimes referred to as the Lagrangian and is written L. This Lagrangian differs from that which appears in the Lagrangian formulation of classical mechanics which we define in (D.38) of Appendix D. Here we show the connection. In (2.2) the function H is defined by, H( r,λ) ≡ f(r) + Σ i=1Cλiai(r) . (2.2) We require that d H = 0 ( 2 . 1 1 ) from which we conclude in (2.4) that a i(r) = 0 and 14 ∂H ∂xk = ∂f(r) ∂xk + Σi=1C λi ∂ai(r) ∂xk = 0 k = 1,2,....N . (2.12) The N coordinates x k are the components of r = (x1, x2, ....xN). In Appendix D the object that most closely corresponds to H is the action S, where S = ∫t1 t2 dt L(qk(t), q•k(t), t ) . (2.13) where L is the classical Lagrangian. The N coordinates q k(t) are the generalized coordinates (f unctions) of the Lagrange problem. We require, similar to (2.11), that δS = 0 . ( 2 . 1 4 ) This produces a set of modified Euler-Lagrange equations (D.55) which closely ressemble those in (2.12), ( ∂L ∂q k – d dt ∂L ∂q•k ) + Σi=1C λi ∂ai ∂qk = 0 k = 1,2...N (2.15) ∂f ∂x k + Σi=1C λi ∂ai ∂xk = 0 k = 1,2,....N . (2.12) Here the role played by ∂f ∂xk in (2.12) is assumed by ( ∂L ∂qk – d dt ∂L ∂q•k ) in (2.15). Just how this last expression arises from δS = 0 is shown at the end of Appendix D. The main point is that H is not the Lagrangian L, but is rather analogous to the time integral of the Lagrangian L which is the action S of (2.13). In the case of H we do a simple calculus variation dH relative to the coordinates x k, whereas for S we do a functional variation δS relative to the generalized coordinate functions q k(t). It is of course OK to refer to H as "the Lagrangian" as long as one understands this differs from L. 15 3. Three simple examples In this section the Method of Lagra nge Multipliers is used to find stationar y points for simple surfaces subject to simple constraints. Example 1 has r in E2, while Example 2 has r in E3. Note that the surface graphs for these examples, being of the form u = f( r), are in E3 and E4. Each of these examples could be reviewed in perh aps two small paragraphs. Instead, we have chosen to pour out a plethora of descriptive text with draw ings in an attempt to hammer home the basic ideas presented abstractly in the previous sections. The final Example 3 has r in EN and generalizes Examples 1 and 2. 3.1 Example 1: 3-dimensional hemisphere with one simple constraint The situation in E3 Let r = (x,y) in E2 and let u = f(x,y) = ± 4 - x2- y2 represent the surface of a sphere (radius R = 2) in E3 centered at the origin. We consider only the upper half of this surface, and we wish to find (x,y) that maximizes f. If there are no constraints, then (1.16) says ∇f = 0. The only place on our u = f surface having ∇f = 0 is the north pole of the sphere. This is just a regular "critical point" where ∂xf = 0 and ∂ yf = 0, so we are happy with this interpretati on of (1.24b) with no constraints, C = 0. We now add a constraint a(x,y) = 0 where a(x,y) = y- 1. This constraint is the line y=1 in the x-y plane E2. We can extrude this line into a plane y = 1 in E3. The hemispherical surface is a 2D surface in E3, and the extruded constraint plane is also a 2D surface in E3. These 2D surfaces intersect in a 1D surface which is a curve. There is hopefully some point on this curve which corresponds to a stationary point for f. Below is a picture of the sphere. The intersection of the upper spherical surface with the plane y = 1 is shown as a red curve (a half circle ). Due to the constraint, we cannot get to the north pole so the maximum value of f is some value less than that u = 2 at the north pole. The stationary point of this constrained problem will be at point A, and point B is not a stationary point . (3.1.1) 16 We pause to see how this example fits into the ge neral formalism. The space shown above in which we are plotting the graph u = f(x,y) is EN+1 whereas the space which holds r = (x,y) is EN, and N = 2. This r- space E2 is the equatorial plane in (3.1.1) and it is in this plane that we work when doing the Lagrange Multiplier analysis. What is potentially misleading is the notion of surface and gradient. A sphere of radius r is described by R( r) = r where R( r) = x2+y2+u2 . In spherical coordinates one easily computes that ∇R = (∂rR) r^ = r^, and this is indeed normal to the sphere at every po int on it, as predicted by Fact 2 (1.2) . Similarly, the y=1 plane constraint function a(x,y) = y-1 = 0 has ∇a = y^ and this is everywhere normal to the constraint plane. The potential confusion is that the Lagrange Multiplier Show does not play out on the stage shown in Fig (3.1.1). It plays out, as just noted, in the u = 0 plane of Fig (3.1.1) which is E2. The gradients of interest are 2D gradients in this plane, not the 3D gradients of surfaces in (3.1.1). Here then is plane of interest, oriented with the x axis down, (3.1.2) We wish to maximize the function f(x,y) = 4 - x2- y2 subject to the constraint y = 1. In (3.1.1) the spherical surface is u2+x2+y2 = R2 = 4, so f(x,y) is then height u. So the function 4 - x2- y2 = 4 - r2 (r is now the 2D r) is the height of the hemisphere in Fig (3.1.1) above the point (x,y). If we consider f(x,y) = K for a set of K values, we get a set of level curves in (3.1.2) which are circles with r = 4-K2 . On each circle in (3.1.2), the height of the sphere ly ing over the surface in (3.1.1) is constant -- that is why they are called level curves. Think of these as lines of latitude projected onto the u=0 equatorial plane. In this problem there is only C = 1 constraint. In the general formalism we would refer to this constraint and its Lagrange multiplier as a 1 and λ1, but here we shall just call them a and λ. The constraint is a(x,y) = 0 with constraint functi on a(x,y) = y-1, and the constraint y = 1 is shown as the red line in (3.1.2). We have, 17 f = 4 - x2- y2 = 4 - r2 a = y - 1 ( 3 . 1 . 3 ) so ∇f = ∇4 - r2 = [-r/ 4-r2 ] r^ // 2D gradients ∇a = ∇ ( y-1) = y^ . ( 3 . 1 . 4 ) For any point r in the disk of (3.1.2), ∇f therefore points toward the disk center, and this is then the uphill direction for the hemisphere in (3.1.1). In Cartesian coordinates, ∇f = ∇4 - x2- y2 = - (x/f) x^ - (y/f) y^ = - x 4 - x2- y2 x^ - y 4 - x2- y2 y^ . (3.1.5) Now for just one constraint (1.24b) reads, ∇f(r) = – λ∇a(r) . ( 3 . 1 . 6 ) How do we make geometric sen se of this equation? If r is stationary point, we should have df = 0 for any legal small displacement d r starting at this point r. Eq. (3.1.6) dotted into d r then says, 0 = df = ∇ f(r) • dr = – λ∇a(r) • dr . (3.1.7) If we displace d r in any direction on the constraint surface (which from Fact 2 has normal ∇a(r)), then ∇a(r) • dr = 0 and we find df = 0, so we are thus at a stationary point of f. For other values of position r we will find either that df > 0 or df < 0 if we move d r along the constraint surface. If we are looking for a maximum of f, for example, then it is to our advantage to move a little d r in a direction for which df > 0. Doing this repeatedly with a computer program (e.g. GNU Octave freeware) to find a maximum (minimum) of f is called the method of gradient ascent (descent) . In our example with only one constraint, (3.1.6) sa ys that we have reached a stationary point when ∇f and ∇a are collinear. Let's draw onto (3.1.2) the points a and b which lie under points A and B in (3.1.1), and we include a point b' which is a mirror image of b : (3.1.8) 18 At each of the points a,b,b' we show ∇f in black and ∇a in red. At point a we have reached the closest distance to the center that is allowed for points on th e red constraint line, so this will correspond to the maximum value of f, which is point A in (3.1.1). At point a we see that in fact the two gradients are collinear as required by (3.1.6) at a stationary point. Suppose we are at point b. We may compute df for a small displ acement upwards (minus x direction) df( b) = ∇ f • dr = [ - (x/f) x^ - (y/f) y^ ] • |dx|(- x^) = |dx| (x/f) > 0 since x>0 and f>0 . (3.1.9) Since df > 0, it is to our advantage in finding max f to move upwards in (3.1.8). Conversely, suppose we are instead at point b'. Then moving toward a gives, df( b') = ∇ f • dr = [ - (x/f) x^ - (y/f) y^ ] • |dx|( x^) = |dx| (-x/f) > 0 since x<0 and f>0 (3.1.10) and again we find that df > 0. From either starting position b or b', moving toward point a is a win. In this Example ∇ a = ∇ (y-1) = y^ spans the 1 dimensional perp space at any point r on the constraint surface A( r) = 0, that is, y = 1. Similarly, x^ spans the 1 dimensional tangent space at any r on A( r) = 0. Finding the stationary point It is an easy matter to compute the so lution value of the Lagrange Multiplier λ 1 and r. Inserting (3.1.4) into (3.1.6) gives ∇f(r) = – λ∇a(r) ⇒ - x 4 - x2- y2 x^ - y 4 - x2- y2 y^ = -λ y^ . (3.1.11) This says x = 0, and then using the constraint y = 1, - 1 4 - 1 y^ = -λ y^ ⇒ λ = 1/ 3 . (3.1.12) We have thus found the following stationary point r and Lagrange multiplier λ for Example 1: r = (0,1) λ = 1/ 3 . ( 3 . 1 . 1 3 ) Therefore the height of point A in Fig (3.1.1) is u = 4 - x2- y2 = 4 - 02- 12 = 3 . Now consider, setting (x,y) = (x 1,x2), f = 4 - x12- x22 ∂if = (1/2) f-1 (-2xi) = - xi/f ∂ 2 if = - [f * 1 - x i∂if] / f2 = - [f - x i(-xi/f)] / f2 = - (1/f) - (x i/f)2 < 0 (3.1.14) 19 Thus ∂2 1f < 0 and ∂2 2f < 0 for all points (x 1,x2) on the equatorial disk under the hemisphere, as is intuitively obvious. Thus these two inequalities apply at the obtained stationary point r, and this confirms that this stationary point is a maximum of the function f. A look at H(x,y, λ) . Theorem 2 of (2.3) claims that the unconstrained "Lagrangian" H(x,y, λ) = f(x,y) + λ(y-1) should have a stationary point at the point (x,y, λ) = (0,1,1/ 3 ) shown in (3.1.13). Consider first, ∂H/∂λ = (y-1) = 0 at the point (x,y, λ) = (0,1,1/ 3 ) // λ partial is 0 at stationary point ∂2H/∂λ2 = 0 // since H is linear in λ // no λ curvature anywhere (3.1.15) Next, ∂xH = ∂xf + λ∂xa = ∂xf = -x/f = 0 at (0,1,1/ 3 ) // x partial is 0 at stationary point ∂yH = ∂yf + λ∂ya = -y/f + λ = -1/ 3 +1/ 3 = 0 // y partial is 0 at stationary point (3.1.16) as expected. Just for fun, here is a plot of the function H with λ set to the solution value λ = 1/ 3 , H(x,y, λ=1/ 3 ) = 4 - x2- y2 + (1/ 3 ) (y-1) . (3.1.17) Now u = (1/ 3 ) (y-1) appearing in the second term of (3.1.17) is a plane sloping up to the right in Fig (3.1.1). This is different from the plane y = 1 which is the vertical constraint plane in that figure. In (3.1.17) we are adding a spherical surface to a plane sloping up to the right, and the result is an ellipsoidal-like surface (really a quartic surface) which should have a stationary point at r = (0,1) : (3.1.18) 20 One can see from these plots that in fact r = (0,1) is a stationary point of H(x,y, λ=1/ 3 ) and is in fact a maximum. In this example it just happens that the stationary point is a maximum for both f and H. More complicated Examples of the same type If we had a surface in Fig (3.1.1) more complicated than a sphere, and a constraint more complicated than y = 1, the nature of the interpretation of (3.1.2) do es not change. For example, if Fig (3.1.1) showed an oblong bowl surface u = f(x,y) whose level curves were those of (1.4) (u = K values increasing toward the outside) and if the constraint a(x,y) = 0 were some arbitrary constraint curve (shown in red), we would have this picture: (3.1.19) The black arrows are the 2D "up-hill" bowl gradients ∇f (larger toward the outside) while the red arrows ∇a are normals to the constraint surface ( ∇a now varies along that surface). The solution point is a where ∇f and ∇a are collinear. At point b the gradient arrows do not line up, and there is advantage in this case for decreasing f by moving toward point a. One can see from the "topo lines" that a is indeed the point on the constraint curve that has the least value of u. Notice that saying the two arrows in E2 are collinear is to say they are linearly dependent. 3.2 Example 2: 4-dimensional hemisphere with two simple constraints In order to better dem onstrate the interpretation of (1.24b) concerning the Lagrange multiplier gradients, we upgrade Example 1 and then allow two constraints. We take the surface of Fig (3.1.1) to be the "upper half" of a 4D sphere of radius 2, u = f(x,y,z) = 4 - x2- y2- z2 . We must imagine a 4D drawing of this sphere which is the new Fig (3.1.1) with a vertical u axis and three "horizontal" axes x,y,z. We cannot draw this picture, but luckily the Lagrange Multiplier Show for this problem plays out in E3 where we can draw pictures. Now with r = (x,y,z) and r = x2+y2+z2 the level surfaces (really surfaces this time) are given by u = 4 - r2 = K for a set of constant K values. Since then r = 4-K2 , we see that these level surfaces are in fact a set of concentric spheres of radius 4-K2 ( so we shall restrict to |K| ≤ 2 ). We take as our two constraints the equations y = 1 and x = 1, just to keep it simple. Thus 21 f(x,y,z) = 4 - x2- y2- z2 = 4 - r2 a1(x,y,z) = y-1 a2( x , y , z ) = x - 1 . ( 3 . 2 . 1 ) The Lagrange method gradients of interest are now 3D gradients in E3 space, so ∇f = ∇ (4 - r2 ) = [-r/ 4-r2 ] r^ = - (x/f) x^ - (y/f) y^ - (z/f) z^ points to sphere center ∇a1 = ∇(y-1) = y^ ∇a2 = ∇(x-1) = x^ . ( 3 . 2 . 2 ) Equation (1.24b) reads ∇f(r) + λ 1∇a1(r) + λ2 ∇a2(r) = 0 . (3.2.3) This says that the three gradients must be linearly de pendent, which means they must be coplanar! Before solving the problem, we would like to draw a picture in E3 corresponding to Fig (3.1.2) in E2 for Example 1. Although one could in fact draw such a pict ure with some effort, we shall decline this task and instead draw two z = constant slices of the desired picture. On the left below is a slice at z = 0, while on the right is a slice at z = 1 (these are not precision drawings) : ( 3 . 2 . 4 ) In these pictures, we are viewing the two red constraint planes y = 1 and x = 1 edge on. The surface which satisfies both constraints is a line normal to the plane of paper which is th e intersection of the two constraint planes. The circles on the left are slices of the family of concentric spheres taken at the equator z = 0. At z = +1 since our slicing plane is closer to the north pole (z=2), some of the inner spheres no longer cut the plane z = 1, and those that still do have smaller diameter circles. 22 Consider now the point r = b shown as a black dot on the right. The two red arrows are ∇a1 = y^ and ∇a2= x^, each normal to its constraint plane. The black arrow dips into the plane of paper because ∇f = [-r/ 4-r2 ] r^ which points toward the sphe re center. The three arrows ∇f,∇a1,∇a2 are therefore not coplanar, so they are linearly independent. That mean s that equation (3.2.3) cannot exist for this point r. This black point b is located outside the r = 7 units sphere (in our crude picture). Consider next the point r = a shown as a black dot on the left. Since this is the equatorial slice, the black arrow ∇f lies in the plane of paper. This means the three vectors ∇f,∇a1,∇a2 are coplanar in the z=0 plane so they are linearly dependent . For this point a , the equation (3.2.3) can and does exist, and therefore point a is the problem solution. This point is located inside the r = 6 unit sphere. Once again, at an extremum point we should have df = 0. Moving an amount d r which is consistent with both constraints (meaning d r is in the z direction) we then have from (3.2.3), df = ∇f • dr = – λ1[∇a1(r)• dr] – λ2 [∇a2(r)• dr] = –λ1 [0] - λ2[0] = 0 (3.2.5) and sure enough, df = 0. At point b one finds for a d r pointed down toward the equatorial plane, df = ∇f • dr = [ - (x/f) x^ - (y/f) y^ - (z/f) z^ ] • (-|dr| z^ ) = (z/f)dr > 0 (3.2.6) and so it is advantageous to move d r = |dr| z^ toward a and thereby increase f, so b is not an extremum. Evaluating (3.2.6) at z = 0 for point a again shows df = 0. Finally, we solve the problem. Inserting (3.2.2) into (3.2.3) gives ∇f(r) = – λ1∇a1(r) – λ2 ∇a2(r) ⇒ [- x 4 - x2- y2- z2 x^ - y 4 - x2- y2- z2 y^ - z 4 - x2- y2- z2 z^ ] = – λ1y^ – λ2 x^ . (3.2.7) We see at once that z = 0 and then the above becomes x 4 - x2- y2 = λ2 y 4 - x2- y2 = λ1 . (3.2.8) But the constraints say x = 1 and y = 1 so, 1 4 - 12- 12 = λ2 1 4 - 12- 12 = λ1 ⇒ λ1 = λ2 = 1/ 2 . (3.2.9) Therefore the solution to Example 2 is this: r = (1,1,0) λ 1 = 1/ 2 λ2 = 1/ 2 . (3.2.10) The value of u at this point is u = 4 - x2- y2- z2 = 4 - 12- 12- 02 = 2 . It is simple matter conceptually to generalize our Example 2 to more complicated functions f,a 1,a2. For some arbitrary function f(r ) there will be a set of 2D surfaces in E3 which are the level surfaces on which 23 f(r) = K for various values of K. These level surfaces th en replace the set of concentric spheres of our simple case. We then imagine replacing our simple planar constraint functions a 1(x,y,z) and a 2(x,y,z) with general functions which then result in arbitrarily shaped 2D constraint surfaces a 1(x,y,z) = 0 and a2(x,y,z) = 0 in E3. The gradients ∇a1 and ∇a2 vary over their respective surfaces (being normal vectors). The intersection of these two 2D constraint surfaces in E3 will be a curve in E3 (1D surface in E3). Potential solution points r must lie on this curve. At each point r on this curve, the gradients ∇a1(r) and ∇a2(r) define a local plane. At a stationary point one will find that the gradient ∇f(r) lies in the same plane as ∇a1(r) and ∇a2(r) so the three gradients are linearly dependent and then ∇f(r) + λ1∇a1(r) + λ2 ∇a2(r) = 0 for some λ1 and λ2. At such a point r, if one considers any d r which lies on both constraint surfaces, one will then find that df = ∇f •dr = 0 so the point r is therefore a stationary point. The situation is described by a 3D version of drawing (3.2.4). Note that the intersection of the two constraint su rfaces might result in multiple curves, and each must then be considered. For example, the intersection of two thin (prolate) ellipsoids might be two closed curves. If the intersection of the constraint surfaces is null ( 1st ellipsoid inside 2nd ), the problem has no solutions. 3.3 Example 3: N-dimensional hemisphere of radius R in EN with C simple constraints We now generalize the above two examples. Consider an upper-half hypersphere which is an N-1 dimensional surface embedded in E N. The sphere has radius R and there are C constraints x i = 1 for i = 1,2..C. Proceed as in the previous examples where now f( r) is defined on EN and the hypersphere surface "graph" is given by u = f( r) in EN+1. r = (x1,x2....xN) f(r ) = f(x1,x2....xN) = + R2- Σi=1Nxi2 // u = f( r) is upper-half hypersphere in EN ai = xi - 1 i = 1,2..C . // C constraints are x i = 1 for first C coordinates (3.3.1) The gradients of interest (N-dimensional) are ∇f = - Σi=1N(xi/f)x^i ∇ai = x^i i = 1,2..C . (3.3.2) Insert these into (1.24b) , ∇f(r) = - Σ i=1Cλi∇ai(r) (1.24b) so - Σi=1N(xi/f)x^i = - Σi=1Cλi x^i or Σ i=1C(xi/f)x^i + Σi=C+1N(xi/f)x^i = Σi=1Cλi x^i or - Σ i=1C [ (xi/f) - λi] x^i + Σi=C+1N (xi/f) x^i = 0 . (3.3.3) Solve to get 24 λ i = xi/f i = 1,2..C x i = 0 i = C+1,....N . (3.3.4) Since x i = 1 for i = 1 to C one has λi = 1/f where f = R2- Σi=1Cxi2 = R2- C . ( 3 . 3 . 5 ) The stationary point is then determined by λi = 1/ R2- C i = 1,2...C // the C Lagrange multipliers are all the same r = (1,1...1,0,0...0) . // C ones fo llowed by N-C zeros (3.3.6) There is no stationary point if C > R2 since then the λ i become imaginary. Finally, we verify that the results of Examples 1 and 2 are recovered: Example 1: N = 2, C = 1, R = 2, λ 1 = 1/ 3 , r = (1,0) Example 2: N = 3, C = 2, R = 2, λ i = 1/ 2 , r = (1,1,0) (3.3.7) 3.4 Example 1 Revisited: What is the half-width of the peak? The question pursued here seems of small interest, but it serves as a model for the situation which arises in Section 5.3 where we are interested in the width ΔNi of the peak of the Boltzmann microstate distribution Ω in (5.3.16). Recall from (3.1.1) where u = f(x,y) = 4-x2-y2 , (3.1.1) 25 The upper spherical surface is a graph in E3 of u = f(x,y) = 4-x2-y2 . The variables x,y are not independent. When the constraint y = 1 is applied we get u = f(x,1) = 3-x2 ≡ F(x) which is the equation of the red curve shown in the picture. Here we have invented a new function name F(x) to represent the function f(x,y) after y has been eliminated by the constraint. One is left with only one independent variable which is x. We could then ask about the half-width of the "p eak" which is the red curve. As an estimate, we model the red curve as a parabola using a Taylor expansion around the peak point A, F(x+Δx) ≈ F(x) + (∂ F/∂x) (Δx) + (1/2) ( ∂ 2F/∂x2) (Δx)2 (3.4.1) Maple computes the two derivatives and then evaluates them at point A where x = 0: (3.4.2) As expected, ( ∂F/∂x) = 0 since A is a stationary point for the function F(x) of one independent coordinate x. Using ( ∂ 2F/∂x2) = -1/ 3 we rearrange (3.4.1) to get ΔF = F(x+Δx) - F(x) = -(2 3 )-1(Δx)2 . ( 3 . 4 . 3 ) To get the half-width of the peak we set (the half width at half maximum). |ΔF / F ( 0 ) | = 1 / 2 ( 3 . 4 . 4 ) so that | -(2 3 )-1(Δx)2 / 3 | = 1/2 ⇒ (Δx) ≈ 3 = 1.73 . (3.4.5) 26 The following plot shows the parabolic approximatio n in blue and marks the estimated half width (3.4.6) As a percentage of the peak value, the half width here is 100%, so it is a fairly large half-width. In the Boltzmann problem of Section 5, we shall start with a function f(N 1,N2.....Nm) of m variables which are not independent. We use up the two probl em constraints to eliminate the variables N 1 and N2 and then define a function F(N 3,N4.....NM) ≡ f(N1(N3,... Nm), N2(N3,... Nm), N3, .....Nm). In the space of the independent variables N 3, N4....Nm we use a vector version of the above procedure to estimate the half width of the peak of the function F in any of the variables N i, as illustrated in (5.3.18). For large N i values, we find that the peak is very narrow compared to its height, unlike the peak shown in (3.4.6). 27 4. Example 4: Extremal distances between an ellipse and a circle Many examples of the use of the Method of Lagrange Multipliers can be found in texts and on the web. See for example Trench or the wiki Lagrange Multipli er page. Example 4 is a bit more complex than the toy examples in texts and in our Section 3, and we obtain a numerical solution after carrying out the analytic steps of the Lagrange Multiplier method. Example 4 has two constraints so C = 2 and there are f our coordinates so N = 4. It happens that these four coordinates are not those of a single point in E 4 but represent two points in E2. We adopt several new notations in this section : 1. The constraint functions will be called a and b instead of a 1 and a2, but the multipliers are still λ1,λ2. 2. The partial derivatives of a function F(x 1, x2 ....) shall be represented as F i : Fi ≡ ∂iF(x1, x2 ....) ≡ ∂iF(x1, x2 ....) ∂xi which is a notation promoted by Buck with generalization F ijk.. = ∂i∂j∂k ... F(x1, x2 ....). 3. The components of the vector r will be called (x,y,x',y') instead of (x 1,x2,x3,x4) . Consider then these two curves in E2: x2/A2 + y2/B2 = 1 ellipse centered at the origin, semimajor/minor axes A and B (x'-α) 2 + (y'-β )2 = R2 . circle of radius R centered at ( α,β) (4.1) Problem : What is the minimum and maximum distance between these two curves? The distance squared is given by f = (x-x') 2 + (y-y')2 = f(x,y,x',y') N = 4 coordinates (4.2) while the constraint functions are a(x,y,x',y') = x 2/A2 + y2/B2 - 1 a(x,y,x',y') = 0 b(x,y,x',y') = (x'-α ) 2 + (y'-β)2 - R2 b(x,y,x',y') = 0 . (4.3) The various derivatives are, 28 f1 = ∂xf = 2(x-x') a 1 = ∂xa = 2x/A2 b1 = ∂xb = 0 f2 = ∂yf = 2(y-y') a 2 = ∂ya = 2y/B2 b2 = ∂yb = 0 f3 = ∂x'f = -2(x-x') a 3 = ∂x'a = 0 b 3 = ∂x'b = 2(x'-α) f4 = ∂y'f = -2(y-y') a 4 = ∂y'a = 0 b 4 = ∂y'b = 2(y'-β ) . (4.4) The H function (2.2) is then, H(x,y,x',y', λ 1,λ2) = f + λ1a + λ2 b = (x-x') 2 + (y-y')2 + λ1 [x2/A2 + y2/B2 - 1] + λ2 [(x'-α)2 + (y'-β )2 - R2] . (4.5) Compute the six H i derivatives as prescribed in (2.4), H 1 = f1 + λ1a1 + λ2 b1 = 2(x-x') + λ 1 2x/A2 H2 = f2 + λ1a2 + λ2 b2 = 2(y-y') + λ 1 2y/B2 H3 = f3 + λ1a3 + λ2 b3 = -2(x-x') + λ22(x'-α) H4 = f4 + λ1a4 + λ2 b4 = -2(y-y') + λ 22(y'-β) H 5 = a H6 = b . ( 4 . 6 ) Set these derivatives to 0 as in (2.4) to find this set of six equations (x-x') + λ 1 x/A2 = 0 1 (y-y') + λ 1 y/B2 = 0 2 - (x-x') + λ2 (x'-α) = 0 3 - (y-y') + λ 2 (y'-β) = 0 4 x 2/A2 + y2/B2 - 1 = 0 5 (x'-α)2 + (y'-β)2 - R2 = 0 6 ( 4 . 7 ) where the 6 unknowns are x,y,x',y',λ 1,λ2 . Notice that a possible λ i solution is λ1 = λ2 = 0 in which case one gets x = x' and y = y'. Certainly this is an extremum situation since distance squared is f = (x-x')2 + (y-y')2 = 0, a minimum. Insertion of x' = x and y' = y into 5 and 6 above and el imination of x yields (after some algebra) an equation for y of the form k1y4 + k2y3 + k3 y2 + k4 y + k4 = 0 . x = ± A 1-(y/B )2 (4.8) If the ellipse and circle intersect, one will find that the above equation has 1,2,3 or 4 real roots which then represent the intersection point(s) of the two curves. If the curves don't intersect, the solution y values will be complex so there are no physical solutions for λ1= λ2 = 0 . 29 So how does one go about solving the set of 6 equations shown in (4.7)? A possible first step is to find a viable set of λi. From a version of (2.9) using the first and third rows of (2.8) one may write, - ⎝⎛ ⎠⎞f1 f3 = ⎝⎛ ⎠⎞ a1 b1 a3 b3 ⎝⎛ ⎠⎞λ1 λ2 ≡ M ⎝⎛ ⎠⎞λ1 λ2 (2.9) where M = ⎝⎛ ⎠⎞ a1 b1 a3 b3 = ⎝⎛ ⎠⎞ a1 0 0 b3 // b 1 and a3 vanish from (4.4) (4.9) so that M -1 = ⎝⎛ ⎠⎞ b3 0 0 a1 / det(M) = 1 a1b3 ⎝⎛ ⎠⎞ b3 0 0 a1 . ( 4 . 1 0 ) Therefore, - ⎝⎛ ⎠⎞λ1 λ2 = M-1 ⎝⎛ ⎠⎞f1 f3 = 1 a1b3 ⎝⎛ ⎠⎞ b3 0 0 a1 ⎝⎛ ⎠⎞f1 f3 ( 4 . 1 1 ) and thus -λ 1 = 1 a1b3 b3f1 - λ2 = 1 a1b3 a1f3 ⇒ λ1 = -f1/a1 = -2(x-x')/[2x/A2] = - A2(x-x')/x λ 2 = - f3/b3 = 2(x-x')/[2(x'-α)] = (x-x')/(x'-α ) . We then have a candidate solution set { λi}. The first four equations of (4.7) become (x-x') + [- A 2(x-x')/x] x/A2 = 0 1 (y-y') + [- A2(x-x')/x] y/B2 = 0 2 - (x-x') + [ (x-x')/(x'- α)] (x'-α) = 0 3 - (y-y') +2[ (x-x')/(x'- α)] (y'-β) = 0 4 or (x-x') - (x-x') = 0 1 (y-y') - (A/B) 2(x-x')y/x = 0 2 - (x-x') + (x-x') = 0 3 - (y-y') + (y'- β)(x-x')/(x'- α) = 0 4 . (4.12) Two of these equations are identities, while the other two are 30 (y-y') - (A/B)2(x-x')y/x = 0 2 - (y-y') + (y'- β)(x-x')/(x'- α) = 0 4 or (y-y')x - (A/B) 2(x-x')y = 0 2 - (y-y')(x'-α ) + (y'-β)(x-x') = 0 4 . (4.13) Including the last two equations of (4.7), one h as four equations in four unknowns x,y,x',y' : (y-y')x - (A/B)2(x-x')y = 0 2 - (y-y')(x'-α ) + (y'-β)(x-x') = 0 4 x 2/A2 + y2/B2 - 1 = 0 5 (x'-α)2 + (y'-β)2 - R2 = 0 . 6 (4.14) This is a system of four 2nd degree polynomial equa tions in four variables. When we ask Maple to analytically solve these equations for the four unknowns x,y,x',y', it succeeds but the results are very messy and involve roots of fourth degree polynomials. In this situation, the analytic solution exists, but is so complicated it hardly seems very useful. We selecte d this example just to show that even a relatively simple problem can be nearly intractable analytically. We turn then to a numerical solution for the particular ellipse + circle case shown in this scaled figure where the background boxes are 1 unit squares, (4.15) The circle has radius R = 2 and is centered at ( α,β) = (3,1). The semi ellipse axes are A = 8 and B = 4. We enter into Maple the 6 original equations (4.7) in the 6 unknowns x,y,x',y', λ1,λ2 (xp = x', L1 = λ1, etc) : 31 ( 4 . 1 6 ) We next set in the specific parameters for the figure shown above ( 4 . 1 7 ) Without given a starting point, Maple finds the following solution which we approximate as λ 1 = -2.8 λ2 = .32 (x,y) = (3.66, 3.56) (x', y') = (3.50, 2.94) . (4.18a) Plotting this on our figure, 32 (4.18b) we see that this solution corresponds to the minimu m distance between the two curves. Maple computes this minimum distance to be (x-x')2+(y-y')2 = .6406470846 units. Another solution is found by giving Maple a search starting point (-8,0) on the ellipse and (5,1) on the circle : which we approximate as λ 1 = -104 λ2 = -6.5 (x,y) = (-7.99, -.22) (x', y') = (4.99, 1.22) . (4.19a) Plotting this solution, (4.19b) we see that this solution corresponds to the maximu m distance between the two curves. Maple computes this maximum distance to be 13.05541338 units. A third candidate solution is found by giving Maple a starting point (-8,0) on the ellipse and (1,1) on the circle: 33 which we approximate as λ1 = -72 λ2 = 4.5 (x,y) = (-7.99, -.22) (x', y') = (1.01, 0.78) . (4.20a) Plotting this solution, (4.20b) one sees that this solution represents a "stationary poin t" but is not in fact a solution to the problem. A slight d r = (dx,dy,dx',dy') variation along both constraint curves from this solution point gives df = 0 where f is the squared distance (4.2). Finally, we find this fourth solution, which we approximate as λ 1 = -20.2 λ2 = -2.3 (x,y) = (3.66, 3.56) (x',y') = (2.50, -.94) . (4.21a) Plotting this solution, (4.21b) we find a situation similar to the third candidate so lution -- a stationary point which is neither a global maximum nor a global minimum. 34 5. Example 5: The Boltzmann factor in Statistical Mechanics In the following phy sics example the role of r =(x1,x2....xm) is played by N = (N1,N2...Nm). The function f(r) to be maximized will be Ω(N) = g1N1 N1! g2N2 N2! ....gmNm Nm! . The two constraints are a( N) = ΣiNi - M = 0 and b( N) = ΣiεiNi-U = 0, all as described below. 5.1 Statement of The Boltzmann Extremum Problem I magine a physical system which has m state levels into which a set of M identical particles can be placed. At a given energy level εi, the number of available states is gi, known as the degeneracy at energy εi. A particular state of the entire system of particles is characterized by the number of particles Ni in each of the states, so system state = (N 1,N2...Nm) = N , a vector. The total number of particles in the system is fixed at M = ΣiNi, and the total energy is U = Σ iεiNi. The system is isolated so its initial energy U does not change, nor can particles be created or destroyed, so U and M are both fixed constants. The system of particles (perhaps a box of gas atoms) is assumed to be at thermal equilibrium so all its macroscopic characteristics (like pressure and temperature) are stable. Imagine a billion of these particle systems (an "ensemble"). Each system settles into its own partition set {Ni} = N. If M is very large, one will find that those billion N vectors are all very similar. This is so because certain {N i} partitions are statistically favored over ot her partitions just because there are a lot more "states" available to a system with those favored {N i} partitions (hence "statistical" mechanics). So, how many states are available ("accessible") to such a system characterized by vector N which represents some partitioning of the particles {N i}? Placing N identical particles in the g degenerate states of a system at energy level ε is like placing N balls in a set of g boxes. The number of different "ways" of placing N indistinguishable balls into g boxes (bins) is a fascinating elementary problem and the answer is ⎝⎛ ⎠⎞ N+g-1 N, a binomial coefficient. Footnote : One arrives at this answer by thinking of N balls and g-1 "partitions" between boxes, all of which have to be laid out in a row, such as | * * * | * | | * for the case N = 5 and g-1 = 4. Here the leftmost bin is empty, as is the second bin from the right. The asterisks are the identical balls. Each layout of the 9 objects is characterized by pick ing a committee of 5 balls from the 9 positions (or a committee of 4 partitions). In this case the number of unique layouts is ⎝⎛ ⎠⎞ 9 5 = 9! 5!4! = ⎝⎛ ⎠⎞ 9 4 . (5.1.1) The count ⎝⎛ ⎠⎞ N+g-1 N allows multiple balls (particles) to be in any given box (state). This count is then associated with so-called Bose-Einstein statistics -- multiple bosons are allowed in the same state. Bosons are particles with integral spin such as photons (S=1) and ground-state helium atoms (S=0). If at most one ball is allowed in an y given box, then one must have N ≤ g and the count of ways is just the number of ways to pick a committee of N boxes (one ba ll into each box) from the set of g boxes, and this 35 count is ⎝⎛ ⎠⎞ g N . This count is then associated with so-called Fermi-Dirac statistics -- only one fermion is allowed in a state. Fermions are particles having half-integral spin. Because electrons (S=1/2) are fermions, atomic energy levels get "filled up" as N increases, allowing there to be some small number of electrons at the highest energy level ("valence electrons") which cause an atom of atomic number N (an "element") to have certain chemical properties. If electrons were bosons, they would all to go into the lowest or "ground" atomic state and the world would be a very different place. In other words, the fact that electrons are fermions is the reason the periodic table of elements exists. However, we shall assume that g i >> Ni (known as the Boltzmann limit), and this allows a simplification of the state count. Both the Bose and Fermi state counts reduce to the same count in this limit: ways = ⎝⎛ ⎠⎞ N+g - 1 N ≈ ⎝⎛ ⎠⎞ g N = g! N! (g-N)! = g(g-1)..(g-N+1) N! ≈ gN N! . (5.1.2) So in our particle problem, the number of ways of having N 1 particles in states with energy ε1 is g1N N1! . For each of these ways, there are g2N N2! ways of putting N 2 particles in energy level ε2. And so on. Thus, the total number of available states ("microstates") for a partition {N i} of the M particles is given by Ω(N1,N2...Nm) = g1N1 N1! g2N2 N2! .... gmNm Nm! . // see for example Zemansky Eq. (10-4) (5.1.3) Since the particles are identical, there is only one distinct way to partition the set of N particles into sets having counts N 1, N2....Nm (prior to putting the particles of each set into that set's states as done above). For example, put the N particles in a straight line, a nd draw m-1 divider lines as needed to get the desired partition counts. So one could add an overall factor of "1" to (5.1.3) to denote this one partition. Here then is our initial extremum problem : Find N which maximizes Ω(N) = g1N1 N1! g2N2 N2! ....gmNm Nm! subject to these two constraints: Σ iNi = M and ΣiεiNi = U . (5.1.4) When M is a very large number (like Avogadr o's number), it turns out (coming soon) that Ω(N) has a very strong and sharp maximum at a certain N which is the solution value N for this extremum problem. Thus, when one examines an ensemble of such systems, this solution N is the overwhelmingly likely partitioning of the particles into {N i}. Our task is then to solve this problem for N. This vector N = (N 1,N2...Nm) is the vector r = (x1,x2...xN) of our general Lagrange Multiplier presentation in Sections 1 and 2, so in this application the number of components in the vector is N = m, and we don't want to confuse this previous use of symbol N with N or N i of the present problem! 36 One might observe that the x i in the general analysis were reals, whereas the N i here are integers. We can dispense with this issue by simply writing N i! = Γ (Ni+1) which then serves to interpolate N i between its integer values. Another approach is to replace N i everywhere in the above problem statement by ni ≡ Ni/M. Since M and the solution N i values are very large integers, one can regard this n i as essentially a continuous real (albeit rational) variable. It should be noted that the value of U is restricted to a certain range, Mε min ≤ U ≤ Mεmax . ( 5 . 1 . 5 ) The limits represent the two extreme possible partitions {N i} of the particles (all in lowest energy state, or all in highest energy state). It is useful to define the average energy of a particle, u ≡ U/M = average energy of a particle in the system (5.1.6) and then (5.1.5) states the obvious fact that ε min ≤ u ≤ εmax . ( 5 . 1 . 7 ) 5.2 Solution of The Boltzmann Extremum Problem Maxim izing the state count Ω is the same as maximizing f ≡ lnΩ, and from (5.1.3), f(N 1, N2...Nm) ≡ lnΩ(N1,N2...Nm) = Σi [ Ni ln gi ] - Σi [ ln Ni! ] (5.2.1) where Σi means Σi=1m. Assuming all the N i are large numbers, we may use Stirling's formula, x! ≈ 2πx xx e-x ⇒ lnx! ≈ ln 2π + (1/2) lnx + xlnx - x ≈ xlnx - x, x >> 1 . (5.2.2) Then ln N i! ≈ NilnNi- Ni so that, f = ln Ω = Σ i [ Ni lngi ] - Σi [ ln Ni! ] ≈ Σi [ Ni lngi ] - Σi [NilnNi- Ni] = Σi [(1+lngi)Ni - NilnNi ] = Σi Ni [1 + ln(g i/Ni)] . (5.2.3) Here then is a restatement of the extremum problem given in (5.1.4): Find N which maximizes f( N) = Σ i [(1+ln g i)Ni - NilnNi ] subject to these two constraints: a ( N) = ΣiNi - M = 0 b ( N) = ΣiεiNi-U = 0 . (5.2.4) 37 Application of the Method of Lagrange Multipliers Following now the prescription of Section 2, we write our "Lagrangian" H of (2.2) as H( N,λ) ≡ f(N) + λ1a(N) + λ2 b(N) = Σ i [(1+ln g i)Ni - NilnNi ] + λ1 [ ΣiNi - M] + λ2 [ΣiεiNi-U] . (5.2.5) The next step is to compute th e derivatives with respect to N i and λi and then set those derivatives to 0. We know that the two equations ∂ H/∂λi= 0 just replicate our two constraints, so it is the other m derivatives which are of interest. We again use the s ubscript notation for partial derivatives introduced at the start of Section 4, so now F i ≡ ∂F(N)/∂Ni . We need these three contributions to H i, f = Σ i [(1+lngi)Ni - NilnNi ] ⇒ fi = ∂f/∂Ni = (1+lng i) - Ni (1/Ni) - lnNi = lngi - lnNi = ln(gi/Ni) a = Σ iNi - M ⇒ ai = ∂a/∂Ni = 1 b = Σ iεiNi-U ⇒ bi = ∂b/∂Ni = εi . i = 1,2...m (5.2.6) Now set the H i partials to 0 : 0 = H i(N,λ) = fi + λ1ai + λ2 bi = ln(g i/Ni) + λ1*1 + λ2*εi i = 1,2...m or ln(N i/gi) = λ1 + λ2εi . i = 1,2...m (5.2.7) Here then are the equations we need to solve, ln(N i/gi) = λ1 + λ2εi i = 1,2...m // m equations Σ iNi = M / / a ( N) = 0 Σ iεiNi = U / / b ( N) = 0 . (5.2.8) This is a system of m+2 equati ons in m+2 unknowns which are the N i, λ1 and λ2. To solve these equations, we first solve the first set of m equations for N i, leaving λ1 and λ2 as unknowns: 38 ln(N i/gi) = λ1+ λ2εi or Ni = gi eλ1 eλ2εi . ( 5 . 2 . 9 ) Already a major result has emerged from our Lagrange Multiplier analysis! The population N i of states at energy level εi depends exponentially on εi. We intuitively expect higher energy levels to be less populated than lower ones, so we expect λ2 to be negative. For this reas on, and following tradition, we set λ2 = - β where β > 0. And while we're at it, we can set eλ1 = A to simplify notation. So β ≡ - λ 2 // derived Lagrange multipliers A ≡ eλ1 ( 5 . 2 . 1 0 ) and then (5.2.9) reads, Ni = A gi e-βεi . // = M Z gi e-βεi, see below (5.2.11) So of we can somehow determine A and β, the problem is solved and the solution N = (N1,N2...Nm) is known. When a problem has a solution of the form (5.2.11) it is said to have Maxwell-Boltzmann statistics . The notion of absolute temperature Suppose εi > εj and consider from (5.2.11) that N i/Nj = (gi/gj) e-β(εi-εj). As β→+∞, Ni/Nj→ 0 and there are no particles in the upper state i co mpared to the lower state j. Basically, as β→+∞ all boson particles will inhabit only the lowest energy state (the "ground state") of a system. This is a situation one associates with absolute zero temperature T = 0. The system is "frozen" into its ground state. For fermions, the particles fill up the energy levels from the bottom up to some level, and above that level N i = 0 (frozen atoms, white dwarfs, neutron stars, etc.). On the other hand, as β → 0, we find that N i/Nj = gi/gj and particles are then distributed statistically acc ording to state degeneracy counts without regard to the energy level difference between the states. This is the situation at an infinitely high temperature T → ∞. So a simple way to define absolute temperature is T ≡ C/β where C is any positive constant. Measurements show that the T=0 situation exists when T C = - 273.15 , where T C is temperature in Centigrade units. In order to have the size of one de gree of absolute temperature units T be the same as that for T C, one must write T = T C + 273.15. The units of T are called K after Lord Kelvin, 1824-1907, aka William Thomson. Thus, the freezing point of water 0oC corresponds to T = 273.15K. Room temperature T F = 70oF = 21oC corresponds to T = 294K. In Section 5.4 below we shall show that the mean energy of a particle of ideal gas is given by u = (3/2)(1/β). Writing this using T ≡ C/β gives u = (3/2)(1/C)T. In effect, one can take a box of ideal gas at temperature T and measure the mean energy u, and this then determines the constant C. In practice this measurement of C is done indirectly, for example by measuring the speed of sound in argon gas, and C is now measured to about 7 decimal points of accuracy. Historically one deals with the inverse constant k = 39 1/C, and that constant k is known as Boltzmann's constant . The upshot of this discussion is that for systems in thermal equilibrium, one can associate the derived Lagrange multiplier β = - λ2 with absolute temperature T according to : β = 1/(kT) k = 1.380648 x 10 −23 J/K . // Joules/Kelvin (5.2.12) Notice from (5.2.11) or u = (3/2)(1/ β) that the units of β must be inverse energy, so kT has energy units. The factor e -βεi in (5.2.11) is referred to as the Boltzmann factor . The quantity S ≡ k lnΩ = kf is the entropy of a system, where k makes another appearance. These matters are discussed in any thermodynamics or statistical mechanics text. See for example Zemansky p 258-272 and equation (10-17). How the derived Lagrange multiplier values β and A are determined. From (5.2.11) the probability of a particle having energy εi is given by pi = Ni M = A M gi e-βεi . ( 5 . 2 . 1 3 ) Since Σipi = 1 one finds that 1 = A M [Σigi e-βεi] which then determines the derived Lagrange multiplier constant A (expressed in terms of β), A = M Σigi e-βεi . // = M Z , Z ≡ Σigie-βεi (5.2.14) The denominator in (5.2.14) is called the partition function Z for the system, so then A = M/Z. Using (5.2.11) for N i, the two constraint equations shown in (5.2.4) [ ΣiNi= M and ΣiNiεi= U] read A Σigie-βεi = M // N i = A gi e-βεi A Σiεig e-βεi = U . (5.2.15) Dividing these two equations and using (5.1.6) that u = U/M, one gets this self-consistent result, u = U M = Σigiεie-βεi Σigie-βεi = (A/M)Σigiεie-βεi (A/M)Σigie-βεi = Σipiεi 1 = Σipiεi = <εi> . (5.2.16) To find the derived Lagrange multiplier constant β = -λ2, write out the second equation of (5.2.15), Σigiεi e-βεi = U/A = (U/M)(M/A) = (u) ( Σigi e-βεi) or Σ i gi [εi-u] e-βεi = 0 . ( 5 . 2 . 1 7 ) 40 This is the equation one must solve for β. Recall that ε i, gi, and u are all known quantities. For the lowest energy states one will have εi < u ( recall that u = < εj> ) so [εi - u] will be negative. For high energy states [εi - u] will be positive. It is possible then to adjust β to bring these positive and negative contributions into balance to generate a 0 sum. One can see from (5.2.17) that β is a function of mean temperature u. For general values of ε i, equation (5.2.17) is a transcendent al equation which must be solved numerically for β. If it happens that all the ε i are ratios of integers, the equation can be written as a polynomial equation, but since those integers are likely to be large, the order of this polynomial is large and again one must do a numerical solution. Notice that the solution to our particle system problem has the same general form regardless of the number m of energy levels ε i. If we take m very large and make the energy levels εi be very closely spaced, they approach a continuum of energy values and the conclusions apply to a classical system. For small finite m, the discrete energy levels indicate a quantum mechanical system. Summary We have used the method of Lagrange multipliers to solve the Boltzmann extremum problem stated in either (5.1.4) or (5.2.4). Since the problem has two constraints, there are two Lagrange multipliers λ 1 and λ2 which we replaced with the numbers A ≡ eλ1 a n d β ≡ - λ2. We showed how β is determined by numerically solving (5.2.17), and then A is determined by (5.2.14). Then the solution (extremum) vector N = {Ni} is given by (5.2.11) which says N i = A gi e-βεi . What we have not shown is that this extremum solu tion is a maximum. Nor have we shown the fact that for large N i the microstate count Ω(N) has a very sharp peak at the solution value N, which then statistically forces a system to assume a value N very close to this solution N . We shall deal with these issues in the following section. 5.3 More details of the solution Recall fro m (5.2.3) the function for which we seek a stationary point, f = ln Ω = Σi=1m Ni [1 + lng i - lnNi] . (5.2.3) (5.3.1) The constraints Σ i=1mNi= M and Σ i=1Niεi= U can be regarded as two equation in the two unknowns N 1 and N2 , N 1+ N2 = M - Σ i=3mNi ε 1N1+ ε2N2 = U - Σi=3mεiNi , (5.3.2) which are easily solved to obtain, 41 N1 = + (ε2- ε1)-1[ ε2M - U + Σ i=3m (εi- ε2)Ni] = N 1(N3,... Nm) N2 = - (ε2- ε1)-1[ ε1M - U + Σ i=3m (εi - ε1)Ni] = N 2(N3,... Nm) . (5.3.3) Notice that, ∂N1 ∂Ni = + (ε2- ε1)-1 [ (εi- ε2) ] = εi-ε2 ε2-ε1 i = 3,4..m ∂N2 ∂Ni = - (ε2- ε1)-1 [ (εi- ε1) ] = - εi-ε1 ε2-ε1 i = 3,4..m . (5.3.4) Define function F of m-2 independent variables of N 3...Nm as follows : F(N 3,N4.....Nm) ≡ f(N1(N3,... Nm), N2(N3,... Nm), N3, .....Nm) . (5.3.5) F has the same value as f but has a differe nt functional form. We rewrite (5.3.1) as F = N 1[ 1 + ln(g 1) - lnN1] + N2 [ 1 + ln(g 2) - lnN2] + Σj=3m Nj [1 + ln(g j) - lnNj] (5.3.6) where recall that N 1 = N1(N3,... Nm) and N2 = N2(N3,... Nm). First Derivative Compute ∂ iF from (5.3.6) as follows, ∂F ∂Ni = ∂N1 ∂Ni [ 1 + ln(g 1) - lnN1] + N1 [- 1 N1 ] ∂N1 ∂Ni + ∂N2 ∂Ni [ 1 + ln(g 2) - lnN2] + N2 [- 1 N2 ] ∂N2 ∂Ni + Σj=3m { δi,j [ 1 + ln(g j) - lnNj] + Nj [ -1/Nj] δi,j } = ∂N 1 ∂Ni [ ln(g1) - lnN1] + ∂N2 ∂Ni [ ln(g2) - lnN2] + [ ln(g i) - lnNi] = εi-ε2 ε2-ε1 [- ln(N1/g1)] + (- εi-ε1 ε2-ε1 ) [- ln(N 2/g2)] + [- ln(N i/gi) ] // (5.3.4) = - εi-ε2 ε2-ε1 ln(N1/g1)] + εi-ε1 ε2-ε1 ln(N2/g2) - ln(Ni/gi) i = 3,4 .... m . (5.3.7) This result is valid for any values of variables N 3, N4....Nm. In order to evaluate ∂F/∂Ni at the stationary point of interest, recall from (5.2.9) that at the stationary point, ln(N i/gi) = λ1+ λ2εi i = 1,2...m . (5.2.9) 42 Inserting this into (5.3.7) one finds, adding factor 1 = ε2-ε1 ε2-ε1 to the last term, ∂F ∂N i |stat = - εi-ε2 ε2-ε1 (λ1 + λ2ε1) + εi-ε1 ε2-ε1 (λ1 + λ2ε2) - ε2-ε1 ε2-ε1 (λ1 + λ2εi) = 1 ε2-ε1 { - (εi- ε2) (λ1 + λ2ε1) + (εi- ε1) (λ1 + λ2ε2) - (ε2- ε1)(λ1 + λ2εi)} . (5.3.8) Inside the Curly Bracket all terms cancel, as Maple verifies with the final result ∂F ∂Ni |stat = 0 i = 3,4...m . (5.3.9) This is the expected result since the solution {N i} should be a regular stationary point of the unconstrained function F(N 3,N4.....NM). Second Derivative Start with (5.3.7) written this way, ∂F ∂N i = - εi-ε2 ε2-ε1 [ lnN1- lng1] + εi-ε1 ε2-ε1 [ lnN2- lng2] - [ lnN i- lngi] . (5.3.10) Then ∂ 2F ∂(Ni)2 = - εi-ε2 ε2-ε1 1 N1 ∂N1 ∂Ni + εi-ε1 ε2-ε1 1 N2 ∂N2 ∂Ni - 1 Ni i = 3,4...m = - ε i-ε2 ε2-ε1 1 N1 εi-ε2 ε2-ε1 + εi-ε1 ε2-ε1 1 N2 (-εi-ε1 ε2-ε1 ) - 1 Ni // (5.3.4) = - [ ( ε i-ε2 ε2-ε1 )2 1 N1 + ( εi-ε1 ε2-ε1 )2 1 N2 + 1 Ni ] = ∂i2F . (5.3.11) This says that the curvature of F in all independent directions N i for i = 3,4...m is negative and this is true for all vectors N = {Ni} for i = 3,4.... This is of course then true for the solution vector N. But for the 43 solution vector N we also have ∂f/∂Ni = 0 from (5.3.9) and therefore the solution is a maximum of F and therefore of Ω. Note that ∂2 iF is independent of the degeneracy g i (in our g i >> Ni limit of interest) . Finally, we may estimate the width of the peak of F and then of Ω. We follow the procedure outlined earlier in Section 3.4. Taking a variation only in a particular N i ( in the "i" direction), we approximate F( N + ΔNii^) ≈ F(N) + (∂ iF) ΔNi + (1/2) (∂i2F) (ΔNi)2 i = 3,4...m . (5.3.12) But at the solution point (∂F/∂Ni) = 0 as in (5.3.9) so ΔF = F( N + ΔNii^) - F(N) ≈ (1/2) (∂i2F) (ΔNi)2 . Thus |ΔNi| ≈ 2 |ΔF| |∂i2F| i = 3,4...m . (5.3.13) Suppose we are interested in finding ΔNi such that the Ω(N) drops to half its peak value. Then Ωhalf/Ωpeak = 1/2 ⇒ ln Ωhalf - ln Ωpeak = ln(1/2) = - ln(2) = F half - Fpeak |ΔF| = Fpeak - Fhalf = ln(2) = 0.7 . (5.3.14) Consider this toy example where Ω is taken to be a Gaussian function, (5.3.15) On the left one can see that the half point of Ω is about 0.7 the half point of ln Ω. The plot on the right has N = 300 and is more representative of what happe ns in statistical mechanics. The peak of Ω is usually a huge number (example in the next section) and ln Ω is barely visible if Ω and lnΩ are plotted on the same scale. But of course (5.3.14) is still valid. Using |ΔF| = 0.7 in (5.3.13) along with (5.3.11) for the curvature, one finds 44 |ΔNi| ≈ 1.4 |∂i2F| = 1.4 [ ( εi-ε2 ε2-ε1 )2 1 N1 + ( εi-ε1 ε2-ε1 )2 1 N2 + 1 Ni ] = 1.4 Ni ( εi-ε2 ε2-ε1 )2 (Ni/N1) + ( εi-ε1 ε2-ε1 )2 (Ni/N2) + 1 . (5.3.16) The fractional width is then |ΔN i| Ni ≈ 1.2 Ni 1 ( εi-ε2 ε2-ε1 )2 (Ni/N1) + ( εi-ε1 ε2-ε1 )2 (Ni/N2) + 1 . (5.3.17) To clarify the meaning of this fractional width, let us temporarily indicate by {N ¯k} the components of the solution vector N ¯ at which Ω(N1,N2 ...Nm) has its maximum value. Since we have taken N 1 and N2 as dependent variables (functions of all the other N i as in (5.3.3)), we can write Ω(N1,N2 ...Nm) = Ω(N3,N4....Nm) where Ω is used only to distinguish the functional form. Here N 3....Nm are all independent variables. Then consider Ω(N¯3,N¯4 ...Ni.....N¯m) where we treat N i as the only variable and we allow it to vary through N ¯i. Here is a sample plot illustrating the peak in Ω(N¯3,N¯4 ...Ni.....N¯m) versus N i , (5.3.18) Although the value at the peak might be a very large number (like 101500), the "narrowness" of the peak is really related to the ratio | ΔNi|/ N¯i. If N¯i ~ 1020, as might be typical for a macroscopic system, then |ΔNi|/ N¯i <~ 10-10 and the peak would be extremely narrow and sharp. Since this is true for any N i one selects, one would in this large-N ¯i case conclude that in an ensemble of a billion systems, there is only a miniscule variation in the solution vectors N ¯ for the different systems. This is so because for any system in the ensemble, the components N i of the solution vector N ¯ are extremely likely to be within a few ΔNi of the solution value N ¯i. The underlying assumption here is that a system is most likely to be found where its available microstate count is very large. 45 In the above analysis we took N 1 and N2 to be the two dependent variab les. Had we instead selected N a and Nb as the dependent variables (those remaining then being independent), we would get these versions of equations above (italics) by taking 1 →a and 2 →b : Na = + (εb- εa)-1[ εbM - U + Σ i≠a,b (εi- εb)Ni] = N a(Ni≠a,b) Nb = - (εb- εa)-1[ εaM - U + Σ i≠a,b (εi - εa)Ni] = N b(Ni≠a,b) (5.3.3) ∂Na ∂Ni = εi-εb εb-εa i = 1..m but i ≠ a,b ∂Nb ∂Ni = - εi-εa εb-εa i = 1..m but i ≠ a,b (5.3.4) ∂F ∂Ni |stat = 0 i = 1..m but i ≠ a,b (5.3.8) ∂2F ∂(Ni)2 = - [ ( εi-εb εb-εa )2 1 Na + ( εi-εa εb-εa )2 1 Nb + 1 Ni ] = ∂i2F (5.3.11) |ΔNi| Ni ≈ 1.2 Ni 1 ( εi-εb εb-εa )2 (Ni/Na) + ( εi-εa εb-εa )2 (Ni/Nb) + 1 . (5.3.17) (5.3.19) 46 5.4 Example: N particles in 3 states : a numerical example For this example we choose three energ y levels (m = 3) so the microstate count Ω is given by Ω(N1,N2,N3) = g1N1 N1! g2N2 N2! g3N3 N3! (5.1.3) (5.4.1) and correspondingly (using the Stirling approximation for N i!), f = ln Ω = Σ i=13 Ni [1 + ln(g i/Ni)] . (5.2.3) (5.4.2) Assume the three energy levels are ordered in this manner, where ε 1 is the lowest energy state, ε 1 < ε2 < ε3 and recall u = U/M and ε1 ≤ u ≤ ε3 . (5.4.3) We regard both N 1 and N2 as functions of N 3 where, according to (5.3.3), N1 = + (ε2 - ε1)-1[ ε2M - U + ( ε3 - ε2)N3] = N1(N3) N2 = - (ε2 - ε1)-1[ ε1M - U + (ε 3 - ε1)N3] = N2(N3) . (5.3.3) (5.4.4) It seems clear that one must have 0 ≤ N 1 ≤ M 0 ≤ N2 ≤ M . ( 5 . 4 . 5 ) Using the expressions (5.4.4) in these inequalities, and making use of (5.4.3), one finds that all four inequalities can be summarized as just two inequalities which we write as, N3min ≤ N3 ≤ N3max N 3min = max(u-ε2 ε3- ε2 M, 0) N 3max = u-ε1 ε3- ε1 M . (5.4.6) The above fact is tedious to derive and we state it only because it is used in the code below for the cosmetic purpose of restricting the N 3 range of some graphs. Recall from (5.2.10) that the de rived Lagrange multipliers are β ≡ - λ2 and A ≡ eλ1 . To find β we must solve (5.2.17), Σi=13 [ εi - u ] gi e-βεi = 0 . (5.2.17) (5.4.7) Then the other derived multiplier A is given by (5.2.14) 47 A = M Σigi e-βεi = M Z where Z ≡ Σigie-βεi . (5.2.14) (5.4.8) The components N i of the solution vector N are then given by (5.2.11), Ni = A gi e-βεi . i = 1,2,3 (5.2.11) (5.4.9) The curvature of f [= F as in (5.3.5)] at its peak is shown in (5.3.11), ∂2F ∂2N3 = - [ ( εi-ε2 ε2-ε1 )2 1 N1 + ( εi-ε1 ε2-ε1 )2 1 N2 + 1 N3 ] . (5.3.11) (5.4.10) The half-width of the Ω peak is then given by (5.3.16), ΔN3 ≈ 1.4 N3 ( εi-ε2 ε2-ε1 )2 (N3/N1) + ( εi-ε1 ε2-ε1 )2 (N3/N2) + 1 . (5.3.16) (5.4.11) Here then is our Maple code. We first enter Ω, (5.4.12) The quantities N 1 and N2 are then replaced as in (5.4.4) [ we use e i in place of εi] , (5.4.13) For a specific example we assume these values, ε 3 = 3 M = 1000 = number of particles ε2 = 2 u = 1.5 = average energy ε1 = 1 g 1 = g2 = g3 = 5000 = degeneracy . (5.4.14) Enter values and use them to compute N 3min and N3max and then to restate N 1(N3) and N2(N3), 48 (5.4.15) For this example one then has N 1 = 500 + N 3 ⇒ dN 1 = +1 dN 3 N2 = 500 - 2N 3 ⇒ dN 2 = -2 dN 3 ( 5 . 4 . 1 6 ) which we plot as follows with N 1 in red and N 2 in black, (5.4.17) We then enter the degeneracies and take a look at Ω in its numerical form, . (5.4.18) As is typical in statistical mechanics, at its peak the value of Ω is a very large number. Even with the relatively small number of particles M = 1000 and dege neracy g = 5000, the peak value is (see below) 49 Ωpeak = 5.71 x 101519 . ( 5 . 4 . 1 9 ) Maple is uncomfortable plotting numbers larger than about 1040 so we pre-scale Ω down by 101519 to make the scaled Ω be Maple-digestible. This scaling process in turn creates numbers smaller than 10-40 which are also rejected by the Maple plotter, so we use a Heaviside function to pin small numbers to 0 unless they are greater than the arbitrary va lue .01. Here then is a plot of scaled Ω(N1(N3),N2(N3),N3), ( 5 . 4 . 2 0 ) which may be compared with our generic plot of Fig (5.3.18). We broke the display command into two pieces: first create xx then plot xx. Replacing : with ; after the plot command gives one a list of numbers to be plotted, and from this list one can determine a reasonable scale factor. Once that is found, the plotting coordinates need not be displayed. The curve Ω(N 3) has the typical Gaussian (normal) shape which no doubt the energetic reader could s how results from the central limit theorem. A plot of f = ln( Ω) has a much smoother and broader shape, reminiscent of Fig (5.3.15), ( 5 . 4 . 2 1 ) 50 Recall from (5.4.10) that ∂2F ∂2N3 is everywhere negative, consistent with the cupping down seen in this plot. Next, we have Maple solve (5.4.7) for the Lagrange multiplier β, (5.4.22) With w = e -β the above "eq" equation (divided by 2500w) is just quadratic 3w2+w-1 = 0 which has solutions w ± = ( -1 ± 13)/6. Since β = -lnw, the w - solution gives the complex (and therefore illegal) value .26 - iπ shown above in (5.4.22), while the w + solution is β = -ln[( -1 + 13)/6] = .834. From this β value we compute the other Lagrange parameter A from (5.4.8), (5.4.23) Once these parameters are known, we use (5 .4.9) to display the solution vector N = {Ni}. (5.4.24) Since we set u < ε 2, we see here a normal distribution where high er energy states have fewer particles. Our final task is to use (5.4.11) determine the half-width of the Ω peak, (5.4.25) which says Δ N3 ≈ 7.5. This width and the solution N 3 = 116.2 can be verified from this blowup plot of the Ω peak, 51 ( 5 . 4 . 2 6 ) Notice that our Maple code does not use the Stirling approximation for N 3! in the Ω of (5.4.18) but provides a continuous interpolation of N 3! as Γ (N3+1). All the plots above should perhaps be indicated with dotted lines to emphasize the fact that the N 3 axis is really discrete integers. The peak values of Ω and f = ln(Ω) are found to be, ( 5 . 4 . 2 7 ) To see more directly where this huge Ω number is coming from, we have Maple compute the three factors in the expression (5.4.1) for Ω using the approximate solution values for N 1, N2, N3 shown in (5.4.20): Ω(N1,N2,N3) = g1N1 N1! g2N2 N2! g3N3 N3! (5.4.1) ( 5 . 4 . 2 8 ) 52 Consider now the following schematic diagram: (5.4.29) The red dots show the solution values for the integers (N 1,N2,N3). These would be the values obtained by averaging the values measured in 1 billi on systems in an ensemble. If we take N 3 to be the independent variable, it has a probability distribution (bell curve) shown in (5.4.20). For any point on this curve, N1(N3) and N2(N3) will also be displaced from their red-dot va lues. We roughly show this variation with little bell curves in the figure. For example, a partic ular system taken from the ensemble might have (N1,N2,N3) values indicated by the black dots in the figur e. The displacements shown are roughly these, dN3 = -50 ⇒ dN1 = -50 and dN 2 = +100 . // see (5.4.16) (5.4.30) In the above example we used u = 1.5 < ε2 = 2. If we instead use u = 2.5 > ε2 = 2, we find that the state populations are as shown above but in reverse order. Such an inverted situation corresponds to a negative absolute temperature T ( β = -.834 = 1/(kT) ) and cannot therefor e arise in thermal equilibrium. It does arise in a laser medium which is driven by some power source ("population inversion"). Geometric Interpretation The above example is similar to our Example 2 of Section 3.2 involving a hypersphere. In both problems there are three components, here N = (N 1,N2,N3) and there r = (x,y,z). In both problems there are two constraints, here ΣNi= N and ΣεiNi= U and there x=1 and y=1. In both problems we can imagine plotting the function of interest in E4, here f = Σi=13 Ni [1 + ln(g i/Ni)] and there f = 4 - x2- y2- z2 . In both problems the "intersection constraint surface" of Section 1 has dimension N-C = 3-2 = 1. For the hypersphere problem, that intersection constraint surface is the line which is the intersection of the x=1 and y=1 planes drawn in E3. In our current problem the intersection constraint surface is the line which is the intersection of the two planes N 1+N2+N3 = N and ε1N1+ε2N2+ε3N3 = U in E3. The first plane has 53 normal (1,1,1) and its closest approach to the or igin is distance N. The second plane has normal ( ε1,ε2,ε3) and its closest approach to the origin is distance U. Before discussing our current problem, it is useful to review the hypersphere Example 2. Recall that in Example 2 the "level surfaces" in E 3 are determined by K = f( r) = 4 - x2- y2- z2 for various values of K, and these surfaces are a set of concentric spheres. Here is some crude Maple code to display three of these spherical level surfaces : (5.4.31) We have restricted the z range to obtain a cutaway view of the concentric spheres. In this picture the constraints are the planes x = 1 and y = 1 which intersect in a line parallel to the z axis. The solution must lie on this line, and there is one point on the line where the gradients to the planes x=1 and y=0 are coplanar with the gradient (normal) of one of the con centric spheres. As we saw, and as is clear here, that point will lie at the z = 0 plane and has the value r = (1,1,0). This point corresp onds to the smallest radius concentric sphere touched by the constraint line and thus gives a maximum for f( r) = 4 - x2- y2- z2 subject to the constraints x=1 and y = 1. In the current problem, the level surfaces are determined by K = f(N ) = Σi=13 Ni [1 + ln(g i/Ni)] (5.4.32) for various values of K. For more familiarity, write this with (N 1,N2,N3) = (x,y,z) to get K = f ( r) = x[1+ln(g 1/x)] + y[1+ln(g 2/y)] + z[1+ln(g 3/z)] = (1+ lng 1) x + (1+ lng 2)y + (1+ lng 3)z - (xlnx+ylny+zlnz) . (5.4.33) 54 The two constraint planes are then x 1+x2+x3 = N and ε1x1+ε2x2+ε3x3 = U . (5.4.34) Using gi = 5000 as in our current example, we have Mapl e plot some level surfaces . First, we construct the function f( r) as follows : Then we set the g i values and restate f(x,y,z), (5.4.35) Finally we display level surfaces f(x,y,z ) = K for K = 1000, 1500, 2000 .... 8000, (5.4.36) Each level surface has the shape of a smooth scallop shell. The two constraint planes x 1+x2+x3 = N and ε1x1+ε2x2+ε3x3 = U intersect in a skewed straight line in this figure (the "intersecting constraint surface"). There is only one point on this line wh ere the constraint plane normals (1,1,1) and ( ε1,ε2,ε3) are coplanar with the normal to one of the scallop surfaces. For ( ε1,ε2,ε3) = (1,2,3) we found that the point is roughly (x,y,z) = (N 1,N2,N3) = (616,238,116) and this is the point which maximizes the function 55 f(N) = Σi=13 Ni [1 + ln(g i/Ni)] subject to the constraints ΣNi= 1000 and ΣεiNi= 1000*1.5. Eq. (5.4.27) shows that f max ≈ 3500. 5.5 Example: The Maxwell-Boltzmann Distribution The Sy stem In Sections 5.1-5.4 we dealt with a system having a finite number of discrete energy levels εi with degeneracy g i. In the current section we shall consider a s ituation which has an infinite number of energy levels which form a continuum whic h has no upper limit. However, we sh all treat this continuum as if it were a finely grained set of discrete energy levels so we can use the methods deve loped in Section 5.1. In the end, we shall apply quantum mechanics to show th at in fact the energy levels really are discrete. The system of interest is a small box of helium gas n ear room temperature. Our "particles" then are the individual gas atoms. The gas is treated as an "ideal gas" : atoms are indistinguishable (identical); there are no at-distance interactions between the atoms; th e atoms are small compared to the distances between them; collisions between atoms and at the walls of the box are "elastic"; no energy is stored in any kind of rotational or vibrational modes of the atoms. The at oms have "spin" 0 so they are bosons and therefore any number of them can be placed into any given energy state, as in our previous examples. The atoms are in a stable state of "thermal equilibri um" at temperature T. The box is sealed and insulated to prevent the transfer of particles or energy in or out. The box contains some total large fixed number of atoms M (to be computed below) and has some total energy U = Mu where u is the average energy of a particle. Connection to Sections 5.1 and 5.2 An atom of helium gas has kinetic energy ε i = (1/2)mv2 ≡ ε(v) . (5.5.1) This is the total energy of a particle having speed v = | v| and mass m. There is no potential energy because the particles don't interact and because we ignore grav ity which is only a miniscule effect for a small box of gas. In classical m echanics, this energy ε(v) is continuous because the speed v is a continuous variable. So ε(ν) = (1/2)mv2 represents the continuum of "energy le vels" for atoms in a box of helium. Let dv be some small but finite speed difference. The num ber of atoms with speed in this range is taken to be N(v)dv, and the count of available states in the ra nge we write as g(v)dv. The amount of energy in the atoms with speed between v and v+dv is ε(v)N(v)dv. We shall determine N(v) and g(v) below. To get a number to represent the degeneracy g i, we hop into velocity space where the axes are v x,vy,vz and we note that the volume in this space co rresponding to v in the range (v,v+dv) is 4 πv2dv, since this is the volume of a shell of radius v and thickness dv. As an artificial (at this point) construct, we shall 56 assume that velocity space is not continuous but is a 3D lattice of allowed values and that the volume of a lattice cube is V. Therefore, we write (at least for now), gi = [4πv2dv]/V = g(v)dv where g(v) = (4 πv2/V) . (5.5.2) All the g i states in this thin shell have energy εi = ε(v). Note that volume V is a tiny assumed cell size in velocity space, it is not the volume of our box of gas. The number N i of atoms having energy εi is now, N i = N ( v ) d v , ( 5 . 5 . 3 ) and then g i/Ni = [g(v) dv]/[N(v)dv] = g(v)/N(v) . (5.5.4) We have the same two constraints as in S ection 5.1, but here they take the form M = Σ iNi = ∫0 ∞ N(v)dv U = Σ iεiNi = ∫0 ∞ ε(v) N(v)dv . (5.5.5) As noted, the distribution of atoms into the energy states is described by N(v). The elastic collisions between atoms and the walls do not affect N(v) but me rely redirect an atom's velocity vector without changing the speed v. However, elastic collisions be tween atoms can and do alter velocities. Before a collision two atoms might have velocities v1 and v2, and after v'1 and v'2. But since the collisions are elastic one has E 1 + E2 = E'1 + E'2, so the total energy U is not affected (nor is M affected). It is true that a collision rearranges the two atoms' positions in the en ergy level diagram since in general one won't have E1 = E1' and E2 = E'2. But for every such rearrangement, the exact reverse rearrangement is equally likely in some other collision in the box (which contains perhaps 1020 atoms), so overall the inter-atomic collisions don't affect N(v). Recall from (5.1.3) that we had this microstate count, assuming g i >> Ni. Ω(N) = g1N1 N1! g2N2 N2! ....gmNm Nm! = Πi=1m giNi Ni! . (5.1.3) What happens to this expression for Ω when there are an infinite number of energy states which are spaced closely together? We might write g iNi = (g(v)dv)N(v)dv Ni! = [N(v)dv]! (5.5.6) 57 and one gets a product of an infinite number of terms, Ω(N) = Πn=0∞ (g(ndv))N(ndv)dv (N(ndv)dv)! . (5.5.7) Just as an infinite sum Σ nan can converge only if a n → 0, an infinite product Πnan can only converge if an→ 1. There is a developed theory of such infinite products, but we leave this topic to the interested reader. Here are a few examples of infinite products, https://en.wikipedia.or g/wiki/Infinite_ product ( 5 . 5 . 8 ) and one sees that a n → 1 in all cases. The object f = ln( Ω) shown in (5.2.3) is a little more tractable since it is merely an infinite sum, f = Σ i Ni [1 + ln(g i/Ni) ] (5.2.3) = Σn=0∞ N(ndv)dv [ 1 + ln{g(ndv)/N(ndv)}] . (5.5.9) We shall avoid dealing directly with Ω or f = ln( Ω) and instead just make use of results derived in Sections 5.1 and 5.2 for a finite number of energy levels. The Maxwell-Boltzmann speed distribution Our main result of interest is (5.2.11) which says that the components N i of the vector N which causes Ω to have a maximum are given by, Ni = A gi e-βεi , (5.2.11) where A ≡ e λ1 and β ≡ - λ2 are the derived Lagrange multipliers. Using (5.5.1,2,3) for εi, gi and Ni we write the above as {N(v)dv} = A {[4 πv2dv]/V} {exp(- βmv2/2) } or N(v) = 4 π(A/V) v 2 exp(-βmv2/ 2 ) . ( 5 . 5 . 1 0 ) 58 We thus arrive (with very little work) at the fam ous Maxwell-Boltzmann speed distribution. For now we ignore the constant 4 π(A/V) -- it will be dealt with later. The general shape of this distribution is as follows : (5.5.11) The quadratic left end of the plot arises from v2 where exp ≈ 1, while the right end comes from the decaying exponential which here is just exp(-v2). Mean Values Knowing the form (5.5.10) of N(v) allows us to compute various mean values. The first item of interest is the mean energy of a particle in the gas, <ε> = u = ∫0 ∞ (1/2)mv2 N(v)dv ∫0 ∞ N(v)dv = (m/2) ∫0 ∞ v4 exp(-βmv2/2)dv ∫0 ∞ v2 exp(-βmv2/2)dv . We can make use of the following integral ∫0 ∞ vn exp(-av2)dv = (1/2) Γ[(n+1)/2] a-(n+1)/2 // Spiegel 15.77 (5.5.12) with a = βm/2 to obtain <ε> = u = (m/2) Γ(5/2)(β m/2)-5/2 Γ(3/2)(β m/2)-3/2 = (m/2) (3/4)π (2/4)π (2/βm) = (3/2)(1/ β) = (3/2)kT , (5.5.13) a result we quoted earlier. This is an example of the "equipartition theorem" which says that each independent quadratic term in the energy expression e nds up getting (1/2)kT of energy in equilibrium. In our case we have three terms since ε(v) = (1/2)mv x2 + (1/2)mv y2 + (1/2)mv z2. 59 Next we calculate the mean speed v of a particle, vmean = <v> = ∫0 ∞ v N(v)dv ∫0 ∞ N(v)dv = ∫0 ∞ v3 exp(-βmv2/2)dv ∫0 ∞ v2 exp(-βmv2/2)dv = Γ(2)(β m/2)-2 Γ(3/2)(β m/2)-3/2 = 1 π /2 (βm/2)-1/2 = (2/π ) 2/βm = 8/βπm = 8kT/πm . (5.5.14) The root mean square speed is given by vrms ≡ <v2> = <(1/2)mv2> (1/2)m = 2 m <ε> = 2 m 3 2 kT = 3kT m . (5.5.15) Finally, the peak speed (most likely speed) occurs where the distribution slope is zero : ( 5 . 5 . 1 6 ) and the meaningful root is the second item, so vpeak = 2kT m . ( 5 . 5 . 1 7 ) We now summarize our conclusions about mean values, <ε> = u = (3/2)kT vpeak = 2kT m = 1.41 kT m = 1.00 v peak vmean = 8kT πm = 1.60 kT m = 1.13 v peak vrms = 3kT m = 1.73 kT m = 1.23 v peak ( 5 . 5 . 1 8 ) The three speed values are not very far apart. 60 Determination of the deri ved Lagrange Multiplier A Recall from (5.5.10) that. N(v) = 4 π(A/V) v 2 exp(-βmv2/ 2 ) . (5.5.10) The total number constraint requires that M = ∫0 ∞ N(v)dv = 4 π(A/V) ∫0 ∞ v2 exp(-βmv2/2) = 4π(A/V) (1/2)( π/2)(β m/2)-3/2 = 4π (A/V)(π /4) (β m/2)-3/2 . (5.5.19) Therefore 4π(A/V) = M (4/ π) (βm/2)3/2 and then N(v) = (4M/ π) (βm/2)3/2 v2 exp(-βmv2/ 2 ) ( 5 . 5 . 2 0 ) which agrees with Zemansky p. 159 Eq. (6-25). Notice that the lattice cell size V in velocity space does not appear in this result, so it applies whether V is finite or V →0 which is the classical limit. In our study of the discrete-state Boltzmann problem in Section 5.2 it was noted that N is the "solution vector" which has "components" N i = Agie-βεi. In the continuum problem the "solution vector" is the entire continuous function N(v) while the "components" are N(v)dv for specific values of v. The other constraint U = ∫0 ∞ ε(v) N(v)dv tells us what we already know from (5.5.13), U = ∫0 ∞ ε(v) N(v)dv = < ε> ∫0 ∞ N(v)dv = < ε>M = uM = (3/2) kT M . (5.5.21) We would like now to verify that N i << gi since (5.2.11) that N i = A g i e-βεi depends on this assumption. From (5.5.4) we must therefore show that N(v) << g(v). But so far we have g(v) = (4 πv2/V) from (5.5.2) where V was our assumed lattice cube size in v-space. Taking V →0 results in g(v) = ∞ and then certainly N(v) << g(v). We now visit the quantu m mechanics department to obtain the correct value for V. Quantum Mechanics and Determination of V To determine the quantum theory value of V in (5.5.2 ), we consider the problem of a single theoretical point particle in a cubic box of edge L. The Schrodinger equation says H^ψ(r) = εψ (r) for a single particle 61 wavefunction ψ, where the Hamiltonian operator H^ is just p ^2/2m. Here p^ = (-h⁄ /i)∇ (a key assumption of quantum theory), so p ^2 = - h⁄ 2∇2. The constant h⁄ is h/(2π) where h is Planck's constant, h = 6.62607 x 10−34 J s h⁄ =1.0545718 x 10−34 J s h⁄ = h/(2π) // Joule-sec . The Schrodinger equation then reads -(1/2m) h ⁄ 2∇2ψ = εψ . We seek a solution to this equation which vanishes at the 6 walls of a cubic b ox of edge L. This is so because ψ must be continuous and ψ = 0 everywhere outside the box since there is no probability | ψ|2 that the particle is outside the box. That solution is ψnxnynz(x,y,z) = K sin( πnxx/L)sin(πnyy/L)sin(πnzz/L) n i = 1,2,3... . (5.5.22) The constant K is determined by requiring that the pr obability of the particle being in the box is 1, ∫0 L dx ∫0 L dy ∫0 L dz |ψnxnynz(x,y,z)|2 = 1 . (5.5.23) Since ∫0 L dx sin2(πnxx/L) = (L/2) the above says K2(L/2)3 = 1 so K = (2/L)3/2. ( 5 . 5 . 2 4 ) In (5.5.22) setting n x = 0 results in ψ = 0 which cannot be normalized to 1 and is a non-solution. A solution with n x = -2 is minus of the solution with n x = +2 and so these two solutions are really the same solution. Solutions must be linearly independent to be separately counted. That is why we have written ni = 1,2,3... . The Schrodinger equation then says p^2 2m ψ = εψ ⇒ - h⁄ 2 2m ∇2ψ = εψ ⇒ - h⁄ 2∇2ψ = 2mεψ ⇒ - h⁄ 2 ( [ -(πnx/L)2 - (πny/L)2 - (πnz/L)2 ] )ψ = 2mεψ // using (5.5.22) for ψ or (πh ⁄ /L)2 [ nx2+ny2+nz2] = 2mε or ε = (2m) -1 (πh⁄ /L)2 [ nx2+ny2+nz2] n i = 1,2,3... . (5.5.25) We arrive at the interesting conclusion that the ener gy spectrum is not continuous but is discrete. It is quantized. Setting ε = (1/2)mv2 we find that v2 = (πh⁄ /mL)2 [ nx2+ny2+nz2] ( 5 . 5 . 2 6 ) so the speed variable is also discrete, not continuous. 62 If we now think of velocity space with v = (vx,vy,vz) then the above says vx2 = (πh⁄ /mL)2nx2 v y2 = (πh⁄ /mL)2ny2 v z2 = (πh⁄ /mL)2nz2 (5.5.27) so then v = (v x,vy,vz) = (πh ⁄ /mL)(nx,ny,nz) . ( 5 . 5 . 2 8 ) In the space (n x,ny,nz) the lattice cube size is 1. Therefore, in v-space this cube size is (π h ⁄ /mL)3. We have now found a value for the parameter V in (5.5.2) : gi = [4πv2dv]/V = g(v)dv where g(v) = (4 πv2/V) . (5.5.2) V = ( πh ⁄ /mL)3 . ( 5 . 5 . 2 9 ) But now we make a correction. Since v = (πh ⁄ /mL)(nx,ny,nz) and since n i = 1,2,3..., when we count states we should only include the first octant shell of v-space where v x,vy,vz are all positive. Thus, we correct (5.5.2) as follows: gi = [(1/8)4 πv2dv]/V = g(v)dv where g(v) = (1/8)(4 πv2/V) . (5.5.2)corr V = ( πh ⁄ /mL)3 . ( 5 . 5 . 3 0 ) The degeneracy function is then g(v) = πv2/2V = (1/2) πv2 (πh ⁄ /mL)-3 = (1/2)π v2 (h/2mL)-3 = (1/2)πv2(2mL/h)3 = 4 π(mL/h)3v2 . ( 5 . 5 . 3 1 ) As a check on this result, integrating it from 0 to v' gives (4π/3)(mL/h) 3(v'2)3/2 = (4π/3)(mL/h)3(2ε'/m)3/2 = (4π/3)(L/h)3(2ε'm)3/2 and this agrees with Zemansky p. 272 Problem (10-2) (answers on p 645). At this point we have determined that, 63 N(v) = (4M/ π) (βm/2)3/2 v2 exp(-βmv2/2), β = 1/(kT) (5.5.20) g(v) = 4 π(mL/h)3v2 (5.5.31) (5.5.32) where m = m He is the mass of a helium atom. Numeric Values for N(v) and g(v) and g(v)/N(v) We now seek numeric values for N(v) and g(v) for a 1 cm3 box of helium at room temperature. Before entering constants, we create expressions for N, g, g/N and v m (vmean from (5.5.14)) as follows, where the Maple v -> construct is used to define a function of v : (5.5.33) Detail: The strange "unapply" command causes g_over_N(v) to be a function of v, namely g(v)/N(v). If we don't do it this way, g_over_N(0) reports a divide by 0 error at v = 0 since N(0) = 0. Notice that the g/N ratio takes its smallest value at v = 0 where exp( βm Hev2/2) = 1. Since we want to show that g/N >> 1, showing this at v = 0 guarantees that g/N >> 1 for all v. Our first task is to compute the number M of helium atoms in the box. This can be found from the ideal gas law PV = NkT which applied to our situation says PL 3 = MkT. We display units as if they were Maple variables, which allows a "dimension check" for all our results. 64 (5.5.34) So there are M ≈ 2.5 x 1019 helium atoms in our 1 cm3 box. Next compute the derived Lagrange multiplier value β = -λ 2 = 1/(kT), and take a look at the basic energy scale for a helium atom kT: (5.5.35) Since it takes 19.8 eV to excite helium out of its elect ronic ground state, and since (3/2)kT = .038 eV, the elastic collision aspect of our ideal gas assumption is well justified at room temperature. Next, set the mass m He of a helium atom, (5.5.36) 65 Recall that g i = g(v)dv and N i = N(v)dv so we expect g and N to ha ve units of inverse velocity. Here then are the numeric values for v m, N(vm) and g(v m) : (5.5.37) The mean helium atoms are moving along at a respectable clip, 1243 m/sec (2,780 mph). Of course they don't go far between collisions. The ratio g(v m)/N(vm) is thus about 106, which justifies the assumption that g(v) << N(v) and g i << N, at least at the mean velocity v m. As noted above, the lower bound for this ratio occurs at v = 0 where we have, (5.5.38) and even here we have g(v) >> N(v). Here is plot of g(v)/N(v) from v = 0 up to 4 times the mean velocity, (5.5.39) 66 One can view the first octant of the velocity sphere as consisting of a large number of tiny cubes (size V) which are the available quantum states. Only a tiny frac tion of these cubes are occupied by particles. For example, in the shell of the sphere at speed v = 2200 m/sec, the above graph shows that only about 1 in 10 7 cubes is occupied, all the rest are empty. The o ccupancy here is much sparser than in our 3-level example of Section 5.4 where the ε 2 state had N 2 = 268 out of g 2 = 5000 states occupied. The ensemble The smooth N(v) distribution shows the ensemble average for N i = N(v)dv at various speeds. If one were to examine a particular system in an ensemble of one billion 1 cm3 boxes of helium, and if one were to divide the distribution into a set of vertical strips each having a large number of atoms, one would measure a set of N i values that do not quite match the distribu tion curve. For each such strip there is a "bell curve" of the type shown in Fig (5.3.18) which has a fractional half-width on the order of 1/ Ns where N s is the (large) number of atoms in a strip. We illustrate this idea in the following drawing, (5.5.40) which is analogous to Fig (5.4.29) for the 3-level system of Section 5.4. One could display on this picture a set of red dots right on the curve, and black dots within the bell curve at each v (we drew one pair of such dots only), as discussed relative to (5.4.29). In the current situation, if one were to dramatically increase the number of strips to drive N s down to a much lower number, the bell curve widths would increase dramatically, and a pattern of black dots if connected by lines would have a very noise-like appearance, hardly resembling the red curve. Summary Applying the discrete-energy-level result N i = Agie-βεi developed in (5.2.11) (using the method of Lagrange multipliers) we have obtained a set of resu lts concerning the nature of a small box of helium atoms at room temperature. The distribution of particle speeds has the Maxwell-Boltzmann shape shown in Fig (5.5.11) and the corresponding expression for the speed distribution N(v) is given in (5.5.32). The 67 distribution of energy level degeneraci es is given by g(v) also in (5.5.32). We used quantum theory only to determine the tiny lattice size V in v-space which cau ses the speed and energy of a helium atom to be quantized. For this specific example we showed that g i >> Ni thus justifying the use of (5.2.11). We showed along the way that the average helium atom en ergy is (3/2)kT and hinted how one might derive the equipartition theorem of which this is an example. Expressions were found for the mean, rms and peak speeds of the N(v) distribution. 68 Section 6. Matrix proof of the Method of Lagrange Multipliers In Section 2 it was shown that the Method of Lagrange Multipliers (Theorem 2) follows directly from Theorem 1, and is really just a restatement of Theorem 1 which we quote: Theorem 1: A point r is a "stationary point for f( r) subject to constraints a i(r) = 0" ⇔ (1.24) (a) a i(r) = 0 for i = 1,2...C (point r must satisfy all the constraints ) (b) There exist C constants λi such that ∇f(r) + Σi=1C λi ∇ai(r) = 0 . Our proof of Theorem 1 was entirely geometric in natu re. We talked about surfaces, normals to surfaces, tangent spaces, perp spaces, the intersection constraint surface, and so on. In this section, we shall rederive Theorem 1 using a certain "R matrix", though some geometry still appears. The reader should regard this simply as an alternative derivation of Theorem 1. As in Section 1, a certain amount of background mate rial is necessary to prov ide a framework for the proof. We define S = C+1 so the number of constraint functions is C = S-1, and we rename the constraint functions to be a,b,c....q instead of a 1, a2, a3....aC. The purpose of this renami ng is to avoid having double subscripts on a constraint function, one identifying it and one indicating a partial derivative. 6.1. The R matrix and its Rank We shall operate in N dime nsional Euclidean space EN with coordinates r = (x1, x2 ... xN). We seek stationary point s of a real function f( r) subject to S-1 constraints a( r) = 0, b(r ) = 0, ... q( r) = 0 where we use q generically to represent the last c onstraint function. For example, if there are two constraints, then S = 3 and the two constraint equations are a( r) = 0 and b( r) = 0. So here is the problem: find the stationary points of f( r) subject to constraints a( r)=0, b( r)=0, c( r)=0, .... q( r)=0 . (6.1.1) 1 2 3 S-1 Construct the following matrix of partial derivatives, where f i(r) ≡ ∂if(r) = ∂ f(r)/∂xi. That is, the subscript on f indicates which of the arguments x 1,x2, x3....xN the partial is with respect to. So here is that matrix, which we shall call R : (6.1.2) For example, f 2 ≡ ∂f(x1,x2, x3....xN)/∂x2. On the right we abbreviate each row of functions in bold , indicating a row vector like f = ∇f. 69 The R matrix has S rows and N columns. We shall only be interested in N ≥ S which is the same as saying N ≥ C+1 or N > C. Recall from the discussion below (2.8) that if N ≤ C, the problem is overconstrained and one cannot seek stationary points where df = 0. If N = S, the R matrix is square, otherwise for N > S it has more columns than rows, so it has a horizontal band shape. The maximum possible rank of the above matrix is S since this is the smaller of the number of rows and columns (see below). S would then be the "full rank" of R. Since all the functions are functions of r , we really have R( r). As the point r is varied, the elements of the matrix R( r) generally change. If S = N, the matrix is square and we can talk about the determinant det(R). In this case there is only one SxS "submatrix" and it is the entire matrix R. If S < N, the matrix is wider than it is tall. In this case, we can talk about square submatrices within R which are of shape SxS. For example, for N = 6 and S = 4 we have three constraints and the R matrix is this, (6.1.3) Here we have outlined two 4x4 submatrices in red and blue. The columns of a submatrix need not be contiguous, so the three green boxes indicate anothe r 4x4 submatrix. In this example, the number of submatrices is given by the binomial factor (6,4) -- pick a committee of 4 columns from 6 candidates. So we have shown only 3 out of (6,4) = 6!/(4!2!) = 15 possible 4x4 submatrices. We can indicate a submatrix using this notation: red = (f 1,f2,f3,f4) blue = (f 2,f3,f4,f5) green = (f 1,f2,f4,f6) . (6.1.4) That is to say, we abbreviate the submatrix by sta ting only the values in the first row of the submatrix. Since our main interest will be whether or not the determinants of these 4x4 submatrices vanish, the ordering of the columns within a submatrix is not important. For the general R matrix shown in (6.1.2) with N columns and S rows, there are (N,S) possible submatrices of size SxS. A particular submatrix is then indicated by a sequence of S subscripted f values. We now quote a few facts from linear algebr a concerning the rank of a matrix: Definitions : A minor of matrix A is the determinant of any square submatrix of A obtained by crossing out some rows and columns The rank r of an m x n matrix A is the dimension of the largest non- vanishing minor within A. [Shilov 1.92 ] Thus, r ≤ min(m,n). (6.1.5) 70 Definitions : [Shilov 2.21] Consider an m x n matrix A. Let r1, r2, ...rN be an arbitrary selection of N rows (row vectors) from the matrix A. This set of N rows is linearly dependent if one can find a set of constants k i at least one of which is non-zero such that k1 r1 + k2 r2 + .... + k N rN = 0 . ( 6 . 1 . 6 ) There could be several such equations, or there could be only one. This generally means that at least one row can be expressed as a linear combination of the other rows. A legal such equation would be r2 = 0 in which case k 2 = 1. One could then say that r2 is a linear combination of the other rows with coefficients all zero. In any event, if one or more rows are all zeros, a set of N rows are linearly dependent. If no such equation (6.1.6) exists, then the set of N rows is linearly independent . For example, the set of vectors r1 = (1,0) and r2 = (0,1) is linearly independent. The above discussion also applies to columns. Take "row" → "column" and ri → ci . Facts : The rank r as defined in (6.1.5) has these prope rties for any m x n matrix A: (6.1.7) (a) the rank r equals the number of linearly independent columns of A [Shilov 3.12 ] (b) the rank r equals the number of linearly independent rows of A [Shilov 3.13 ] In Appendix A below we provide our own proof of the claims of (6.1.7). With the above as background, the following theorem is fundamental to our matrix proof of Theorem 1 : Theorem 3 . Point r is a "stationary point of f( r) subject to the constraints a=0,b=0...q=0" ⇔ rank[R( r) ] < S . ( 6 . 1 . 8 ) But, according to (6.1.7b), rank[R( r)] < S ⇔ the rows of R shown in (6.1.2) are linearly dependent since there are S rows in R. Therefore rank[R( r)] < S ⇔ ∇f(r) + λ 1∇a(r) + λ2∇b(r) + ... + λ S-1∇q(r) = 0 for some λi. Thus Theorem 3 ⇔ Theorem 1 of (1.24). But we know from Section 2 that Theorem 1 ⇔ Theorem 2, so proving Theorem 3 provides an alternate proof of the Method of Lagrange Multipliers. According to the rank definition (6.1.5), rank(R) < S ⇔ all SxS submatrices must have zero determinant (6.1.9) If, for example, the red submatrix in Fig 1.3 had a non-zero determinant, then rank(R) = 4 = S. The point of Theorem 3 is that the solution point r of the constrained stationary point problem must be a point where the R( r) matrix drops below full rank. We shall prove Theorem 3 in the ⇒ direction in Section 6.2 by showing that, if r is a stationary point, then all t hose SxS submatrices discussed above have zero determinant and thus rank(R) < S. If one is searching for candidate solutions r to the stationary point problem, one can restrict one's search to points r for which rank[R( r)] < S. The proof of Theorem 3 in the ⇐ direction is fairly simple. If rank(R) < S, then the rows of R are linearly dependent by (6.1.7) and then ∇f(r) + λ 1∇a(r) + λ2∇b(r) + ... + λS-1∇q(r) = 0. This implies that df( r) = ∇ f(r) • dr = – λ1[∇a(r) • dr] – λ2 [∇b(r) • dr] – ...... – λS-1 [∇q(r) • dr] (6.1.10) 71 for any d r. Restricting to d r which don't violate any of the constrai nts, we know that each term on the right vanishes. For example, ∇a(r) • dr = 0 because ∇a(r) is normal to the surface a(r ) = 0, see (1.2). Since all the terms on the right above then vanish, we get df( r) = 0 and r is a stationary point of f by our definition (1.18) of a constrained stationary point. Comment on the R matrix. Imagine a general transformation from x-space x = (x1, x2, ....xN) to u-space with coordinates u = (u1, u2, ... uS). One might write such a transformation as u = F(x) where F:EN→ ES. Although this transformation in general is non-linear, in a tiny neighborhood of point x in x-space (and the corresponding point u in u-space), the transformation will be linear if certain conditions are met. That local linear relation is then described by d u = Rd x where R is a matrix (N co lumns, S rows) whose matrix elements are R ij = ∂ui/∂xj. If we think of f = u 1(x1, x2, ....xN) R 11 = ∂u1/∂x1 = ∂f/∂x1 = f1 R 12 = f2 etc. a = u 2(x1, x2, ....xN) R 21 = ∂u2/∂x1 = ∂a/∂x1 = a1 R 22 = a2 etc b = u 3(x1, x2, ....xN) ... q = u S(x1, x2, ....xN) ( 6 . 1 . 1 1 ) then the matrix shown in (6.1.2) is exactly the R matrix for this general transformation u = F(x). If S < N, then the R matrix is not square, the relation du = Rd x is a "projection" of a larger space into a smaller one, and the equation therefore cannot be inverted. This notation "R" is used extensively in our Tensor Analysis document, see in particular Chapter 2, though only invertible transformations F:EN→ EN are considered there. Sometimes R is referred to as "the differential" of a transformation. 6.2. Proof of Theorem 3 ( ⇒) Theore m 3. Point r is a "stationary point of f( r) subject to the constraints a=0,b=0...q=0" ⇔ rank[R( r) ] < S . (6.1.8) (a) Proof for N = 6 and S = 4 The more general the proof is made, the less clear it becomes due to the notational baggage that must be added. Therefore we shall prove Theorem 3 for the specific case shown above in Fig (6.1.3) where N = 6 and S = 4 (3 constraints). The reader should have no trouble following each step for this sample case. When we reach the proof conclusion for the sample case, we start over again doing the general case and compare each step to the steps of this sample case. Although a lot of e quations appear below (and a lot of column inches are consumed), everything is really quite simple. So, for our sample case we start off with this R ma trix, again with bolded row vectors on the right, f 1 f2 f3 f4 f5 f6 f a 1 a2 a3 a4 a5 a6 a R = b 1 b2 b3 b4 b5 b6 = b (6.2.1) c 1 c2 c3 c4 c5 c6 c 72 We want then to find an stationary point of f(x 1, x2, x3, x4, x5, x6) subject to these constraints: a(x 1, x2, x3, x4, x5, x6) = 0 b(x1, x2, x3, x4, x5, x6) = 0 c(x1, x2, x3, x4, x5, x6) = 0 . ( 6 . 2 . 2 ) In practice, a person solving this problem might try to eliminate variables. For example, somehow with enough labor one should be able to "use up" the three c onstraint equations to eliminate three of the six variables x i. There are of course (6,3) ways to select three variables for elimination. For starters, imagine that we have used up the 3 constraints to eliminate variables x 1, x2 and x3. Then we write a(x 1, x2, x3, x4, x5, x6) = 0 x 1 = X1(x4, x5, x6) b(x1, x2, x3, x4, x5, x6) = 0 ⇒ x2 = X2(x4, x5, x6) c(x1, x2, x3, x4, x5, x6) = 0 x 3 = X3(x4, x5, x6) (6.2.3) where X 1, X2 and X3 are three resulting functions. This is a "theoretical" elimination, one does not actually have to do the process, one need only underst and that in principle it could be done and the functions Xi could in principle be found. See Appendi x C for an elaboration of this point. Now rewrite the function f and the three constraint f unctions in this manner, substituting for example the function X1 for x1 everywhere it appears, f(X1(x4, x5, x6), X2(x4, x5, x6), X3(x4, x5, x6), x4, x5, x6) ≡ F(x4, x5, x6) a(X1(x4, x5, x6), X2(x4, x5, x6), X3(x4, x5, x6), x4, x5, x6) ≡ A(x4, x5, x6) b(X1(x4, x5, x6), X2(x4, x5, x6), X3(x4, x5, x6), x4, x5, x6) ≡ B(x4, x5, x6) c(X1(x4, x5, x6), X2(x4, x5, x6), X3(x4, x5, x6), x4, x5, x6) ≡ C(x4, x5, x6) . (6.2.4) In this way we have defined four new functions F,A,B,C each of the 3 variables shown (those that were not eliminated). Next, compute the first partial derivatives of F using the chain rule. For example ∂F ∂x4 = ∂f ∂x1 ∂X1 ∂x4 + ∂f ∂x2 ∂X2 ∂x4 + ∂f ∂x3 ∂X3 ∂x4 + ∂f ∂x4 . (6.2.5a) In our compact notation where F i means ∂F/∂xi this can be restated as F4 = f1X1 4 + f2X2 4 + f3X3 4 + f4 . (6.2.5b) Doing this for each of the non-eliminated variables one gets, F4 = f1X1 4 + f2X2 4 + f3X3 4 + f4 F5 = f1X1 5 + f2X2 5 + f3X3 5 + f5 F6 = f1X1 6 + f2X2 6 + f3X3 6 + f6 . (6.2.6) 73 Notice that the subscripts on the leftmost three f i correspond to those of the eliminated variables i = 1,2,3. Now, since we are seeking an extremal point of F(x 4, x5, x6) with no constraints at all (the three initial constraints have all been absorbed in the process of eliminating three variables), we require that (x 4, x5, x6) be a "critical point" for the functio n F, which means we require that F 4 = F5 = F6 = 0. Then the above equations become, f1X1 4 + f2X2 4 + f3X3 4 + f4 = 0 f1X1 5 + f2X2 5 + f3X3 5 + f5 = 0 f1X1 6 + f2X2 6 + f3X3 6 + f6 = 0 . (6.2.7) Now apply the same process to the three constraint functions A,B,C. First compute the derivatives, and then set the derivatives to zero. One does this for A(x 4, x5, x6), for example, because the constraint condition a = 0 says that that A(x 4, x5, x6) = 0 = a constant, for all values of x 4, x5, x6 , so all first derivatives must vanish. Doing this for f unction A then gives these three equations, a1X1 4 + a2X2 4 + a3X3 4 + a4 = 0 a1X1 5 + a2X2 5 + a3X3 5 + a5 = 0 a1X1 6 + a2X2 6 + a3X3 6 + a6 = 0 . ( 6 . 2 . 8 ) We have just taken the previous equations and replaced f → a everywhere. Then do this for B and C. One ends up then with this set of 12 equations: f 1X1 4 + f2X2 4 + f3X3 4 + f4 = 0 f1X1 5 + f2X2 5 + f3X3 5 + f5 = 0 f1X1 6 + f2X2 6 + f3X3 6 + f6 = 0 a 1X1 4 + a2X2 4 + a3X3 4 + a4 = 0 a1X1 5 + a2X2 5 + a3X3 5 + a5 = 0 a1X1 6 + a2X2 6 + a3X3 6 + a6 = 0 b 1X1 4 + b2X2 4 + b3X3 4 + b4 = 0 b1X1 5 + b2X2 5 + b3X3 5 + b5 = 0 b1X1 6 + b2X2 6 + b3X3 6 + b6 = 0 c 1X1 4 + c2X2 4 + c3X3 4 + c4 = 0 c1X1 5 + c2X2 5 + c3X3 5 + c5 = 0 c1X1 6 + c2X2 6 + c3X3 6 + c6 = 0 . ( 6 . 2 . 9 ) Here there are S = 4 equation groups, and each group h as N-S+1 = 6-4+1 = 3 equations (this last 3 is the number of non-eliminated variables). Next reorder these equations into three groups of four, for example taking the first equation in each group above to make the first new group below, 74 f1X1 4 + f2X2 4 + f3X3 4 + f4 = 0 a1X1 4 + a2X2 4 + a3X3 4 + a4 = 0 b1X1 4 + b2X2 4 + b3X3 4 + b4 = 0 c1X1 4 + c2X2 4 + c3X3 4 + c4 = 0 f 1X1 5 + f2X2 5 + f3X3 5 + f5 = 0 a1X1 5 + a2X2 5 + a3X3 5 + a5 = 0 b1X1 5 + b2X2 5 + b3X3 5 + b5 = 0 c1X1 5 + c2X2 5 + c3X3 5 + c5 = 0 f 1X1 6 + f2X2 6 + f3X3 6 + f6 = 0 a1X1 6 + a2X2 6 + a3X3 6 + a6 = 0 b1X1 6 + b2X2 6 + b3X3 6 + b6 = 0 c1X1 6 + c2X2 6 + c3X3 6 + c6 = 0 . (6.2.10) Now there are N-S+1 = 3 equation groups , and each group has S = 4 equations. Next, write each of these three equation sets as a matrix equation, ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞f1 f2 f3 f4 a1 a2 a3 a4 b1 b2 b3 b4 c1 c2 c3 c4 ⎝⎜⎛ ⎠⎟⎞ X1 4 X2 4 X3 4 1 = (f 1,f2, f3, f4) ⎝⎜⎛ ⎠⎟⎞ X1 4 X2 4 X3 4 1 ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞f1 f2 f3 f5 a1 a2 a3 a5 b1 b2 b3 b5 c1 c2 c3 c5 ⎝⎜⎛ ⎠⎟⎞ X1 5 X2 5 X3 5 1 = (f 1,f2, f3, f5) ⎝⎜⎛ ⎠⎟⎞ X1 5 X2 5 X3 5 1 ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞f1 f2 f3 f6 a1 a2 a3 a6 b1 b2 b3 b6 c1 c2 c3 c6 ⎝⎜⎛ ⎠⎟⎞ X1 6 X2 6 X3 6 1 = (f 1,f2, f3, f6) ⎝⎜⎛ ⎠⎟⎞ X1 6 X2 6 X3 6 1 . (6.2.11) On the right we show our abbreviated notation for a matrix, just showing the first row. Now we claim that the three matr ices shown must have zero determinant! Suppose the first matrix equation had a non-zero determinant. It could then be inverted (see (B.4.8)) to give ⎝⎜⎛ ⎠⎟⎞ X1 4 X2 4 X3 4 1 = (f 1,f2, f3, f4)-1 ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 . (6.2.12) 75 But this produces a blatant contradiction that 1 = 0 (not to mention the other rows), and therefore we must have det (f 1,f2, f3, f4) = 0 det (f 1,f2, f3, f5) = 0 det (f 1,f2, f3, f6) = 0 or det (f 1,f2, f3, fm) = 0 m ≠ 1,2,3 . (6.2.13) Therefore, we have shown that three of the SxS = 4x 4 submatrices of Fig (1.3) vanish. But there are (6,4) = 15 such submatrices, and in order to show that rank(R) < S, we have to show that all 15 SxS determinants vanish. But that can be demonstrated by starting over several times and each time c hoosing to eliminate three different variables. Notice in the det = 0 equations above that the first three f i values correspond to the three eliminated variables x 1, x2 and x3, Had we instead chosen to eliminate x 1, x2 and x4, we would have obtained, det (f 1,f2, f4, f3) = 0 det (f 1,f2, f4, f5) = 0 det (f 1,f2, f4, f6) = 0 or det (f 1,f2, f4, fm) = 0 m ≠ 1,2,4 . (6.2.14) More generally, had we chosen to eliminate variables x i, xj and xk ( i ≠ j ≠ k), we would have obtained det (f i,fj, fk, fm) = 0 m ≠ i,j,k . (6.2.15) As i,j,k range over all possible values 1,2,3,4,5,6, we clearly hit all possible 4x4 subdeterminants, and thus we have shown that they all vanish, and therefore rank(R) < S, QED. Lest it be overlooked, one must keep in mind that r = (x1, x2, x3, x4, x5, x6) must be a stationary point of f( r) subject to the constraints. And then the c onclusion is that, for such a solution point r, one has rank[R( r)] < S. To summarize, one theoretically eliminates S-1 of the N variables to create the functions F,A,B,C of (6.2.4). One then looks fo r the normal (unconstrained) critical points of F where all partials vanish. Since the constraints all have the form A = 0, B = 0, C = 0, their partials vanish as well. Writing all these equations in matrix form (for various sets of eliminat ed variables) leads to the conclusion that all the SxS minors of matrix R vanish so that rank(R) < S. (b) Proof for general S ≤ N To generalize now t o general N and S, we restate some of the previous equations (adding a prime to the equation number) using a more generalized notation. The reader w ill note how even our simple notation becomes rather unpleasant. 76 We start off with this R matrix which has N columns and S ≤ N rows. The S-1 constraint functions are a,b,....q. Remember that f 2 means ∂ f/∂x2. f 1 f2 f3 ... f N f a 1 a2 a3 ... a N a R = b 1 b2 b3 ... b N = b (6.2.1)' .... ... q 1 q2 q3 ... q N q Here are the S-1 constraint equations, a(x 1, x2, ...xN) = 0 // S-1 constraint equations b(x1, x2, ...xN) = 0 ... q(x 1, x2, ...xN) = 0 . ( 6 . 2 . 2 ) ' Eliminate (theoretically) the first S-1 variables x 1,x2....xS-1 using the S-1 constraint equations: a(x 1, x2, ...xN) = 0 x 1 = X1(xS, xS+1 ... xN) b(x1, x2, ...xN) = 0 ⇒ x2 = X2(xS, xS+1 ... xN) ... ... q(x 1, x2, ...xN) = 0 x S-1 = XS-1(xS, xS+1 ... xN) . (6.2.3)' Substitute to get the F,A,B...Q equations (a set of S equations), f(X 1(xS, xS+1 ... xN), X2(xS, xS+1 ... xN), .. XS-1(xS, xS+1 ... xN) , xS, xS+1 ... xN) ≡ F(xS, xS+1 ... xN) a(X1(xS, xS+1 ... xN), X2(xS, xS+1 ... xN), .. XS-1(xS, xS+1 ... xN) , xS, xS+1 ... xN) ≡ A(xS, xS+1 ... xN) b(X1(xS, xS+1 ... xN), X2(xS, xS+1 ... xN), .. XS-1(xS, xS+1 ... xN) , xS, xS+1 ... xN) ≡ B(xS, xS+1 ... xN) ... q(X 1(xS, xS+1 ... xN), X2(xS, xS+1 ... xN), .. XS-1(xS, xS+1 ... xN) , xS, xS+1 ... xN) ≡ Q(xS, xS+1 ... xN) . ( 6 . 2 . 4 ) ' Now take first derivatives as shown in (6.2.5) to get this generalized version of (6.2.6), FS = f1X1 S + f2X2 S + f3X3 S + ... + f S-1XS-1 S + fS FS+1 = f1X1 S+1 + f2X2 S+1 + f3X3 S+1 + ... + f S-1XS-1 S+1 + fS+1 FS+2 = f1X1 S+2 + f2X2 S+2 + f3X3 S+2 + ... + f S-1XS-1 S+2 + fS+2 ... F N = f1X1 N + f2X2 N + f3X3 N + ... + f S-1XS-1 N + fN . // N-S+1 equations (6.2.6)' Requiring (x S, xS+1 ... xN) to be a "critical point", we set all the these first derivatives to 0 to get, 77 f1X1 S + f2X2 S + f3X3 S + ... + f S-1XS-1 S + fS = 0 f1X1 S+1 + f2X2 S+1 + f3X3 S+1 + ... + f S-1XS-1 S+1 + fS+1 = 0 f1X1 S+2 + f2X2 S+2 + f3X3 S+2 + ... + f S-1XS-1 S+2 + fS+2 = 0 ... f1X1 N + f2X2 N + f3X3 N + ... + f S-1XS-1 N + fN = 0 . // N-S+1 equations (6.2.7)' A similar set of equations is obtained for a,b,c...q. Just replace f →a, then f →b and so on. We shall not write out all these equations as we did in (6.2.9). We have then S sets of equations, each set containing N-S+1 equations. As in the sample case, we then re order the equations to get first this group, f 1X1 S + f2X2 S + f3X3 S + ... + f S-1XS-1 S + fS = 0 a1X1 S + a2X2 S + a3X3 S + ... + a S-1XS-1 S + aS = 0 b1X1 S + b2X2 S + b3X3 S + ... + b S-1XS-1 S + bS = 0 .... q 1X1 S + q2X2 S + q3X3 S + ... + q S-1XS-1 S + qS = 0 . // S equations (6.2.10)' S The next group has S →S+1 in all subscript positions. The last group has S →N, and here is that last group, f 1X1 N + f2X2 N + f3X3 N + ... + f S-1XS-1 N + fN = 0 a1X1 N + a2X2 N + a3X3 N + ... + a S-1XS-1 N + aN = 0 b1X1 N + b2X2 N + b3X3 N + ... + b S-1XS-1 N + bN = 0 .... q 1X1 N + q2X2 N + q3X3 N + ... + q S-1XS-1 N + qN = 0 . // S equations (6.2.10)' N So there are now N-S+1 groups of equations and each group has S equations each having S terms. The next task is to write each set of S equations as a matrix equation. Here we do that just using the abbreviated notation for the matrix. Ther e are then N-S+1 matrix equations, ⎝⎜⎜⎛ ⎠⎟⎟⎞ 0 0 ... 0 0 = (f 1,f2, f3, ... fS-1, fS) ⎝⎜⎜⎛ ⎠⎟⎟⎞ X1 S X2 S ... XS-1 S 1 ⎝⎜⎜⎛ ⎠⎟⎟⎞ 0 0 ... 0 0 = (f 1,f2, f3, ... fS-1, fS+1) ⎝⎜⎜⎛ ⎠⎟⎟⎞ X1 S+1 X2 S+1 ... XS-1 S+1 1 •••• ⎝⎜⎜⎛ ⎠⎟⎟⎞ 0 0 ... 0 0 = (f 1,f2, f3, ... fS-1, fN) ⎝⎜⎜⎛ ⎠⎟⎟⎞ X1 N X2 N ... XS-1 N 1 . (6.2.11)' 78 We then argue as before that, to avoid the contradic tion 0 = 1, the determinants of these SxS submatrices must be zero, det(f 1,f2, f3, ... fS-1, fS) = 0 det(f 1,f2, f3, ... fS-1, fS+1) = 0 ... det(f 1,f2, f3, ... fS-1, fN) = 0 or det(f 1,f2, f3, ... fS-1, fm) = 0 m ≠ 1,2,3...S-1 . (6.2.13)' Finally, if we start off by eliminating the S-1 variables x i, xj, xk.... xr where i ≠ j ≠k....≠ r, we would find that det(f i,fj, fk, ... fr, fm) = 0 m ≠ i,j,k...r . (6.2.15)' Since this includes any possible SxS subdeterminant of the R matrix, we have then shown that all (N,S) SxS subdeterminants vanish, and therefore rank[R( r)] < S when r is a stationary point of f(r ) subject to the constraints. 79 Appendix A: Matrix rank equals the number of independent columns and rows This Append ix proves our claim of (6.1.7) that, for a general m x n matrix M, the number of linearly independent columns and the number of linearly inde pendent rows are the same, and both numbers are equal to the rank of the matrix. The reader will appr eciate that this fact is non-obvious. Shilov proves it but the proof is a bit spread out over several sections. Here we provide a self-contained alternate proof. We first quote a series of square-matrix theorems th at are proved from scratch in our Appendix B. These are doubtless very familiar to the reader, but we just want to get them stated. The derivations of these theorems may be less familiar. Theorem 1 : det(M T) = det(M). Switching rows with column s does not change a determinant. (B.1.10) Theorem 2: det(M) can be represented in these two ways, where ε is the permutation tensor (B.2.8) : det(M) = Σ a1a2... an εa1a2... an M1a1M2a2 ...... Mnan (B.2.10a) det(M) = Σ a1a2... an εa1a2... an Ma11Ma22 ...... Mann . (B.2.10b) Theorem 3: For a square matrix M, adding a multiple of one row to another does not change det(M). The same is true for adding a multiple of one column to another. (B.2.12) Theorem 4: Swapping two rows (columns) of square matrix M causes det(M) → - det(M). (B.2.13) Corollary 4: If two rows (columns) of square matrix M are the same, then det(M) = 0. (B.2.14) Theorem 5 : det(AB) = det(A)det(B) (B.2.15) Corollary 5 : det(ABC) = det(A)det(B)det(C) and so on. (B.2.16) Theorem 6 (Cofactor Expansions): det(M) = Σ n Msn cof(Msn) = Σn (rs)n cof(Msn) work across row s s = 1,2....n (B.3.16a) det(M) = Σn Mns cof(Mns) = Σn (cs)n cof(Mns) work down column s s = 1,2....n (B.3.16b) Theorem 7: M-1 = CT det(M) = [cof(M)]T det(M) = cof(MT) det(M) for square matrix M. (B.4.8) Theorem 8 (Cramer's Rule): y = M x ⇒ xs = det(M[ cs→ y]) det(M) (B.4.14) To this list, we shall now add four new Theorems which will be proven right here. The definition of linear independence and linear dependence is given in the discussion surrounding (6.1.6) and won't be repeated here. 80 Theorem 9: Columns of square matrix M are linearly dependent ⇔ det(M) = 0 (A.1) Contrapositives : Columns of square matrix M are linearly independent ⇔ det(M) ≠ 0 If we can prove Theorem 9 as stated, then the theorem is also true fo r rows. The reason is that swapping rows and columns corresponds to M ↔ MT and Theorem 1 says det(MT) = det(M). Proof of ⇒: According to (6.1.6) and nearby discussion, our premise of linear dependence says that either one or more columns of M are zero, or at least one column can be written as a linear combination of the other columns, so perhaps cj = Σi≠jkici . If a column is zero, det(M) = 0 from Theorem 6 going down that column. If c j = Σi≠jkici, one can add - Σi≠jkici to cj without changing det(M) according to Theorem 3 . But this makes the new column j vanish, so again det(M) = 0. Proof of ⇐ : Here we shall prove the contrapositive instead of the claim: claim: det(M) = 0 ⇒ the columns of M are linearly dependent contrapositive: the columns of M are linearly independent ⇒ det(M) ≠ 0 If the column vectors c i of matrix M are linearly independent, then (6.1.6) says Σxici = 0 ⇒ x = 0. Therefore we know that ΣxiMji = 0 ⇒ xj = 0 which is the same as Mx = 0 ⇒ x = 0. Thinking of M:En → En, we claim that mapping M x = y is one-to-one. Certainly each x goes into a single y, Could a y map back into two different values of x, call them x and x' with x ≠ x' ? Then M x = y and M x' = y . Subtract to get M( x-x') = 0. But since M z = 0 ⇒ z = 0 (as just shown), we find x = x' . Thus, the mapping Mx = y really is one-to-one, and that means it is invertible and the inverse is unique. From Theorem 7 the inverse is in fact given by M-1 = cof(MT)/det(M). For M-1 to exist, one must have det(M) ≠ 0. Another proof is very simple but perhaps less c onvincing. One can show (see our Tensor Analysis document in Refs) that if a set of column vectors ci of matrix M spans an n-piped in En, the "volume" of that n-piped is given by V = |det(M)|. If those vectors are linearly independent, then V ≠ 0. This is obvious in 2D (volume = area) and 3D but less obvious for general En. For example, in 3D if a third vector lies in the plane of the other two (is dependent), there is no volume. Theorem 10 (mini-columns theorem). Let R be a subset of N rows of nxm matrix A, and let S be a subset of columns {c i} in A. The matrix elements of A included in the intersection of these two sets form a set of "mini-columns" { Ci}. This set of mini-columns forms a matrix B within A. Matrix B could be contiguous, or it could be non-contig uous in one or both directions. Th e claim and its contrapositive are: (a) { Ci} linearly independent in B ⇒ {ci} linearly independent in A (b) {c i} linearly dependent in A ⇒ {Ci} linearly dependent in B (A.2) Below we shall prove (b) which then implies (a). A picture is worth a thousand words. Here matrix B happens to be contiguous in both directions : 81 (A.3) Proof of (b): If the set { ci} is linearly dependent in A, then from (6.1.6) one can write Σi∈S kici = 0 with k ≠ 0. This is a set of m scalar equations (one for each row in A) which contains as a subset the set of N equations Σi∈S kiCi = 0 (one for each row in B). This last equation then says that the { Ci} are linearly dependent in B. There is of course a similar "mini rows theorem". Theorem 11 . If all kxk minors in k columns of matrix A vanish, the k columns are linearly dependent. ( A . 4 ) Contrapositive : If k columns are linearly independent, they must contain at least one non-zero kxk minor. Once proven for columns, replacing A →AT gives the same theorem for rows. Proof : Assume matrix A has m rows and n columns. In order for kxk minors to exist, one must have k ≤ min(m,n). Gather up the k columns and make them be the leftmost k columns of a new matrix B. Since these columns have m elements, and since k ≤ m, add m-k arbitrary new columns to the right of the k columns so that matrix B is then a square matrix. On the left below we show the k columns taken from matrix A in red, and then the arbitrary added columns are shown in blue. (A.5) Consider now the process of computing det(B) by going down the rightmost column using the standard cofactor sum formula of Theorem 6 . This det(B) is a linear combinati on of (m-1)x(m-1) minors all in the left m-1 columns. Consider one of these minors. If we evaluate it using the same cofactor formula (of one less dimension), it will be a linear combination of (m-2)x(m-2) minors in the left m-2 columns. We keep going until we are evaluating a set of kxk cofactors in the left k columns. But these all vanish by the theorem premise. Thus, reversing this logic we conclude that det(B) = 0. The picture on the right above shows one minor at each level of the descent just described (red dot goes with red minor, etc). Now suppose it were possible that the k red columns were independent. Since the added blue columns are arbitrary, we could certainly find a set of added blue columns so that all m columns were independent. 82 But then we would have det(B) ≠ 0 by Theorem 9 . But we just showed that det(B) = 0, so it must not be possible to have the k red columns be i ndependent, so they must be dependent. We are finally ready for the Main Act. Recall first that : Definition : The rank r of an m x n matrix A is the dimension of the largest non-vanishing minor within A. [Shilov 1.92 ] Thus, r ≤ min(m,n). (6.1.5) (A.6) Theorem 12 . For a general nxm matrix A, (A.7) rank(A) = r ⇔ A has exactly r linearly independent columns Once this theorem is proved for columns, it is also true for rows since rank(A T) = rank(A) by (A.6). Proof of ⇒ : If rank(A) = r, A must have at least one non-vanishing rxr minor. Think of this minor being a set of mini-columns as in Theorem 10. Since then det(minor) ≠ 0, this set of mini-columns is linearly independent according to Theorem 9 . Thus, the set of corresponding full columns (those that pass down through the minor) are linearly independent from Theorem 10 (a). So matrix A has at least these r independent columns. Now let k = r+1. Since rank(A) = r we know that all kxk minors in A vanish. By Theorem 11, any set of k columns must be linearly dependent. Since k = r+1, any r+1 columns are linearly dependent, so the number of independent columns is just r, the number of columns passing down through the minor. Proof of ⇐ : If there are r linearly independent columns, th en at least one rxr minor within those columns must be non-zero (contrapositive of Theorem 11). Suppose there were an (r+k)x(r+k) non-vanishing minor. Then the r+k mini-columns of that minor would be independent, and thus by Theorem 10 there would be r+k independent full columns. But the th eorem premise says there are only r independent columns, so all minors larger than rxr must vanish. Therefore rank(A) = r. Sometimes the number of linearly independent rows of a matrix is called the row rank , and the number of independent columns is called the column rank . Theorem 12 and its row version show that column rank = row rank = rank = dimension of largest non-vanishing minor (A.8) Basis Minors and Basis Columns. If rank(A) = r, we know there must exist at least one non-vanishing rxr minor in A. Shilov refers to such a minor as a basis minor and the columns passing through this minor are called basis columns. We showed in Theorem 10 that the set of such basis columns is linearly independent (which is why they are called basis columns). Obviously any of these columns can be written as a linear combination of th e basis columns, such as c2 = Σi=1k kici = 1 c2. If k = r+1, all kxk minors vanish and by Theorem 11 all sets of k = r+1 columns are linearly dependent. Thus, every non-basis column is a linear combination of the basis columns. We have thus proved Shilov's Basis Minor Theorem which we quote from section 1.93 of his book (p 25) (A.9) In his proof that column rank = rank, Shilov uses the above theorem as a starting point. 83 Appendix B: Determinants An alterative title for this Appendix m ight be "More th an you really want to k now about determinants". We find it convenient to gather these Facts all in on e place. Many of these Facts are used in Appendix A where we prove other Facts about matrix rank which in turn are used in Section 6 to support our alternate derivation of the Method of Lagrange Multipliers. It is hoped that the encapsulated presentation below might be useful to the reader for other applications of determinants. Of special interest to us is that permutations form a group which im plies certain "rearrangement theorems" which in turn provide quick proofs of many well-known and some lesser-known Facts about determinants. B.1 Definition of the determinant First, we define z 0 to be the n-component vector of increasing integers 1 to n, z0 ≡ ⎝⎜⎛ ⎠⎟⎞ 1 2 ... n . ( B . 1 . 1 ) Let a be some permutation (reordering) of these integers, so write a ≡ ⎝⎜⎛ ⎠⎟⎞ a1 a2 ... an = A ⎝⎜⎛ ⎠⎟⎞ 1 2 ... n = A z 0 ( B . 1 . 2 ) where A is an nxn matrix which has 1's in th e right places to create this permutation vector a. For example, a = Az o ↔ ⎝⎜⎛ ⎠⎟⎞ 1 3 2 = ⎝⎜⎛ ⎠⎟⎞ 1 0 0 0 0 1 0 1 0 ⎝⎜⎛ ⎠⎟⎞ 1 2 3 . The vector a can be obtained from the vector z 0 by making S a pairwise swaps of elements of z0. Although S a is not unique, the number (-1)Sa is unique. For example, to get from (1,2,3) to (1,3,2) one could swap the second pair so S a = 1, but one could then swap the first pair twice and then S a = 3. This number (-1)Sa is of course ±1 and we shall call it the parity of the permutation a, Parity( a) ≡ (-1)Sa . ( B . 1 . 3 ) Fact: Parity(C a) = Parity(CA z 0) = (-1)Sc+Sa = Parity(c ) Parity( a) = Parity(C-1a) (B.1.4) Proof : Doing permutation CA involves first doing A with its S a swaps, and then doing C with its S c swaps, for a total of S a+Sc swaps. The inverse permutation C-1 obviously involves the same number of swaps as the permutation C. 84 Moving now toward the definition of the determinant of a matrix M, define ( Π means "product") Π(a,b; M) ≡ Mab ≡ Ma1b1Ma2b2 ...... Manbn = product of n factors . (B.1.5) The notation Π(a,b; M) is easier to deal with than M ab, but we shall use both these notations. Suppose c = Cz 0 is some arbitrary permutation of z0. Then: Fact: Π(Ca,Cb; M) = Π(a,b; M ) ( B . 1 . 6 ) Proof : Applying the same permutation to both the a k and bk indices of M a1b1Ma2b2 ...... Manbn merely reorders the terms in the product but the product stays the same. Fact : Π(a,b; MT) = Π(b,a; M ) ( B . 1 . 7 ) Proof : Π(a,b; MT) = MT a1b1MT a2b2 .... MT anbn = Mb1a1Mb2a2 .... Mbnan = Π(b,a; M) . We shall now define the determinant of an n x n matrix M in the following admittedly obscure manner (later we will show that it reduces to more familiar forms), det(M) ≡ 1 n! Σa Σb (-1)Sa+Sb Π(a,b; M) = 1 n! Σa,b (-1)Sa+Sb Mab = 1 n! Σa Σb (-1)Sa+Sb Ma1b1Ma2b2 ...... Manbn . (B.1.8) Here Σa means the sum over all permutations a of z0. In (B.1.8) the columns and rows of M are on a completely equal footing. Notice that a and b are in effect dummy summation indices. If we do a↔b, the expression for det(M) is unchanged. Thus one can rewrite (B.1.8) as det(M) = 1 n! Σb Σa (-1)Sb+Sa Π(b,a; M) . (B.1.9) We are now ready for our first determinant theorem: Theorem 1 : det(M T) = det(M). Switching rows with columns does not change a determinant. (B.1.10) Proof : det(MT) = 1 n! Σb Σa (-1)Sb+Sa Π(b,a; MT) // (B.1.9) applied to MT = 1 n! Σb Σa (-1)Sb+Sa Π(a,b; M) // (B.1.7) applied to MT = det(M) // (B.1.8) 85 B.2 The permutation group, the permutation tensor ε, and expansions for det(M) Definition. A set of ele ments {g i} form a group G if : • gigj is also in the group (closure) • gi-1 exists for each g i in G, where g i-1 is also in G (inverse) • (gigj)gk = gi(gjgk) (associative) (B.2.1) Comment : It is implicit in the above definition that th e group has some "operation" which gives meaning to gigj, which we shall just think of as "multiplica tion". When group elements are represented by matrices, that operation is multiplication of those ma trices, and that will apply to the permutation group below. Fact : giG = G (the rearrangement theorem) (B.2.2) This says that multiplication of all the elements of a group by an element g i in the group creates a reordering of the group elements. Here G is the set of group elements. Proof . Consider g iG = gi [g1, g2....gn] = [gig1, gig2....gign] = set of n elements. Unless two elements are the same, this must exhaust the entire group. How do we know that g ig1 and gig2 might not be the same? Apply g i-1 from the left and that would say g 1 = g2 which is not the case. Fact: Σ g f(g) = Σ g f(g1g) if Σg runs over the entire group G (B.2.3) Proof: In Σg f(g1g), as g runs over G, the argument g' ≡ g1g runs over G by the above rearrangement theorem. Thus, the sum Σg f(g1g) is just a reordering of the terms in the sum Σg f(g). It is easy to show that the set of permutations of z 0 forms a group and therefore the above facts can be used. For example, the product of two permutations is a permutation, and ever y permutation clearly has an inverse, and (AB)C = A(BC) for any matrices. Fact (B.2.3) can be written ΣB f(Bzo) = ΣB f(CB z0) // b = Bz0 Here the sum Σ B is over all permutation matrices which correspond to permutations b of z0, and C is some arbitrary permutation matrix associated with some permutation c of z0. An equivalent way of stating the above involves a direct summation Σb over all permutation vectors b of z0 , Σ b f(b) = Σb f(Cb) . // permutation sum rearrangement theorem (B.2.4) Throwing in the arbitrary permutation C merely causes a reordering of the sum. This is an extremely powerful and useful fact. Consider then : 86 det(M) ≡ 1 n! Σa [Σb (-1)Sa+Sb Π(a,b; M) ] // (B.1.8) = 1 n! Σa [Σb (-1)Sa+Sb Π(A-1a,A-1b; M)] // (B.1.6) with C = A-1; A-1a = A-1Az0 = z0 = 1 n! Σa [ Σb Parity(A-1b) Π(z0,A-1b; M) ] // (B.1.4) with C = A-1 and a = b = 1 n! Σa [ Σb Parity( b) Π(z0,b; M) ] // (B.2.4) with C = A-1 (key step) = Σb Parity( b) Π(z0,b; M) // n! identical terms in Σa = Σa Parity( a) Π(z0,a; M) // rename dummy sum variable = Σa Parity( a) M1a1M2a2 ...... Mnan . // (B.1.5) (B.2.5) Since det(MT) = det(M) from Theorem 1 (B.1.10), we can also write this result as det(M) = Σa Parity( a) Ma11Ma22 ...... Mann ( B . 2 . 6 ) Proof : det(M) = det(MT) = Σa Parity( a) MT 1a1MT 2a2 ...... MT nan = Σa Parity( a) Ma11Ma22 ...... Mann Here then are the last two results: ( recall that S a is the number of pairwise swaps to get from z0 to a ) det(M) ≡ Σa Parity( a) M1a1M2a2 ...... Mnan Parity( a) = (-1)Sa (B.2.7a) det(M) ≡ Σa Parity( a) Ma11Ma22 ...... Mann Parity( a) = (-1)Sa . (B.2.7b) Definition : The permutation tensor εijk.. (n subscripts) : • ε12...n = 1 • for any index swap, ε changes sign: ε..i..j.. = - ε..j..i.. • therefore, if any two indices are the same, ε = 0 ε..i..i.. = - ε..i..i.. = 0 (B.2.8) Fact : If a is a permutation of z 0, then εa1a2... an = Parity( a) = (-1)Sa ≡ εa . (B.2.9) Proof : Since indices a i represent a permutation of z0 , it takes S a swaps to get from εa1a2... an to ε123...n by the definition of ε. In dense notation, one could say εa = (-1)Sa εz0 . Using (B.2.9) in (B.2.7) then gives these two classic det(M) expressions : Theorem 2: det(M) can be represented in these two ways: det(M) = Σ a1a2... an εa1a2... an M1a1M2a2 ...... Mnan (B.2.10a) det(M) = Σ a1a2... an εa1a2... an Ma11Ma22 ...... Mann . (B.2.10b) 87 In dense notation, one could write the above as ( a is now a vector), det(M) = Σ aεa Mz0a ( B . 2 . 1 0 a ) ' det(M) = Σ aεaMaz0 . (B.2.10b)' In most applications, the following simpler notation suffices: det(M) = Σ abc..q εabc..q M1aM2bM3c....Mnq ( B . 2 . 1 1 a ) det(M) = Σ abc..q εabc..q Ma1Mb2Mc3....Mqn . (B.2.11b) Theorem 3: For a square matrix M, adding a multiple of one row to another does not change det(M). The same is true for adding a multiple of one column to another. (B.2.12) Proof : (rows) Suppose we replace r3 → r3 + α r2. This says M 3i → M3i + α M2i . Eq (B.2.11a) says : det(M') = Σ abc..q εabc..q M1aM2b(M3c + α M2c) ....Mnq = det(M) + α Σ abc..q εabc..q M1aM2bM2c ....Mnq . Since M 1aM2bM2c ....Mnq is symmetric under b ↔c while εabc..q is anti-symmetric, the extra α term vanishes. That is to say, sum = Σbc AbcSbc = Σcb AcbScb = Σcb (-Abc)(Sbc) = - sum = 0. Using (B.2.11b) this argument shows that c3 → c3 + α c2 similarly does not alter det(M). Theorem 4: Swapping two rows (columns) of square matrix M causes det(M) → - det(M). (B.2.13) Proof : Let's swap rows 1 and 3 in M to get M'. Then from (B.2.11a), det(M') = Σ abc..q εabc..q M3aM2bM1c....Mnq = Σ abc..q [- εcba..q] M3aM2bM1c....Mnq // swap a↔c on ε = - Σ abc..q εabc..q M3cM2bM1a....Mnq // dummy rename a ↔c = - Σ abc..q εabc..q M1aM2bM3c....Mnq // reorder product = - det(M) . Using Theorem 1 that det(A) = det(A T), we then know that det(M'T) = - det(MT), so swapping columns 1 and 3 negates det(M). An obvious corollary follows if the two swapped columns have identical data : Corollary 4: If two rows (columns) of square matrix A are the same, then det(M) = 0. (B.2.14) 88 Here is one more well-known de terminant theorem which again demonstrates the power of the rearrangement theorem (B.2.4). Theorem 5 : det(AB) = det(A)det(B) (B.2.15) Proof : Since we already use A in a = Az0 , we shall prove det(XY) = det(X)det(Y) to avoid overloading symbols. Note that Parity( a) Parity( b) = (-1)Sa(-1)Sb = (-1)Sa+Sb = Parity(A b) from (B.1.4). Then : det(X)det(Y) = [ Σ a Parity( a) Π(z0,a; X)] [ Σb Parity(b) Π(z0,b; Y)] // (B.2.5) twice = Σa [ Σb Parity( a) Parity( b) Π(z0,a; X) Π(z0,b; Y) ] // (B.1.4) used next line = Σ a [ Σb Parity(A b) Π(z0,a; X) Π(Az0,Ab; Y) ] // (B.1.6) with C = A = Σ a [ Σb Parity( Ab) Π(z0,a; X) Π(a,Ab; Y) ] // A z0 = a = Σ a [ Σb Parity( b) Π(z0,a; X) Π(a,b; Y) ] // (B.2.4) with C = A = Σ b Parity( b) Σa Π(z0,a; X) Π(a,b; Y) // move Σa = Σ b Parity( b) Σa Π(z0,b; XY) // see below = d e t ( X Y ) . / / ( B . 2 . 5 ) The idea here is to get a in the right place on both Π's so that Σ a Π(z0,a; X) Π(a,b; Y) = Σ a1a2... an X1a1X2a2 ...... Xnan Ya1b1Ya2b2 ...... Yanbn = Σ a1a2... an (X1a1Ya1b1)(X2a21Ya2b2).....(XnanYanbn) = (XY) 1b1 (XY)2b2 ... (XY) nbn = Π(z0,b; XY) which in dense notation one would write as Σ a Xz0aYab = (XY) z0b . Corollary 5 : det(ABC) = det(A)det(B)det(C) and so on. (B.2.16) Proof : det(ABC) = det(A[BC]) = det( A)det(BC) = det(A)det(B)det(C) . 89 It is not hard to slightly generalize the results of Theorem 2 to obtain : Theorem 2A: If b is a permutation of z 0, then det(M) can be represented in these two ways: det(M) = εb1b2... bn Σa1a2... an εa1a2... an Mb1a1Mb2a2 ...... Mbnan (B.2.17a) det(M) = εb1b2... bn Σa1a2... an εa1a2... an Ma1b1Ma2b2 ...... Manbn . (B.2.17b) In dense notation, one could write the above as ( a and b are now vectors), det(M) = ε bΣaεaMba ( B . 2 . 1 8 a ) det(M) = ε bΣaεaMab . (B.2.18b) Since ε b2 = 1, one can move εb to the left side of these equations, so in the dense notation, ε bdet(M) = Σ aεaMba ( B . 2 . 1 9 a ) εbdet(M) = Σ aεaMab . (B.2.19b) Equations (B.2.19) are valid even if b is not a permutation of z 0. In this case one has 0 = 0. For the special case b = z0 Theorem 2A reduces to Theorem 2 since εz0 = ε12...n = 1. Proof of Theorem 2A: We shall prove (B.2.17a) and then (B.2.17b) follows by taking M →MT and using det(MT) = det(M). det(M) = Σ a Parity( a) Π(z0,a; M) // line 6 of (B.2.5) = Σ a Parity( a) Π(Bz0,Ba; M) // (B.1.6) with C = B = Σ a Parity( Ba) Parity( b) Π(Bz0,Ba; M) // Parity(B a) Parity( b) = Parity( a) = Σ a Parity( a) Parity(b) Π(Bz0,a; M) // (B.2.4) that Σa f(Ba) = Σa f(a) = Parity( b) Σa Parity( a) Π(b,a; M) // b = Bz0 = εb1b2... bn Σa1a2... an εa1a2... an Mb1a1Mb2a2 .... Mbnan // (B.1.5) and (B.2.9) = εb Σa εa Mba . // previous line in dense notation (B.2.20) Here it is assumed that B and hence b is associated with a permutation of z 0, In this case, ( εb)2 = 1 so one can move εb to the left side to get ε bdet(M) = Σ a εa Mba . ( B . 2 . 2 1 ) 90 In this form, the equation is valid whether or not b is a permutation of z0 . If in b = ( b1,b2...bn) the bi are arbitrary elements of the set {1,2,...n} where two or more b i are the same (that is, b is not a permutation of z0), then Σa εa1a2... an Mb1a1Mb2a2 .... Mbnan = 0 by symmetry. For example, if b 1 = b2 then S ≡ Σa εa1a2... an Mb1a1Mb1a2 .... Mbnan = Σa εa2a1... an Mb1a2Mb1a1 .... Mbnan // rename dummy indices a 1 ↔ a2 = Σa [-εa1a2... an ] Mb1a1Mb1a2 .... Mbna // (B.2.8) for ε and slide M b1a2 to the right = - Σaεa1a2... anMb1a1Mb1a2 .... Mbna = - S ⇒ S = 0 . (B.2.22) So in this case the right side of (B.2.21) is zero. But from (B.2.8) εb = 0 so the left side is also zero. Reader Exercise : Prove that ε1i2... indet(M) = Σ A parity(A) M 1,A(i1)M2,A(i2)...... Mn,A(in) . (B.2.23) Proof in dense notation: εidet(M) = εiΣa εa Mz0,a // this is εi times usual expansion (B.2.7a) = parity(I) Σ A parity(A) M z0,Az0 // a = Az0 = parity(I) Σ A parity(AI) M z0,AIz0 // rearrangement theorem ΣAf(A) = ΣA f(AI) = ΣA parity(A) M z0,Ai // i = Iz0 QED 91 B.3 Minors, Cofactors, and the Cofactor Expansions of det(M) Definition : The minor of a matrix element M rs is the determinant of the submatrix of M obtained by crossing out the rth row and the sth column. It is a challenge, however, to write this out in symbols. Let's start with a more detailed version of (B.2.11a), where the a i are column summation indices, det(M) = Σ a1a2a3a4a5...an εa1a2a3a4a5...an M1a1M2a2M3a3M4a4M5a5 ...... Mnan . (B.3.1) We shall make the following conjecture for the form of minor(M 23) minor(M 32) = (-1)3-2 Σa1a2a4a5...an εa1a22a4a5...an M1a1M2a2M4a4M5a5 ...... Mnan (B.3.2) Compared to det(M) shown in (B.3.1), we have made these changes : • Removed the factor M 3a3 ( since row-3 matrix elements cannot appear in minor(M 32) ) • Removed the sum over a 3. • Replaced a 3 by the number 2 on the ε tensor. • added a sign factor (-1)3-2. Notice that since there is a 2 on the ε tensor, any time a summation index a i = 2 there is no contribution since then the ε tensor has two indices the same, so in effect the value 2 has been removed from all the residual summations a i. That is good, since column 2 is supposedly "crossed out" in minor(M 32). The factor (-1)3-2 is added so that the "diagonal term" in the minor will be positive. This sign (-1)3-2 gets used up if we slide the "2" to its natural position (position 2) on the ε tensor, minor(M 32) = Σa1a2a4a5...an εa12a2a4a5...an M1a1M2a2M4a4M5a5 ...... Mnan . (B.3.3) The diagonal term in minor(M 32) is then positive, matching the mechanical crossing-out method, ε 12345..n M11M22M44M55 ...... Mnn = + M11M22M44M55 ...... Mnn . (B.3.4) A more compact notation for (B.3.2) would be minor(M 32) = (-1)3-2 Σai,i≠3 εai,a3=2 Πi≠3 (Mi,ai) . (B.3.5) Starting over, we could show similarly that minor(M 42) = (-1)4-2 Σai,i≠4 εai,a4=2 Πi≠4 (Mi,ai) . (B.3.6) We now have a sign factor (-1)4-2 because the "2" on ε has to be slid 2 positions to get to its natural location (position 4 on ε). More generally we can say minor(M s2) = (-1)s-2 Σai,i≠s εai,as=2 Πi≠s (Mi,ai) (B.3.7) 92 where the slide is now s-2 places. Still more generally we find, replacing 2 by r, minor(M sr) = (-1)s±r Σai,i≠s εai,as=r Πi≠s (Mi,ai) . (B.3.8) This then is our "best form" expression for a minor of matrix element M sr. Either sign will do, since (-1)s-r = (-1)s-r(-1)2r = (-1)s+r since (-1)2r = 1 . Choosing the + sign in (B.3.8) and putting (-1)s+r on the left side we get, (-1)s+r minor(M sr) = Σai,i≠s εai,as=r Πi≠s (Mi,ai) . (B.3.9) The left side here is called the cofactor of Msr, written cof(M sr). Thus we have shown that cof(M sr) ≡ (-1)s+r minor(M sr) = Σai,i≠s εai,as=r Πi≠s (Mi,ai) . (B.3.10) Fact : Neither minor(M sr) nor cof(M sr) are functions of the M sr matrix element of M! In (B.3.10) this is so because: (1) row s is excluded in Π i≠s(Mi,ai); (2) ai = r is excluded by the factor εai,as=r . More intuitively, this is so because to get minor(M sr) we "cross out" row s and column r. Thus, ∂cof(Msr) ∂Msr = ∂minor(M sr) ∂Msr = 0 for any pair r,s in 1,2...n . (B.3.11) Finally, suppose we replace r with an integer which we call a s. Doing this causes εai,as=r → εai,as=as = εai = εa1a2a4a5...an = the normal ε tensor form . (B.3.12) Then (B.3.10) becomes the following, cof(M sas) ≡ (-1)s+as minor(M sas) = Σai,i≠s εai Πi≠s (Mi,ai) . (B.3.13) Now start again with det(M) of (B.3.1) and rewrite it in our compact notation: det(M) = Σ a1a2a3a4a5...an εa1a2a3a4a5...an M1a1M2a2M3a3M4a4M5a5 ...... Mnan = Σ ai εai Πi(Mi,ai) // next, extract the a s sum and its M sas factor : = Σas Msas [Σai,i≠s εai Πi≠s (Mi,ai)] // next use (B.3.12) to get = Σas Msas [cof(Msas)] // next, change summation index from a s to n = Σn Msn cof(Msn) . // valid for any s in 1,2...n (B.3.14) 93 This is the cofactor expansion of det(M) where one "works across row s" and n is a column index. Another form of this expansion is det(M) = det(M T) = Σn MT sn cof(MT sn) = Σ n Mns cof(Mns) ( B . 3 . 1 5 ) and in this form one "works down column s" wher e n is now a row index. We have just proven : Theorem 6 (Cofactor Expansions): det(M) = Σ n Msn cof(Msn) = Σn (rs)n cof(Msn) work across row s s = 1,2....n (B.3.16a) det(M) = Σn Mns cof(Mns) = Σn (cs)n cof(Mns) work down column s s = 1,2....n (B.3.16b) The term "cofactor" presumably arose because each sum term consists of a "factor" M ij and a "co-factor" which cooperates with the factor to make the su m, as coauthors cooperate to write a book. A common notation is to define Mij ≡ cof(Mij) and then the cofactor expansions become det(M) = Σ n MsnMsn det(M) = Σn MnsMns . Although compact, this is not so handy for tensor notation with its up and down indices (see Lucht Ref). B.4 Expressions for the inverse matrix M-1 and Cramer's Rule Consider the expansion (B.3.16a) working across row s, det(M) = Σ n Msn cof(Msn) . ( B . 4 . 1 ) Suppose we replace the elements of row s (M sn) with the elements of some other row r in M (M rn). In doing so we have created a new matrix, call it M'. Notice that cof(M' sn) = cof(M sn) because , although going from M to M' we have altere d row s, we have not altered cof(M sn) since this depends only on the entries in all the rows other than row s (think "crossing out" row s for minor(M sn) ). See Fact (B.3.11) and text above. Therefore we find, det(M') = Σn M'sn cof(M'sn) = Σn Mrn cof(Msn) . (B.4.2) But by Corollary 4 (B.2.14) det(M') = 0 si nce two rows of M' are the same. Thus 0 = Σn Mrn cof(Msn) . r ≠ s ( B . 4 . 3 ) Combining this with (B.4.1) gives, Σ n Mrn cof(Msn) = det(M) δr,s . (B.4.4) 94 Now for clarity define a cofactor matrix C in this manner C sn ≡ cof(Msn) . ( B . 4 . 5 ) In matrix notation, we could define the matrix cof(M) to be matrix C and then C = cof(M) ⇒ Csn = [cof(M)] sn = cof(M sn) . (B.4.6) Now (B.4.4) says Σ nMrn Csn = det(M)δ r,s or Σ nMrn CT ns = det(M) δr,s and finally in matrix notation, MCT = det(M) 1 or M [ CT det(M) ] = 1 . (B.4.7) Therefore we find this classic square matrix inversion formula, Theorem 7: M -1 = CT det(M) = [cof(M)]T det(M) = cof(MT) det(M) . (B.4.8) To verify the last equality in (B.4.8) consider : = [cof(M)] T ns = [cof(M)] sn // meaning of transpose = cof(M sn) // (B.4.6) = cof(M T ns) = [cof(MT)]ns // (B.4.6) applied to MT and therefore we have this matrix identity, [cof(M)] T = [cof(MT) ] . ( B . 4 . 9 ) Using (B.4.8) it is trivial to solve a non-singular (detM ≠0) system of linear equations y = M x : y = M x ⇒ x = M-1 y = 1 det(M) cof(MT) y = 1 det(M) [cof(M)]T y . (B.4.10) 95 In components, x s = 1 det(M) Σn ([cof(M)]T)sn yn = 1 det(M) Σn yn [cof(M)] ns = 1 det(M) Σn yn cof(Mns) . (B.4.11) Now recall the cofactor expansion (B.3.16b), det(M) = Σ n Mns cof(Mns) = Σn (cs)n cof(Mns) . (B.4.12) If one replaces column c s in M by y one gets det(M[ c s→ y]) = Σn yn cof(Mns) // (B.4.12) = det(M) x s . // (B.4.11) (B.4.13) Solving this for x s we obtain another classic result, Theorem 8 (Cramer's Rule, 1750): // Gabriel Cramer (1704-1752), Swiss, made it to age 47 y = M x ⇒ x s = det(M[ cs→ y]) det(M) . ( B . 4 . 1 4 ) Comment: Names of Matrices Historically the matrix CT = [cof(M)]T = [cof(MT)] was called the adjoint of matrix M and was often denoted by M^ and then (B.4.8) reads M-1 = M^ /|M|. But this then conflicted with another usage, namely the matrix M† ≡ (MT)* = (M*)T being the adjoint of M (the conjugate transpose of M). When M† = M, M is said to be self-adjoint and has significance in quantum (matrix) mechanics: the quantum operators of all physical observables are self-adjoint and thus have real eigenvalues. Nowadays the matrix [cof(M)] T is called the classical adjoint of M, while M† is then the Hermitian adjoint (or Hermitian conjugate ) of M and is sometimes written MH. If M† = M then M is said to be Hermitian. Historically when [cof(M)]T was the adjoint of M, M† was called the associate of M. The classical adjoint sometimes appears as the adjugate or adj(M). Transpose matrices MT are often denoted M~ . On a related matter, some older authors used the word minor to refer to what we now call a cofactor. See for example Morse & Feshbach page 509. These authors never use the word cofactor. 96 Appendix C: The Intersection of Constraint Surfaces In the proof of Theorem 3 in Section 6.2 (a) we rely on our ability, at least in theory, to eliminate x 1, x2, x3 from the three constraint equations shown in (6.2.3). a(x1, x2, x3, x4, x5, x6) = 0 x 1 = X1(x4, x5, x6) b(x1, x2, x3, x4, x5, x6) = 0 ⇒ x2 = X2(x4, x5, x6) c(x1, x2, x3, x4, x5, x6) = 0 x 3 = X3(x4, x5, x6) . (6.2.3) How does one know this is possible? The equations mi ght be very complicated, some of the coordinates might not appear in some of the equations, and so on. We look first at three examples in E 3 then in Example 4 we return to our question. Example 1 : Consider these two constraint equations in E3, [this is the same example as appears in (1.9)] a(x 1,x2,x3) = x12 + x22 + x32 - 22 = 0 b(x1,x2,x3) = (x1-2)2 + x22 + x32 - 22 = 0 . ( C . 1 ) Each constraint represents a sphere of radius 2 and so is a 2-dimensional smooth surface in E 3. The first sphere is centered at the origin, while the second has its center at (2,0,0). Subtracting the first equation from the second gives -4x 1+4 = 0 so x 1= 1, then inserting this into both equations gives 1 + x 22 + x32 - 22 = 0 1 + x 22 + x32 - 22 = 0 ⇒ x22 + x32 = 3 and x 1 = 1 . (C.2) The intersection of the two constraint surfaces is a circle, a 1-dimensional smooth surface in E 3. Given a(x 1,x2,x3) = 0 and b(x 1,x2,x3) = 0, is it possible to eliminate x 1 and x2 from the two constraint equations? The answer is yes, as follows, x1 = X1(x3) = 1 x2 = X2(x3) = ± 3 - x32 . ( C . 3 ) The following schematic drawing show s the intersection surface in red, (C.4) 97 Notice in this example that ther e are two possible solutions for x 2. One can regard the intersection surface (the circle x 22 + x32 = 3) as having two pieces (front half, back half) and the ± sign selects one of these pieces. It is easy to imagine more complicat ed examples for two constraints in E3 where perhaps one or both of the constraint surfaces contain disjoint pieces (e.g. a hyperboloid), and where the intersection surface also has several disjoint pieces (ellipsoid intersecting a hyperboloid) . Example 2 : Consider these two constraint equations in E3 : a(x1,x2,x3) = x12 + x22 + x32 - 22 = 0 b(x1,x2,x3) = x 22 + x32 - 12 = 0 . ( C . 5 ) The second equation is missing the x 1 coordinate. The first equation describes the same origin-centered radius 2 sphere of the previous example. The second e quation is that of a circle of radius 1, but as a function of three variables it is in fact a cylinder whose axis is the x 1 axis. If x 2 and x3 lie on the circle, any value of x 1 satisfies the second equation. Missing coordina tes result in surfaces which are "extruded" in the dimensions of the missing coordinates. Subtracting the second equation from the first gives x 12 - 3 = 0 so x1 = ± 3 , then inserting this into both equations gives 3 + x 22 + x32 - 22 = 0 x 22 + x32 - 12 = 0 ⇒ x22 + x32 = 1 and x 1 = ± 3 . (C.6) Given a(x 1,x2,x3) = 0 and b(x 1,x2,x3) = 0, is it possible to eliminate x 1 and x2 from the two constraint equations? The answer is yes, as follo ws (the signs are independent), x1 = X1(x3) = ± 3 x2 = X2(x3) = ± 1 - x32 . ( C . 7 ) The following schematic drawing show s the intersection surface in red, (C.8) In this example there are two possible solutions for x 1 and two for x 2 so overall there are four solutions. The intersection surface consists of the two circ les each having a front half and a back half. 98 Example 3 : Consider these two constraint equations in E3 : a(x 1,x2,x3) = x12 + x22 + x32 - 22 = 0 b(x1,x2,x3) = x12 + x22 - 12 = 0 . ( C . 9 ) Now x 3 is missing from the second equation instead of x 1. Subtracting the second equation from the first gives x32 - 3 = 0 so x 3 = ± 3 , then inserting this into both equations gives x 12 + x22 + 3 - 22 = 0 x12 + x22 - 12 = 0 ⇒ x12 + x22 = 1 and x 3 = ± 3 . (C.10) Given a(x 1,x2,x3) = 0 and b(x 1,x2,x3) = 0, is it possible to eliminate x 1 and x2 from the two constraint equations? The answer is yes, as follows, where α is an arbitrary real parameter in [-1,1], x1 = X1(x3 = ± 3 ) = α -1 ≤ α ≤ 1 x2 = X2(x3 = ± 3 ) = ± 1 - α2 . // the two ± signs are independent (C.11) In the previous two examples the argument x 3 was a free parameter whos e variation mapped out the smooth constraint intersection surface piece(s). In Example 3 x 3 is fixed at ± 3 and a new parameter α must be introduced to map out the intersection surf aces. There are again four solution half-circles. The following schematic drawing show s the intersection surface in red, (C.12) Example 4 : Now consider these three constraint equations in E6 : a(x 1, x2, x3, x4, x5, x6) = 0 b(x1, x2, x3, x4, x5, x6) = 0 c(x1, x2, x3, x4, x5, x6) = 0 . (6.2.2) (C.13) We assume that each equation describes a smooth surface of dimension 5 in E 6 (each is a hypersurface) and each surface may have several pieces. In a problem with constraints, one assumes that there is some 99 non-null surface of intersection of all the constraint surfaces. This intersection surface may have multiple pieces and each piece (barring degenerate cases) is itself a smooth surface of dimension 3 in E6. A candidate solution point r must lie on one of these intersection surface pieces. In general, each constraint knocks down the dimension of the intersection surface by one degree of freedom. So let us assume that x 4, x5, x6 are the last three components of some point(s) on the overall intersection surface. Each such point of course has some x 1, x2, x3 components. If x 4, x5, x6 are varied slightly, x 1, x2, x3 will also vary slightly, and this is what we mean by writing the functions x 1 = X1(x4, x5, x6) x2 = X2(x4, x5, x6) x3 = X3(x4, x5, x6) . (6.2.3) (C.14) If the intersection surface has multip le pieces which have the same x 4, x5, x6 value, then any of the three functions Xn above may be multi-valued, as occurred in our earlier examples. Conver sely, if we select a triplet x4, x5, x6 for which there are no points on the intersection surface, the solution x 1, x2, x3 does not exist. We don't care about such points in E6. So this then is what we mean in Section 6.2 (a) when we say above (6.2.3) that one can use the three constraint equations to eliminate the three variables x 1, x2, x3. In the concluding equations (6.2.11) one can regard the functions X1, X2, X3 as being any of the multi-valued functions just discussed if the intersection surface has multiple pieces. For example, we had in (6.2.11) that ⎝⎜⎛ ⎠⎟⎞ 0 0 0 0 = ⎝⎜⎛ ⎠⎟⎞f1 f2 f3 f4 a1 a2 a3 a4 b1 b2 b3 b4 c1 c2 c3 c4 ⎝⎜⎛ ⎠⎟⎞ X1 4 X2 4 X3 4 1 part of (6.2.11) (C.15) and here X1 4 refers to ( ∂/∂x4)X1(x4, x5, x6) evaluated at coordinates x 4, x5, x6 for some candidate solution point r on the intersection surface. The fact that the matrix in (C.15) must have zero determinant is not affected by the possible existence of multiple solution X1, X2, X3 functions. Notice that the derivatives appearing in the matr ix are also unaffected by the possibility of multi- valued functions for X1, X2, X3. In the same manner as above, we can "eliminate" any triplet of variables x i, xj, xk ( i ≠ j ≠ k) from the three constraint equations (C.13). In the general case where there are S-1 constraints and r has N components x 1,x2...xN with N > S, the dimensionality of the intersection surface is N - (S-1) in EN. The conclusions of the Theorem 3 proof outlined in Section 6.2 (b) are similarly not affected by the possible existence of an intersection surface having multiple pieces and the functions Xn possibly having multiple values. 100 Footnote: Why is a(x 1, x2, x3, ...xN) = 0 an (N-1)-dimensional surface in EN ? (C.16) This fact probably seems obvious to the reader and is certainly obvious in the first three examples above. Here we attempt a simple "engineering" explanation. We assume that a constraint equation is a "smooth" equa tion, meaning it is continuous and differentiable in all its arguments except perhaps at certain isolated locations which can be handled on an ad hoc basis either by fiat or by taking limits. Consider that, since a( r) = constant (namely 0), da = 0 so 0 = da = Σ i=1N aidxi , // a i = ai(r) = ai(x1, x2, x3, ...xN) (C.17) where a i is usual means ∂ a/∂xi. Let r be some point for which a( r) = 0. At this point we shall assume that there exists at least one a s(r) which is non-zero. If all ai(r) = 0, then from (1.2) that r is a point on a surface which has a null normal vector, which is not possible, so such an r value could not lie on a(r ) = 0. Then (C.17) may be written dxs = (1/as) Σi≠s aidxi . ( C . 1 8 ) Create a coordinate system whose origin is located at this point r for which a( r) = 0, and whose axes are aligned with those of EN. Now imagine an arbitrary tin y displacement of all the dx i other than dx s. For this set of N-1 dx i there is only one possible dx s which causes the point r + dr to satisfy the equation a(r+dr) = a( r) + da = 0 + 0 = 0, where d r ≡ (dx1, dx2.....dxN). That one possible dx s value is that given by (C.18). Imagine repeating this process for a continuum of values for the dx i other than dx s, and in each case we obtain the unique solution dx s from (C.18). We can think of the x s axis as being "vertical" and all the other axes being "horizontal" inasmuch as they are all perpendicular to the x s axis. The vectors d r generated in this manner comprise a tiny patch of a "plane" of di mension N-1 (a hyperplane) in EN which contains the point r. This is so because Σi=1N aidxi = a • dr = 0 is the equation of a plane in EN passing through our origin at r with a being that plane's normal vector. More generally, n^ • r = d is the equation of a hyperplane in En having a unit normal n^ and whose closest approach to the origin is distance d. Thus at a point r satisfying a( r) = 0 we have constructed a tiny neighborhood of nearby points r which also satisfy a(r ) = 0. This neighborhood is a patch of a hyperplane in EN which is certainly a piece of "surface" in EN having dimension N-1. By repeating this process using a mesh of points ri satisfying a(ri) = 0, we then map out a triangulated polyhedral surface of dimension N-1 in EN. In the limit the mesh size goes to 0, we arrive at a smooth surface of dimension N-1 in EN and that surface is a(r ) = 0. The tiny planar patch at any point on the surface is part of the "tangent plane" to the surface at that point. The constraint a( r) = 0 is assumed to be "locally smooth" in the immediate region of any point r on the operational constraint surface, so a small local planar region (an open set) can be constructed around any such point. This is the basic idea of a surface bein g a manifold M, and the set of vectors in the tangent plane of a point r on M comprise the "tangent space" of M at r, usually denoted by T rM. ( C . 1 9 ) 101 Appendix D: Lagrange Multipliers in Classical Mechanics Here is an outline of this A ppendix showing the equation number nearest the start of each subsection: Constraints and virtual displacements (D.1) Assumption that constraints do no virtual work (D.11) Comment on internal forces (D.13) First Application of Lagrange Multipliers: the forces of constraint (D.13) Summary of the above Lagrange Multiplier Application (D.25) Generalized coordinates, generalized forces, and Lagrange's Equations (D.27) Footnote: Derivation of the result (D.31) used above (D.40) Second Application of Lagrange Multipliers: non-holonomic constraints ( D.49) The Case of C holonomic constraints treated as if they were non-holonomic (D.55) Footnote: Carry out the δS = 0 functional variation of the action (D.56) Our main focus in this document is the u se of Lagrange Multipliers to find candidates r for which a scalar function f(r ) is stationary, df( r) = 0, subject to constraints. In this Appendix, however, we consider applications of Lagrange Multipliers which relate to classical mechanics for a system of N particles in the presence of constraints. In this case the function f( r) with r lying in EN is replaced by a functional f(φ) where the φ lie in a space of functions. The goal is to find the functions φ which create a stationary point of the functional, and one writes this as δf(φ) = 0. A certain amount of background information is required to put these applications into context, which we now present. Constraints and virtual displacements Let rk(t) be the position of the kth particle of an N particle system in 3 dimensional space. At any fixed time t, if the particles are unencumbered with constraints, the variables ( r1, r2..... rN) can be regarded as a single point in the space E3N and the "configuration space" of the system indicated by ( r1, r2..... rN) is all of this space E3N. All points in E3N are "legal" points for the system. A holonomic constraint has the form a( r1, r2..... rN, t ) = 0 . ( D . 1 ) At a given time t, this equation represents a surface of dimension 3N-1 in E3N. If there are s such holonomic constraints ai(r1, r2..... rN, t) = 0 i = 1,2...s (D.2) then at time t the system vector ( r1, r2..... rN) is constrained to lie on a surface of dimension 3N-s in E3N, as discussed in Appendix C, so each extra constraint lowers the operating surface dimension by 1. One can think of the term holonomic as meaning that the c onstraints are functions of "whole" coordinates like r2 as opposed to differential coordinates like d r2. The constraints which appear in (1.7) are holonomic constraints of the form a(x 1, x2, ... xn) = 0. 102 A non-holonomic constraint is one that is not holonomic. An inequality constraint like a( r1, r2..... rN, t) < 0 which might bound the particles to one side of a surface is an example of a nonholonomic constraint, but our interest lies more in nonholonomic constraints of the form Σk=1N bk• drk + btd t = 0 ( D . 3 ) in which differentials rather than "whol e" coordinates appear. The "coefficients" bk and bt can in general be functions of ( r1, r2..... rk, t). If it happens that such a constraint can be integrated by some hook or crook to give an equation of the form a( r1, r2..... rk, t) = 0, then it is regarded as being a holonomic constraint, so a non-holonomic constraint must be non-integrable. If there are m such non-holonomic constraints imposed on a system, one could write them as Σk=1N Aik • drk + Aitdt = 0 i = 1,2...m o r ( D . 4 ) Σk=1N Aik • r•k + Ait = 0 i = 1,2...m . Here bolded Aik represents 3 matrices each of which is m x N (rows x columns). Since these constraints are non-integrable, one cannot really described them in terms of surfaces in E3N. If one differentiates the holonomic constraints (D.2) one gets Σk=1N [∇(k)ai] • drk + [∂tai] dt = 0 i = 1,2...s o r ( D . 5 ) Σk=1N [∇(k)ai] • r•k + [∂tai] = 0 i = 1,2...s where ∇(k)≡ (∂/∂(rk)1, ∂/∂(rk)2, ∂/∂(rk)3) is a gradient with respect to the components of rk. These equations have the same form as those in (D.4), but these are holonomic because they are integrable to give ai(r1, r2..... rN, t) = 0. For convenience below, we shall now combine (D.4 ) and (D.5) into a single set of equations, Σk=1N Bik • drk + Bitdt = 0 i = 1,2....C // C = s+m o r ( D . 6 ) Σk=1N Bik • r•k + Bit = 0 i = 1,2....C , where B ik = Aik B it = Ait for i = 1,2...s non-holonomic B ik = ∇(k)ai Bit = ∂tai for i = s+1, s+2....C holonomic . (D.7) In general, the functions a i, Bik and Bit are all functions of ( r1, r2..... rk, t). 103 In the above we have followed the notation of Ray and Shamanna (2006) and we continue our path in the manner they propose. Suppose the set of N differentials {dr k} satisfies all of equations (D.6) for some time interval dt. Such a set {d rk} is then called a set of allowed displacements . Imagine that {dr' k} is some other distinct allowed displacement set. Elements of the difference set defined as {δrk} = {d rk} - {d r'k} k = 1,2...N (D.8) are then called virtual displacements . This is how Ray and Shamanna define the vague term virtual displacement which appears in many texts (about which they have much to say in their Section I A). It follows then that for any virtual displacement set { δrk} equations (D.6) take this form Σk=1N Bik • δrk = 0 i = 1,2....C (D.9) simply because the last term in (D.6) cancels out when one writes δrk = drk - dr'k. Notice that the δrk are not arbitrary differential displacements. They are differences of allowed displacements. Meanwhile, Newton's Law for a system of N part icles subject to constraints has the form m r••k = Fk + Rk k = 1,2...N (D.10) where Rk is the sum of all constraint forces applied to particle k and Fk is the sum of all other forces applied to particle k. Assumption that constraints do no virtual work Next, assume that ( Ray and Shamanna and others refer to this as being an "ideal constraint" situation), 0 = Σ k=1N Rk • δrk . ( D . 1 1 ) In Goldstein this equation corresponds to setting the sec ond term in (1-40) to zero. The main justification for making this assumption is that it is valid for pa rticles which form a rigi d body and for many other constraint situations. There seems to be no general justification of this equation for all constraint situations, and it is certainly not va lid if friction is present. The equation can be interpreted as saying that the total virtual work done by all the constraint forces on all the particles of the system vanishes. As a first simple example, consider two partic les 1 and 2 independently sliding down a frictionless static ramp. In this case for particle 1 the constraint force R1 is the normal force of the ramp acting on the particle, and this is cl early perpendicular to dr 1 for any allowed displacement set {dr k}, and thus R1 is also perpendicular to any δr1 which is part of any virtual displacement set { δrk}. In this case we have R1• δr1 = 0 and R2• δr2 = 0 so the terms in (D.11) are separately zero. A second example is more interesting. Particles 1 and 2 of arbitrary masses accelerate down (or up) a one dimensional ramp but are glued to the two ends of a massless stick so they comprise a rigid body. The upper mass is attached to a string which provides a constant tension T. Here is a picture (for m 2 = 2m1), 104 ( D . 1 2 ) Particle 2 exerts a constraint force T12 on particle 1 through the stick, and Particle 1 exerts a constraint force T21 on particle 2. Since the stick is massless, T21 = -T12 (T12+T21= mstick a). All allowed displacement sets {d r1, dr2} are of the form {d r, dr}, and therefore all virtual displacement sets are of the form { δr1, δr2} = {δr, δr}. In (D.12) we have drawn the two normal forces N i and the two total constraint forces Ri. It is no longer true that R1• δr1 = 0 and R2• δr2 = 0 as the picture shows. However, it is true that Σ k=12 Rk • δrk = [ Σk=12 Rk ] • δr = 0 because R1 and R2 have equal and opposite components along the ramp, so equation (D.11) is valid for this example. Now suppose the ramp is translating in some direction not along the ramp. The displacements d r1= dr2 are then no longer along the ramp so that Σk=12 Rk • drk ≠ 0. However, any virtual displacement set would still have δr1 = δr2 along the ramp and then again Σk=12 Rk • δrk = 0. Remember that a virtual displacement set is the difference of any two allowed di splacement sets, and in this subtraction the effect of the ramp translation cancels out. So equation (D.11) is still valid. Other examples are given in Ra y and Shamanna Section III. Reader Exercise: Derive the equations shown in (D. 12). Show that things makes sense at θ = 0, π/2, and also for cases m 1→0 and m 2→0 . (Use F1 = m1g + R1+ T, F2 = m2g + R2, a1 = a2.) Comment on internal forces In (D.10) we have classified the forces acting on particle k into two groups and we write Rk + Fk, where Rk are "constraint" forces and the F k are any "other" forces. Goldstein refers to the Fk forces as "applied" forces in (1-39), while Ray and Shamanna call them "ext ernal" in Sec IV. Consider a system consisting of N masses all connected by a network of massless sp rings. The force acting on particle k due to all the springs to which it is attached would be best classified in the Fk "other" category. Were we to classify this force into the Rk group, we could not use the no-constraint-work property (D.11) since spring forces do work as the springs compress and expand. An example of such a system is the Solar System where the spring forces are represented by gravity. It would seem a misnomer to classify the internal forces in this case as "applied" or "external". However, if all th e springs are slowly adju sted until they become infinitely stiff, in that limit it would be better to classify the internal spring forces into the Rk constraint forces group, since these springs do no work and (D.11) applies. The system is then a "rigid body" and in this case, with the intern al forces classified as R k constraint forces, one benefits from the many simplifications which arise for rigid body mechanics. The internal forces in our spring example fall into 105 the category where Fij = - Fji and Fij is in the direction ri- rk . When this is not the case, the rigid body limit does not apply, but the internal forces can still be classified as "constraint" plus "other". First Application of Lagrange Multip liers: the forces of constraint Recall the C = m+s constraint equations (D.9), Σk=1N Bik • δrk = 0 i = 1,2...C . // C ≡ m+s = total number of constraints (D.9) Multiply each of these C equations by a Lagrange Multiplier λ i and then add the equations to get 0 = Σi=1C λi [Σk=1N Bik • δrk] = 0 = Σ k=1N [ Σi=1C λi Bik ] • δrk . ( D . 1 3 ) Now subtract equation (D.13) from (D.11) [ 0 = Σ k=1N Rk • δrk] to get 0 = Σk=1N [ Rk - Σi=1C λi Bik] • δrk . ( D . 1 4 ) We know that the N vectors in the set { δrk} are not linearly independent because they must satisfy the equations (D.9), so we cannot claim from (D.14) that [ Rk - Σi=1C λi Bik] = 0 for each k. But now write out equation (D.14) in full detail, showing all vector components, 0 = Σk=1N Σj=13 [ (Rk)j - Σi=1C λi (Bik)j] (δrk)j . ( D . 1 5 ) Next, relabel the 3N terms in this sum us ing a single index n defined in this way (k,j) = (1,1), (1,2), (1,3), (2,1), (2 ,2), (2,3),..........(N,1), (N,2), (N,3) n = 1, 2, 3, 4, 5, 6 , 3N-2, 3N-1, 3N s o ( D . 1 6 ) n = 3(k-1)+j k = 1+ Int[(n-1)/3] j = 1 + Rem[(n-1)/3] . Then define (here the n value is imp lied by the k,j labels as per above) x n ≡ (rk)j Rn ≡ (Rk)j B in ≡ (Bik)j . ( D . 1 7 ) Equation (D.15) can then be expressed as 0 = Σ n=13N [ Rn - Σi=1C λi Bin] δxn ( D . 1 8 ) and (D.9) becomes 0 = Σ k=1N Σj=13 (Bik)j (δrk)j or 0 = Σ n=13N Bin δxn i = 1,2...C . (D.19) 106 Once again, we cannot set the square brackets in (D.18) to 0 because the 3N coordinates δxn are not linearly independent due to the C equations (D.19). At any given time t, due to these equations (D.19), there must exist at least one subset of the 3N variables δx n which are independent and the remaining C δxn are dependent. Assume that the dependent variables have indices n = n 1, n2,...nC. Consider then the following set of C equations involving the dependent- term square brackets in (D.18), [ R n - Σi=1C λi Bin] = 0 n = n 1, n2,...nC . ( D . 2 0 ) This represents C linear equations in the C parameters λi so, when all is said and done, there must exist a set of values { λi} which make these C equa tions be true. We assume these values will be λ i(t) and this is justified below when the full set of equations is considered. Now the terms in (D.18) indicated by n = n 1, n2,...nC all vanish due to (D.20), so (D.18) can be written, 0 = Σn≠n1,n2,....nC [ Rn - Σi=1C λi Bin] δxn . ( D . 2 1 ) But all the δ xn appearing in (D.21) are independent, so the square brackets here must also vanish. We then end up with all the square brackets in (D.18) being 0 (assuming the correct solution values of λi), [ R n - Σi=1C λi Bin] = 0 n = 1,2....3N or R n = Σi=1C λi Bin n = 1,2....3N or ( Rk)j = Σi=1C λi (Bik)j k = 1,2..N j = 1,2,3 or Rk = Σi=1C λi Bik k = 1,2...N . (D.22) Thus the constraint forces are certain linear combinations of the B ik constraint functions shown in (D.7). Meanwhile, Newton's Law (D.10) says that m kak = Fk + Rk where the Fk are the non-constraint forces for the problem. Inserting (D.22) gives mkr••k = Fk + Σi=1C λiBik . Bik = ⎩⎪⎨⎪⎧ Aik i = 1 to m (non-holonomic) [∇(k)ai] i = m+1 to C (holonomic) (D.23) The problem is then summarized in the following set of equations from (D.23) and (D.6), mkv•k(t) = Fk(r1, r2..... rN, t) + Σi=1C λi(t)Bik(r1, r2..... rN, t) k = 1..N 3N equations r•k(t) = vk(t) k = 1..N 3N equations Σk=1N Bik(r1, r2..... rN, t) • r•k(t) + Bit(r1, r2..... rN, t) = 0 i = 1,2....C . C equations ( D . 2 4 ) 107 This can be regarded as a set of 6N+C scalar fi rst-order ODEs in the single variable t for the 6N+C unknown scalar functions rk(t), vk(t) and λi(t). Although the first 3N equations happen to be linear in the λi(t), the equations are in general non-linear in the rk(t) so this is then a non-linear s ystem of ODE' s. Nevertheless, a well-known existence theorem (e.g., Ince Secti on 3.3) states that, provided the Fk, Bik and Bit functions are continuous and diff erentiable in all their arguments, a unique solution exists for any set of initial values rk(0), vk(0) and λ i(0). As (D.22) shows, the λi(0) are related to the constraint forces at time 0. The uniqueness of the solution agrees with our intuitio n that a physical system will evolve in a unique manner. Finding an analytic form for the solution might not be possible, but a numeric form can always be obtained. Summary of the above Lagrange Multiplier Application We have used a set of Lagrange Multipliers as an aid to solving a mechanics problem involving N particles with C constraints. As already noted, the above presentation is based on Section IV of Ray and Shamanna. The steps were as follows, considered more generically: 1. One has a sum Σn=1NRnδxn = 0 . ( D . 2 5 ) 2. One also has a set of sums Σn=1NBinδxn = 0 for i = 1,2...C with C < N (constraints in the above case). This implies that one can find C of the δ xn which are linearly dependent on the other δxn. 3. If λi are C arbitrary factors, certainly Σi=1Cλi(Σn=1NBinδxn) = 0 based on 2. These λi are "Undetermined Lagrange Multipliers". 4. Subtract sum 3 from sum 1 to get the new sum Σ n=1N[ Rn - Σi=1Cλi Bin] δxn = 0. 5. The C λi are selected so that [..] = 0 in the C dependent δxn terms of this sum. The remaining sum has only independent δxn terms so [..] = 0 for all those terms as well. Then all [...] = 0 in the sum of step 4. 6. The conclusion is that [ R n - Σiλi Bin] = 0 for all n, assuming the solution λi values of step 5. In contrast with our main topic of finding stationary points df = 0 of a function f( r), in the Lagrange Multiplier application just presente d there is no explicit function f( r) we are making stationary. On the other hand, one can regard (D.11) that ΣkRk • δrk = 0 as the statement δ W = 0 where δW = ΣkRk • δrk is the total differential "virtual work" done by the constr aint forces during a differential virtual displacement {δrk}. Using Newton's law (D.10) this equation can then be written 0 = δW = Σk[mkak - Fk] • δrk // D'Alembert (D.26) 108 in which form it is known as D'Alembert's Principle (1743, Goldstein (1-42)). The analog of f( r) is then an action-type integral W of the differential virtual work which integral one renders stationary to get δW = 0. The notion of "the action" is discussed at the end of this appendix. We now present another application of the Lagrange Multiplier technique outlined in (D.25). As before, a certain amount of background is needed befo re the application can be demonstrated. Generalized coordinates, generalized forces, and Lagrange's Equations Usually the Cartesian components of the particle positions r k are not the most convenient variables to use in solving a problem. Angles are often more useful. With C holonomic constraints one can replace the 3N variables {x n} in (D.17) with some set of 3N-C independent variables {q n} where x n = xn(q1, q2, ....q3N-C, t) n = 1,2..3N (D.27) and one can then rework the above presentation in these new generalized coordinates q j which then won't always have the dimensions of length. For example, one then has, for virtual displacements, δxn = Σj=13N-C(∂xn/∂qj) δqj n = 1,2..3N or δr k = Σj=13N-C(∂rk/∂qj) δqj k = 1,2..N // as in Goldstein (1-44) . (D.28) The usual presentation of Lagrangian dynamics with holonomic constraints avoids dealing with the forces of constraint and also avoids the need for using a set of Lagrange Multipliers. That presentation appears in Section V of the Ray and Shamanna and Goldstein p 17 (and many other places) and proceeds like this: 1. Start with Σ kRk • δrk = 0 in (D.11) which we write in the form of D'Alembert's Principle (D.26) mentioned above, 0 = Σ k=1N [mkak - Fk] • δrk (D.26) where Fk is the total non-constraint force acting on particle k. 2. Install (D.28) into D'Alembert's Principle to get 0 = Σ k=1N [mkak - Fk] • Σj=13N-C(∂rk ∂qj )δqj = Σj=13N-C { Σk=1Nmkak • ∂rk ∂qj - Σk=1N Fk • (∂rk ∂qj ) } δ qj = Σj=13N-C { Σk=1Nmkak • ∂rk ∂qj - Qj } δqj ( D . 2 9 ) where Q j ≡ Σk=1N Fk • (∂rk ∂qj ) . ( D . 3 0 ) 109 The Qj defined above are the generalized forces associated with the generalized coordinates q j . 3. Make use of the following non-obvious result (derived in Footnote below starting at (D.40)), Σk=1N mkak • ∂rk ∂qj = d dt (∂T ∂q•j ) - ∂T ∂qj ( D . 3 1 ) where T is the total kinetic energy of the N particles, T ≡ (1/2) Σ k=1N mk (vk • vk) . ( D . 3 2 ) Then (D.29) reads, 0 = Σ j=13N-C { d dt (∂T ∂q•j ) - ∂T ∂qj - Qj } δqj . ( D . 3 3 ) Since the δqj are independent variables, {..} = 0 and we end up with one form of "Lagrange's Equations" d dt (∂T ∂q•j ) - ∂T ∂qj = Qj j = 1,2....3N-C . // as in Goldstein (1-50) (D.34) If the "other" forces can be derived from a potential V(r1, r2, ..... rN) such that Fk = - ∇(k)V(r1, r2, ..... rN) k = 1,2..N // no t or r•i dependence in V (D.35) then from (D.30), Q j ≡ Σk=1N Fk • (∂rk ∂qj ) = - Σk=1N ( ∂V ∂rk ) • (∂rk ∂qj ) = - ∂V ∂qj j = 1,2....3N-C (D.36) where V(q i) ≡ V(rk(qi)) is a function of the generalized coordinates q i. Since V is not a function of the q•i we know that ∂V ∂q•j = 0. Using this fact, and (D.36) that Q j = – ∂V ∂qj , (D.34) can be written d dt (∂(T-V) ∂q•j ) - ∂(T-V) ∂qj = 0 . j = 1,2....3N-C (D.37) The last step is to define the classical Lagrangian L ≡ T - V ( D . 3 8 ) to obtain 110 d dt (∂L ∂q•j ) - ∂L ∂qj = 0 j = 1,2....3N-C . // as in Goldstein (1-53) (D.39) which is the more conventional fo rm of "Lagrange's Equations" w ith C holonomic constraints. The equations of (D.39) are usually called the Euler-Lagrange Equations . One can define generalized (canonical, conjugate) momenta pj ≡ ∂L ∂q•j so the Lagrange equations (D.39) are just p•j = ∂L ∂qj or p• = ∇(q)L . If V = V(q i), then p j = ∂T ∂q•j and then if T = 1 2 Σkmkq•k2 one has the familiar result p j = mjq•j. When T = T(q•i) one finds Q j ≡ - ∂V ∂qj = ∂L ∂qj = p•j so p• = Q, which is the generalized Newton's Law in terms of generalized coordinates and generalized forces. Footnote: Derivation of the result (D.31) used above Define the 3N-C independent generalized coordinates q i as in (D.27) according to r k = rk(q1, q2....q3N-C, t) ≡ rk(qi,t ) . ( D . 4 0 ) Compute the total time derivative of r k , vk = r•k = drk dt = Σj ∂rk ∂qj ∂qj ∂t + ∂rk ∂t = Σj ∂rk ∂qj q•j + ∂rk ∂t . // = vk(qi, q•i, t) (D.41) Then based on this result comput e these two partial derivatives of v k , ∂vk ∂qi = Σj ∂2rk ∂qi∂qj q•j + ∂2rk ∂qi∂t ( D . 4 2 ) ∂v k ∂q•i = ∂rk ∂qi . ( D . 4 3 ) Next, the total kinetic energy of all k particles is given by (D.32), T ≡ (1/2) Σ k mk (vk • vk) (D.32) so that ∂T ∂q• i = Σk mk vk • ∂vk ∂q•i = Σk mk vk • ∂rk ∂qi . // using (D.43) (D.44) 111 Apply d/dt to the above, with ak = r•k : d dt (∂T ∂q•i ) = Σk mk ak • ∂rk ∂qi + Σk mk vk • d dt (∂rk ∂qi ) . (D.45) Compute the total time derivative appearing in the second term, d dt (∂rk ∂qi ) = Σj ∂2rk ∂qj∂qi q•j + ∂2rk ∂t ∂qi . ( D . 4 6 ) Meanwhile, ∂T ∂q i = Σk mk vk • ∂vk ∂qi = Σk mk vk • [ Σj ∂2rk ∂qi∂qj q•j + ∂2rk ∂qi∂t ] // (D.42) = Σ k mk vk • d dt (∂rk ∂qi ) . // (D.46) (D.47) Subtract (D.47) from (D.45) to get d dt (∂T ∂q•i ) - ∂T ∂qi = Σk mk ak • ∂rk ∂qi ( D . 4 8 ) which is the result quoted above in (D.31). Second Application of Lagrange Mu ltipliers: non-holonomic constraints This is our final Lagrange Multiplie r example, and it exactly follows the generic prescription of (D.25). The general topic is well described in Goldstein Chapter 2. As the starting point we invoke Hamilton's Principle which is this : δS = 0 where S = ∫t1 t2 dt L(qn(t), q•n(t), t ) n = 1,2....M (D.49) where there are M independent generalized coordinates q n. The notation is a shorthand to imply that L is 112 a function of all the qn(t) and all the q•n(t) functions. One might better write L({q n(t)}, {q•n(t)}, t ) where the {..} indicate sets. L = T - V is the Lagrangian (D.38) with T(q•i) and V(q i) being the kinetic and potential energies associated with the M generalized coordinates. We shall assume that M = N-D where there were D holonomic constraints on N initial problem coordinates. As shown in (D.61) below ( or see Goldstein p 41 (2-20') based on p 37 (2-15) ) the equation δS = 0 results in the following sum being 0, 1. Σn=1M [ ∫t1 t2 dt { d dt (∂L ∂q•n ) - ∂L ∂qn } ] δqn(t) = 0 . (D.50) The number 1. on the left refers to the first step outlined in (D.25). In this equation the δ qn(t) are M independent functions of t ("path variations"), restricted only by δqn(t1) = δ qn(t2) = 0. In order for (D.50) to be true for M arbitrary independent functions δqn(t), one must have {} = 0, d dt (∂L ∂q•n ) - ∂L ∂qn = 0 n = 1,2....M (D.51) which are the Euler-Lagrange Equations (D.39). Thus it is that the Lagrange Equations can be derived from Hamilton's Principle in the case of only holonomic constraints. But suppose that in addition to the D holonomic constraints there are C non-holonomic constraints, so only M-C of the δ q n are independent. In this case the above conclusion (D.51) is incorrect. 2. The C non-holonomic constraints result in C sums Σ n Ani δqn = 0 i = 1,2...C analogous to (D.9) above. 3. If λ i are C arbitrary factors, certainly Σiλi(ΣnAinδqn) = 0 based on 2. These λi are "Undetermined Lagrange Multipliers". Then ∫t1 t2 dt ( ΣiλiΣnAinδqn ) = 0 as well. 4. Subtract sum 3 from sum 1 to get the new sum Σ n=1M [ ∫t1 t2 dt { d dt (∂L ∂q•n ) - ∂L ∂qn - Σi=1C λiAin } ] δ qn(t) = 0 . // as in Goldstein (2-26) (D.52) 5. The C λi are selected so that {..} = 0 in the C dependent δ qn(t) terms of this sum. The remaining sum has only independent δ qn(t) terms so {..} = 0 for all those terms as well. Then all {...} = 0 in the sum of 4. 6. The conclusion is that d dt (∂L ∂q•n ) - ∂L ∂qn = Σi=1CλiAin . n = 1,2....M . // as in Goldstein (2-30) (D.53) 113 This then is the corrected form of the Euler-L agrange Equations in the presence of C non-holonomic constraints. Recalling from (D.7) that Bik = Aik for non-holonomic constraints, the second line of (D.22) with B in = Ain reads R n = Σi=1CλiAin . Thus one can interpret the sum on the right of (D.53) as being the force of constraint R n arising from the C non-holonomic constraints. The forces of the D holonomic constraints are already incorporated into the left side of (D.53) in the reduction of independent coordinates q n from N to M = N-D. The integral S shown in (D.49) is called "the acti on" and Hamilton's Principle is a classical mechanics incarnation of the Principle of Least Action . In relativistic quantum field theory the action is a spacetime integral of the Lagrangian density L(φn(xμ), ∂μφn(xμ)) where the coordinates q n of L are replaced by fields φn(xμ), and t is replaced by xμ, a point in spacetime. In this case the Lagrange Equations (D.39) take the form (the second line is for comparison), ∂ ∂x μ ( ∂L ∂(∂μφn) ) - ∂L ∂φn = 0 where ∂μ ≡ ∂ ∂xμ // Bjorken and Drell (11.30) (D.54) d dt (∂L ∂q•n ) - ∂L ∂qn = 0 t → xμ q n(t) → φn(xμ) q•n(t) → ∂μφn(xμ) (D.39) The Case of C holonomic constraints treated as if they were non-holonomic Suppose there are N generalized coordinates q i with C holonomic constraints and 0 non-holonomic constraints. As noted in (D.6) and (D.7), the holonom ic constraints can be written in differential form so they look just like non-holonomic constraints. For generalized coordinates q n these two equations would appear as Σn=1N Bindqn + Bitdt = 0 i = 1,2....C // C = s+m o r ( D . 6 ) ' Σn=1N Bin q•n + Bit = 0 i = 1,2....C where B in = Ain B it = Ait for i = 1,2...s non-holonomic B in = ∂ai/∂qn Bit = ∂tai for i = s+1, s+2....C holonomic . (D.7)' Setting s = 0 we can carry out the procedure of the previous section with A ik replaced by ∂ai/∂qk to obtain this version of (D.53) re ndering the action S stationary, d dt (∂L ∂q•n ) - ∂L ∂qn = Σi=1Cλi ∂ai ∂qn n = 1,2....N (D.55) where the Lagrange multipliers λi are chosen to make this equation be valid for a set of C dependent q n coordinates, as in Step 5 above. Only N-C of the coordinates q n are independent but the equations are 114 written for all N coordinates. Here the functions a i(q1, q2, .....qn, t) = 0 are the C holonomic constraint equations. Footnote: Carry out the δS = 0 functional variation of the action . The variation in each function q n(t) is parameterized in terms of an arbitrary independent function ηn(t) (but which vanishes at both time integra tion endpoints) and one scalar parameter α : qn(t,α) = qn(t,0) + α ηn(t) δqn(t) = α ηn(t) ∂qn ∂α = ηn(t) ηn(t1) = ηn(t2) = 0 q•n(t,α) = q•n(t,0) + α η•n(t) δq•n(t) = α η•n(t ∂q•n ∂α = η•n(t) . (D.56) Then S(α) ≡ ∫t1 t2 dt L(qn(t,α), q•n(t,α), t ) // action (D.49) ∂S ∂α = ∫t1 t2 dt Σn ( ∂L ∂qn ∂qn ∂α + ∂L ∂q•n ∂q•n ∂α ) = ∫t1 t2 dt Σn ( ∂L ∂qn ηn + ∂L ∂q•n dηn dt ) . (D.57) But for the second term do a parts integration, ∫t1 t2 dt ( ∂L ∂q•n dηn dt ) = [ ∂L ∂q•n ηn(t) ]t2 t1 - ∫t1 t2 dt ( d dt ∂L ∂q•n ) ηn(t) . (D.58) The "parts" terms vanish since ηn vanishes at both endpoints, so (D.57) becomes, ∂S ∂α = ∫t1 t2 dt Σn ( ∂L ∂qn - d dt ∂L ∂q•n ) ηn(t) (D.59) so δS ≡ ∂S ∂α dα = { ∫t1 t2 dt Σn ( ∂L ∂qn - d dt ∂L ∂q•n ) ηn(t) } dα. (α/α) = { ∫t1 t2 dt Σn ( ∂L ∂qn - d dt ∂L ∂q•n ) δqn(t) } dα/α . (D.60) Then, 115 δS = 0 ⇔ ∫t1 t2 dt Σn [ ( ∂L ∂qn - d dt ∂L ∂q•n ) ηn(t) ] = 0 . o r ( D . 6 1 ) δS = 0 ⇔ ∫t1 t2 dt Σn [ ( ∂L ∂qn - d dt ∂L ∂q•n ) δqn(t) ] = 0 . But the functions δqn(t,α) = αηn(t) are arbitrary and independent functions of t, so the integral can only vanish if the function in the parentheses vanishes for each term in the sum on n, so δS = 0 ⇔ ( ∂L ∂qn - d dt ∂L ∂q•n ) = 0, n = 1,2...m (D.62) which are the Euler-Lagrange Equations (D.39) and (D.51). For more detail see Goldstein Chapter 2. 116 References The method of Lagrange multipliers often makes cam eo appearances in textbooks on the calculus of multiple variables. There are also several pages and pdf's on the web, search on "Lagrange Multipliers". Whole books on the use the Lagrange multipliers for sp ecialized applications may be found at the usual book vendor websites. The links below we re last verified 25 Oct 2016. J.D. Bjorken and S.D. Drell, Relativistic Quantum Fields (McGraw-Hill, New York, 1965). R.C. Buck with E.F. Buck, Advanced Calculus , 2nd Ed. (McGraw-Hill, New York, 1965). The discussion of extremum problems with constraints occupies pp 359-363. In the 3rd Edition (McGraw-Hill, New York, 1978), reprinted by (Waveland Press, Long Grove IL, 2003), this discussion has moved to pages 536-540. This excellent and now classic book was first published in 1956. H. Goldstein, Classical Mechanics (Addison-Wesley, Boston, 1950), another classic. There is a 3rd edition 2001 by Goldstein, Safco and Poole, but our page and equation references are to the 1950 edition. E.L. Ince, Ordinary Differential Equations (Longmans, Green & Co. (Ltd.), London, 1927). This classic text became a Dover paperback in 19 56 and is available online for $4. S. Jensen, An Introduction to Lagrange Multipliers, http://www.slimy.com/~steuard/teaching/tutorials/Lagrange.html ( includes video lectures) P. Lucht, Tensor Analysis and Curvilinear Coordinates (2016), http://user.xmission.com/~rimrock . If not there, search on . P.M. Morse and H. Feshbach, Methods of Theoretical Physics ( McGraw-Hill, New York, 1953). S. Ray and J. Shamanna, "On Virtual Displaceme nt and Virtual Work in Lagrangian Dynamics", European Journal of Physics, 27 (2006) 311-329. Also https://arxiv.org/abs/physics/0510204v2 . G.E. Shilov, Linear Algebra (Dover Publications, 1977). M.R. Spiegel, Mathematical Handbook of Formulas and Tables ( Schaum's Outline Series, McGraw-Hill, New York, 1968). This is now at 4th Ed (2012) with two coauthors and new reference numbers. W. F. Trench, The Method of Lagrange Multipliers (2013), http://digitalcommons.trinity.edu/mono/7/ , Supplement 2 (31p). His Theorem 1 (Eq. 6) is basically the same as our Theorem 1 (1.24). Wiki on Lagrange Multipliers: https://en.wikipedia.org/wiki/Lagrange_multiplier M.W. Zemansky, Heat and Thermodynamics , 5th Ed. (McGraw-Hill, New York, 1968). Our page references are to this edition, but there is a 7th Ed . (1997) with coauthor R.H Dittman and an 8th Ed (2011) in the Special Indian Edition (SIE) series.