plane pendulum

The Plane Pendulum Problem in Classical Mechanics PhL 4.10.08 Comment: I did this analysis as an "exercise" after reading Goldstein through Chapter 9. In my current pass through his book (my first complete pass, really), I did not do any of his suggested problems, so this is an attempt to at least do something. I take this relatively simple plane pendulum problem and look at it in each of the formalisms that Goldstein has presented. The main motivation was to see it in the Hamilton-Jacobi theory and in particular the action-angle variables. I got sidetracked, however, in solving the actual problem which, for large angles, was much harder than I thought it would be, involving the Jacobi elliptic integrals and functions. Since I knew nothing about this Jacobi world, I had to go off and study it, and that is in a separate writeup called "elliptic functions". When I finished all this stuff, I took a quick look at what the quantum solution is to this same problem. The notes below are summarized in a separate "executive summary" document which goes much lighter on the technical details. I would recommend that to the casual reader. 1. Setup . 2 2. Compute the Hamiltonian 3 3. Write down Hamilton's Equations of motion 3 5. Initial conditions 4 7. Now let's try this method 4 8. Take small angle limit of the exact solution 7 9. Examine the exact solution 8 10. Question: How can we make our solution show "rotation"? 10 11. Comment on the arcsin function. 13 12. Phase Space. 16 13. Hamilton-Jacobi theory applied to this problem. 17 14. Action-angle variables applied to this problem. 19 15. What are the classical waves doing in this problem? 23 16. What is the solution to the quantum plane pendulum problem? 24 Contents of this very long document: 1. Setup for the pendulum problem 2. Compute Lagrangian and Hamiltonian and canonical momenta. 3. Hamilton's Equations of Motion, ω2 = g/l. Mass does not affect anything. 4. Statement of the key equation of motion: = – ω2 sinθ, 2nd order non-linear (in θ) ODE. But the equation is trivially integrable twice. 5. Statement of the initial conditions:, start at bottom with = v, then C = ½ v2 – ω2 6. General method for solving such a non-linear ODE, called an "autonomous" ODE = f(θ) You just integrate twice. 7. Apply this method to the case f(θ) = – ω2 sinθ Implications of initial conditions, set the two integration constants. Look up integral in GR, find a useful form for the "libration" solution. Statement of problem solution as sin(θ/2) = sin(θ0/2) sn(ωt; k= sin(θ0/2)), web confirms. Here, θ0 is max angle during swinging, sn is a certain Jacobi elliptic function, sine-like. 8. Taking the small oscillation limit of the general result, agrees with simple result 9. Examine the exact libration solution, plot for different values of θ0. Comment on the period as a function of θ0: T = 4K(k), k= sin(θ0/2) = v/(2ω) Here, K(k) is the complete elliptic integral of the first kind. 10. Looking for the "rotation" type solution Use a different GR integral, then realize this is just a "transformation" on the libration solution. The rotational solution is sin[θ/2] = sn(vt/2; k=2ω/v) Take the large v limit of this solution and get θ(t) = v t. 11. Comment on the arcsin function and its infinite number of "branches" for plotting How to make Maple plot the true arcsin(x). How to make Maple plot the true θ = arcsin(sinθ) as an ongoing straight line How to make Maple plot the rotation solution, example shown. 12. Plot of the phase space orbits picture, see both librations and rotations. 13. Hamilton-Jacobi theory applied to this problem. Write the HJE for function W(θ,E). Find that W is in fact an E -- second kind incomplete elliptic integral. Examine t + β1 = ∂W/∂α1, use known solution to verify this form is correct 14. Action-angle variables applied to pendulum Compute Jθ and wθ. Use the known general result that νθ = constant to solve the entire problem! Summary of "this little section", use HJ to derive the final result instead of confirm it. 15. What are the classical waves doing in this problem? Not much here. 16. Solve the same problem in quantum mechanics, get Mathieu functions, not Jacobi. The ODE is completely different, but problem is typical of a QM eigen problem. ____________________________________________________________________________ 1. Setup . Here is a picture: 2. Compute the Hamiltonian as a function of its correct arguments: When the pendulum is at rest straight down, I would like to have energy E = 0. Then E > 0 for all possible motions. This means we have V(θ) = mgh' where h' = l-h and picture shows h = lcosθ. Therefore V(θ) = mg[l-lcosθ] = mgl(1-cosθ) ok H = T + V = p2/2m + mgl(1-cosθ) L = T - V = p2/2m – mgl(1-cosθ) ok Since r = l = constant, there is really only one coordinate in this problem, which is θ. vθ = l L = ½ ml22 – mgl(1-cosθ) pθ = ∂L/∂ = ml2 T = ½ mv2 = ½ m (l)2 = ½ m (pθ/ml)2 = pθ2/(2ml2) H(qi, pi) = H(θ,pθ) = pθ2/(2ml2) + mgl(1-cosθ) = Hamiltonian. 3. Write down Hamilton's Equations of motion which are = ∂H/∂p and = – ∂H/∂q. Can remember the sign from free particle where H = p2/2m so = p/m = V. In our case then = ∂H/∂pθ = pθ/(ml2) θ = – ∂H/∂θ = – mglsinθ 4. Let's try to solve the right equation with θ as variable of interest: ml2 = – mgl sinθ or l = – gsinθ = – (g/l) sinθ = – ω2 sinθ where we have defined ω2 = g/l. For small angles, get = –ω2θ and get traditional result of θ = A sin(ωt + B). Note 1.27.17: Belendez et all agree with the above basic equation, which they write as so we are "OK to here". At least I have the correct starting equation. 5. Initial conditions. Suppose our start condition is θ=0 and = v. Then we have B = 0 and (0) = ωAcos(ωt) t=0 = ωA = v, so solution is , and ωt = arcsin[ ωθ/v ] and θ(t) = (v/ω)sin(ωt). So our assumption should be that (v/ω)<< 1 for this approximation to work. [ small osc only ] In general, what do I know about = – ω2sinx as a second order ODE ? Not obvious how to solve this thing. Let's ask Maple. Maple shows a mess. Web has help: 6. General method of solving this kind of non-linear ODE: [ It is an "autonomous" equation because the variable t does not appear explicitly ] d2x/dt2 = f(x(t)) Let w(x(t)) = dx(t)/dt. Then d2x/dt2 = dw/dt = dw/dx * dx/dt = w dw/dx. Then have w dw/dx = f(x) => wdw = f(x) dx => ½ w2 = f(x')dx' + C ok // note that w(x(t)) = dx(t)/dt = v(t), and this will be in our application ok Then we have (dx/dt)2 = 2 [ f(x')dx' + C ] which says dx/dt = ± = ± ≡ ± h(x) ok h(x) ≡ ok and then dx/h(x) = ± dt ok so get ± t = (1/) dx" [1/] + C' ≡ g(x) ok But this only tells you ± t = g(x). To find x(t), you need then to solve ±t – g(x) = 0. In other words, you need the "roots" x(t) of the equation ±t – g(x) = 0. This is what Maple is trying to say. ok 7. Now let's try this method for f(x) = –ω2sinx. We find // = – ω2 sinθ ok ½ w2 = f(x')dx' + C ok ½ w2 = f(x')dx' + C = -ω2sinx'dx' + C = ω2cosx + C ok So end up with ± t = (1/) dx" [1/] + C' ≡ g(x) ok Initial Conditions: Suppose we "launch" our pendulum from θ(0) = 0 but with some velocity (0) = v. We had at one point, where x = θ and w = : // unusual initial condition. ½ w2 = f(x')dx' + C = ω2cosx + C ok At t=0 the LHS is ½ v2 while the RHS is ω2+ C. Therefore, ½ v2 = ω2 + C C = ½ v2 – ω2 = ½ [ v2- 2ω2] ok The energy E of this pendulum is E = T + V = ½ m (l)2 + mgl(1-cosθ) and E(t=0) = ½ ml2v2 ok so ½ v2 = E/(ml2) and then C = E/(ml2) - ω2. ok Now assume that our initial v can kick the pendulum up to angle θ0. How high does the pendulum swing for libration solutions? Set E(0) = ½ ml2v2 = mgl(1-cosθ0) ok which says ½ v2 = (g/l) (1-cosθ0) and find that v2 = 2ω2(1-cosθ0) so that ok C = ½ [ v2- 2ω2] = ½ [2ω2(1-cosθ0)- 2ω2] = ω2(1- cosθ0 - 1) = -ω2cosθ0 ok Meanwhile, when t = 0, we know θ = 0 (which is x above), so we know C' = 0 [ assuming the lower endpoint of the integral below is 0 ] . Our solution is then ± t = (1/) dx" [1/] ok ±ωt = (1/) dθ [1/] where v2 = 2ω2(1-cosθ0) ok But (1-cosθ0) = 2 sin2(θ0/2) so v2 = 4ω2 sin2(θ0/2) and so v = 2ω sin (θ0/2) ok This integral fits the mold of GR page 154 #4 which says dθ [1/] = 2/ F[ x/2, ] a = -cosθ0 ok Here it is from GR7 and set b = 1: But = sin(θ0/2) and r = = / [sin(θ0/2) ] = csc(θ0/2) ok so our result is this dθ [1/] = 2/( sin(θ0/2) ) F[ x/2, ] ok or (1/) dθ [1/] = csc(θ0/2) F[ x/2, csc(θ0/2) ] ok But I don't like this form because the second argument k > 1 ! So let's try the other form that GR give in the same place which is this dθ [1/] = F(γ, /) γ = sin-1 ok = F(sin-1 , sin(θ0/2)) ok Our problem's solution is then this ±ωt = (1/) dθ [1/] = (1/) F(sin-1 , sin(θ0/2)) ok = F[ sin-1() , sin(θ0/2)] = F(φ,k) = u // This is the GR7 F ok Recall that if someone hands you F(φ,k) = u and asks you to find φ(u), the answer is sin(φ) = sn(u). ok We can restate this as φ = am(u) and sn(u) = sin(φ) How does GR7 support these claims? First, here is their definition of F(φ,k) from page 860: Then many pages down they define sn this way which I can write as u = F(sin-1(sn(u)),k) Now set u = sn-1x to get sn-1(x) = F(sin-1(x),k) Now let sin-1(x) = φ to get sn-1(sinφ) = F(φ,k) Now someone hands you an equation u = F(φ,k). You can conclude then that u = sn-1(sinφ) or sn(u) = sin(φ) What about the am(u) function? Here is what GR7 say So this says u = F(φ,k) = F(am(u),k) and φ = am(u). So both claims are OK Now notice that the second argument k =sin(θ0/2) has our desired property that k ≤ 1. Next, we know that φ = am(u) = sin-1() ok sn(u; k= sin(θ0/2)) = sinφ = ranges from 0 to 1 during motion ok Inserting ±ωt for u, we find that sin(θ/2) = sin(θ0/2) sn(±ωt; k= sin(θ0/2)) ok v2 = 2ω2(1-cosθ0) v = 2ω sin (θ0/2) ok θ(0)=0 and (0) = v ok ω2 = g/l ok Now which sign of ±ωt should I take ? Expansion in GR7 says that sn(x) is odd in its argument. So I guess either sign is fine, so I take the + sign. I have confirmation from the web that this is the correct solution of the plane pendulum, see last equation of this paragraph taken from http://www.mai.liu.se/~halun/complex/elliptic/ : New site: http://users.mai.liu.se/hanlu09/complex/elliptic/ I go to that site again and blow up and it does confirm my solution for θ(0) = 0. Comments: The above solution is for θ(0) = 0 and θ0 as the max. Same source above says Thus K(k) is a quarter period. So if you time shift my solution by a quarter period, but starting it at the max value, you will get sin(θ/2) = sin(θ0/2) sn(ωt; k= sin(θ0/2)) // my original solution θ/2 = sin-1 [ sin(θ0/2) sn(ωt; k= sin(θ0/2)) ] = sin-1 [ k sn(ωt; k) ] k = sin(θ0/2) Here is a plot: Replace t by t' = t - K(k). Then at t' = 0 you will have t = K(t) = one quarter cycle in! Then sin(θ/2) = sin(θ0/2) sn(ω[K(k)+t']; k= sin(θ0/2)) // my original solution 1. This source is a paper just on this subject and their solution is this, whereas my solution is sin(θ/2) = sin(θ0/2) sn(ωt; k= sin(θ0/2)) we seem to disagree on both arguments of the sn function! Let's find another source. 2. Here is another solution from another paper At least we agree on the second argument, but notice that vertical bar notation 8. Take small angle limit of the exact solution. Now what can we say about sn(ωt;k) for small k? Before going on, we need to learn something about the famous K and K', to wit, K = K(k) ≈ π/2 for small k GP page 905 K' = K'(k) = K() → ∞ as k → 0 because K(1) diverges. Therefore, in the small k regime, we know that the "nome" q = exp(-πK'/K) ≈ exp(-2K') = very small << 1. Therefore we can use the very first term in the GR page 911 expansion for snu which becomes sn(u; k) ≈ 2π/(kK) /(1-q) sin(πu/2K) = (4/k) sin(u). But in the k→ 0 limit, we need to study /k ≈ exp(- K')/k. This is a 0/0 limit. We know that K' = K'(k) = K(k') where k' = so we are near k' = 1. Use the expansion #2 shown page 905 GR [ item 8.113.3 ]. Swap k and k' to get K'(k) = K(k') = ln(4/k) + [...]k2 + .... near k = 0 Thus we approximate K(k') ≈ ln(4/k) for k near 0 and we then have limk→0 [/k] = limk→0 [exp{- K'}/k] = exp{-ln(4/k)}/k = (k/4)/k = 1/4 Thus our we get desired result that (4/k)→ 1 and then our limit is sn(u; k) ≈ sin(u) as k→0 // but of course I knew that so that our problem solution becomes sin(θ/2) = k sn(ωt; k= sin(θ0/2)) ≈ k sin(ωt) = sin(θ0/2) sin(ωt) and for the small angles as we are assuming, this says θ = θo sin(ωt) θ(0) = 0 // this seems correct !!! = ω θo cos(ωt) (0) = ωθ0 v2 = 2ω2(1-cosθ0) = ω2 sin2(θ0/2)/4 v = (ω/2) sin(θ0/2) ≈ ω θ0 and we finally recover our small angle result. [that only took me a whole day and a half! ] In the above, I wanted to see an expansion so I could in theory look at other terms in the approximation sn(u; k) ≈ sin(u). The series I used is sort of a Fourier harmonic expansion. 9. Examine the exact solution. We have confirmed that this is the complete pendulum solution: sin(θ/2) = sin(θ0/2) sn(ωt; k=sin(θ0/2)) where we use the notation sn(u,k) given on GR page 904. We can also express this as ωt = F[ sin-1() , sin(θ0/2)] = F(φ,k) Thus we could say that we first "relate" φ to our physical angle θ, and then sinφ = and F(φ,k) = ωt k = sin(θ0/2) But consider the first form above which says this: sin(θ/2) = sin(θ0/2) sn(ωt; k) We know that sn(ωt; k) is "sinusoidal" in a general sense and goes between the limits of +1 and -1. We know that for small k at least, as can approximate sn(ωt; k) ≈ sin(k'ωt) So what exactly is the above saying? It says θ goes back and forth between the two limits θ0 in some sine like fashion. θ(t) = 2 arcsin [sin(θ0/2) sn(ωt; k) ] This is the exact solution for any k. It is clear what happens for k < 1, meaning sin(θ0/2) < 1 meaning θ0/2 < π/2 meaning θo < π. We can have Maple plot this for some values of k. Well, here is a 3D plot where we vary both t and k plot3d(2*arcsin( sin(theta0/2) * JacobiSN(t,k)), t = 0..20, k=0..0.9, grid = [40,40]); For small k we get the near curves which are very sine like, but as k increases, the tops of the sines flatten off as the pendulum starts to "hold" at the top of its swing! Also, we clearly see the pendulum period increasing according to K(k) as plotted above. Here is a 2D plot for k = .99 plot(2*arcsin( sin(theta0/2) * JacobiSN(t,.99)), t = 0..20); Fascinating! Fact: The period of the librating pendulum is determined by ωT = 4*K(k) where k = sin(θ0/2) so the swinging frequency is given by ω' = 2π/T = [2π/4K(k)] ω. For example, if you have small amplitude, then K(k) ≈ K(0) = π/2 as for sin, ω' = ω. But suppose the pendulum goes half way up to θ0 = π/2. Then we have k = sin(π/4) = 1/ and then K(1/) = 1.854074677 so [2π/4K(k)] = .847 and our swinging frequency is about 15% slower than that of the small oscillations. By increasing θ0, we can get the frequency to be arbitrarily low (in theory). I just made a simple pendulum and measured 86% for the .847 above. Hard to get it to swing full up where I set it up, and there is air friction for the big swing since more velocity, but close enough!!! 10. Question: How can we make our solution show "rotation"? Our libration solution was this: θ(t) = 2 arcsin [k sn(ωt; k) ] v = 2ω sin (θ0/2) = 2ωk k = sin (θ0/2) sin(θ/2) = sin(θ0/2) sn(ωt; k) But let's go back and resolve this problem for general constant C. We find ½ 2 = ω2cosθ + C so = = ω = dθ/dt => dθ/ = ± ω dt and as usual the ± just means our solution can run either way in time, so take the +. Then ω t = dθ/ where I continue to assume that at t=0, θ = 0 which is fine. Repeating what we said above: This integral fits the mold of GR page 154 #4 which says dθ [1/] = 2/ F[ x/2, ] a =C/ω2 valid a > 1 hence C > ω2 AND valid 0 ≤ x ≤ π Now in the libration case, we had C = -ω2cosθ0 which was always negative, so we had a = -cosθ0 and then we had k = = cscθ0 and I complained that this was > 1 and so I used a different GR integral form. But now we have C = -ω2 + ½ v2 and we will be looking beyond the critical point which we can see is v = ω. Beyond this point C will be positive and the nature of the solution is going to change! The energy of this pendulum is always E = ½ ml2v2 since V = 0 at the bottom. We find that k = = 2ω/v which is rather convenient. Now our integral says valid only if C > ω2 which means -ω2 + ½ v2 > ω2 or ½ v2 > 2ω2 or v > 2ω. This means that k = 2ω/v ≤ 1 which is what we expect. So here is the solution I think I have found using this GR integral: ω t = dθ/ = 2/ F[ θ/2, ] = k F(θ/2,k) k = 2ω/v Thus we have for our rotational solution ωt/k = (v/2)t = F(θ/2,k) k = 2ω/v valid only for 0 ≤ θ ≤ π Notice that our libration solution was this: ωt = F[ sin-1() , sin(θ0/2)] = F(φ,k) Recall that if someone hands you F(φ,k) = u and asks you to find φ(u), the answer is sin(φ) = sn(u). We can restate this as φ = am(u) and sn(u) = sin(φ) sin(θ/2) = sin(θ0/2) sn(ωt; k) θ(t) = 2 arcsin [k sn(ωt; k) ] v = 2ω sin (θ0/2) = 2ωk k = sin (θ0/2) = v/(2ω) In our rotation situation we have this for our solution: ωt/k = (v/2)t = F(θ/2,k) k = (2ω/v) valid only for 0 ≤ θ ≤ π Recall that if someone hands you F(φ,k) = u and asks you to find φ(u), the answer is sin(φ) = sn(u). We can restate this as φ = am(u) and sn(u) = sin(φ) In our present situation we have φ = θ/2 and so we know that sin(θ/2) = sn((v/2)t). So here is our answer: sin[θ/2] = sn(vt/2; k=2ω/v) valid only for 0 ≤ θ ≤ π For very large v, we get k ≈ 0 so we get sin[θ(t)/2] = sin(vt/2) The conclusion is then θ(t)/2 = vt/2 => θ(t) = v t = v This is the limit really of g = 0 and ω = 0, a hard knock. The initial energy will be as before which was E = ½ ml2v2. All this energy goes into rotational motion, E = ½ ω'2 I where I is the moment of inertia which is ml2 so we get ½ ml2v2 = ½ ω'2 * ml2 so ω'2 = v2 and ω' = v. This is pretty obviously the right limit for a very hard kick. Furthermore, we know that the period of the above sn function is 4K(k). How much time goes by for this function sn(vt/2; k=2ω/v) to repeat itself ? Δ(vt/2) = 4K(k) = (v/2)Δt => Δt = (8/v)K(k) k = 2ω/v During this time Δt, the function sin(θ/2) makes a complete cycle as well, so this means that the variable θ/2 has gone 2π, which means θ has gone 4π. But this would be 2 cycles of θ or two "periods of the rotational motion", so we would say 2Trot = Δt = (8/v)K(k) and we then get this result for the period of our rotational motion Trot = 4K(k)/v k = 2ω/v and the frequency and angular frequency are νrot = 1/Trot = v/(4K(k)) ωrot = 2π νrot = v (π/2)/K(k) Let's check this in the high energy limit. In this case, k→0 and K = π/2 so we get ωrot = v = (0) which is what we expect. Review of the two solution regimes: [v = 2ω sin (θ0/2), E = ½ ml2v2 ] libration sin(θ/2) = sin(θ0/2) sn(ωt; k) k= sin(θ0/2) = v/(2ω) ωt = F[ sin-1() ,k] period = 4K(k) rotation sin(θ/2) = sn(vt/2; k) k = csc(θ0/2)= (2ω/v) ωt/k = (v/2)t = F(θ/2,k) period = 4K(k) One is appropriate for small v, the other for large v. The solutions are the same and have the same period when k = 1 which is the critical point between the two solution modes, corresponding to θ0 = π. The solutions are not the same! Perhaps they are related by a transformation of the sn function, however. Yes! Write the first as sin(θ/2) = k*sn(ωt; k) k = v/2ω Use this transformation rule sn(ku,1/k) = k sn(u,k) using u = ωt. The LHS is then sn(v/2ω*ωt, 2ω/(v)) = sn(vt/2, 2ω/v) which agrees! 11. Comment on the arcsin function. Maple provides a single-valued "function" like so: It cannot continue this because then we don't really have a function with a single value at each point. But we can define a more generalized function like so. We would like to define another function to plot at the same time as the above, which will be a continuation. Use these variable names x = sin(θ) θ = arcsin(x) For θ in range (-π/2,π/2) we can use the Maple function arcsin. n=0 For θ in range (π/2,3π/2) we need to instead use θ = - arcsin(x) + π n=1 For θ in range (3π/2,5π/2) we need to instead use θ = + arcsin(x) + 2π n=2 In general then , for θ in range (-π/2 + nπ, π/2 + nπ) θ = + (-1)n arcsin(x) + nπ So here is a set of functions θn(x) = (-1)n arcsin(x) + nπ n = any integer OK, here is how to make Maple plot a bunch of these functions together: plot([seq((-1)^n*arcsin(x)+n*Pi, n=-5..5)], x=-1..1, color = [seq(red,n=-5..5)]); We are just doing a plot of multiple functions and we enumerate them using the seq( fn, n=1..3) type construction. Something inside a [...] in Maple is called a list, whereas {...} denotes a set. Here we just want lists of things. If we don't specify colors as shown, the various curves come out different colors, but I want them all the same color. Here is what the above command yields: which is exactly what I want. Comment on plotting u = sin-1[ w=sinθ ] What we want of course is a plot that shows u = θ, but this is tricky to obtain because the sin-1 function acts on a domain of w = (-1,1) and a range of u = (-π/2,π/2). We want it to work out of this range, so we have to do a lot of work. Here is my solution to this problem to plot two adjacent regions: > K := Pi/2: > q := Pi/2: > f1 := arcsin(sin(t)): > f2 := arcsin(sin(t-2*K))+2*q: > plot(piecewise(t>=-K and t<K,f1,t>K and t<=3*K, f2, 0),t=-K..3*K); I have to make a separate piecewise function for each of the two ranges of t. Quantity K is the quarter period of the periodic function which is the argument of the arcsin, while q is the quarter period of the sin function. Using this idea, we can examine our "rotation problem" solution sin[θ/2] = sn(vt/2; k=2ω/v) θ = 2*arcsin[sn(vt/2; k=2ω/v)] I will select v = 2 to make the scale simple, so we get θ = 2*arcsin[sn(t; k= ω] We know that K(k) is the quarter period of sn, and we have an extra factor of 2* out front, so here is a way to get Maple to plot this thing: > k := .999: > K := EllipticK(k); K := 4.495596396 > f1 := 2*arcsin(JacobiSN(t,k)): > f2 := 2*( arcsin(JacobiSN(t-2*K,k))+2*q): > q := Pi/2: > plot(piecewise(t>=-K and t<K,f1,t>K and t<=3*K, f2, 0),t=-K..3*K); So finally we are getting the solution θ(t) for our "rotation problem". We have chosen k = .999 to emphasize slow motion at the top of the pendulum. It just barely reaches the top and then falls ahead to continue for another loop, so the angle θ as presented here just increases all the time. The thing never swings backwards of course (which would cause me more plotting problems). So in reality, what we do is compute our solution in the range -K to K, then we just stamp that picture down again to the right and shifted up π. So I think this is a complete understanding of our rotation solution which is θ(t) = 2*arcsin[sn(vt/2; k=2ω/v)] 12. Phase Space. I am happy now with the solution I have found. I would now like to generate that famous phase space picture I have seen in several places. From above we had ½ 2 = ω2cosθ + C C = E/(ml2) - ω2. ω2 = g/l pθ = ml2 From this I think I find pθ2 = 2ml2E - (2ml2ω)2 sin2(θ/2) = [2m2gl3] ( E' – 2sin2θ/2 ) = [2m2gl3] ( {E'– 1}+cosθ ) where E is the pendulum energy. This of course has the form pθ2 = E' - 2 sin2(θ/2) This means that, horizontally, what you see in the range -π,π for θ will just repeat in windows to the left and right. Because there are ± square roots, the picture is reflected in the θ axis, so I will only plot it above this axis. If we just name E' to be E, and set A = 1, we get pθ = + If E < 1, we get one kind of curve, and it E > 1 we get another. In this plot, I will run E in steps of 0.1 from .1 up to 3.0: plot([seq((E/10 - 2*sin(theta/2)^2)^(1/2),E = 1..30)],theta = -Pi..Pi,color =[seq(red,E = 1..30)], numpoints=1000); So this is the famous picture! See Marion page 164 or Goldstein page 290. My 10th curve from the inside has E=1 and you see it is the critical curve as shown better in Marion. As you run the pendulum with different E, you will loop on the different phase space orbits shown, or you will do rotation. The larger orbits have longer periods as we know already, and as we will probably find soon from Hamilton-Jacobi theory if we ever get to it with this pendulum which has consumed multiple days of time so far. 13. Hamilton-Jacobi theory applied to this problem. From Section 2 above we know that H(qi, pi) = H(θ,pθ) = pθ2/(2ml2) + mgl(1-cosθ) = Hamiltonian. pθ = ml2 so the HJE says this ( using the fixed energy E = α1 situation) H(θ, pθ= ∂W/∂θ,pθ = αθ) = α1 = E which we write out as (∂W/∂θ)2/(2ml2) + mgl(1-cosθ) = E On scratch paper I fiddle this to get dW = A dθ A = ml2ω E' = E/mgl W = A ∫ dθ ∂W/∂α1 = ( 1/ω) ∫dθ 1/ " second integral" To do the second integral we go to our same GR integrals on p 154 dθ [1/] = F(sin-1) , sin(θ0/2)) // libration form = csc(θ0/2) F[ θ/2, csc(θ0/2) ] // rotation form where for now θ0 represents the unknown constant of integration. Then according to page 287 we should have t + β1 = ∂W/∂α1 = ( 1/ω) F(sin-1) , sin(θ0/2)) = (1/ω) F(sin-1) , sin(θ0/2)) = (1/ω) csc(θ0/2) F[ θ/2, csc(θ0/2) ] k = csc(θ0/2) = 2ω/v which we recognize as our libration solution from above, with β1 = 0 for our initial conditions. Just for fun, what is the W function itself? W = A ∫ dθ A = ml2ω We can use GR page 156 with a = E'-1 and b = 1 to get and a+b = E' W(θ, E) = A 2 E(θ/2, r= = csc(θ0/2) ) // see below on last part = 2 E(θ/2, csc(θ0/2)) = B E(θ/2, csc(θ0/2)) B = 2 which is a second-kind elliptic integral function. Now we should be able to differentiate this to get our previous result. In this set of equations, we shall use k = csc(θ0/2): E = ½ ml2v2 where v = 2ω sin (θ0/2) => E = 2mgl sin2(θ0/2) so that sin2(θ0/2) = E/(2mgl) = E'/2 = k-2 => 2mgl = Ek2 and E = 2mglk-2 => 0 = E2kdk + dEk2 => dk/dE = -k/2E and k = csc(θ0/2) = Then we get ∂W(θ, E)/∂E = B[ ∂E(θ/2, k)/∂k * ∂k/∂E + (1/(2)) E(θ/2,k) ] k = csc(θ0/2) For the derivative, we consult GR p age 907 8.123.3 which says ∂E/∂k = (E-F)/k so we get ∂W(θ, E)/∂E = B[ ∂E(θ/2, k)/∂k * ∂k/∂E + (1/(2)) E(θ/2,k) ] = B[ [(E(θ/2,k)-F(θ/2,k))/k] * [-k/2E] + (1/(2)) E(θ/2,k) ] = B[ - (1/(2)[(E(θ/2,k)-F(θ/2,k))] + (1/(2)) E(θ/2,k)] = B(1/(2) F(θ/2,k) = (k/ω) F(θ/2, k = csc(θ0/2)) which agrees exactly with our previous result. Just more exercise in use of elliptic formulas. 14. Action-angle variables applied to this problem. Goldstein 9-32 page 290 gives the left equation here, and from above I keep going: pθ = ∂W/∂θ = A A = ml2ω E' = E/mgl But inside the square root we have E' - (1-cosθ) = E' - 2sin2(θ/2) => pθ = ∂W/∂θ = A We have seen in Section 12 that this provides the characteristic phase space graphic. In our current notation for A and E' (different from what I used earlier), the critical energy is E' = 2. Notice that A2E' = 2ml2E and reproduces our first term in our previous pθ2 expression. Similarly, notice that 2A2 = (2ml2ω)2 which reproduces our second coefficient, so our result here is exactly the same as the previous result in Section 12 above. Now let's compute the "action" variable Jθ [ remember: this is the area inside the phase space loop ] Jθ = dθ (∂W/∂θ) = dθ [A ] A = ml2ω E' = E/mgl We know the curve is symmetric, so the area is 2 times the integral 0 to π (see G p 290) Jθ = 2A dθ = 4A E(π/2, k ) = 4A E(k) k = = csc(θ0/2) where E is the complete second kind elliptic function. Next we want to find the "angle" variable which is w ≡ ∂W/∂Jθ = ∂W/∂E' * ∂E'/∂Jθ = ∂W/∂E' / ∂Jθ /∂E' Here is what we know so far: W(θ, E) = 2A E(θ/2,k ) k = = csc(θ0/2) A = ml2ω E' = E/mgl Jθ = 4A E(k) = 2W(π/2, E) So we can compute the required derivatives as follows: ∂W/∂E' = 2A [ ∂E(θ/2,k )/∂k * ∂k/dE' + (1/2) E(θ/2,k ) ] But we know that k = so k2E' = 2 so ∂k/dE' = - k/(2E'). Also we look up GR page 197 which tells us that ∂E(θ/2,k )/∂k = [(E-F)/k] so we get ∂W/∂E' = 2A [[(E-F)/k] * -k/(2E') + (1/2) E(θ/2,k ) ] = 2A(1/2) [[- (E-F)] + E(θ/2,k ) ] = (A/) F(θ/2,k) = (A/) F(θ/2,k) The other derivative is the same as what we just did except it has θ = π and factor 2, so get ∂Jθ /∂E' = (2A/) F(π/2,k) = (2A/) K(k) Thus it appears that our angle variable is this: [ all this time k = = csc(θ0/2) ] wθ = (A/) F(θ/2,k) / (2A/) K(k) = ½ F(θ/2,k)/ K(k) = ½ F(θ/2,k) / F(π/2,k) which is mildly interesting. But now it is supposed to be true that νθ = θ In our rotation solution above, we found that ωt/k = (v/2)t = F(θ/2,k) k = (2ω/v) = csc(θ0/2) so let's now use this in our last result wθ = ½ F(θ/2,k)/ K(k) = ½ (v/2)t /K(k) = (v/4) t /K(k) Then we find that νθ = θ = (v/4)/K(k) = νrot so we are consistent with our previous calculation in the rotation section above. Summary of this little section: The W function for this problem is a second-kind elliptic function, W(θ,E) = A dθ = 2A E(θ/2,k ) where A = ml2ω E' = E/mgl k = k = = csc(θ0/2) The k value here is appropriate for the "rotation" regime of the pendulum where E' > 2 so the integrand is always real and k is always < 1. You could express this in a form appropriate for libration by using this identity which I derived from the GR page 908 transformation table, E(φ,k) = k F(ksinφ, 1/k) + (1-k2)/k * F(ksinφ, 1/k) with φ = θ/2 but the other form is much simpler so we stick with it. We can compute the action and the angle variables and we find that Jθ = 2A = 4A E(π/2, k ) = 4A E(k) k = wθ = ½ F(θ/2,k)/ K(k) = ½ F(θ/2,k) / F(π/2,k) Remember that we did a canonical transformation from θ, pθ to wθ, Jθ which is a particular case of doing the HJ transformation with a W. Notice that Jθ is a constant since E' is a constant. Instead of what we did above, at this point we could use the HJ knowledge that νθ = θ = constant, and then we would say [ this is going to lead us to a solution of the problem! ] νθ = θ = (1/2K(k) ) ∂F(φ,k)/∂φ |φ=θ/2 * ∂φ/∂t|φ=θ/2 But from the definition of the integral F, we know that ∂F(φ,k)/∂φ |φ=θ/2 = 1/Δφ|φ=θ/2 = 1/( ) ∂φ/∂t|φ=θ/2 = ½ (t) to get that νθ = θ = (1/2K(k) ) 1/ * ½ (t) which we could then solve for (t) to get (t) = 4 νθ K(k) => dθ 1/ = 4 νθ K(k) t Change variables to φ = θ/2 in the integral so it becomes 2 dφ 1/ = 2 F(θ/2,k) GR page 904 Then we have as our problem's solution [2 νθ K(k)] t = F(θ/2,k) Looking up 5 equations, we could look at t=0 where we have θ=0 so νθ = (1/2K(k) ) * 1 * ½ (0) = v/(4K(k)) which tells us that [2 νθ K(k)] = v/2 so our problem solution is then vt/2 = F(θ/2,k) which we can compare to our known rotational solution which is ωt/k = (v/2)t = F(θ/2,k) k = (2ω/v) = csc(θ0/2) which is the same! So in this way, we have used the action-angle variables version of the Hamilton-Jacobi method to find the solution to our pendulum problem. We never had a non-linear differential equation which we had to solve. This is really a different method of solution! Comment on the rotational case. At the top of the rotation we have θ = π. Our motion equation says (t) = 4 νθ K(k) = v Velocity at the top is then top = v k' . Ttop = ½ m l2 2 = ½ m l2 v2 (1-k2) k = 2ω/v Vtop = 2mgl Tbottom = ½ m l2 v2 Vtop = 0 Ttop/Vtop = (1-k2)/k2 = cot2α if k = sinα Here is a possible interpretation: top/bot = k' = cosα Conclusion: during rotation mode, I can find no "angle" in the picture which has any correspondence with k. I was wondering about my idea of k = csc(θ0/2). This is just not something useful. 15. What are the classical waves doing in this problem? The action S is following surfaces such that W(θ,z,E) = Et as t increases. So we need to solve this for φ to see what the surfaces look like. W(θ,E) = 2A E(θ/2,k ) But what is the inverse of the E function? This is not obvious from a quick look, it is not sn or cn or dn, but of course it is "well defined" and here then are our surfaces: 2A E(θ/2,k ) = Et As t increases, we get some θ = f(t,E). This is a flat surface at angle θ(t) which stretches flat into the z direction. As t increases, these surfaces move sort of perpendicular to the pendulum mass, just as you would expect. Show them with your right hand. We could compute the wave velocity u, and we know the wavelength is λ = h/mv. The waves don't keep up with the mass. Easier to study this for a free particle as done in Goldstein, I just thought I would mention it here. 16. What is the solution to the quantum plane pendulum problem? It would have bound states for low E. Not clear how you would tie a string onto a non-localized mass, so I don't know how you would even set this problem up. I guess the string follows the same ψ as the mass. Setup. By setting up in terms only of angle θ as in the classical case, we are saying "there is a string". The TI SE Says this: (Laplacian in cylindricals in effect has that 1/r2 out frong. (1/l2) ∂2ψ(θ)/∂θ2 + A-1 (E-V)ψ(θ) = 0 A-1 = 2m/2 V(θ) = mgl(1-cosθ) (A/l2) ∂2ψ(θ)/∂θ2 + (E-V)ψ(θ) = 0 Let θ = x for the moment and write this as - (A/l2)ψxx + mgl(1-cosx)ψ = Eψ E' = E/mgl Define B = A/mgl3 = 2/(2m2gl3) and use same E' = E/mgl as in the classical stuff above, so then -Bψxx + (1-cosx)ψ = E'ψ Notice that B and E' are dimensionless. The angular momentum of the pendulum going at its basic ω rate would be L = rp = rmv = lmlω = l2 m so that L2 = m2l3g which is what appears in the denominator of B, and of course also has dimensions of L. This is linear in ψ, so a linear ODE. Compare to θtt= -ω2sinθ we had before which was "auto". Change from θ to ψ and t to x to get ψxx = -ω2 sinψ which is "autonomous". So now we have a completely different equation. Write our new equation as -Bψθθ = [ (E'-1) + cosθ] ψ Maybe this helps to write it using x = t: -B = [ (E'-1) + cos(t)] ψ(t) = f(t)ψ So how do you go about solving a "new" ODE that you are not familiar with? Could always do power series then try to find the series. Then end up with mysterious recursion relation. How about our old friend Maple! Well, Maple cannot do it. You can turn on the info level and see it try. Here is the Maple output: So here we are again with B + [ (E'-1) + cos(t)] ψ(t) = 0 B = A/mgl3 = 2/(2m2gl3) E' = E/mgl I have no idea what to do! I look up "SE" and pendulum web, was led to this fact: That was a good lead, GR has a whole section on this guy starting page 991. Page 993 shows "negative q" which is our case. In fact, here is what GR show ψtt + (a + 2q cos(2t))ψ = 0 If I define t' = 2t so that ψt = 2ψt' etc, then we have 4 ψt't'(t') + (a + 2q cos(t'))ψ(t') = 0 Define a' and q' as a/4 and q/4 so we then have ψt't'(t') + (a' + 2q' cos(t'))ψ(t') = 0 then set t' back to t, arbitrary ψtt(t) + (a' + 2q' cos(t))ψ(t) = 0 But this is our pendulum SE where we have a' = (E' - 1)/B and 2q' = 1/B. Therefore, here is our claim. "The pendulum SE B ψtt + [ (E'-1) + cos(t)] ψ(t) = 0 is really the Mathieu equation, ψtt + (a + 2q cos(2t))ψ = 0 with these parameters a = (E' - 1)/B and 2q = 1/B and B = A/mgl3 = 2/(2m2gl3) E' = E/mgl so we have dimensionless q and a as follows, q = m2gl3/ 2 = 1/(2B) a = 2q(E' - 1) end of claim. " Then on GR page 993 we see solution functions called ce and se. It is noted that the periodic solutions exist for certain eigenvalues of a. Each solution is presented as a Fourier sum of cosines or sines, hence the names ce and se. There is a label 2n and 2n+1 which is the usual kind of thing. Since I have the TISE with t really θ, periodic would mean periodic in θ, not time. Meanwhile, here is what the potential V(θ) looks like: I am not used to this kind of potential. Since it is periodic, I guess you would expect the Mathieu solution to be periodic. But in our physical system, we really have this But we can think of θ going up forever in the rotation solution, and then the first does apply! The second applies for the libration solutions of lower energy. That would be bound states in this well. The rotation solution I guess would be a wave packet perhaps moving to the right through the wells! Very good. A&S also have a section on Mathieu functions page 722. The parameters are a and q, and page 724 shows a plot of the eigenvalues ar as functions of q for the various solutions "r". " The French mathematician E.L. Mathieu investigated this equation in 1868 while seeking a description of the vibrations of an elliptical membrane (drumheads)." I finally found something on the web that at least addresses this problem, a PDF, from which I quote In the above, AMF refers to "angular Mathieu functions", and the text continues: I think you need to have eigenvalues in order for the solutions to have the right periodicity to match the π periodicity of the potential in this problem. That sets the eigenenergies, and then you can look at the lowest states, here are some graphs from the same paper Each graph is shown for several values of parameter q. The angle η = θ/2 is zero at the TOP of the pendulum (my classical analysis had it at the bottom.). So the bottom is θ = π or η = π/2. The top left graph shows the ground state ce0 which has probability peaking at the bottom position. Then the next state up has none at the bottom, and some on the two sides. As usual, the number of nodes increases with the index. The se functions are also solutions, and I don't understand why things are doubled. The se are odd if you go to the - angle region and the ce are even, but not sure we care about that region. I suppose all these solutions exist and are different. The AS graph shows a0 and b1 are very close. "for cem the eigenvalues are usually denoted as am(q), whereas for sem they are represented as bm(q)." Summary: The quantum plane pendulum can be studied using the TISE with 2 in cylindrical coordinates which just means 1D with an extra 1/l2. The differential equation which results is completely different from that which arises in the classical plane pendulum problem: θtt = -ω2 sinθ // classical solutions are Elliptic functions ψtt + (a' + 2q' cos(t))ψ = 0 // quantum solutions are Mathieu functions As expected, the TDSE has eigenstates with the "usual shapes" in the θ space. We would then construct the classical solutions by taking wave packets formed from these states. At low energy, these packets would bounce back and forth, emulating the classical librating pendulum. At higher energy, these packets would skid along the potential to ever-increasing θ, emulating the rotational solutions. So yes, this is a pretty good problem to compare classical and quantum solutions. It corresponds to a real physical situation, unlike particles in boxes.