Advanced Mechanics, Lecture Notes - Eric Poisson, Department of Physics University of Guelph, Lecture notes of Physics

Download Advanced Mechanics, Lecture Notes - Eric Poisson, Department of Physics University of Guelph and more Lecture notes Physics in PDF only on Docsity! Advanced mechanics PHYS*3400 Lecture notes (January 2008) -15 -10 -5 0 5 10 15 -3 -2 -1 0 1 2 3 p θ pppp Eric Poisson Department of Physics University of Guelph Contents iii 3.5.2 Hamilton-Jacobi equation 148 3.5.3 Case study: Linear pendulum 149 3.6 Problems 152 A Term project: Motion around a black hole 155 A.1 Equations of motion 155 A.2 Circular orbits 156 A.3 Eccentric orbits 156 A.4 Numerical integration of the equations of motion 157 Chapter 1 Newtonian mechanics 1.1 Reference frames An important aspect of the fundamental law of Newtonian mechanics, F = ma, (1.1.1) is that it is formulated in a reference frame which is either at rest or moving with a uniform velocity (the velocity must be constant both in magnitude and in direc- tion). Such frames are called inertial frames. A reference frame is a set of three axes attached to a point O called the origin. The position of the origin in space is arbitrary, but some specific choices are sometimes convenient. For example, when describing a system of N bodies it is usually a good idea to place the origin at the centre of mass (which will be introduced below). The origin of an inertial frame is either fixed or moving uniformly relative to another inertial frame. The orienta- tion of the axes is also arbitrary, but some specific choices can again simplify the description. For example, when studying the motion of a particle in a gravitational field it is convenient to align one of the coordinate axes with the direction of the gravitational force. The coordinate axes define a set of basis vectors x̂, ŷ, and ẑ. (These are some- times denoted i, j, and k.) These vectors point in the directions of increasing x, y, and z, respectively, and they all have a unit norm: x̂ · x̂ = ŷ · ŷ = ẑ · ẑ = 1; this property is indicated by the “hat” notation. Relative to a choice of origin O, a particle has a position vector r(t) at time t. This is decomposed in the basis as r(t) = x(t)x̂ + y(t)ŷ + z(t)ẑ. (1.1.2) The functions x(t), y(t), and z(t) are the particle’s coordinates relative to the refer- ence frame. The coordinates change as t varies, and the particle traces a trajectory in three-dimensional space. The central goal of Newtonian mechanics is to deter- mine this trajectory, assuming that the force F acting on the particle is known at all times. The particle’s velocity vector is v(t) = dr dt = ẋ(t)x̂ + ẏ(t)ŷ + ż(t)ẑ, (1.1.3) where we have introduced the notation ẋ = dx/dt = vx; we shall also use ṙ = dr/dt as an alternative notation for the vector v. The particle’s momentum vector is defined by p = mv, (1.1.4) where m is the particle’s mass. The particle’s acceleration vector is a(t) = dv dt = ẍ(t)x̂ + ÿ(t)ŷ + z̈(t)ẑ, (1.1.5) 1 2 Newtonian mechanics b r r S S0 1 1 0 Figure 1.1: Two reference frames, S0 and S1, separated by a displacement b. with the notation ẍ = d2x/dt2 = v̇x = ax. Newton’s equation, ma = F , has the mathematical structure of a system of second-order differential equations for the coordinates x(t), y(t), and z(t). To describe the particle’s trajectory, knowing the force, it is necessary to integrate these differential equations. Suppose that we have two reference frames, S0 and S1, separated by a dis- placement b (see Fig. 1.1). Relative to S1 the position vector of a particle is r1; relative to S0 it is r0. The transformation between the two position vectors is clearly r0 = b + r1, or r1 = r0 − b. (1.1.6) Suppose now that S1 moves relative to S0, so that the vector b depends on time. Since the position vectors also depend on time, Eq. (1.1.6) should be written as r1(t) = r0(t)− b(t). Taking a time derivative produces the transformation between the velocity vectors: v1 = v0 − ḃ. (1.1.7) Taking a second time derivative gives us the transformation between the acceleration vectors: a1 = a0 − b̈. (1.1.8) If S0 is an inertial frame, then the equations of motion for the particle as viewed in S0 are ma0 = F . In S1 the equations are instead ma1 = F − mb̈. (1.1.9) We see that Newton’s equation is preserved only if b̈ = 0, that is, if ḃ is a constant vector. In this case S1 moves relative to S0 with a constant velocity, and it is also an inertial frame. When, however, S1 is not inertial, the equations of motion do not take the Newtonian form. We have instead Eq. (1.1.9), which can be rewritten as ma1 = F + Ffictitious, with Ffictitious = −mb̈. The second term on the right can be thought of as a fictitious force that arises from the fact that the reference frame is not inertial. A well-known example is the centrifugal force, which arises in a rotating (and therefore non-inertial) frame of reference. We now consider a situation in which S1 and S0 are both inertial. We assume, in fact, that they share a common origin O, but that they differ in the orientation of the coordinate axes. A concrete example (see Fig. 1.2) is one in which S1 is obtained from S0 by a rotation around the z axis. In this case the basis vectors x̂1 and ŷ1 differ in direction from x̂0 and ŷ0. Similarly, the particle’s coordinates x1(t) 1.2 Alternative coordinate systems 5 (∂r/∂φ) = r2 sin2 φ+r2 cos2 φ = r2, and to get a unit vector we must divide ∂r/∂φ by r. We conclude that φ̂ = 1 r ∂r ∂φ = − sin φ x̂ + cos φ ŷ (1.2.8) is the desired basis vector. Exercise 1.2. Check that r̂ · φ̂ = 0. Let us now work out the components of the vectors r, v, and a in the basis (r̂, φ̂). According to Eqs. (1.2.2) and (1.2.7) we have r · r̂ = [ (r cos φ)x̂ + (r sinφ)ŷ ] · [ cos φ x̂ + sin φ ŷ ] = r cos2 φ + r sin2 φ = r. Similarly, Eqs. (1.2.2) and (1.2.8) give r · φ̂ = [ (r cos φ)x̂ + (r sin φ)ŷ ] · [ − sin φ x̂ + cos φ ŷ ] = −r sin φ cos φ + r sinφ cos φ = 0. From these results we infer that r = r r̂, (1.2.9) and this expression should not come as a surprise, given the meaning of the quanti- ties involved. Proceeding similarly with the vectors v and a, we find that they are decomposed as v = ṙ r̂ + rφ̇ φ̂ (1.2.10) and a = ( r̈ − rφ̇2 ) r̂ + 1 r d dt ( r2φ̇ ) φ̂ (1.2.11) in the new basis. As we have pointed out, the components of r in the polar basis are obvious, and the components of v also can be understood easily: The radial component of the velocity vector must clearly be vr = ṙ, and the tangential compo- nent must be vφ = rφ̇ because the factor of r converts the angular velocity φ̇ into a linear velocity. The components of the acceleration vector are not so easy to interpret. It is im- portant to notice that the radial component of the acceleration vector is not simply ar = r̈, and the angular component is not simply aφ = φ̈. It is a general observation that the components of the acceleration vector are not simple in nonCartesian coor- dinate systems. It should be observed that the radial component of the acceleration vector contains both a radial part r̈ and a centrifugal part −rφ̇2 = −v2φ/r. Exercise 1.3. Verify by explicit calculation that Eqs. (1.2.10) and (1.2.11) are correct. Suppose now that the force F has been resolved in the polar basis (r̂, φ̂). We have F = Frr̂ + Fφφ̂, (1.2.12) and Newton’s law F = ma breaks down into two separate equations, the radial component r̈ − rφ̇2 = Fr m (1.2.13) 6 Newtonian mechanics and the angular component d dt ( r2φ̇ ) = rFφ m . (1.2.14) These are the equations of motion for a particle subjected to a force F , expressed in polar coordinates (r, φ). When, for example, Fφ = 0 and the force is purely radial, then according to Eq. (1.2.14), r2φ̇ = rvφ is a constant of the motion. When, in addition, r̈ = 0 and the particle travels on a circle r = constant, then Eq. (1.2.13) reduces to rφ̇2 = v2φ/r = −Fr/m; this is the familiar equality between the centrifugal acceleration v2φ/r and (minus) the radial component of the force (divided by the mass). Exercise 1.4. Consider the spherical coordinates (r, θ, φ) defined by x = r sin θ cos φ, y = r sin θ sin φ, and z = r cos θ. Show that in this alternative coordinate system, the basis vectors are given by r̂ = ∂r ∂r = sin θ cos φ x̂ + sin θ sin φŷ + cos θ ẑ, θ̂ = 1 r ∂r ∂θ = cos θ cos φ x̂ + cos θ sin φŷ − sin θ ẑ, φ̂ = 1 r sin θ ∂r ∂φ = − sin φ x̂ + cos φŷ. Verify that these vectors are all orthogonal to each other. 1.3 Mechanics of a single body In this section we explore some consequences of the law F = ma when it applies to a single particle. 1.3.1 Line integrals We begin with a review of some relevant mathematics. Let A be a vector field in three-dimensional space. (A vector field is a vector that is defined in a region of space and which may vary from position to position in that region.) Let C be a curve in three-dimensional space, and let ds be the displacement vector along the curve. The displacement vector is defined so that ds is everywhere tangent to the curve, and such that its norm ds = |ds| is equal to the distance between two neighbouring points on the curve; the total length of the curve is the integral ∫ C ds. Now introduce ∫ 2 1 A · ds, the line integral of the vector field A between point 1 and point 2 on the curve C. Such integrals occur often in physics. In the present context the force F will play the role of the vector field A, and the particle’s trajectory will play the role of the curve C; we then have ds = dr = vdt and the line integral will be the work done by the force as the particle moves from point 1 to point 2. It is a fundamental theorem of vector calculus that if a line integral between two fixed points in space does not depend on the curve joining the points, then the vector field A must be the gradient ∇f of some scalar function f . This theorem is essentially a consequence of the identity ∫ 2 1 ∇f · ds = ∫ 2 1 df ds ds = f(2) − f(1) independently of the curve, 1.3 Mechanics of a single body 7 x = 1x = −1 x = −1 x = 1 y A A y θ Figure 1.3: Line integrals of a vector field A. which is a generalization of the statement ∫ b a (df/dx) dx = f(b)−f(a) from ordinary calculus. Another way of presenting this result is to say that if A = ∇f , then ∮ A · ds = 0 for any closed curve C in three-dimensional space. This last statement follows because if the curve C is closed, point 2 is identified with point 1, and ∮ ∇f · ds = f(1) − f(1) = 0. To illustrate these notions let us work through a concrete example. Consider the vector field A = (x, y) in two-dimensional space. We wish first to evaluate the line integral of A along the x axis, from x = −1 to x = +1 (see Fig. 1.3). The safest way to proceed is to first obtain a parametric description of the curve C, which in this case is the line segment that links the points x = ∓1. We may describe this curve in the following way: x(u) = −1 + 2u, y(u) = 0, where the parameter u is restricted to the interval 0 ≤ u ≤ 1. (The choice of param- eterization is arbitrary; we might just as well have chosen x as the parameter, but it is generally a good idea to keep the parameter distinct from the coordinates.) From these equations it follows that the displacement vector on C has the components dx = 2 du and dy = 0, so that ds = (2 du, 0). The vector field evaluated on C is A = (−1 + 2u, 0), and we have A · ds = 2(−1 + 2u) du. The line integral is then ∫ C A · ds = ∫ 1 0 2(−1 + 2u) du. Evaluating this ordinary integral is straightforward, and the result is zero. We therefore have ∫ C A · ds = 0 for this choice of curve linking the points (x = −1, y = 0) and (x = 1, y = 0). Let us now evaluate the line integral of A along a different curve C ′ which joins the same two endpoints (refer again to Fig. 1.3); we choose for C ′ a semi-circle of unit radius, which we describe by the parametric relations x(θ) = − cos θ, y(θ) = sin θ, 10 Newtonian mechanics This is the same statement as in Eq. (1.3.7), and we have established the work- energy theorem. In very many situations the line integral ∫ 2 1 F · dr is actually independent of the trajectory adopted by the particle to go from point 1 to point 2. In these situations we must have that F is the gradient of some scalar function f(r). We write f = −V , inserting a minus sign for reasons of convention, and express the force as F = −∇V (r). (1.3.8) The scalar function V is known as the potential energy of the particle. When F is expressed as in Eq. (1.3.8) the line integral of Eq. (1.3.5) becomes W12 = − ∫ 2 1 ∇V · dr = − [ V (2) − V (1) ] , and this is clearly independent of the particle’s trajectory: The total work done is equal to the difference V (1) − V (2) no matter how the particle moves from 1 to 2. Equation (1.3.7) then becomes V (1)−V (2) = T (2)−T (1), or T (1)+V (1) = T (2)+ V (2). This tells us that the quantity T + V stays constant as the particle moves from point 1 to point 2. We therefore have obtained the statement of conservation of total mechanical energy E = T + V = 1 2 mv2 + V (r) (1.3.9) for a particle moving under the action of a force F that derives from a potential V . We can verify directly from Eq. (1.3.9) that the total energy is a constant of the motion. We have dE dt = 1 2 m dv2 dt + dV dt . As we have seen, dv2 dt = 2 dv dt · v. The potential energy V depends on time only through the changing position of the particle: V = V (r(t)) = V (x(t), y(t), z(t)). We therefore have dV dt = ∂V ∂x dx dt + ∂V ∂y dy dt + ∂V ∂z dz dt = ∇V · v. All of this gives dE dt = ma · v + ∇V · v = F · v − F · v = 0, as expected. An example of a force that derives from a potential is gravity: The force Fgravity = mg = mg(0, 0,−1) (1.3.10) is the negative gradient of Vgravity = mgz. (1.3.11) We have indicated that the vector g points in the negative z direction (down, that is); its magnitude is the gravitational acceleration g ≃ 9.8 m/s2. The total 1.3 Mechanics of a single body 11 mechanical energy E is conserved when a particle moves under the action of the gravitational force. An example of a force that does not derive from a potential is the frictional force Ffriction = −kv, (1.3.12) where k > 0 is the coefficient of friction; this force acts in the direction opposite to the particle’s motion and exerts a drag. It is indeed easy to see that Ffriction cannot be expressed as the gradient of a function of r. (The expression Vfriction = kv · r might seem to work, but this potential depends on both r and v, and this is not allowed.) This implies that in the presence of a frictional force, the total mechanical energy of a particle is not conserved. The reason is that the friction produces heat, which is rapidly dissipated away; because this heat comes at the expense of the particle’s mechanical energy, E cannot be conserved. Energy conservation as a whole, of course, applies: the amount by which E decreases matches the amount of heat dissipated into the environment. It is important to understand that the work-energy theorem of Eq. (1.3.7) is always true, whether or not the force F derives from a potential. But whether E is conserved or not depends on this last property: When F = −∇V we have dE/dt = 0 and the total mechanical energy is conserved; but E is not in general conserved when the force does not derive from a potential. 1.3.5 Case study #1: Particle in a gravitational field To illustrate the formalism presented in the preceding subsections we now review the problem of determining the motion of a particle in a gravitational field. The force is given by Eq. (1.3.10), F = mg = mg(0, 0,−1), and the potential by Eq. (1.3.11), V = mgz. The equations of motion are ẍ = 0, ÿ = 0, z̈ = −g. (1.3.13) These are easily integrated: x(t) = x(0) + vx(0)t, y(t) = y(0) + vy(0)t, z(t) = z(0) + vz(0)t − 1 2 gt2. (1.3.14) These equations describe parabolic motion. Here x(0), y(0), z(0) are the positions at time t = 0, and vx(0), vy(0), and vz(0) are the components of the velocity vector at t = 0; these quantities are the initial conditions that must be specified in order for the motion to be uniquely known at all times. The velocity vector at time t is obtained by differentiating Eqs. (1.3.14); we get vx(t) = vx(0), vy(t) = vy(0), vz(t) = vz(0) − gt. (1.3.15) With Eqs. (1.3.14) and (1.3.15) we have sufficient information to compute the total mechanical energy E = T + V of the particle. After some simple algebra we obtain E = 1 2 m [ vx(0) 2 + vy(0) 2 + vz(0) 2 ] + mgz(0) (1.3.16) for all times t; this is clearly a constant of the motion. Exercise 1.6. Verify that Eqs. (1.3.14) really give the solution to the equations of motion r̈ = g. Then compute E and make sure that your result agrees with Eq. (1.3.16). 12 Newtonian mechanics 1.3.6 Case study #2: Particle in a gravitational field subjected to air resistance We now suppose that the particle is subjected to both a gravitational force mg and a frictional force −kv supplied by the ambient air. For convenience we set k = m/τ , thereby defining the quantity τ , and the total applied force is F = m(g − v/τ). (1.3.17) The equations of motion are ma = F , or a = g − v/τ , or again r̈ + ṙ/τ = g. (1.3.18) We assume that the particle is released from a height h with a zero initial velocity. The initial conditions are therefore z(0) = h and ż(0) = 0. We assume also, for simplicity, that there is no motion in the x and y directions. The only relevant component of Eq. (1.3.18) is therefore v̇ + v/τ = −g, (1.3.19) where we have set v = ż. To arrive at Eq. (1.3.19) we have used the fact that g = g(0, 0,−1). Our task is to solve the first-order differential equation of Eq. (1.3.19). We use the method of variation of parameters. Suppose first that g = 0. In this case the equation becomes dv/dt = −v/τ or dv/v = −dt/τ . This is easily integrated, and we get ln(v/c) = −t/τ , or v = c e−t/τ . This is the solution for g = 0, and the constant of integration c is the solution’s parameter. To handle the case g 6= 0 we allow c to depend on time — we vary the parameter — and we substitute the trial solution v(t) = c(t)e−t/τ into Eq. (1.3.19). We have v̇ = ċe−t/τ − v/τ and −g = v̇ + v/τ = ċe−t/τ . The differential equation for c(t) is therefore ċ = −get/τ , so that c(t) = −gτet/τ + c0, where c0 is a true constant of integration. The result for v(t) is then v(t) = −gτ + c0e−t/τ . To determine c0 we invoke the initial condition v(0) = 0. Because v(0) = −gτ + c0 we have that c0 = gτ . Our final answer is therefore v(t) = −gτ [ 1 − e−t/τ ] . (1.3.20) This is ż, the z component of the particle’s velocity vector. Integrating Eq. (1.3.20) gives z(t), the position of the particle as a function of time. Exercise 1.7. Integrate Eq. (1.3.20) and obtain z(t). Make sure to impose the initial condition z(0) = h. Equation (1.3.20) simplifies when t is much smaller than τ = m/k. At such early times, when t/τ ≪ 1, the exponential is well approximated by e−t/τ ≃ 1− t/τ and Eq. (1.3.20) becomes v(t) ≃ −gt, 1.3 Mechanics of a single body 15 the net torque acting on it. Work through the details and verify that this equation does indeed lead to Eq. (1.3.24). This method of derivation does not require the new basis of unit vectors; all calculations can be carried out in the Cartesian basis. The second-order differential equation of Eq. (1.3.24) determines the motion of the pendulum. It can immediately be integrated once with respect to time. The trick is to multiply Eq. (1.3.24) by θ̇; this gives θ̈θ̇ + (ω2 sin θ)θ̇ = 0. Now note that θ̈θ̇ = 1 2 d dt θ̇2 and (sin θ)θ̇ = − d dt cos θ. We therefore have d dt ( 1 2 θ̇2 − ω2 cos θ ) = 0, or 1 2 θ̇2 − ω2 cos θ = ε = constant. (1.3.26) This is a first-order differential equation for θ(t). It seems intuitively plausible that the conserved quantity ε should have some- thing to do with the pendulum’s total energy E. This is indeed the case. The kinetic energy is T = 12m(ẋ 2 + ż2) = 12mℓ 2θ̇2, according to our previous results. The po- tential energy associated with the gravitational force is V = −mgz = −mgℓ cos θ = −mℓ2ω2 cos θ, where we have used Eq. (1.3.25). The potential energy associated with the rod’s tension is zero: The tension always acts in the rod’s direction, which is always perpendicular to the direction of the motion; the tension does no work on the pendulum. We finally have E = T + V = mℓ2( 12 θ̇ 2 − ω2 cos θ), or E = mℓ2ε. (1.3.27) We shall call ε the pendulum’s reduced energy. Similarly, we shall call 12 θ̇ 2 the reduced kinetic energy and ν(θ) ≡ −ω2 cos θ the reduced potential energy. The qualitative features of the pendulum’s motion can be understood without further calculation, purely on the basis of the following graphical construction. We draw an energy diagram, a plot of the reduced potential energy ν(θ) = −ω2 cos θ as a function of θ, together with the constant value of the reduced energy ε (see Fig. 1.6). According to Eq. (1.3.26), which we rewrite as 1 2 θ̇2 = ε − ν(θ), ν(θ) = −ω2 cos θ, (1.3.28) the difference between ε and ν(θ) is equal to the reduced kinetic energy 12 θ̇ 2. For motion to take place this difference must be positive, and a quick examination of the diagram reveals immediately the regions for which ε − ν(θ) ≤ 0. Motion is possible within these regions, and impossible outside. For example, when ε < ω2 we see that the motion of the pendulum takes place between the two well-defined limits θ = ±θ0; motion is impossible beyond these points. This situation corresponds to ordinary pendulum motion: The weight oscil- lates back and forth around the horizontal axis (θ = 0), with an amplitude θ0. The diagram reveals that the angular velocity |θ̇| is maximum when the weight crosses 16 Newtonian mechanics -1 -0.5 0 0.5 1 1.5 -4 -2 0 2 4 θ Figure 1.6: Energy diagram for the pendulum. The difference between the line ε = constant and the curve ν(θ) = −ω2 cos θ is the reduced kinetic energy 1 2 θ̇2, which must be positive for motion to take place. The lower value of ε is such that ε < ω2. The higher value is such that ε > ω2. In the plot ω2 is set equal to 1. θ = 0, and that the pendulum comes to a momentary rest (θ̇ = 0) when θ = ±θ0. This amplitude is determined by setting θ̇ = 0 in Eq. (1.3.28); we have ε = ν(θ0) = −ω2 cos θ0. (1.3.29) This equation can be solved for θ0 whenever ε < ω 2; there are no solutions otherwise. When ε > ω2 the diagram reveals that there are no intersections between the line ε = constant and the curve ν(θ). There are no points at which 12 θ̇ 2 = 0, θ is allowed to increase without bound, and the motion is not limited. This high-energy situation corresponds to the weight doing complete revolutions around the pivot point. Points in the energy diagram at which the line ε = constant meets the curve ν(θ) are called turning points. At these points the reduced kinetic energy 12 θ̇ 2 drops to zero and θ̇ changes sign, either from the positive to the negative (if θ was increasing toward θ0), or from the negative to the positive (if θ was decreasing toward −θ0). These are the points at which the pendulum reaches its maximum angle and turns around. Combining Eqs. (1.3.28) and (1.3.29) gives 1 2 θ̇2 = ω2(cos θ − cos θ0), (1.3.30) and this is a first-order differential equation for θ(t). This equation, unfortunately, cannot be solved in closed form, unless θ0 is assumed to be very small (we shall deal separately with this simple case at the end of this subsection). The best we can do is to express t in terms of an integral involving θ. First we take the square root of Eq. (1.3.30), θ̇ = ± √ 2ω √ cos θ − cos θ0, 1.3 Mechanics of a single body 17 -3 -2 -1 0 1 2 3 0 0.5 1 1.5 2 2.5 3 sw in g an gl e time Figure 1.7: Motion of a pendulum with four different amplitudes. We plot the swing angle θ (in radians) as a function of time. The unit of time is 2π/ω. Notice that the period of oscillation increases with the amplitude. Notice also that for high amplitude, the curve differs significantly from a sinusoid. and we solve for dt. After integration we get t = ± 1√ 2ω ∫ dθ√ cos θ − cos θ0 + constant. (1.3.31) This integral must be evaluated numerically, and the result t(θ) must be inverted to give θ(t); the inversion must also be done numerically. To obtain these details requires some labour, and this will not be pursued here. The results of a numerical integration are presented in Fig. 1.7. The motion of the pendulum is clearly periodic, and Eq. (1.3.31) allows us to calculate the period P , the time required for the pendulum to complete a full cycle of oscillation (θ going from −θ0 to +θ0 and then back to −θ0.) This is twice the time required to go from −θ0 to +θ0, or four times the time required to go from θ = 0 to θ = θ0. So the period is given by P = 4√ 2ω ∫ θ0 0 dθ√ cos θ − cos θ0 . To put this integral in standard form we change the variable of integration to z = sin 12θ sin 12θ0 and introduce the parameter s = sin 12θ0. (1.3.32) Simple manipulations reveal that dz dθ = √ 1 − s2z2 2s , √ cos θ − cos θ0 = √ 2s √ 1 − z2, 20 Newtonian mechanics 1.4.2 Centre of mass The centre of mass of a system of N bodies is at a position R which is defined by R = 1 M ∑ A mArA, (1.4.4) where M = ∑ A mA (1.4.5) is the total mass of the system. The centre of mass moves in accordance with Newton’s law, which implies MR̈ = ∑ A mAaA = ∑ A FA = ∑ A,B,A 6=B FAB , where we have used Eq. (1.4.2). In the last line we sum over both A and B (both from 1 to N), but we make sure to exclude all terms for which A = B. Let us examine the double sum in the special case of three particles. We have ∑ A,B,A 6=B FAB = N ∑ A=1 N ∑ B=1 FAB = N ∑ A=1 ( FA1 + FA2 + FA3 ) = ( F21 + F31 ) + ( F12 + F32 ) + ( F13 + F23 ) = ( F21 + F12 ) + ( F31 + F13 ) + ( F32 + F23 ) = 0. The double sum vanishes by virtue of Newton’s third law, and this property remains true for arbitrary values of N . We therefore have R̈ = 0, ⇒ R(t) = R(0) + Ṙ(0)t. (1.4.6) The centre of mass moves with a uniform velocity, and it therefore defines the origin of another inertial frame. It is usually convenient to shift the origin of the reference frame to the centre of mass, by defining new positions vectors r′A(t) according to r′A = rA − R. (1.4.7) It should be kept in mind that the centre of mass defines the origin of an inertial frame only when there are no external forces acting on the particles. When external forces are present each particle moves according to mAaA = F internal A + F external A , where the first term represents the internally-produced force acting on A, and the second term represents the external force. It is then easy to show that the centre of mass will move according to MR̈ = ∑ A F external A ; it is accelerated by the net sum of all the external forces. Exercise 1.12. Prove the preceding statement. 1.4 Mechanics of a system of bodies 21 1.4.3 Total linear and angular momenta The total linear momentum of the system of N bodies is defined by P = ∑ A pA = ∑ A mAvA, (1.4.8) where pA are the individual momenta. We have P = d dt ∑ A mArA, or, according to Eq. (1.4.4), P = MṘ. (1.4.9) The total momentum therefore follows the motion of the centre of mass. Because Ṙ(t) = Ṙ(0) according to Eq. (1.4.6), we have the important statement that the total linear momentum is a constant vector. If the origin of the inertial frame is at the centre of mass, then R = 0 and Ṙ = 0; this means that P = 0. In this centre-of-mass frame, the total momentum of the system of particles is zero. The total angular momentum of the system is L = ∑ A rA × pA = ∑ A mArA × vA. (1.4.10) Its rate of change is calculated as L̇ = ∑ A mA ( vA × vA + rA × aA ) = ∑ A rA × FA = ∑ A,B,A 6=B rA × FAB , where we have again involved Eq. (1.4.2). Let us examine the double sum for the special case of three particles. We have ∑ A,B,A 6=B rA × FAB = ∑ A ( rA × FA1 + rA × FA2 + rA × FA3 ) = ( r2 × F21 + r3 × F31 ) + ( r1 × F12 + r3 × F32 ) + ( r1 × F13 + r2 × F23 ) = ( r1 − r2 ) × F12 + ( r1 − r3 ) × F13 + ( r2 − r3 ) × F23, where we have used Eq. (1.4.3). The vector r1−r2 is directed from body 2 to body 1. In most circumstances the force F12 also is directed from body 2 to body 1 (or in the opposite direction). Under these conditions the vector product (r1−r2)×F12 is zero, and this is true for all other pairs of bodies. The double sum is therefore zero. These considerations generalize to an arbitrary number of bodies, and we conclude that L̇ = 0 (1.4.11) whenever the force FAB points in the direction of the relative separation rA − rB . Under these conditions we have conservation of the system’s total angular 22 Newtonian mechanics momentum. Exercise 1.13. Calculate dP /dt and dL/dt when there are also external forces acting on the particles. Let us express the position vector of body A as in Eq. (1.4.7), rA = R + r ′ A, (1.4.12) where r′A is its position relative to the centre of mass. We write, similarly, vA = Ṙ + v ′ A. (1.4.13) We make these substitutions into Eq. (1.4.10), and get L = ∑ A mA ( R + r′A ) × ( Ṙ + v′A ) = ∑ A mA ( R × Ṙ + R × v′A + r′A × Ṙ + r′A × v′A ) = ( R × Ṙ ) ∑ A mA + R × ∑ A mAv ′ A − Ṙ × ∑ A mAr ′ A + ∑ A mAr ′ A × v′A. This mess simplifies. For the first term on the right-hand side we have ∑ A mA = M , the total mass of the system. In the second term we recognize that ∑ A mAv ′ A is the system’s total momentum as measured in the centre-of-mass frame; this is zero. The third term vanishes also, and we finally have L = MR × Ṙ + ∑ A mAr ′ A × v′A. (1.4.14) In this expression, the first term represents the angular momentum of the centre of mass, while the second term is the total angular momentum of the system of particles relative to the centre of mass. When the origin of the inertial frame is placed at the centre of mass, we have R = 0 and the first term disappears. In general, we see that L depends on the choice of origin. 1.4.4 Conservation of energy The presentation here parallels closely our discussion of Sec. 1.3.4 on energy con- servation for a single particle. The notation of this section, however, will be slightly more cumbersome, because we now have to keep track of many particles. We begin by calculating the total work done on all the particles as they move from a configuration labeled 1 to another configuration labeled 2. (This means that in the interval of time over which we follow the particles, each moves from a point 1 to a point 2 on its trajectory.) This is W12 = ∑ A ∫ 2 1 FA · drA = ∑ A ∫ 2 1 FA · vA dt, (1.4.15) where drA = vA dt is the displacement vector on the trajectory of particle A. Substituting the equations of motion (1.4.1) gives W12 = ∑ A ∫ 2 1 mA dvA dt · vA dt. 1.4 Mechanics of a system of bodies 25 But dr1 − dr2 = d(r1 − r2) = dr12, so this can be expressed as W12 = ∫ 2 1 ( F12 · dr12 + F13 · dr13 + F23 · dr23 ) . At this stage of the derivation we incorporate the fact that the mutual forces are derived from a potential. As we have seen, F12 = −∇1V12, where ∇1 = (∂/∂x1, ∂/∂y1, ∂/∂z1). But since V12 depends on (x1, y1, z1) only through its de- pendence on (x12, y12, z12) (where, for example, x12 = x1 − x2), the force can also be expressed as F12 = −∇12V12, where ∇12 is the gradient operator with respect to r12 = (x12, y12, z12), ∇12 = ( ∂ ∂x12 , ∂ ∂y12 , ∂ ∂z12 ) . This is possible because ∂x12/∂x1 = 1, and so on. So we now have W12 = ∫ 2 1 ( −∇12V12 · dr12 − ∇13V13 · dr13 − ∇23V23 · dr23 ) . Each integral can be evaluated (refer back to Sec. 1.3.1), giving W12 = − [ V12(2) − V12(1) ] − [ V13(2) − V13(1) ] − [ V23(2) − V23(1) ] ≡ − [ V12 + V13 + V23 ]2 1 . Since V21 = V12 and so on, we may write this as W12 = − 1 2 [ V12 + V13 + V21 + V23 + V31 + V32 ]2 1 , where we now sum over all possible pairs of indices, provided that each index is not repeated. Generalizing to an arbitrary number of particles, this is W12 = − [ 1 2 ∑ A,B,A 6=B VAB ]2 1 . We define the total potential energy of the system to be V = 1 2 ∑ A,B,A 6=B VAB . (1.4.23) With this notation our previous result is W12 = −[V (2) − V (1)], and Eq. (1.4.17) becomes −V (2) + V (1) = T (2) − T (1) or T (1) + V (1) = T (2) + V (2). We have finally established that the total mechanical energy of the system, E = T + V = ∑ A 1 2 mAv 2 A + 1 2 ∑ A,B,A 6=B VAB(rAB), (1.4.24) stays unchanged as the particles move from configuration 1 to configuration 2. We recall that the mutual potentials VAB are assumed to depend on rAB = |rA − rB | only; the mutual forces are then given by Eqs. (1.4.21) and (1.4.22). This is the statement of energy conservation for a system of particles. Exercise 1.15. Starting from the definition of Eq. (1.4.24), prove directly that dE/dt = 0. 26 Newtonian mechanics 1.5 Kepler’s problem To give concreteness to the formal developments of the preceding section we exam- ine, in this section, the specific situation of two bodies subjected to their mutual gravitational forces. This could be the Earth-Moon system, or the Sun-Jupiter sys- tem, or again a binary system of two main-sequence stars. Our goal is to determine the motion of the two bodies, that is, to find a solution to Kepler’s problem. 1.5.1 Gravitational force The force acting on body 1 due to the gravity of body 2 has a magnitude Gm1m2/r 2 12, where G is Newton’s gravitational constant, m1 the mass of body 1, m2 the mass of body 2, and r12 is the distance between the two bodies. The force is directed along the vector r2 − r1, which points from body 1 to body 2. Introducing the notation r = r1 − r2, r = |r1 − r2| ≡ r12, (1.5.1) we write F12 = −Gm1m2 r r3 . (1.5.2) The force acting on body 2 due to the gravity of body 1 is F21 = Gm1m2 r r3 , (1.5.3) and it is directed along r1 − r2, which points from body 2 to body 1. These forces can be derived from a mutual potential V12 = − Gm1m2 r . (1.5.4) This means that the force of Eq. (1.5.2) is given by F12 = −∇1V12, (1.5.5) where ∇1 is the gradient operator with respect to the coordinates r1 = (x1, y1, z1) of body 1. Similarly, the force of Eq. (1.5.3) can be expressed as F21 = −∇2V12, (1.5.6) where ∇2 is the gradient operator with respect to the coordinates r2 = (x2, y2, z2) of body 2. To verify these statements, let us calculate, say, the z component of F21. We have F21,z = − ∂V12 ∂z2 = −dV12 dr ∂r ∂z2 . The first factor is dV12 dr = Gm1m2 r2 , and to calculate the second factor we start with r2 = (x1 − x2)2 + (y1 − y2)2 + (z1 − z2)2 and differentiate both sides with respect to z2. This gives 2r ∂r ∂z2 = −2(z1 − z2) or ∂r ∂z2 = −z1 − z2 r . 1.5 Kepler’s problem 27 So finally, F21,z = Gm1m2 z1 − z2 r3 , and this is clearly compatible with Eq. (1.5.3). Similar calculations would return all other components of F21 and all components of F12, and Eqs. (1.5.5) and (1.5.6) would be fully verified. According to Eq. (1.4.23), the total potential energy of the two-body system is V = 1 2 ∑ A,B,A 6=B VAB = 1 2 (V12 + V21), or V = V12. (1.5.7) This result will allow us, in the following subsections, to omit the label “12” from the mutual potential; we shall write, simply, V12 = V = −Gm1m2/r. 1.5.2 Equations of motion Newton’s equations for the two bodies are m1r̈1 = F12 = −Gm1m2r/r3 and m2r̈2 = F21 = Gm1m2r/r 3. Simplifying, we arrive at r̈1 = −Gm2 r r3 (1.5.8) and r̈2 = Gm1 r r3 , (1.5.9) where, we recall, r = r1 − r2 and r = |r|. The position vectors r1 and r2 can be expressed in terms of R, the position of the centre of mass, and r, the relative position. We have, according to Eq. (1.4.4), MR = m1r1 + m2r2, where M = m1 + m2 is the total mass. Simple algebra gives r1 = R + m2 M r (1.5.10) and r2 = R − m1 M r. (1.5.11) The motion of the centre of mass is determined by the equation MR̈ = m1r̈1 + m2r̈2 = −Gm1m2r/r3 + Gm1m2r/r3 = 0. As we had discovered in Sec. 1.4.2, the centre of mass moves uniformly: R(t) = R(0) + Ṙ(0)t. (1.5.12) The motion of the relative position, on the other hand, is determined by the equation r̈ = r̈1 − r̈2 = −Gm2r/r3 − Gm1r/r3, or r̈ = −GM r r3 , M = m1 + m2. (1.5.13) Exercise 1.16. Verify Eqs. (1.5.10) and (1.5.11). The centre of mass defines the origin of an inertial frame, and the mathematical description of the two-body system is simplest in this reference frame. We shall therefore set R = 0, which brings Eqs. (1.5.10) and (1.5.11) to the simpler form r1 = m2 M r, r2 = − m1 M r. (1.5.14) 30 Newtonian mechanics φr d rφd Figure 1.10: The position vector moves by an angle dφ during a time dt. 1.5.5 Kepler’s second law Kepler’s second law states that the position vector of a planet orbiting the Sun sweeps out equal areas in equal times. The statement generalizes to any two massive bodies, and in this case the position vector refers specifically to the relative separation r(t) between the two bodies. This law comes as an immediate consequence of Eq. (1.5.26), which implies the conservation statement r2φ̇ = constant. Let us first show that this constant is in fact h, the magnitude of the constant vector h defined by Eq. (1.5.15). According to Eq. (1.5.20), this is given by h = xẏ − yẋ. Making use of Eqs. (1.5.21), we write this as h = (r cos φ)(ṙ sin φ + rφ̇ cos φ) − (r sin φ)(ṙ cos φ − rφ̇ sin φ). Simplifying this we obtain h = r2φ̇ = constant, (1.5.27) the expected result. The fact that r2φ̇ is conserved gives us our statement of the second law. Consider Fig. 1.10. During an interval dt of time the position vector moves by an angle dφ and sweeps out an area dA. To a good approximation the area is shaped as a triangle and we have dA = 12r(r dφ) = 1 2r 2 dφ. The rate at which the position vector sweeps out this area is therefore dA dt = 1 2 r2φ̇ = 1 2 h. (1.5.28) This is a constant, and we have the mathematical statement of Kepler’s second law. 1.5.6 Conservation of energy With the substitution φ̇ = h/r2 obtained from Eq. (1.5.27), Eq. (1.5.25) becomes r̈ + GM r2 − h 2 r3 = 0. (1.5.29) This equation can immediately be integrated by multiplying all members by ṙ. (Recall that we used the same trick back in Sec. 1.3.7.) We have r̈ṙ = d dt ( 1 2 ṙ2 ) , 1.5 Kepler’s problem 31 GM r2 ṙ = d dt ( −GM r ) , −h 2 r3 ṙ = d dt ( h2 2r2 ) , and the preceding equation becomes d dt ( 1 2 ṙ2 − GM r + h2 2r2 ) = 0. This implies 12 ṙ 2 − GM/r + h2/(2r2) = ε, where ε is the constant of integration. We shall write this result in the form 1 2 ṙ2 + ν(r) = ε, (1.5.30) with ν(r) = −GM r + h2 2r2 . (1.5.31) The first term on the left of Eq. (1.5.30) can be thought of as a reduced kinetic energy for the radial component of the motion. The second term is a reduced effective potential for this motion, and the constant ε is a reduced total energy. (Recall that we introduced this terminology back in Sec. 1.3.7; by “reduced” we mean a rescaled version of the usual quantities.) The reduced energy ε is directly related to the system’s true total energy E. Let us calculate it. The system’s total kinetic energy is T = 12m1|ṙ1|2 + 12m2|ṙ2|2, and according to Eqs. (1.5.4) and (1.5.7), the potential energy is V = −Gm1m2/r. The total energy is therefore E = 1 2 m1|ṙ1|2 + 1 2 m2|ṙ2|2 − Gm1m2 r . After involving Eq. (1.5.14) and cleaning up the algebra, this becomes E = m1m2 M ( 1 2 |ṙ|2 − GM r ) . (1.5.32) Now Eq. (1.5.23) states ṙ = v = ṙr̂ + rφ̇φ̂, so that |ṙ|2 = ṙ2 + r2φ̇2 = ṙ2 + h2/r2. So E = m1m2 M ( 1 2 ṙ2 − GM r + h2 2r2 ) , and comparing this with Eqs. (1.5.30) and (1.5.31) yields E = m1m2 M ε. (1.5.33) As promised, ε is a rescaled version of the total energy E. Recall that back in Sec. 1.5.3 we had similarly obtained L = (m1m2/M)h. Exercise 1.18. Verify Eq. (1.5.32) and check the algebra leading to Eq. (1.5.33). 1.5.7 Qualitative description of the orbital motion The equations of motion have been reduced to the first-order form of Eqs. (1.5.27) and (1.5.30), with the effective potential ν(r) given by Eq. (1.5.31). The equation φ̇ = h/r2 informs us that φ increases monotonically with t: If h is positive φ̇ is 32 Newtonian mechanics -1 -0.5 0 0.5 1 0 2 4 6 8 10 12 r Figure 1.11: Energy diagram for the radial component of the motion. When ε < 0 the motion takes place between two turning points; this is elliptical motion. When ε > 0 the motion takes place on the right of a single turning point; this is hyperbolic motion. When ε = 0 the motion is parabolic. always greater than zero and φ(t) is an increasing function; if h is negative φ̇ is always smaller than zero and φ(t) is then a decreasing function. (The case h = 0 will be dealt with separately later.) The equation 1 2 ṙ2 + ν(r) = ε governs the radial component of the motion. As in Sec. 1.3.7 we will describe this qualitatively by constructing an energy diagram, a plot of the effective potential ν(r) together with the constant value of the reduced energy ε. The energy diagram is shown in Fig. 1.11. Recall the two main features of such diagrams: (i) The difference between ε and ν(r) represents 12 ṙ 2, the reduced kinetic energy, which must be positive; and (ii) points on the diagram for which ν(r) = ε represent turning points of the motion, at which ṙ changes sign, either from positive to negative, or from negative to positive. Because the effective potential can be negative, it is possible for ε to be either negative or positive. The nature of the motion depends sensitively on this sign. When ε < 0 the motion takes place between two turning points at r = rmin and r = rmax. The motion is bounded, and as we shall see, the orbit possesses an elliptical shape. When ε is equal to the minimum value of the effective potential, ε = νmin < 0, motion can take place only at r = r0, the radius at which the minimum occurs, which is defined by ν(r0) = νmin. The orbit is then circular, because ṙ is always zero and r can therefore never change with time. When ε > 0 the motion takes place only to the right of a single turning point at r = rmin. The motion proceeds from r = ∞ (where ν = 0 and 12 ṙ2 = ε) down to r = rmin (where ṙ changes sign from negative to positive), and then back to r = ∞. The particle traces a hyperbola in the orbital plane, and the motion is said to be hyperbolic. 1.5 Kepler’s problem 35 The final result for r(φ) is r(φ) = p 1 + e cos φ , (1.5.40) in which we have set φ0 = 0 to simplify the expression. This involves two constants: We have p, which plays the role of average radius and is known officially as the semilatus rectum, and we have e, which measures the range over which r varies and is known as the eccentricity. We have seen that p is related to the reduced angular momentum h by h = √ GMp. (1.5.41) The eccentricity, on the other hand, can be related to the reduced energy ε; as we shall calculate in a following paragraph, ε = −GM 2p ( 1 − e2 ) . (1.5.42) This equation is valid for e < 1, which means that ε < 0, and it is valid also for e ≥ 1, which means that ε ≥ 0. We have just observed that ε < 0 when e < 1. This is the case of bound motion, which takes place between two turning points at r = rmin = p/(1 + e) and r = rmax = p/(1 − e). As we see from Eq. (1.5.40), the motion proceeds from r = rmin (known as the orbit’s pericentre) when φ = 0, to r = rmax (known as the orbit’s apocentre) when φ = π, and then back to r = rmin when φ = 2π. When e < 1 the equation r = p/(1 + e cos φ) describes an ellipse. The maximum length of this ellipse is rmin + rmax = p/(1 + e) + p/(1 − e) = 2p/(1 − e2). Half of this is the ellipse’s semi-major axis, a = p 1 − e2 . (1.5.43) These statements give rise to Kepler’s first law: A body moving under the grav- itational influence of another body follows an elliptical orbit when the motion is bounded. When e = 0 we have that Eq. (1.5.40) reduces to r(φ) = p, and the el- lipse has become a circle. In this case we have p ≡ r0, Eq. (1.5.41) becomes identical to Eq. (1.5.34), and Eq. (1.5.42) becomes Eq. (1.5.35). Exercise 1.19. Look up a reference book on elementary geometry and review the properties of ellipses. Answer the following questions: (1) Is Eq. (1.5.40) really the equation of an ellipse? (2) Where are the two foci of the ellipse? (3) What is the semi-minor axis b of the ellipse? A good way to answer some of these questions is to show that Eq. (1.5.40) is equivalent to the usual description of an ellipse via the equation X2/a2 + Y 2/b2 = 1. But be careful! In this equation X is not equal to r cos φ and Y is perhaps not equal to r sin φ, because the two coordinate systems do not share the same origin. We have also observed that ε > 0 when e > 1. This is the case of unbound motion, which takes place to the right of a single turning point at r = rmin = p/(1+e). In this case the equation r = p/(1+e cos φ) describes a hyperbola, and we have hyperbolic motion. In the special case e = 1 we have ε = 0, and the equation r = p/(1+cos φ) describes a parabola; in this special case we have parabolic motion. You may have learned that an ellipse, a hyperbola, and a parabola are all special cases of a general family of curves called conic sections. The conic sections are all described by the parametric equation r(φ) = p/(1 + e cos φ). Let us now return to the derivation of Eq. (1.5.42). We go back to Eqs. (1.5.30) and (1.5.31), ε = 1 2 ṙ2 − GM r + h2 2r2 , 36 Newtonian mechanics and we compute each member of the right-hand side. To evaluate ṙ we start with Eq. (1.5.40) and get ṙ = ep sin φ (1 + e cos φ)2 φ̇ = e p r2φ̇ sin φ = e p h sin φ. This gives us ε = 1 2 e2 p2 h2 sin2 φ − GM p (1 + e cos φ) + h2 2p2 (1 + e cos φ)2, and replacing h2 by GMp in this equation yields ε = GM 2p [ e2 sin2 φ − 2(1 + e cos φ) + (1 + e cos φ)2 ] . After simplification the expression within the square brackets becomes e2 sin2 φ − 1 + e2 cos2 φ = −1 + e2, and we arrive at ε = −GM 2p (1 − e2), the same statement as in Eq. (1.5.42). 1.5.10 Motion in time Now that r(φ) is known we must relate φ to the time t in order to have a complete description of the motion. The relevant equations are φ̇ = h/r2, h = √ GMp, and r = p/(1 + e cos φ). Putting this all together, we obtain dφ dt = √ GM p3 (1 + e cos φ)2. (1.5.44) This is the differential equation that must be solved in order to obtain φ(t). Unless e is very small, in which case approximate analytical results can be obtained, this equation must be integrated numerically. Results of numerical integrations are displayed in Fig. 1.12. The integral form of Eq. (1.5.44) is t = √ p3 GM ∫ dφ (1 + e cos φ)2 + constant. (1.5.45) This indefinite integral cannot be evaluated in closed form, but it provides a nice way of calculating the orbital period P of bound orbits (e < 1). Because this is equal to the time required for φ to advance by 2π, or twice the time required for φ to advance by π, we have P = 2 √ p3 GM ∫ π 0 dφ (1 + e cos φ)2 . This definite integral can be evaluated, and the result is π/(1−e2)3/2. We therefore have P = 2π √ [p/(1 − e2)]3 GM . We obtain a cleaner form of this result by involving Eq. (1.5.43). In terms of the semi-major axis a = p/(1 − e2), the orbital period is P = 2π √ a3 GM . (1.5.46) We have that P 2 ∝ a3, and this is the general statement of Kepler’s third law. 1.5 Kepler’s problem 37 -5 0 5 10 15 20 25 0 0.5 1 1.5 2 2.5 3 ra di us , r ad ia l v el oc ity , a ng ul ar v el oc ity time/(orbital period) radius radial velocity angular velocity Figure 1.12: Numerical integration of the equations of motion for an orbit with eccen- tricity e = 0.5. The blue curve shows the angular velocity φ̇ as a function of time, the green curve shows the radial velocity ṙ as a function of time, and the red curve shows the radial position r as a function of time. The time variable is scaled by the orbital period P , and three complete orbital cycles are displayed. Notice that the motion starts at the pericentre with maximum angular velocity and zero radial velocity. 40 Newtonian mechanics result to y0 to get y1, and repeat the procedure to get y2, y3, and so on. Euler’s method is very simple and economical, but because its error term is of order ∆2, it is not very accurate. With a little cleverness, however, it is possible to improve the accuracy of the method so that its error term becomes of order ∆3 ≪ ∆2. One way of achieving this would be to use Eq. (1.6.3) instead of its cruder version. The price to pay would be the need to evaluate f ′(y), the derivative of the function with respect to y. This may not be practical in some circumstances, and there is an alternative method. Consider evaluating the function f not at y = y0, but at the midpoint between y0 and y1 ≃ y0 + f(y0)∆: f(y0) → f(y0 + 12f0∆), where we use the notation f0 = f(y0). By Taylor expansion we have f(y0 + 1 2f0∆) = f(y0) + f ′(y0)( 1 2f0∆) + O(∆ 2) = f(y0) + 1 2 ff ′(y0)∆ + O(∆ 2), and this shows that Eq. (1.6.3) is equivalent to y1 = y0 + f(y0 + 1 2f0∆)∆ + O(∆ 3). (1.6.4) This approximation for y1 has an error term of order ∆ 3, and it is obtained simply by evaluating the function f at the midpoint; the value of its derivative is not needed. By being increasingly clever it is possible to decrease further the size of the error term. The fourth-order Runge-Kutta method consists of the following recipe. Suppose that the differential equation has been integrated up to x = xn, and that we wish to proceed to the next grid point, at x = xn+1. We have therefore obtained yn and we wish to calculate yn+1. First we compute the auxiliary quantities k1 = f(yn)∆, k2 = f(yn + 1 2k1)∆, k3 = f(yn + 1 2k2)∆, k4 = f(yn + k3)∆, and next we approximate yn+1 by yn+1 = yn + 1 6 k1 + 1 3 k2 + 1 3 k3 + 1 6 k4 + O(∆ 5). (1.6.5) As indicated, the judicious choice of coefficients in front of k1, k2, k3, and k4 ensures that the error term is now of order ∆5, and therefore very small. The Runge-Kutta method is easy to implement, it is accurate, and it is robust; it works well for most functions f(y). The method can also be generalized to handle functions f(x, y) that depend on both variables. The method also generalizes to a set of differential equations dy[i] dx = f [i] ( y[1], y[2], · · · ) , i = 1, 2, · · · (1.6.6) for a set of dependent variables y[i]. In this case the auxiliary quantities k1, k2, k3, and k4 acquire an index [i]; for example we now have k1[i] = f [i] ( yn[1], yn[2], · · · ) ∆ 1.7 Problems 41 and k2[i] = f [i] ( yn[1] + 1 2k1[1], yn[2] + 1 2k1[2], · · ·)∆. This generalization is useful, because it allows us to use the method to integrate second-order differential equations. Consider, for example, the pendulum equation of Eq. (1.3.24), θ̈ + ω2 sin θ = 0. This can be recast as a system of two first-order differential equations. To do this we define y[1] = θ and y[2] = θ̇. When then have the system dy[1] dt = y[2], dy[2] dt = −ω2 sin(y[1]). In this instance we find that f [1](y[1], y[2]) = y[2] and f [2](y[1], y[2]) = −ω2 sin(y[1]). This system of equations can be integrated straightforwardly, and the result for y[1](t) is the numerical approximation to θ(t), the solution to the original second- order equation. 1.7 Problems 1. Let A = (3x2 − 6yz)x̂ + (2y + 3xz)ŷ + (1 − 4xyz2)ẑ. Calculate ∫ C A · ds along the following paths that link the point (0, 0, 0) to the point (1, 1, 1): (a) The curve described by x = u, y = u2, and z = u3, in which the parameter u is restricted to the interval 0 < u < 1. (b) The straight line that joins these points. 2. Let A = (2xy + z3)x̂ + (x2 + 2y)ŷ + (3xz2 − 2)ẑ. Find the function f such that A = ∇f . Then evaluate ∫ C A · ds along any path C that links the point (1,−1, 1) to the point (2, 1, 2). 3. Evaluate ∮ C r · ds along all closed loops C, where r = xx̂ + yŷ + zẑ is the position vector. 4. A projectile is launched with initial speed v0 at an angle α with the horizontal. Calculate: (a) the position vector as a function of time; (b) the time required to reach the highest point; (c) the maximum height reached by the projectile; (d) the time of flight back to the Earth’s surface; (e) the range of the projectile; (f) the angle α which maximizes the range. 5. Suppose that in the preceding problem, the projectile is also subjected to a frictional force equal to −kv, where v is the velocity vector and k a positive constant. Find: (a) the velocity vector as a function of time; 42 Newtonian mechanics (b) the position vector as a function of time; (c) the terminal velocity of the projectile. For ease of notation set k = m/τ . 6. A particle of mass m is traveling in the x direction. At time t = 0 it is located at x = 0 and has a speed v0. The particle is subjected to a frictional force which opposes the motion; its magnitude is equal to βv2, where v = v(t) is the particle’s speed at time t and β is a positive constant. (a) What is the speed of the particle as a function of time? (b) What is the position of the particle as a function of time? 7. A particle of mass m rests on top of a sphere of radius R. The particle is then displaced slightly so that it starts to move down the sphere. (It is assumed that the particle slides down without rolling and without friction.) As it moves down the sphere, the particle makes an angle θ with the vertical direction. At some point the particle loses contact with the surface of the sphere, and it proceeds to fall freely. We are interested in the motion of the particle from the initial moment where it is at rest to the final moment where it leaves the sphere. (a) Derive an equation of motion for θ(t), and find an expression for N , the magnitude of the normal force. (b) At which angle θ does the particle leave the surface of the sphere? (c) What is the speed of the particle when it leaves the surface of the sphere? [Hint: This problem is involved. You may find it useful to resolve the force and acceleration vectors into a basis that consists of r̂, a unit vector that points in the direction normal to the sphere, and θ̂, a unit vector that points in the direction of increasing θ.] 8. The planar pendulum of Sec. 1.2.7 is now subjected to a frictional force Ffriction = −(m/τ)v, where τ is a positive constant. Derive the new equa- tion of motion for the swing angle θ. 9. The equation of motion of the preceding problem reduces to θ̈ + 2γθ̇ + ω2θ = 0 when the oscillations have a very small amplitude; here γ is a positive constant that is related to τ in the preceding problem. Find the general solution to this equation. Assume that ω2 > γ2, so that the oscillations are underdamped. [Be sure that your final expression for θ(t) is a real (not complex) function.] 10. A mass m is allowed to move along the x axis, either in the positive or in the negative direction. It is subjected to a constant force +F when x < 0 and to a constant force −F when x > 0 (here F is positive). (a) Describe the motion qualitatively with the help of an energy diagram. (b) Calculate the period of the motion; express your result in terms of m, F , and the amplitude A of the motion. 11. A cylindrical cork is partially immersed in a liquid of (mass) density ρ. The cork’s axis is oriented with the vertical direction, and the cork floats in the liquid. The cylinder’s cross-sectional area is A, and its mass is m. The cork is gently pushed down into the liquid and then released; it starts to 1.8 Additional problems 45 (b) Determine the integer n. 20. A two-body system moves under the influence of a central force given by f = a r2 + b r3 , where a and b are constants. (a) Show that the shape of the orbit is described by r = p 1 + e cos(kφ) , where p, e, and k are constants. Express p and k in terms of a, b, h2, and µ. (Assume that b < µh2.) (b) Plot the orbit in the x-y plane. Set p = 1, e = 0.6, k = 0.99, and let φ range from 0 to 16π. What is happening to the major axis of the ellipse? 1.8 Additional problems 1. An inclined plane makes an angle α with the horizontal. A projectile is launched from point A at the bottom of the inclined plane. Its initial speed is v0, and its initial velocity vector makes an angle β with the horizontal. The projectile eventually hits the inclined plane at point B. Air resistance is negligible. (a) Calculate the range R of the projectile, the distance between points A and B. Show that it it can be expressed in the form R = R0 sin(β − α) cos β and find an expression for R0. (b) Find the angle βmax which maximizes the range. 2. A particle traveling in the positive x direction is subjected to a force F = kx3. The particle started from an initial position x0 < 0. Draw an energy diagram for this situation and provide a qualitative description of the possible motions. 3. Two bodies of masses m1 and m2 are subjected to a mutual attractive force F12 = −km1m2r, where k is a constant and r = r1−r2 is the relative position vector. (a) Show that the equation of motion for r(t) can be put in the form of an energy equation, 1 2 ṙ2 + ν(r) = ε, and find an expression for ν(r), the effective potential. Draw an energy diagram for this system and give a qualitative description of the possible motions. (b) Prove that r(φ) = r0 √ 2 − e − e cos(2φ) describes the shape of the orbit, and solve for r0 in terms of the constants e, M = m1 + m2, h, and k. 46 Newtonian mechanics 4. The parabolic coordinates u and v are sometimes useful to describe the motion of a particle in a two-dimensional plane. These are related to the Cartesian coordinates x and y by x = uv, y = 1 2 (u2 − v2). (a) Sketch the shapes of the curves u = constant in the x-y plane. (b) Sketch the shapes of the curves v = constant in the x-y plane. (c) Find the unit vectors û and v̂ associated with this coordinate system. Chapter 2 Lagrangian mechanics 2.1 Introduction: From Newton to Lagrange The methods of Newtonian mechanics, based on the vectorial equation F = ma, are very powerful and they can be applied to all mechanical systems. But they lack in efficiency when Cartesian coordinates (x, y, z) do not give the simplest description of a mechanical system. An example is the problem of the pendulum (Sec. 1.3.7), which is best analyzed in terms of the swing angle θ; we have seen that to derive the equation of motion for θ(t) requires somewhat laborious calculations, and the reason is precisely that θ is not a Cartesian coordinate. Another example is Kepler’s problem (Sec. 1.5), which is best analyzed in terms of the polar coordinates (r, φ); again we saw (back in Sec. 1.5.4) that to derive equations of motion for r(t) and φ(t) required some long calculations. To increase the efficiency of the theoretical methods of mechanics, a number of scientists in the centuries following Newton endeavoured to recast the Newtonian laws into a more flexible formulation. The most famous players include Leonhard Euler (1707–1783), Joseph Lagrange (1736–1813), William Rowan Hamilton (1805– 1865), and Carl Gustav Jacobi (1804–1851). Their new techniques proved extremely useful, and they allowed them and others to solve increasingly challenging problems, most notably in the context of celestial mechanics. These new powerful techniques are the topic of this chapter on Lagrangian mechanics, and the following chapter on Hamiltonian mechanics. It is important to point out that the Lagrangian and Hamiltonian formulations of the laws of mechanics are largely restricted to forces that can be derived from a potential. For other problems, such as a particle subjected to air resistance, the new techniques cannot be applied in a very straightforward way, and it is usually best to go back to the old Newtonian methods. In this chapter and the next, we shall consider only forces that can be derived from a potential. The entire content of Lagrangian mechanics is summarized in the following sim- ple recipe: 1. Select generalized coordinates qa to describe the degrees of freedom of a me- chanical system. These coordinates are completely arbitrary. They need not be the original Cartesian coordinates associated with an inertial frame. In- deed, there is no need for the coordinates to even be attached to an inertial frame. The index a = 1, 2, · · · labels each one of the generalized coordinates; there is one coordinate for each degree of freedom. 2. In terms of the generalized coordinates, calculate the system’s total kinetic energy T and total potential energy V . Then form what is known as the Lagrangian function of the system, which is denoted L(qa, q̇a); this depends on the generalized coordinates qa and the generalized velocities q̇a = dqa/dt. 47 50 Lagrangian mechanics y x = −1 x = 1 Figure 2.1: A curve in the x-y plane that links the points (−1, 0) and (1, 0). 2.2 Calculus of variations In this section we introduce the mathematical tools — the calculus of variations — that are required in the derivation of the Euler-Lagrange (EL) equations from the principle of least action. We will look at this issue from a purely mathematical point of view, and return to the physics in the next section. 2.2.1 Curve of maximum area Let us examine the following mathematical problem. We consider the infinite num- ber of curves in the x-y plane that link the point (x = −1, y = 0) to the point (x = +1, y = 0); see Fig. 2.1. Of all these curves we select those that have a total arc length (the total distance traveled along the curve) equal to π. Of all the curves that are left we wish to find the one which maximizes the area under the curve. (Notice that the mathematical problem involves maximization of an area, while the physical problem involves minimization of an action. The mathematical techniques to be developed below work for both cases, maximization and minimization, and they do not care about the identity of the quantity to be extremized.) We describe the family of curves introduced in the previous paragraph by para- metric relations x(s) and y(s), in which the parameter s is the curve’s arc length, calculated from the starting point (−1, 0). Because all the curves within the family have a total arc length of π, the parameter s ranges from 0 to π as each curve runs from (−1, 0) to (+1, 0). We have ds2 = dx2 + dy2, and this relation implies that the functions x(s) and y(s) are not independent of each other. The area under the curve is obtained by integration, A = ∫ y dx, which we write as A = ∫ π 0 y(s) dx ds ds. We can replace the factor dx/ds by √ 1 − y′2, where y′ = dy/ds. This gives us, finally, A = ∫ π 0 y √ 1 − y′2 ds. (2.2.1) We wish to find the function y(s) that produces the largest possible value for A. Once this function is identified, x(s) can be obtained by integrating the equation x′ = √ 1 − y′2. (2.2.2) The maximal curve is then fully determined. 2.2 Calculus of variations 51 f(x) xx _ Figure 2.2: A function with a maximum point at x = x̄. Because this is an extremum point, a displacement around x̄ produces the smallest change in the function. 2.2.2 Extremum of a functional To proceed it is helpful to broaden the scope of the preceding discussion and to examine the general structure of the mathematical problem. We are given a func- tional A[y], a function A of a function y(s), which we wish to maximize, or perhaps minimize, with respect to the choice of path y(s). (In general we say that we wish to find the extremum of the functional, and we shall never need to distinguish between a maximum and a minimum.) The functional has the following structure: A[y] = ∫ s1 s0 G(y, y′) ds; (2.2.3) it is given by an integral over a parameter s of a function G which depends on the path y(s) and its derivative y′(s) = dy/ds. The integral can be evaluated for any choice of trial function ytrial(s), and the result is a number Atrial. We are looking for the function ȳ(s) that produces the largest (or smallest) number. In mathematical terms, we are looking for the extremum of the functional A[y]. The mathematical task of extremizing a function f(x) with respect to its argu- ment x — the argument being a number — is a simple one: We simply calculate the derivative of the function and set the result equal to zero; the solutions to df/dx = 0 are all extremum points (minima and maxima) of the function. To ex- tremize a functional A[y] with respect to a functional argument y(s) is a much more delicate task. How does one do this? Let us examine more closely the straightforward task of finding an extremum of a function f(x). We imagine, for concreteness, that the function has a single maximum at x = x̄; this is represented in Fig. 2.2. We have, of course, f ′(x̄) = 0, with a prime indicating differentiation with respect to x. An important property of x̄ is that it is the point from which the function f(x) changes the least when x is displaced from x̄ to a neighbouring point x̄ + δx. That this is so can easily be seen from the figure, but it is just as easy to prove it mathematically. Let us calculate δf , the change induced in the function when its argument x is moved to a neighbouring point x + δx. By Taylor’s theorem we have δf ≡ f(x + δx) − f(x) = f ′(x)δx + 1 2 f ′′(x)(δx)2 + · · · . 52 Lagrangian mechanics y(s) s s y y 0 0 1 1 s Figure 2.3: A family of paths which leave y = y0 when s = s0 and arrive at y = y1 when s = s1. From this calculation we learn that in general, the change in the function is propor- tional to δx, as we might have expected. But when we let x become an extremum point x̄, we get a different result. In this case we have f ′(x̄) = 0 and the preceding equation becomes δf = 1 2 f ′′(x̄)(δx)2 + · · · . Now the change in the function is proportional to (δx)2, and this is much smaller than what we get in the general case. We have just found that the variation δf is smallest when it is taken at an extremum point. A useful way of characterizing an extremum point is therefore to say that it is a point from which a displacement δx produces a vanishing change δf , to linear order in δx. (The change is not actually zero, but it is of second order in δx, as we have shown.) We shall use the same idea to find the extremum path of a functional. We will look for a path ȳ(s) — analogous to the extremum point x̄ — that has the property that a displacement away from this path produces no change in the functional A[y], to linear order in the displacement δy(s). In other words, if we evaluate the function on the extremum path ȳ(s) and get the number Ā, we will find that if we then evaluate the functional on the displaced path y(s) = ȳ(s) + δy(s), we will still get the number Ā, except for a correction of second order in the displacement; the change δA is zero to first order in δy(s). To flesh this out let us consider all paths y(s) that leave the point y = y0 when s = s0 and arrive at the point y = y1 when s = s1; members of this family of curves are displayed in Fig. 1.3. Out of all these possible paths that link y0 and y1 we wish to find the one which extremizes the functional A[y]. Our strategy will be to assume the existence of an extremum path, which we denote ȳ(s), and which we treat as a reference path. We shall examine what happens to A[y] when we displace the path from y(s) = ȳ(s) to y(s) = ȳ(s) + δy(s). While we shall find that in general, this produces a change δA that is proportional to δy(s), we will instead demand that δA vanish to first order in the displacement; as we shall see, this procedure will permit us to identify the extremum path ȳ(s). To carry out this procedure properly it is important to ensure that all the considered paths begin and end at the same two end points. The reference path ȳ(s) and the displaced paths y(s) = ȳ(s) + δy(s) must all satisfy y(s0) = y0 and y(s1) = y1. This implies that the displacement δy(s), which are completely arbitrary in the interval s0 < s < s1, must satisfy the 2.2 Calculus of variations 55 is obtained by finding a solution ȳ(s) to the EL equation d ds ∂G ∂y′ − ∂G ∂y = 0. This statement is true whether the extremum is a maximum or a minimum, and it is independent of the detailed nature of the function G(y, y′). Any function of the two variables y(s) and y′(s) can thus be substituted inside the functional, and our calculus of variations applies to a very wide range of situations. 2.2.3 Curve of maximum area (continued) The function G that corresponds to our original problem is G(y, y′) = y √ 1 − y′2. (2.2.7) Substitution of this function into the EL equation will produce a second-order dif- ferential equation for y(s). Solving this will give us the curve that maximizes the area. When we substitute Eq. (2.2.7) into Eq. (2.2.6) we must first calculate the derivative of G with respect to y′, treating y as a constant parameter. This is ∂G ∂y′ = −yy′ [ 1 − y′2 ]−1/2 . We next differentiate this with respect to s. Because ∂G/∂y′ ≡ Gy′ depends on s through its dependence on both y and y′, we must apply the chain rule. This gives d ds ∂G ∂y′ = ∂Gy′ ∂y dy ds + ∂Gy′ ∂y′ dy′ ds = ∂Gy′ ∂y y′ + ∂Gy′ ∂y′ y′′. We have ∂Gy′ ∂y = −y′ [ 1 − y′2 ]−1/2 and ∂Gy′ ∂y′ = −y [ 1 − y′2 ]−3/2 , so that d ds ∂G ∂y′ = −y′2 [ 1 − y′2 ]−1/2 − yy′′ [ 1 − y′2 ]−3/2 . The remaining quantity to calculate is ∂G ∂y = [ 1 − y′2 ]1/2 . After cleaning up the algebra we find that the EL equation is yy′′ − y′2 + 1 = 0. (2.2.8) This is a nonlinear, second-order differential equation for the function y(s). Exercise 2.1. Make sure that you can reproduce the computations that lead to Eq. (2.2.8). The general solution to Eq. (2.2.8) is y = 1 c1 sin c1(s + c2), 56 Lagrangian mechanics where c1 and c2 are two constants. That this is indeed a solution can be verified by direct substitution; that this is the general solution can be seen from the fact that it depends on two arbitrary constants, the correct number for a second-order differential equation. These constants are determined by enforcing the boundary conditions y(s = 0) = 0 and y(s = π) = 0, which follow from the requirement that the maximum curve must link the points (−1, 0) and (+1, 0). The first condition gives (1/c1) sin(c1c2) = 0, which implies that c2 = 0. The second condition gives (1/c1) sin(c1π) = 0, which implies that c1 must be an integer, which we call n. We therefore have y(s) = 1 n sinns, y′(s) = cos ns. We may now look for x(s), which is determined by Eq. (2.2.2), x′ = √ 1 − y′2 = √ 1 − cos2 ns = sin ns. This integrates to x = x0−(1/n) cos ns, where x0 is another constant of integration. We must now impose the boundary conditions x(s = 0) = −1 and x(s = π) = +1. The first condition gives x0 − 1/n = −1, so that x0 = −1 + 1/n. The second condition gives −1 + (1 − cos nπ)/n = 1, or cos nπ = 1 − 2n, which implies that n = 1. We therefore have x = − cos s, and the constraint n = 1 also implies y = sin s. Exercise 2.2. Verify that y = (1/c1) sin c1(s + c2) is a solution to Eq. (2.2.8), and verify that the choices c1 = 1 and c2 = 0 are appropriate given the boundary conditions. Our final result is this: The curve that maximizes the area A is described by the parametric relations x̄(s) = − cos s, ȳ(s) = sin s, 0 < s < π. (2.2.9) This is a half-circle of unit radius that links the points (−1, 0) and (+1, 0). The maximum area is then given by Amax = ∫ π 0 ȳ(s) dx̄ ds ds = ∫ π 0 sin2 s ds = π 2 ≃ 1.5708. To test whether this is really a maximum we evaluate A for a different choice of curve, one which consists of two straight segments. The first segment connects the points (−1, 0) and (0, y0), while the second segment connects the points (0, y0) and (1, 0). The length of each segment is ℓ = √ 1 + y20 . Because the total length of the curve must be equal to π, we must set y0 = √ (π/2)2 − 1. The area under this curve is the area of a triangle of base 2 and height y0, so A = 1 2 (2)(y0) = √ (π/2)2 − 1 ≃ 1.2114. This area is indeed smaller than Amax. 2.2.4 Path of minimum length The calculus of variations, introduced in Sec. 2.2.2, can be employed to solve many different problems involving either the maximization or minimization of a functional. A simple example is the problem of finding the curve y(x) that minimizes the distance between two fixed points in the x-y plane. We already know that the answer is a straight line, but it will be comforting to use the calculus to give a mathematical proof of this statement. 2.2 Calculus of variations 57 We shall take the two points to be (0, 0) and (x1, y1), respectively. We want to calculate the distance s measured along the curve y(x), and we want to find the path ȳ(x) that minimizes this distance. The increment of distance ds along the curve is easy enough to calculate; it is given by ds = √ dx2 + dy2 = √ 1 + (dy/dx)2 dx = √ 1 + y′2 dx, where we have set y′ = dy/dx. The total distance along the curve is obtained by integration. We have s = ∫ x1 0 √ 1 + y′2 dx. (2.2.10) This is a functional of the path y(x), and we wish to minimize this functional. So here s plays the role of A[y], and x plays the role of the old parameter s. The function G is given by G(y, y′) = √ 1 + y′2. (2.2.11) Notice that this depends only on y′; there is no explicit dependence on y. The EL equation for this situation is d dx ∂G ∂y′ − ∂G ∂y = 0. Because G does not depend explicitly on y we have that ∂G/∂y = 0. The EL equation implies d dx ∂G ∂y′ = 0, and this states that the quantity ∂G/∂y′ is in fact a constant, independent of x. We shall call this constant c. Calculating ∂G/∂y′ gives y′/ √ 1 + y′2, and we have obtained the statement y′ √ 1 + y′2 = c. This equation can easily be solved for y′, and we get y′ = c√ 1 − c2 ≡ m, where m is a new constant. Integration of this equation is straightforward, and we obtain y(x) = mx + b, where b is a final constant of integration. This is the equation of the straight line, the result we expected. The constants m and b can be determined from the boundary conditions, y(x = 0) = 0 and y(x = x1) = y1. The first condition implies b = 0, while the second condition implies m = y1/x1. The final result is therefore that the path which minimizes the distance between (0, 0) and (x1, y1) is described by ȳ(x) = y1 x1 x. (2.2.12) That this is indeed a minimum, instead of a maximum, is obvious from the fact that the maximum distance between two points is always infinite. 60 Lagrangian mechanics 2 1.5 1 0.5 0 0 1 2 3 4 5 6 z x z Figure 2.6: A cycloid that connects the points (0, 0) and (5, 1). This feature of the minimal slide is surprising. Can we be sure that this slide truly minimizes the time? Would not a straight slide do a better job? To convince ourselves that we do have the minimal slide, let us compare the times required for the particle to go from (0, 0) to (5, 1) when it uses either the cycloid or a straight slide. We shall calculate √ 2gt[x] for each case and compare the answers. For the cycloid we have √ 2gtcycloid = ∫ z1 0 √ 1 + x′2√ z dz. With the change of variables introduced above we have x′ = (dx/dθ)/(dz/dθ) = (1 − cos θ)/ sin θ, so that √ 1 + x′2 = √ sin2 θ + (1 − cos θ)2 sin θ = √ 2(1 − cos θ) sin θ . It follows that √ 2gtcycloid = ∫ θ1 0 √ 2(1 − cos θ) sin θ a sin θ dθ √ a(1 − cos θ) = √ 2a ∫ θ1 0 dθ, or √ 2gtcycloid = √ 2aθ1 ≃ 6.1466, using the numerical values listed previously. The shape of the straight slide is described by x = 5z, which implies that x′ = 5. In this case we have √ 2gtstraight = ∫ 1 0 √ 26√ z dz = 2 √ 26z1/2 ∣ ∣ ∣ 1 0 , or √ 2gtstraight = 2 √ 26 ≃ 10.198, 2.2 Calculus of variations 61 and this is a larger number. We have found that, sure enough, tcycloid < tstraight. The particle spends less time on the cycloid than on the straight slide, in spite of the fact that loses speed on the way up toward (x1, z1). The reason is that it picks up a lot of speed on the way down, and this more than makes up for the loss of speed on the way up. The straight slide just does not measure up. 2.2.6 Multiple paths It is useful, and necessary, to generalize the calculus of variations to functionals A that depend not on one path only, but on a collection of paths. In this final subsection we consider the task of extremizing the multi-path functional A[y1, y2, · · ·] = ∫ s1 s0 G(y1, y ′ 1; y2, y ′ 2; · · ·) ds (2.2.17) with respect to each individual path ya(s); the index a = 1, 2, · · · is used to label each path within the collection. This generalization is straightforward. For each variable ya(s) within the col- lection we select a reference path ȳa(s) and we calculate A[ȳ1, ȳ2, · · ·]. We then displace each path from ȳa(s) to ȳa(s) + δya(s) and calculate the new value A[ȳ1 + δy1, ȳ2 + δy2, · · ·] for the functional. The extremum of A is found by demanding that the variation δA = A[ȳ1 + δy1, ȳ2 + δy2, · · ·]−A[ȳ1, ȳ2, · · ·] vanish to first order in the displacements δya(s). As before we impose that the reference and displaced paths all begin and end at the same end points, ya(s0) and ya(s1). We therefore impose that the variations δya all vanish at the end points, δya(s0) = 0 = δya(s1). The change in functional that occurs when we displace the paths from the ref- erence paths ȳa(s) is δA = ∫ s1 s0 [ G(ȳ1 + δy1, ȳ ′ 1 + δy ′ 1; ȳ2 + δy2, ȳ ′ 2 + δy ′ 2; · · ·) − G(ȳ1, ȳ1′; ȳ2, ȳ′2; · · ·) ] ds. By Taylor’s theorem, G(ȳ1 + δy1, ȳ ′ 1 + δy ′ 1; ȳ2 + δy2, ȳ ′ 2 + δy ′ 2; · · ·) = G(ȳ1, ȳ1′; ȳ2, ȳ′2; · · ·) + ∂G ∂y1 ∣ ∣ ∣ ∣ y1=ȳ1,y′1=ȳ ′ 1 ;y2=ȳ2,y′2=ȳ ′ 2 ;··· δy1 + ∂G ∂y′1 ∣ ∣ ∣ ∣ y1=ȳ1,y′1=ȳ ′ 1 ;y2=ȳ2,y′2=ȳ ′ 2 ;··· δy′1 + ∂G ∂y2 ∣ ∣ ∣ ∣ y1=ȳ1,y′1=ȳ ′ 1 ;y2=ȳ2,y′2=ȳ ′ 2 ;··· δy2 + ∂G ∂y′2 ∣ ∣ ∣ ∣ y1=ȳ1,y′1=ȳ ′ 1 ;y2=ȳ2,y′2=ȳ ′ 2 ;··· δy′2 + · · · . Here G is differentiated with respect to each one of its variables, and the partial derivatives are evaluated on the reference paths; we discard all terms that are not linear in the displacements δya and δy ′ a. We have, in a more compact notation, δA = ∫ s1 s0 ∑ a ( ∂G ∂ya δya + ∂G ∂y′a δy′a ) ds, where we sum over all the variables and omit the warning that all partial derivatives must be evaluated on the reference paths, at ya = ȳa and y ′ a = ȳ ′ a. We now write δy′a ≡ y′a − ȳ′a = d ds ( ya − ȳa ) = d ds δya and express the second term within the integral as ∂G ∂y′a d(δya) = d ( ∂G ∂y′a δya ) − δyad ( ∂G ∂y′a ) . 62 Lagrangian mechanics Integrating this term by parts gives δA = ∑ a ∂G ∂y′a δya ∣ ∣ ∣ ∣ s1 s0 + ∑ a ∫ s1 s0 ( ∂G ∂ya − d ds ∂G ∂y′a ) δya ds. Because the displacements must vanish at the end points, the two boundary terms disappear. And because each displacement δya(s) is independent of any other dis- placement, and the displacements are arbitrary in the interval s0 < s < s1, we conclude that δA = 0 ⇒ d ds ∂G ∂y′a − ∂G ∂ya = 0. (2.2.18) We have one EL equation for each path ya(s). This simple statement provides the desired generalization of the calculus of variations to multi-path functionals. 2.3 Hamilton’s principle of least action In Chapter 1 we saw that Newton’s law, F = ma, can serve as the very foundation of all of mechanics; conservation of momentum, angular momentum, and energy could be derived as a consequence of this dynamical law. In this section we offi- cially replace this old foundation by a new one, which is at once more practical, more powerful, and more easily generalizable to other areas of physics. This new foundation will be Hamilton’s principle of least action; the dynamical law F = ma, and the statements of conservation, will all be derived as consequences of this new principle. The principle of least action states that of all the paths qa(t) that a system of particles could take to go from an initial configuration qa(t0) to a final configuration qa(t1), the paths q̄a(t) that the particles actually take are the ones that minimize the action functional S[qa] = ∫ t1 t0 L(qa, q̇a) dt, (2.3.1) where L = T − V (2.3.2) is the Lagrangian function of the mechanical system. The Lagrangian is the differ- ence between T , the system’s total kinetic energy, and V , the total potential energy. The Lagrangian can be expressed in any system of generalized coordinates qa that conveniently describe the system’s degrees of freedom. Because the Lagrangian is a scalar function (as opposed to a vectorial function), the choice of coordinates is immaterial to the formulation of Hamilton’s principle. In particular, it is not nec- essary to adopt Cartesian coordinates attached to an inertial frame. (Of course, nothing prevents us from making this choice if it is convenient.) To find the paths q̄a(t) that minimize the action functional we follow the tech- niques developed in Sec. 2.2. Here S[qa] is a multi-path functional, and the paths qa(t) play the role of the functions ya(s); the Lagrangian plays the role of the func- tion G, and the parameter is the time t. There is no need to repeat the calculations described in Sec. 2.2.6; the conclusion is δS = 0 ⇒ d dt ∂L ∂q̇a − ∂L ∂qa = 0. (2.3.3) These are the Euler-Lagrange (EL) equations for the mechanical system; there is one EL equation for each degree of freedom. The EL equations, when fully worked out, become a set of second-order differential equations for the paths qa(t). The solutions to these equations, which much be subjected to the boundary conditions at t = t0 and t = t1, are the paths q̄a(t) that minimize the action functional. 2.4 Applications of Lagrangian mechanics 65 The equations of motion (2.4.4)–(2.4.6) could also be derived by resolving New- ton’s equation F = ma in the vectorial basis (ρ̂, φ̂, ẑ). The results would be identical, but the computations would be much more laborious. Exercise 2.6. Challenge yourself: Derive Eqs. (2.4.4)–(2.4.6) the hard way, as described in the previous paragraph. Begin by computing the acceleration vector a in terms of the cylindrical coordinates (ρ, φ, z). Next, find the basis vectors ρ̂, φ̂, and ẑ using the method outlined in Sec. 1.2. Finally, resolve the equation ma + ∇V = 0 in this basis, and use the chain rule to calculate ∂V/∂ρ and ∂V/∂φ in terms of ∂V/∂x and ∂V/∂y. The end result should resemble Eqs. (2.4.4)–(2.4.6). If you are not already, after all this you will be fully convinced of the superiority of the Lagrangian methods! 2.4.2 Equations of motion in spherical coordinates Suppose now that a particle moves in the presence of a potential V that is most simply expressed in terms of spherical coordinates (r, θ, φ). Their relation with the usual Cartesian coordinates (x, y, z) is x = r sin θ cos φ, y = r sin θ sinφ, z = r cos θ. (2.4.7) The use of spherical coordinates would be advantageous, for example, when the potential is axially symmetric, so that it depends only on r and θ, or spherically symmetric, when it depends only on r. From Eq. (2.4.7) we obtain the total differentials dx = (sin θ cos φ) dr + (r cos θ cos φ) dθ − (r sin θ sin φ) dφ, dy = (sin θ sinφ) dr + (r cos θ sinφ) dθ + (r sin θ cos φ) dφ, dz = (cos θ) dr − (r sin θ) dθ. It follows that the squared distance between two neighbouring points is given by ds2 = dr2 + r2 dθ2 + r2 sin2 θ dφ2. (2.4.8) The squared velocity is then v2 = ṙ2 + r2θ̇2 + r2 sin2 θ φ̇2, and the Lagrangian is L(r, ṙ; θ, θ̇;φ, φ̇) = 1 2 m(ṙ2 + r2θ̇2 + r2 sin2 θ φ̇2) − V (r, θ, φ). (2.4.9) Exercise 2.7. Verify Eq. (2.4.8). The equations of motion for the particle are obtained by substituting L into the EL equations for qa = (r, θ, φ). We have, from Eq. (2.4.9), ∂L ∂ṙ = mṙ, so that d dt ∂L ∂ṙ = mr̈. We also have ∂L ∂r = mr(θ̇2 + sin2 θ φ̇2) − ∂V ∂r , and the EL equation for r is mr̈ − mr(θ̇2 + sin2 θφ̇2) + ∂V ∂r = 0. (2.4.10) 66 Lagrangian mechanics Moving on, we now have ∂L ∂θ̇ = mr2θ̇, which implies d dt ∂L ∂θ̇ = m d dt ( r2θ̇ ) . We also have ∂L ∂θ = mr2 sin θ cos θ φ̇2 − ∂V ∂θ , and the EL equation gives m d dt ( r2θ̇ ) − mr2 sin θ cos θ φ̇2 + ∂V ∂θ = 0. (2.4.11) Finally, we have ∂L ∂φ̇ = mr2 sin2 θ φ̇, which implies d dt ∂L ∂φ̇ = m d dt ( r2 sin2 θ φ̇ ) . We also have ∂L ∂φ = −∂V ∂φ , and the EL equation gives m d dt ( r2 sin2 θ φ̇ ) + ∂V ∂φ = 0. (2.4.12) The equations of motion (2.4.10)–(2.4.12) could also be derived by resolving Newton’s equation F = ma in the vectorial basis (r̂, θ̂, φ̂). The results would be identical, but as in the preceding subsection the computations would be much more laborious. Exercise 2.8. Challenge yourself once again: Derive Eqs. (2.4.10)–(2.4.12) the hard way, as described in the previous paragraph. Or finally cry uncle and pledge allegiance to the Lagrangian way of life! 2.4.3 Motion on the surface of a cone As our first real application of the Lagrangian formalism, we consider a particle that is constrained to move on the surface of a cone, subjected to gravity. As shown in Fig. 2.7, the cone has an opening angle of 2α, and it is placed vertically in the gravitational field. The particle is at a distance r(t) from the cone’s apex, and at an angle φ(t) relative to the x axis. Because the particle is confined to the cone’s surface, its angle θ with respect to the z axis is a constant; it is in fact equal to α. The motion of the particle is best described in terms of spherical coordinates (r, θ, φ), with θ restricted at all times to the value α. According to the results of Sec. 2.4.2, its kinetic energy is T = 12m(ṙ 2 + r2 sin2 α φ̇2), and its potential energy is V = mgz = mgr cos α. The Lagrangian is therefore L(r, ṙ;φ, φ̇) = 1 2 m(ṙ2 + r2 sin2 α φ̇2) − mgr cos α. (2.4.13) The equations of motion for r(t) and φ(t) are obtained by substituting this La- grangian into the EL equations. 2.4 Applications of Lagrangian mechanics 67 y x z α φ r Figure 2.7: A particle of mass m moves on the surface of a cone of opening angle 2α. The motion is described by the coordinates r(t) and φ(t). We have ∂L ∂ṙ = mṙ, so that d dt ∂L ∂ṙ = mr̈. We also have ∂L ∂r = mr sin2 α φ̇2 − mg cos α, and the EL equation for r is r̈ − r sin2 α φ̇2 + g cos α = 0. (2.4.14) Moving on, we observe that L is independent of φ, and the fact that ∂L/∂φ = 0 means that the EL equation for φ reduces to d dt ∂L ∂φ̇ = 0. This implies that the quantity ∂L/∂φ̇ is a constant, which we shall call mh. Calcu- lating the partial derivative gives ∂L/∂φ̇ = mr2 sin2 α φ̇, and we finally obtain the statement r2 sin2 α φ̇ = h = constant. (2.4.15) The quantity h is readily interpreted as the z component of the particle’s reduced angular momentum vector, and it is a constant of the motion. Equation (2.4.15) shows that φ̇ is always of the same sign; the angular part of the motion is monotonic. Substituting φ̇ = h/(r2 sin2 α) into Eq. (2.4.14) produces r̈ − h 2 r3 sin2 α + g cos α = 0. This equation can be integrated by using the standard trick of multiplying each term by ṙ (recall that we used this trick back in Sec. 1.5.6). We have r̈ṙ − h 2ṙ r3 sin2 α + gṙ cos α = 0, 70 Lagrangian mechanics θ z φ x y l Figure 2.10: The motion of a spherical pendulum is described in terms of the angles θ(t) and φ(t). Exercise 2.9. Verify the quoted relation between h2 and r±. 2.4.4 Spherical pendulum We now examine the situation of a pendulum which is free to move in all directions about its pivot point. The pendulum has a mass m, a constant length ℓ, and its motion is described in terms of the two angles θ(t) and φ(t), as shown in Fig. 2.10. These coordinates are related to the standard Cartesian coordinates by x = ℓ sin θ cos φ, y = ℓ sin θ sin φ, z = ℓ cos θ. As shown in the figure, the z axis is pointing down, in the direction of the grav- itational acceleration g. It is clear that we are once more dealing with spherical coordinates. This time, however, it is the radial coordinate r that is held fixed to the value ℓ. According to the results of Sec. 2.4.2 the pendulum’s kinetic energy is T = 12mℓ 2(θ̇2 + sin2 θ φ̇2). Its potential energy is V = −mgz = −mgℓ cos θ = −mℓ2ω2 cos θ, where we have re-introduced the quantity ω = √ g/ℓ. (2.4.18) The pendulum’s Lagrangian is L(θ, θ̇;φ, φ̇) = 1 2 mℓ2(θ̇2 + sin2 θ φ̇2) + mℓ2ω2 cos θ. (2.4.19) The equations of motion for θ(t) and φ(t) are obtained by substituting this La- grangian into the EL equations. We compute ∂L ∂θ̇ = mℓ2θ̇, 2.4 Applications of Lagrangian mechanics 71 which implies d dt ∂L ∂θ̇ = mℓ2θ̈. We also have ∂L ∂θ = mℓ2 sin θ cos θ φ̇2 − mℓ2ω2 sin θ, and the EL equation for θ is θ̈ − sin θ cos θ φ̇2 + ω2 sin θ = 0. (2.4.20) Moving on, we observe that L is independent of φ, and the fact that ∂L/∂φ = 0 means that the EL equation for φ reduces to d dt ∂L ∂φ̇ = 0. This implies that the quantity ∂L/∂φ̇ is a constant, which we shall call mℓ2h. Calculating the partial derivative gives ∂L/∂φ̇ = mℓ2 sin2 θ φ̇, and we finally obtain the statement sin2 θ φ̇ = h = constant. (2.4.21) The quantity h is once more interpreted as the z component of the pendulum’s reduced angular momentum vector, and it is a constant of the motion. In the special case h = 0 the pendulum is prevented to move in the φ direction, and Eq. (2.4.20) for θ reduces to θ̈+ω2 sin θ = 0; this is the same equation that was first derived in Sec. 1.3.7, and then again in Sec. 2.1, and which describes the motion of a planar pendulum. In the general case (h 6= 0) we see that φ̇ is always of the same sign, so that φ(t) is a monotonic function of time; this means that the pendulum rotates in a consistent direction around the z axis. With the substitution φ̇ = h/ sin2 θ Eq. (2.4.20) becomes θ̈ − h 2 cos θ sin3 θ + ω2 sin θ = 0. Multiplying each term by θ̇ allows us to integrate this equation. The result is the conservation statement 1 2 θ̇2 + ν(θ) = ε = constant, (2.4.22) where ε is the pendulum’s reduced total mechanical energy, and ν(θ) = h2 2 sin2 θ − ω2 cos θ (2.4.23) is an effective potential for the motion in the θ direction. Equations (2.4.22) and (2.4.23) give rise to the energy diagram of Fig. 2.11. From this diagram we may immediately conclude that the motion takes place between two turning points at θ = θ±; these are determined by the condition ν(θ±) = ε. Exercise 2.10. Verify that Eqs. (2.4.22) and (2.4.23) do indeed follow from the equations of motion. In Fig. 2.12 we present the results of a numerical integration of the equations of motion, which we recast into the first-order form θ̇ = v, v̇ = h2 cos θ sin3 θ − ω2 sin θ, φ̇ = h sin2 θ . 72 Lagrangian mechanics 0 1 2 3 4 5 6 0.5 1 1.5 2 2.5 3 θ Figure 2.11: Energy diagram for the spherical pendulum. The motion always takes place between two turning points at θ = θ±. We start the integration at θ = θ−, setting v = 0 and φ = 0. The constant h can be determined in terms of θ− and θ+ by using the relation ν(θ−) = ν(θ+), which follows from Eq. (2.4.22). The result, after some algebra, is h2 = 2ω2 (sin θ+ sin θ−) 2(cos θ− − cos θ+) (sin θ+ − sin θ−)(sin θ+ + sin θ−) . Exercise 2.11. Verify the quoted relation between h2 and θ±. 2.4.5 Rotating pendulum Another variation on the pendulum theme has the pivot point of a planar pendulum forced to rotate with a constant angular velocity Ω on a circle of radius a. This situation is shown in Fig. 2.13. Once more we describe the motion of the pendulum in terms of the swing angle θ(t), which is defined relative to the vertical direction; this we now associate with the y-direction. The Cartesian coordinates of the pendulum, relative to the pivot point, are xrelative = ℓ sin θ, yrelative = −ℓ cos θ. The Cartesian coordinates of the pivot point are xpivot = a cos Ωt, ypivot = a sin Ωt. The Cartesian coordinates of the pendulum, relative to the origin of the coordinate system, are therefore x = a cos Ωt + ℓ sin θ, y = a sin Ωt − ℓ cos θ. (2.4.24) 2.4 Applications of Lagrangian mechanics 75 We next compute ∂L ∂θ = maℓΩθ̇ cos(θ − Ωt) − mℓ2ω2 sin θ, and substituting all this into the EL equation produces θ̈ + ω2 sin θ − (a/ℓ)Ω2 cos(θ − Ωt) = 0. (2.4.26) This is the equation of motion of a rotating pendulum. Equation (2.4.26) cannot be integrated with the help of the θ̇ trick; this is prevented by the fact that the equation depends explicitly on time through the term in cos(θ − Ωt). As a consequence, the motion cannot be analyzed with the help of an energy diagram; this can be understood from the very fact that the total mechanical energy of the rotating pendulum is not conserved. The only tool that remains at our disposal to analyze the motion is numerical integration, and Fig. 2.14 displays the results. The graphs reveal that when the pendulum is driven at a frequency Ω that is close to its natural frequency ω, the response is more violent: the amplitude of the oscillations is then much larger. This is the phenomenon of resonance. This phenomenon can be illustrated in the context of a simpler model, one which can be solved exactly. We consider a simple harmonic oscillator which is driven by an oscillating external force. The equations of motion for this simplified model is θ̈ + ω2θ = A cos Ωt. (2.4.27) When Ω 6= ω a solution to this equation is θ(t) = A ω2 − Ω2 cos Ωt (Ω 6= ω). (2.4.28) In this situation the pendulum oscillates at the driving frequency Ω, and the oscil- lations have a constant amplitude. Notice, however, that the amplitude grows as Ω approaches the natural frequency ω. The solution of Eq. (2.4.28) is not valid when Ω = ω. In this case we have instead θ(t) = At 2ω sin ωt (Ω = ω). (2.4.29) In this case the oscillations keep growing in amplitude; the simple harmonic oscil- lator is in resonance with the driving force. Exercise 2.13. Verify that Eqs. (2.4.28) and (2.4.29) are solutions to Eq. (2.4.27). A more challenging question: What is the general solution to Eq. (2.4.27) when Ω 6= ω and when Ω = ω? The general solution should be parameterized in terms of the initial angle θ(0) and the initial angular velocity θ̇(0). What choices of initial conditions give rise to Eqs. (2.4.28) and (2.4.29)? When the rotating pendulum is driven at resonance we observe a growth in the amplitude of oscillations, but this growth is bounded; it saturates and the amplitude then starts to decrease. This saturation is produced by nonlinear effects: When the amplitude grows the natural period of the oscillations changes (as we learned back in Sec. 1.3.7) and the pendulum is no longer driven at its natural frequency. As resonance stops the amplitude starts to decrease and the pendulum’s natural frequency returns to its original value. At this stage the conditions are once more suitable for a resonant growth of the amplitude, and the cycle repeats. For certain choices of parameters the driving force can have a dramatic influence on the pendulum. This is illustrated in Fig. 2.15, for which the driving frequency was set to Ω = 0.9ω. Here we see the driving force causing the pendulum to go beyond θ = π, completing one or two revolutions before returning to a short oscillation cycle. 76 Lagrangian mechanics -30 -20 -10 0 10 20 30 0 2 4 6 8 10 sw in g an gl e (d eg re es ) time unperturbed pendulum driven pendulum (low frequency) -30 -20 -10 0 10 20 30 0 2 4 6 8 10 sw in g an gl e (d eg re es ) time unperturbed pendulum driven pendulum (high frequency) -100 -50 0 50 100 0 2 4 6 8 10 sw in g an gl e (d eg re es ) time unperturbed pendulum driven pendulum (resonant frequency) Figure 2.14: The motion of a rotating pendulum. Each graph shows the swing angle θ(t) of the driven pendulum in blue, and the swing angle of a free pendulum in red. In the first graph the pendulum is driven at a low frequency set at Ω/ω = 0.4. In the second graph the pendulum is driven at a high frequency set at Ω/ω = 2.4. In the third graph the pendulum is driven at resonant frequency, so that Ω/ω = 1.0; notice the large amplitude of oscillations in this case. In all cases we have set (a/ℓ)Ω2 = 0.2, and the initial conditions are θ(0) = 0.2 and θ̇(0) = 0.3. 2.4 Applications of Lagrangian mechanics 77 -200 0 200 400 600 800 1000 1200 1400 0 2 4 6 8 10 sw in g an gl e (d eg re es ) time unperturbed pendulum driven pendulum (near-resonant frequency) Figure 2.15: The motion of a rotating pendulum, with Ω/ω = 0.9, (a/ℓ)Ω2 = 0.2, θ(0) = 0.2, and θ̇(0) = 0.3. 2.4.6 Rolling disk As our next application we consider a disk of mass m and radius R that rolls without slipping on an inclined plane of total length ℓ; the plane’s inclination relative to the horizontal is α. As shown in Fig. 2.16, the distance from the top position on the plane to the disk’s centre of mass — its geometric centre — is denoted s, and θ is the angle of a selected point on the disk’s rim relative to an axis perpendicular to the inclined plane. There is both a translational motion of the centre of mass and a rotational motion of the disk in this problem. The disk’s kinetic energy is T = 1 2 mṡ2 + 1 2 Iθ̇2, where I = 12mR 2 is the disk’s moment of inertia. The coordinates s and θ, however, are not independent; they are related by the no-slip condition, which implies s = Rθ. So we have ṡ = Rθ̇ and the kinetic energy becomes T = 1 2 mR2θ̇2 + 1 4 mR2θ̇2 = 3 4 mR2θ̇2. The disk’s potential energy is V = mgz = mg(l − s) sin α = mg(ℓ − Rθ) sinα. The Lagrangian is therefore L(θ, θ̇) = 3 4 mR2θ̇2 − mg(ℓ − Rθ) sinα. (2.4.30) To obtain the disk’s equation of motion we substitute this into the EL equation. We first compute ∂L ∂θ̇ = 3 2 mR2θ̇, which implies d dt ∂L ∂θ̇ = 3 2 mR2θ̈. 80 Lagrangian mechanics and ∂Lrel ∂r = µrφ̇2 − GµM r2 . This gives r̈ − rφ̇2 + GM r2 = 0. (2.4.37) This is the same statement as Eq. (1.5.25). To obtain the equation of motion for φ we observe that since Lrel is independent of φ, we must have that ∂Lrel/∂φ̇ is a constant of the motion. Calling this constant µh and calculating the partial derivative, we get r2φ̇ = h = constant. (2.4.38) This is the same statement as Eq. (1.5.27), and h is identified as the reduced angular momentum of the two-body system. The equations of motion (2.4.37) and (2.4.38) can be analyzed with the same mathematical techniques as those employed in Sec. 1.5. It should be clear that compared with the Newtonian methods of Chapter 1, the Lagrangian methods provide a much simpler way of obtaining these equations. 2.5 Generalized momenta and conservation statements 2.5.1 Conservation of generalized momentum In the applications of Lagrangian mechanics presented in Sec. 2.4 it occurred a number of times that the Lagrangian was independent of one of the generalized coordinates (mostly it was the φ coordinate), and we saw that this fact always translated into the existence of a constant of the motion (which we usually called h). A specific example is the case of a particle moving on the surface of a cone (Sec. 2.4.3), for which the Lagrangian is indeed independent of φ and for which the constant of the motion was h = r2 sin2 α φ̇. A similar situation occurred for the spherical pendulum (Sec. 2.4.4) and for Kepler’s problem (Sec. 2.4.7). It is easy to generalize this discussion and to derive the very useful fact that whenever the Lagrangian does not depend explicitly on one (or more) of the gen- eralized coordinates qa, there exists a corresponding constant of the motion. To establish this statement we shall first introduce the notion of a generalized momen- tum. Consider a Lagrangian L(qa, q̇a) that depends on a number of generalized coor- dinates qa and a number of generalized velocities q̇a. The quantities pa = ∂L ∂q̇a (2.5.1) feature prominently in the EL equations, which can be written in the form ṗa = ∂L ∂qa . (2.5.2) The quantities pa are the generalized momenta of the mechanical system. There is one generalized momentum pa for each generalized coordinate qa. The generalized momenta can represent either a component of the linear-momentum vector or a component of the angular-momentum vector. Generally speaking, when- ever qa represents a linear variable the corresponding pa will be a linear momen- tum; and whenever qa represents an angular variable its corresponding pa will be 2.5 Generalized momenta and conservation statements 81 an angular momentum. Consider, for example, the Lagrangian of a free particle in cylindrical coordinates (Sec. 2.4.1). This is L = 1 2 m(ρ̇2 + ρ2φ̇2 + ż2). The generalized momenta are pρ = ∂L ∂ρ̇ = mρ̇, pφ = ∂L ∂φ̇ = mρ2φ̇, pz = ∂L ∂ż = mż. In the case of pρ and pz we clearly have quantities that represent components of a linear-momentum vector. But the case of pφ is different. Here we have pφ = m(ρ)(ρφ̇), and this clearly represents the component of an angular-momentum vector. Exercise 2.15. Show that pρ = p · ρ̂ and pz = p · ẑ, where p is the particle’s momentum vector. Show, on the other hand, that pφ = L · ẑ, where L is the particle’s angular-momentum vector. Suppose now that a Lagrangian L(q1, q̇1; q2, q̇2, · · ·) happens not to depend ex- plicitly on one of its generalized coordinates, say q∗. Then ∂L ∂q∗ = 0 and it follows from the EL equation for q∗ that dp∗ dt = 0, where p∗ = ∂L/∂q̇∗ is the generalized momentum associated with the coordinate q∗. This equation states that p∗ is a constant of the motion, and we have established the following theorem: Whenever the Lagrangian of a mechanical system does not depend ex- plicitly on a generalized coordinate q∗, the corresponding generalized momentum p∗ = ∂L/∂q̇∗ is a constant of the motion. A coordinate q∗ that does not appear in L is sometimes called a cyclic coordinate. A Lagrangian may contain any number of cyclic coordinates. As an example consider the following Lagrangian, again in cylindrical coordi- nates, L = 1 2 m(ρ̇2 + ρ2φ̇2 + ż2) − V (ρ). Here it is assumed that the potential energy V depends only on ρ; the mechanical system is cylindrically symmetric. This implies that φ and z are cyclic coordinates, and that pφ = mρ 2φ̇ and pz = mż are constants of the motion. This theorem on cyclic coordinates and conserved quantities is extremely im- portant and very useful. To find all the constants of the motion is usually a key step during the integration of the equations of motion, and the theorem provides a very efficient algorithm to identify at least some of them. 82 Lagrangian mechanics 2.5.2 Conservation of energy Conservation of total mechanical energy E is also an important aspect of the motion of a mechanical system and a key to solving the equations of motion. In this subsection we show that energy is conserved whenever the Lagrangian does not depend explicitly on time t. To begin the discussion let us consider a Lagrangian L(qa, q̇a, t) that depends on a number of generalized coordinates qa, a number of generalized velocities q̇a, and let us consider the possibility that it depends also explicitly on time. (An example is the Lagrangian of a rotating pendulum, which was written down in Sec. 2.4.5.) Applying the chain rule, we find that the total time derivative of the Lagrangian is given by dL dt = ∑ a ∂L ∂qa q̇a + ∑ a ∂L ∂q̇a q̈a + ∂L ∂t . The first term accounts for the time dependence contained in each qa(t), the second term for the time dependence contained in each q̇a(t), and the third term accounts for the explicit dependence of the Lagrangian on t. We have defined the generalized momenta pa by pa = ∂L ∂q̇a and the EL equations can be expressed in the form ṗa = ∂L ∂qa . We make these substitutions in the previous equation, and obtain dL dt = ∑ a ( ṗaq̇a + paq̈a ) + ∂L ∂t , or dL dt = d dt ( ∑ a paq̇a ) + ∂L ∂t , which is equivalent to the previous form by virtue of the chain rule. We have obtained the equation d dt ( ∑ a paq̇a − L ) = −∂L ∂t , (2.5.3) and a statement of conservation follows immediately: Whenever L does not depend explicitly on time, so that ∂L/∂t = 0, we have that h(qa, q̇a) ≡ ∑ a paq̇a − L (2.5.4) is a constant of the motion, dh/dt = 0. Surely the function h(qa, q̇a) must have something to do with the system’s total mechanical energy. Let us first figure out the relationship in the context of a sim- ple example. We go back to the Lagrangian of a particle expressed in cylindrical coordinates, L = 1 2 m(ρ̇2 + ρ2φ̇2 + ż2) − V (ρ, φ, z), 2.5 Generalized momenta and conservation statements 85 On the other hand, ∂∆L ∂qa = ∂ ∂qa ( ∑ b ∂f ∂qb q̇b + ∂f ∂t ) = ∑ b ∂2f ∂qa∂qb q̇b + ∂2f ∂qa∂t , and this is equal to the previous result, because the order in which one evaluates second partial derivatives does not matter. We conclude that d dt ∂∆L ∂q̇a = ∂∆L ∂qa and that the change of Lagrangian has no effect on the EL equations. The equations of motion that derive from L′ = L+df/dt are indeed identical to the equations that derive from L. There is a more elegant way to prove this result. In this alternative derivation we appeal directly to Hamilton’s principle. The action S′ = ∫ t1 t0 L′ dt associated with L′, and the action S = ∫ t1 t0 Ldt associated with L, are related by S′ = S + ∫ t1 t0 ∆Ldt = S + ∫ t1 t0 df dt dt = S + f ( qa(t1), t1 ) − f ( qa(t0), t0 ) . The equations of motion are obtained from S′ or S by varying the paths qa(t) and demanding that the variation of the action be zero to first order in the variations δqa(t). The variations, you may recall, must be subjected to the boundary condi- tions δqa(t0) = δqa(t1) = 0; the paths must all begin at the same qa(t0) and end at the same qa(t1). But under these conditions we find that the values f(qa(t0), t0) and f(qa(t1), t1) can never change under a variation of the paths, and we must conclude that δS′ = δS. An extremum of S will also be an extremum of S′, and the equations of motion derived from L and L′ are guaranteed to be the same. While the operation L → L′ = L+df/dt does not affect the equations of motion, it may nevertheless change the expressions for the generalized momenta pa and the total energy h. The new momenta p′a are given by p′a = ∂L′ ∂q̇a = ∂L ∂q̇a + ∂∆L ∂q̇a , or p′a = pa + ∂f ∂qa , (2.5.7) according to our previous computations. The new energy function h′ is given by h′ = ∑ a p ′ aq̇a − L′, so h′ = ∑ a ( pa + ∂f ∂qa ) q̇a − L − ∆L = h + ∑ a ∂f ∂qa q̇a − ∑ b ∂f ∂qb q̇b − ∂f ∂t . 86 Lagrangian mechanics The sums cancel each other out, and we are left with h′ = h − ∂f ∂t . (2.5.8) We find that the expression for the energy is affected only when f depends explicitly on time. 2.6 Charged particle in an electromagnetic field The Lagrangian formulation of mechanics is well suited to mechanical systems for which the forces can all be derived from a potential-energy function V (qa); these forces will depend on the positions qa, but that they might also depend on the velocities q̇a is normally out of the question. There is, however, an important mechanical system for which the forces do depend on velocity: a charged particle moving in the presence of an electromagnetic field. In this case the particle is subjected to the Lorentz force, and the equations of motion are ma = q(E + v × B). (2.6.1) Can this equation be derived on the basis of a Lagrangian? The answer is in the affirmative. An interesting property of this Lagrangian is that it depends on the scalar potential Φ and vector potential A instead of depending on the fields E and B. Recall that the fields can be expressed in terms of the potentials as E = −∂A ∂t − ∇Φ, B = ∇ × A. (2.6.2) The potentials are usually introduced to simplify the structure of Maxwell’s equa- tions. The definition of E implies ∇ × E = − ∂ ∂t ∇ × A − ∇ × (∇Φ); since the curl of a gradient is always zero, this gives ∇ × E = −∂B ∂t , one of the four Maxwell equations. Similarly, the definition of B implies ∇ · B = ∇ · (∇ × A); since the divergence of a curl is always zero, this gives ∇ · B = 0, another one of the Maxwell equations. The remaining two equations can then be recast into equations that Φ and A must satisfy. It is convenient to express the fields in terms of the potentials, but it is important to understand that the potentials do not have direct physical meaning. Indeed, it is even possible to change the potentials by a certain transformation and leave the fields unaffected. This transformation is given by Φ → Φ′ = Φ − ∂f ∂t , A → A′ = A + ∇f, (2.6.3) where f(r, t) is an arbitrary function of position and time. Such a transformation of the potentials is known as a gauge transformation, and its defining property is 2.6 Charged particle in an electromagnetic field 87 that the transformation leaves the fields invariant. Different sets of potentials that are related by a gauge transformation describe the same fields and therefore the same physical situation. Exercise 2.17. Show that the transformation of Eq. (2.6.3) leaves the fields unaffected. That is, show that the transformation produces E → E ′ = E and B → B′ = B. The Lagrangian for a particle of charge q in an electromagnetic field is L = 1 2 mv2 − qΦ + qA · v, (2.6.4) where v2 = v · v. As stated previously, this depends on the potentials Φ and A instead of the fields E and B. Another interesting property of the Lagrangian is that the potential-energy term V = qΦ − qA · v depends on the velocity vector v as well as the position r. The dependence on position, of course, comes from the potentials, which may also depend explicitly on time. Let us verify that the Lagrangian of Eq. (2.6.4) does indeed give rise, via the EL equations, to the Lorentz-force equation of Eq. (2.6.1). It will suffice to verify the x component of the equation, which we write as mẍ = qEx + q(v × B)x = qEx + q(ẏBz − żBy). Similar computations would allow us to verify also the y and z components, but we will not present these here. We begin by presenting the Lagrangian in a more explicit form, as L = 1 2 m(ẋ2 + ẏ2 + ż2) − qΦ + q(ẋAx + ẏAy + żAz). We have ∂L ∂ẋ = mẋ + qAx, and this implies d dt ∂L ∂ẋ = mẍ + q ( ∂Ax ∂x ẋ + ∂Ax ∂y ẏ + ∂Ax ∂z ż + ∂Ax ∂t ) . In this step we took into account the fact that Ax depends on time through its dependence on the coordinates x(t), y(t), and z(t), and also through its own explicit dependence on t; the total time derivative had to be evaluated by using the chain rule. Finally, we have ∂L ∂x = −q ∂Φ ∂x + q ( ẋ ∂Ax ∂x + ẏ ∂Ay ∂x + ż ∂Az ∂x ) . The EL equations are 0 = mẍ+q ( ∂Ax ∂t + ∂Φ ∂x ) +qẋ ( ∂Ax ∂x − ∂Ax ∂x ) −qẏ ( ∂Ay ∂x − ∂Ax ∂y ) +qż ( ∂Ax ∂z − ∂Az ∂x ) , or mẍ = q ( −∂Ax ∂t − ∂Φ ∂x ) + qẏ ( ∂Ay ∂x − ∂Ax ∂y ) − qż ( ∂Ax ∂z − ∂Az ∂x ) . In the first set of brackets we recognize Ex, in the second Bz, and in the third By. We therefore have mẍ = qEx + q(ẏBz − żBy) = q(E + v × B)x, 90 Lagrangian mechanics y y’ z = z’ θ r φ Ωt x’ x P r sinθ Figure 2.17: The rotating frame S of the turntable, as viewed in the inertial frame S′. A point P is referred to the inertial frame by its Cartesian coordinates (x′, y′, z′) or its spherical coordinates (r, θ, φ′ = φ+Ωt). It is referred to the rotating frame by its Cartesian coordinates (x, y, z) or its spherical coordinates (r, θ, φ). and x′ + iy′ = r sin θ(cos φ′ + i sin φ′) = r sin θ eiφ ′ . Then we have x′ + iy′ = r sin θ ei(φ+Ωt) = eiΩtr sin θ eiφ, or x′ + iy′ = eiΩt(x + iy). (2.7.1) When fully expanded, this is x′ = x cos Ωt − y sinΩt, y′ = y cos Ωt + x sin Ωt, z′ = z. (2.7.2) Exercise 2.19. Verify that Eq. (2.7.2) follows from Eq. (2.7.1). Then work out the inverse transformation, (x′, y′, z′) → (x, y, z). A particle moving in the rotating frame S with a position vector r(t) = [x(t), y(t), z(t)] moves in the inertial frame S′ with a position vector r′(t) = [x′(t), y′(t), z′(t)]; these are related by the transformation of Eq. (2.7.2). The components of the velocity vectors are then related by ẋ′ = ẋ cos Ωt − ẏ sinΩt − Ω(x sin Ωt + y cos Ωt), ẏ′ = ẏ cos Ωt + ẋ sinΩt − Ω(y sin Ωt − x cos Ωt), ż′ = ż. After a fairly laborious calculation, we find that the squared velocity, as measured in the inertial frame, is v′2 = ẋ′2 + ẏ′2 + ẏ′2 = ẋ2 + ẏ2 + ż2 − 2Ω(yẋ − xẏ) + Ω2(x2 + y2). (2.7.3) 2.7 Motion in a rotating reference frame 91 The particle’s kinetic energy is then T = 12mv ′2. It contains a contribution from ṙ, the particle’s velocity vector as measured in the rotating frame, and contributions from the rotational motion of the frame (the terms that involve Ω). Exercise 2.20. Verify Eq. (2.7.3). You will save yourself some work if you use the trick of forming complex combinations. The particle’s potential energy can be expressed in terms of the inertial coordi- nates (x′, y′, z′), but after the transformation of Eq. (2.7.2) it becomes a function of the rotating coordinates (x, y, z). Denoting this function V (x, y, z), we find that the particle’s Lagrangian is L = 1 2 m(ẋ2 + ẏ2 + ż2) − mΩ(yẋ − xẏ) + 1 2 mΩ2(x2 + y2) − V (x, y, z). (2.7.4) The equations of motion for the particle are then obtained by substituting this into the EL equations for x, y, and z. The computations that lead to the equations of motion will be left as an exercise for the reader. We find mẍ = −∂V ∂x + 2mΩẏ + mΩ2x, (2.7.5) mÿ = −∂V ∂y − 2mΩẋ + mΩ2y, (2.7.6) mz̈ = −∂V ∂z . (2.7.7) These equations can be expressed in vectorial form if we introduce the angular- velocity vector Ω, defined by Ω = Ωẑ = [0, 0,Ω]. (2.7.8) The vectorial form is mr̈ = Fapplied + FCoriolis + Fcentrifugal, (2.7.9) where Fapplied = −∇V (2.7.10) is the true applied force on the particle, given by the gradient of the potential energy, while FCoriolis = 2mṙ × Ω = [2mΩẏ,−2mΩẋ, 0] (2.7.11) and Fcentrifugal = mΩ × (r × Ω) = [mΩ2x,mΩ2y, 0] (2.7.12) are fictitious forces that arise because the reference frame S is not an inertial frame (refer back to the discussion of Sec. 1.1). The Coriolis force is linear in the angular velocity Ω, and it depends on the particle’s velocity vector ṙ; its effect on the particle depends on its state of motion. The centrifugal force is quadratic in Ω, and it depends only on the position vector r; this is always an outward force that points away from the centre of motion. Exercise 2.21. Verify that Eqs. (2.7.5)–(2.7.7) follow from the Lagrangian of Eq. (2.7.4). Exercise 2.22. Show that Eqs. (2.7.5)–(2.7.7) are equivalent to the vectorial equation 92 Lagrangian mechanics (2.7.9), together with the definitions of Eqs. (2.7.10)–(2.7.12). 2.7.2 Case study #1: Particle attached to a spring As our first application of the rotating-frame formalism we examine a particle at- tached to a linear spring that is free to rotate around the z axis. The particle is thus subjected to the potential V = 12k(x 2 + y2), which we write in the alternative form V = 1 2 mω2(x2 + y2), (2.7.13) in which ω2 ≡ k/m is a stand-in for the spring constant k. For simplicity we assume that the particle is confined to the x-y plane and we set, accordingly, z = ż = 0 in all equations. The equations of motion of Eqs. (2.7.5) and (2.7.6) become ẍ = 2Ωẏ + (Ω2 − ω2)x, ÿ = −2Ωẋ + (Ω2 − ω2)y. (2.7.14) Let us first analyze these equations in the limit of no rotation, Ω = 0. In this case they reduce to ẍ = −ω2x and ÿ = −ω2y, the equations of simple harmonic motion. The general solution to the equations of motion is then xΩ=0(t) = a cos(ωt + α), yΩ=0(t) = b cos(ωt + β). (2.7.15) The four constants a, b, α, and β can be related to the four initial values x(0), ẋ(0), y(0), and ẏ(0). Equations (2.7.5) are parametric equations for the motion of the particle in the x-y plane, and it is easy to show that the trajectory is elliptical. To analyze the equations in the general case (Ω 6= 0) we once more employ the clever trick of forming complex combinations. We introduce ξ = x + iy and we combine the two equations (2.7.14) into a single equation for ξ: ξ̈ = ẍ + iÿ = 2Ω(ẏ − iẋ) + (Ω2 − ω2)(x + iy) = −2iΩ(ẋ + iẏ) + (Ω2 − ω2)(x + iy), or ξ̈ + 2iΩξ̇ − (Ω2 − ω2)ξ = 0. (2.7.16) To find solutions to this equation we use a trial expression of the form ξ = ceiλt, where c and λ are complex constants. Substitution into Eq. (2.7.16) produces a quadratic equation for λ, λ2 + 2Ωλ + (Ω2 − ω2) = 0, which factorizes as (λ + Ω + ω)(λ + Ω − ω) = 0. The solutions, obviously, are λ = −(Ω ± ω), and the general solution for ξ is ξ = c1e −i(Ω+ω)t + c2e −i(Ω−ω)t, or ξ = e−iΩt ( c1e −iωt + c2e iωt ) , where c1 and c2 are complex numbers. To help us understand what we have just found, we observe that if we let Ω go to zero, our general solution for ξ becomes c1e −iωt + c2e iωt. This is the solution in