(*^ ::[ frontEndVersion = "NeXT Mathematica Notebook Front End Version 2.1"; next21StandardFontEncoding; paletteColors = 128; automaticGrouping; currentKernel; fontset = title, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e8, 24, "Times"; ; fontset = subtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, L1, e6, 18, "Times"; ; fontset = subsubtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, L1, e6, 14, "Times"; ; fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, L1, a20, 18, "Times"; ; fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, L1, a15, 14, "Times"; ; fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, L1, a12, 12, "Times"; ; fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 10, "Times"; ; fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L1, 12, "Courier"; ; fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5, 12, "Courier"; ; fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L1, 12, "Courier"; ; fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L1, 12, "Courier"; ; fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L1, 12, "Courier"; ; fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287, L1, 12, "Courier"; ; fontset = name, inactive, noPageBreakInGroup, nohscroll, preserveAspect, M7, italic, B65535, L1, 10, "Times"; ; fontset = header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, L1, 12, "Times"; ; fontset = Left Header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic, L1, 12, "Times"; ; fontset = footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, italic, L1, 12, "Times"; ; fontset = Left Footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, italic, L1, 12, "Times"; ; fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12, "Courier"; ; fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, L1, 12; ] :[font = title; inactive; preserveAspect; startGroup; ] Lagrangian Equations :[font = text; inactive; preserveAspect; ] Jeff Adams, Wolfram Research Inc., October 1992 :[font = subsection; inactive; preserveAspect; startGroup; ] Introduction :[font = text; inactive; preserveAspect; endGroup; ] This notebook illustrates one way of using Mathematica to solve some advanced theoretical mechanics problems using Lagrangian Mechanics. The example problems are taken from Analytical Mechanics, Grant R. Fowles, Fourth Edition, Saunders College Publishing, 1986. ;[s] 5:0,0;43,1;54,2;174,3;194,4;264,-1; 5:1,11,8,Times,0,12,0,0,0;1,10,8,Times,2,12,0,0,0;1,11,8,Times,0,12,0,0,0;1,10,8,Times,2,12,0,0,0;1,11,8,Times,0,12,0,0,0; :[font = subsection; inactive; preserveAspect; startGroup; ] Function :[font = text; inactive; preserveAspect; ] The function LagrangianEquations[T, V, Q, genCoords] takes four arguments: the kinetic energy, T, of the system, the potential energy, V, of the system, the non-conservative forces of the system, Q, and a list of the generalized coordinates of the system. The function returns a list of the lagrangian equations for this system. Note, it is very important to explicitly declare all dependencies of T,V, and Q on the generalized coordinates, as well as the explicit time dependence of the generalized coordinates. :[font = input; preserveAspect; endGroup; ] LagrangianEquations[T_,V_,Q_:0,genCoords_List] := Module[{L=T-V}, (D[D[L,D[#,t]],t] == Q + D[L,#])& /@ genCoords ] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; ] Example 1 : Harmonic Oscillator :[font = input; preserveAspect; ] genCoords = {x[t]}; T = 1/2 m x'[t]^2; V = 1/2 K x[t]^2; Q = 0; :[font = input; preserveAspect; startGroup; ] LagrangianEquations[T, V, Q, genCoords] :[font = output; output; inactive; preserveAspect; endGroup; ] {m*Derivative[2][x][t] == -(K*x[t])} ;[o] {m x''[t] == -(K x[t])} :[font = input; preserveAspect; startGroup; ] DSolve[%,genCoords,t] :[font = output; output; inactive; preserveAspect; endGroup; endGroup; ] {{x[t] -> E^((-I*K^(1/2)*t)/m^(1/2))*C[1] + E^((I*K^(1/2)*t)/m^(1/2))*C[2]}} ;[o] (-I Sqrt[K] t)/Sqrt[m] {{x[t] -> E C[1] + (I Sqrt[K] t)/Sqrt[m] E C[2]}} :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; ] Example 2 : Harmonic Oscillator with velocity damping :[font = input; preserveAspect; ] genCoords = {x[t]}; T = 1/2 m x'[t]^2; V = 1/2 K x[t]^2; Q = -c x'[t]; :[font = input; preserveAspect; startGroup; ] LagrangianEquations[T, V, Q, genCoords] :[font = output; output; inactive; preserveAspect; endGroup; ] {m*Derivative[2][x][t] == -(K*x[t]) - c*Derivative[1][x][t]} ;[o] {m x''[t] == -(K x[t]) - c x'[t]} :[font = input; preserveAspect; startGroup; ] DSolve[%,genCoords,t] :[font = output; output; inactive; preserveAspect; endGroup; endGroup; ] {{x[t] -> C[1]/E^(((c + (c^2 - 4*K*m)^(1/2))*t)/(2*m)) + E^(((-c + (c^2 - 4*K*m)^(1/2))*t)/(2*m))*C[2]}} ;[o] C[1] {{x[t] -> --------------------------------- + 2 ((c + Sqrt[c - 4 K m]) t)/(2 m) E 2 ((-c + Sqrt[c - 4 K m]) t)/(2 m) E C[2]}} :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; ] Example 3 : Single Particle in a Central Field :[font = text; inactive; preserveAspect; ] Using polar coordinates :[font = input; preserveAspect; ] genCoords = {r[t],theta[t]}; T = 1/2 m ( r'[t]^2 + r[t]^2 theta'[t]^2 ); V = v[r[t]]; Q = 0; :[font = input; preserveAspect; startGroup; ] LagrangianEquations[T, V, Q, genCoords] :[font = output; output; inactive; preserveAspect; endGroup; endGroup; ] {m*Derivative[2][r][t] == m*r[t]*Derivative[1][theta][t]^2 - Derivative[1][v][r[t]], 2*m*r[t]*Derivative[1][r][t]*Derivative[1][theta][t] + m*r[t]^2*Derivative[2][theta][t] == 0} ;[o] 2 {m r''[t] == m r[t] theta'[t] - v'[r[t]], 2 2 m r[t] r'[t] theta'[t] + m r[t] theta''[t] == 0} :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; ] Example 4 : Atwood's Machine :[font = text; inactive; preserveAspect; ] A mechanical system known as Atwood's machine consists of two weights of mass m1 and m2 , respectively, connected by a light inextensible cord of length l which passes over a pulley. The system has one degree of freedom. Let x represent the configuration of the system, where x is the vertical distance from the pulley to m1. I is the moment of inertia of the pulley with radius a. :[font = input; preserveAspect; ] genCoords = {x[t]}; T = 1/2 m1 x'[t]^2 + 1/2 m2 x'[t]^2 + 1/2 I x'[t]^2/ a^2; V = -m1 q x[t] - m2 q (l - x[t]); Q = 0; :[font = input; preserveAspect; startGroup; ] LagrangianEquations[T, V, Q, genCoords] :[font = output; output; inactive; preserveAspect; endGroup; ] {(I*Derivative[2][x][t])/a^2 + m1*Derivative[2][x][t] + m2*Derivative[2][x][t] == m1*q - m2*q} ;[o] I x''[t] {-------- + m1 x''[t] + m2 x''[t] == m1 q - m2 q} 2 a :[font = input; preserveAspect; startGroup; ] Solve[%, x''[t]] :[font = output; output; inactive; preserveAspect; endGroup; endGroup; ] {{Derivative[2][x][t] -> -((-(a^2*m1*q) + a^2*m2*q)/(I + a^2*m1 + a^2*m2))}} ;[o] 2 2 -(a m1 q) + a m2 q {{x''[t] -> -(--------------------)}} 2 2 I + a m1 + a m2 :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; ] Example 5 : Double Atwood Machine :[font = text; inactive; preserveAspect; ] Replace one of the weights in the atwood machine by another pulley supporting two weights connected by another cord. The system now has two degrees of freedom. We shall specify the coordinates x1 and x2. For simplicity we shall neglect the masses of the pulleys in this case. :[font = input; preserveAspect; ] genCoords = {x1[t],x2[t]}; T = 1/2 m1 x1'[t]^2 + 1/2 m2 (-x1'[t] + x2'[t])^2 + 1/2 m3 ( -x1'[t] - x2'[t])^2; V = -m1 q x1[t] - m2 q (l1 - x1[t] + x2[t]) - m3 g (l1 - x1[t] + l2 - x2[t]); Q = 0; :[font = input; preserveAspect; startGroup; ] LagrangianEquations[T, V, Q, genCoords] :[font = output; output; inactive; preserveAspect; endGroup; ] {m1*Derivative[2][x1][t] - m3*(-Derivative[2][x1][t] - Derivative[2][x2][t]) - m2*(-Derivative[2][x1][t] + Derivative[2][x2][t]) == -(g*m3) + m1*q - m2*q, -(m3*(-Derivative[2][x1][t] - Derivative[2][x2][t])) + m2*(-Derivative[2][x1][t] + Derivative[2][x2][t]) == -(g*m3) + m2*q} ;[o] {m1 x1''[t] - m3 (-x1''[t] - x2''[t]) - m2 (-x1''[t] + x2''[t]) == -(g m3) + m1 q - m2 q, -(m3 (-x1''[t] - x2''[t])) + m2 (-x1''[t] + x2''[t]) == -(g m3) + m2 q} :[font = input; preserveAspect; startGroup; ] Simplify[ Solve[%, {x1''[t],x2''[t]}] ] :[font = output; output; inactive; preserveAspect; endGroup; endGroup; ] {{Derivative[2][x1][t] -> (-2*g*m2*m3 + m1*m2*q + m1*m3*q - 2*m2*m3*q)/ (m1*m2 + m1*m3 + 4*m2*m3), Derivative[2][x2][t] -> (-(g*m1*m3) - 2*g*m2*m3 + 2*m1*m2*q - m1*m3*q + 2*m2*m3*q)/(m1*m2 + m1*m3 + 4*m2*m3)}} ;[o] -2 g m2 m3 + m1 m2 q + m1 m3 q - 2 m2 m3 q {{x1''[t] -> ------------------------------------------, m1 m2 + m1 m3 + 4 m2 m3 x2''[t] -> -(g m1 m3) - 2 g m2 m3 + 2 m1 m2 q - m1 m3 q + 2 m2 m3 q --------------------------------------------------------}} m1 m2 + m1 m3 + 4 m2 m3 :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; ] Example 6 : Particle Sliding on a Movable Inclined Plane :[font = text; inactive; preserveAspect; ] Consider a particle sliding on a smooth inclined plane which itself, is free to slide on a smooth horizontal surface. There are two degrees of freedom x1 and x2, the horizontal displacement of the plane from some reference point and the displacement of the particle from some reference point on the plane. Using the law of cosines for the velocity of the particle, we have: :[font = input; preserveAspect; ] genCoords = {x1[t],x2[t]}; T = 1/2 m (x1'[t]^2 + x2'[t]^2+2 x1'[t] x2'[t] Cos[theta]) + 1/2 M x1'[t]^2; V = - m q x2[t] Sin[theta] + cons; Q = 0; :[font = input; preserveAspect; startGroup; ] LagrangianEquations[T, V, Q, genCoords] :[font = output; output; inactive; preserveAspect; endGroup; ] {M*Derivative[2][x1][t] + (m*(2*Derivative[2][x1][t] + 2*Cos[theta]*Derivative[2][x2][t]))/2 == 0, (m*(2*Cos[theta]*Derivative[2][x1][t] + 2*Derivative[2][x2][t]))/2 == m*q*Sin[theta]} ;[o] m (2 x1''[t] + 2 Cos[theta] x2''[t]) {M x1''[t] + ------------------------------------ == 0, 2 m (2 Cos[theta] x1''[t] + 2 x2''[t]) ------------------------------------ == m q Sin[theta]} 2 :[font = input; preserveAspect; startGroup; ] Simplify[ Solve[%, {x1''[t],x2''[t]}] ] :[font = output; output; inactive; preserveAspect; endGroup; endGroup; ] {{Derivative[2][x1][t] -> (m*q*Sin[2*theta])/(-m - 2*M + m*Cos[2*theta]), Derivative[2][x2][t] -> (2*(m + M)*q*Sin[theta])/(m + 2*M - m*Cos[2*theta])}} ;[o] m q Sin[2 theta] {{x1''[t] -> -------------------------, -m - 2 M + m Cos[2 theta] 2 (m + M) q Sin[theta] x2''[t] -> ------------------------}} m + 2 M - m Cos[2 theta] :[font = subsection; inactive; Cclosed; preserveAspect; startGroup; ] Example 7 : Spherical pendulum, or bar of soap in a bowl :[font = input; preserveAspect; ] genCoords = {theta[t],phi[t]}; T = 1/2 m l^2 (theta'[t]^2 + Sin[theta[t]]^2 phi'[t]^2); V = m g l (1 - Cos[theta[t]]); Q = 0; :[font = input; preserveAspect; startGroup; ] LagrangianEquations[T, V, Q, genCoords] :[font = output; output; inactive; preserveAspect; endGroup; ] {l^2*m*Derivative[2][theta][t] == -(g*l*m*Sin[theta[t]]) + l^2*m*Cos[theta[t]]*Sin[theta[t]]*Derivative[1][phi][t]^2\ , 2*l^2*m*Cos[theta[t]]*Sin[theta[t]]* Derivative[1][phi][t]*Derivative[1][theta][t] + l^2*m*Sin[theta[t]]^2*Derivative[2][phi][t] == 0} ;[o] 2 {l m theta''[t] == -(g l m Sin[theta[t]]) + 2 2 l m Cos[theta[t]] Sin[theta[t]] phi'[t] , 2 2 l m Cos[theta[t]] Sin[theta[t]] phi'[t] theta'[t] + 2 2 l m Sin[theta[t]] phi''[t] == 0} :[font = text; inactive; preserveAspect; ] Given that phi is an ignorable coordinate then m l^2 Sin[theta[t]]^2 phi'[t] = constant :[font = input; preserveAspect; startGroup; ] % /. phi'[t] -> h1/ (m l^2 Sin[theta[t]]^2 ) :[font = output; output; inactive; preserveAspect; endGroup; ] {l^2*m*Derivative[2][theta][t] == (h1^2*Cot[theta[t]]*Csc[theta[t]]^2)/(l^2*m) - g*l*m*Sin[theta[t]], 2*h1*Cot[theta[t]]*Derivative[1][theta][t] + l^2*m*Sin[theta[t]]^2*Derivative[2][phi][t] == 0} ;[o] 2 2 2 h1 Cot[theta[t]] Csc[theta[t]] {l m theta''[t] == -------------------------------- - 2 l m g l m Sin[theta[t]], 2 h1 Cot[theta[t]] theta'[t] + 2 2 l m Sin[theta[t]] phi''[t] == 0} :[font = input; preserveAspect; startGroup; ] Solve[%,theta''[t]] :[font = output; output; inactive; preserveAspect; endGroup; ] {{Derivative[2][theta][t] -> ((h1^2*Cot[theta[t]]*Csc[theta[t]]^2)/(l^2*m) - g*l*m*Sin[theta[t]])/(l^2*m)}} ;[o] {{theta''[t] -> 2 2 h1 Cot[theta[t]] Csc[theta[t]] -------------------------------- - g l m Sin[theta[t]] 2 l m ------------------------------------------------------}} 2 l m :[font = input; preserveAspect; startGroup; ] thetaEqu = Simplify[ % /. {h1 -> h l^2 m}] :[font = output; output; inactive; preserveAspect; endGroup; ] {{Derivative[2][theta][t] -> h^2*Cot[theta[t]]*Csc[theta[t]]^2 - (g*Sin[theta[t]])/l}} ;[o] {{theta''[t] -> 2 2 g Sin[theta[t]] h Cot[theta[t]] Csc[theta[t]] - ---------------}} l :[font = text; inactive; preserveAspect; endGroup; endGroup; ] In one special case, if phi = constant then phi' is zero and so h=0 and we have the simple pendulum case. ^*)