4. Lagrangians and Hamiltonians

Legal Notice, Disclaimer, and Book Ordering Information:

Null

Setup

\(Clear["Global`*"]; 
Off[General::spell]; 
Off[General::spell1]; 
    \)

Needs["Calculus`VectorAnalysis`"]

<< "D:\Mma\Inits\ch4init.m"

Clear["Global`*"];

4.2 Problems

4.2.2  Orbiting Bodies

Problem 5:  Quadrupole Problem and Perturbative Solutions

\(TextForm
    \`Consider\ the\ orbit\ of\ a\ satellite\ in\ the\ equitorial\ plane\ of\ 
        an\ oblate\ solid\ of\ revolution . 
		The\ gravitational\ 
        potential\ is\ V\  = \ 
      \(-G\)\ M/r\ 
        \((1\  + \ J\_2/2  r\^2\ \((1 - 3  Cos[θ]\^2)\))\)\)

Part (a)   Assume a satellite is in the equitorial orbit plane θ=π/2 and write the equations of motion.  Eliminate the variable ϕ and get a radial equation for the orbit.

\(\(rRule = 
      r -> \((Evaluate[
					
              CoordinatesToCartesian[

                						{r0[#], θ[#], ϕ[#]}, Spherical]]&)\); 
    \)\)

V = 
      -(G m/r0[t]) 
        (1 + (J2/(2 r0[t]^2)) (1 - 3 Cos[θ[t]]^2))

T = 1/2  r'[t] . r'[t]

L0 = T - V

\(L = 
      \((L0 //. rRule // TrigExpand)\)
 //. {Cos[x_]^2 -> 1 - Sin[x]^2} // 
        Simplify\)

θRule = {θ -> ((π/2)&)};

eqMotion = 
      Lag[{r0, θ, ϕ}, L] //. θRule // Simplify

ϕRule = {ϕ -> (Integrate[dϕ[#], #]&)}

ϕEq = eqMotion[[3]] //. ϕRule

\(dϕSol = 
      DSolve[{ϕEq, dϕ[0] == ϵ1/\((m\ r0[0]^2)\)}, dϕ, 
          t]
 // Flatten\)

\(ϕRule2 = {\(ϕ'\)[t] -> dϕ[t], 
          \(\(ϕ'\)'\)[t] -> \(dϕ'\)[t]}
 //. dϕSol\)

eqMotion2 = eqMotion //. ϕRule2

Part (b)  Let u=1/r and derive an equation for u.

uRule = {r0 -> (1/u[ϕ[#]]&)}

\(eq1 = 
      \((eqMotion2[\([1]\)] //. Join[ϕRule2, uRule]
		 // ExpandAll)
        \)\)

\(eq2 = 
      \((\((Solve[eq1, \(\(u'\)'\)[ϕ[t]]] // Flatten)\)

              		 //. {Rule -> Equal} // ExpandAll)\)
 //. {
          ϕ[t] -> ϕ}\)

Part (c)   Consider the elliptical solution that follows from the equation with (Subscript[J, 2])=0.  Find the first order perturbation correction due to the quadrupole term (Subscript[J, 2]) by expanding the equations of motion about (Subscript[J, 2]) and solving.  Keep only the lowest-order terms in ϵ.  Show that the results can be expressed as a precessing ellipse.

initial = {u'[0] == 0, u[0] == u0 (1 + ϵ)}

\(dsol = 
      \((DSolve[{\((eq2 /. {J2 -> 0})\), initial} // Flatten

              			, {u[ϕ]}, ϕ] // Flatten)\) // Simplify\)

\(uSol = 
      Solve[\((\(u[ϕ] /. dsol\) /. {Cos[ϕ] -> 0})\) == u0, u0]
 // 
        Flatten\)

pertRule = {
        u -> ((u0 (1 + ϵ Cos[#]) + J2 u1[#])&)}

eq3 = eq2 /. pertRule

\(eq4 = 
      \((\((\(Solve[eq3, \(\(u1'\)'\)[ϕ]]\)[\([1]\)]
		 //. uSol // 
                    ExpandAll)\)
	 //. {J2 -> 0, ϵ^2 -> 0} // 
              Simplify)\)
 //. {Rule -> Equal} // ExpandAll\)

init1 = {u1'[0] == 0, u1[0] == 0}

\(dsol4 = 
      \(DSolve[Join[eq4, init1]
			, u1[ϕ], ϕ] // Flatten\) // 
        ExpandAll\)

seq5 = u[ϕ] /. pertRule /. dsol4

eq5 = u01 (1 + ϵ1 Cos[(1 - J2 α) ϕ])

seq6 = Series[eq5, {J2, 0, 1}] // Normal

u01Rule = 
      Solve[seq5 == seq6 //. {Cos[__] -> 0, Sin[__] -> 0}, u01] // Flatten
           // ExpandAll

eqϵ1 =  
      seq5 ϵ1^6 == seq6 ϵ1^6 //. {Sin[__] -> 0} //. 
          u01Rule // Simplify

Note:  The polynomial eqϵ1 has 7 roots.  Take the 1st one.

ϵ1Rule=Solve[eqϵ1, ϵ1][[1]]

eqα = 
      seq5 ϵ1^6 == seq6 ϵ1^6 //. {Cos[__] -> 0} // 
        Simplify //. u01Rule //. ϵ1Rule

αRule = 
      (Solve[eqα, α] // Flatten) //. uSol // Simplify

\(\((\({u01, ϵ1, J2\ α}
 //. 
          Join[u01Rule, ϵ1Rule, αRule]\)
	 //. {J2 -> 0})
      \)\)

Part (d)  Use the user-defined function firstOrderPert to find the perturbative solution and show that it agrees with the solution found in Part (c).

Note:  This operation takes about 15 minutes on a 200 MHz Pentium.

\(test1a = 
      firstOrderPert[eq2[\([1]\)], u, {u0 \((1 + ϵ)\), 0}, J2, 
            ϕ]
	\ \ \  //. uSol // Timing\)

test1b = test1a[[2]] // Simplify

Note:  Get rid of the high order terms in ϵ.

test1c=test1b//.{ϵ^2 -> 0}

test1 = test1c // ExpandAll

test1 = test1d // ExpandAll

Note:  The parentheses are important.  Otherwise the uSol substitution doesn't work.

test2=(Series[seq5, {J2, 0, 1}] // Normal)//.uSol//ExpandAll

test1 == test2