(* Content-type: application/mathematica *)
(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)
(* CreatedBy='Mathematica 6.0' *)
(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[ 145, 7]
NotebookDataLength[ 14799, 386]
NotebookOptionsPosition[ 13802, 346]
NotebookOutlinePosition[ 14218, 364]
CellTagsIndexPosition[ 14175, 361]
WindowFrame->Normal*)
(* Beginning of Notebook Content *)
Notebook[{
Cell[CellGroupData[{
Cell["Precession of Planetary Orbits", "Title"],
Cell["Student Projects in AP1603 investigating classical motion.", \
"Subsubtitle"],
Cell["\<\
All projects must have\[Ellipsis]
Name:____________
Email:_____@columbia.edu\
\>", "Subsubtitle"],
Cell[CellGroupData[{
Cell["Introduction", "Section"],
Cell[TextData[{
"During our classroom discussions, we've investigated the motion of planets \
orbiting a massive star. In particular, we've shown how ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" can be used to visualize elliptic orbits, to demonstrate the conservation \
of angular momentum, and to prove Kepler's \"harmonic law\" of planets. \
Kepler's laws result from the inverse-square dependence of the gravitational \
force on the separation between massive bodies. What happens when the force \
does not have a perfect inverse square dependence?\n\nIn this research \
project, you are to investigate the stability of planetary motion when the \
radial-dependence of the gravitational force changes from the ideal ",
Cell[BoxData[
FormBox[
RowBox[{"1", "/",
SuperscriptBox["r", "2"]}], TraditionalForm]]],
". Nature produces these perturbations in two ways. First, the each planet \
accelerates according to a total force: the vector-sum of the forces from the \
Sun (very large) and the other planets (much smaller, but still important). \
The motion of three or more interacting objects can not be integrated \
analytically, and computer simulation of the dynamics of three (or more) \
objects shows the possibility of orbit instability and chaos. Secondly, near \
extremely massive objects, Newton's law of gravity becomes inaccurate, and we \
need to include corrections due to Einstein's theory of general relativity. \
When the gravitational force deviates from a perfect inverse-square \
dependence, elliptical orbits precess.\n\nThe following two sections present \
related suggestions for your study and research using ",
StyleBox["Mathematica",
FontSlant->"Italic"],
". Take one of these as a starting-point, and conduct \"computational \
experiments\" and analysis to explore how elliptical orbits precess. \
Re-examine the notebooks discussed in class, copy relevant expressions, and \
modify them during your investigations. \n\nYour goal is to produce a \
complete ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" notebook containing a summary of your investigations. This will most \
likely be an answer to one or more questions concerning planetary orbits. It \
should contain ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" expressions and graphics which illustrate your solution. Your notebook \
should be concise and easy to read. Try not to include a large number of \
repeated expressions. Instead, generate a table or graphic of your results. \
In all cases, format your notebook and include textual comments and \
descriptions. Your notebook need not be long. It should be interesting and \
arrive at a definite conclusion."
}], "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell["Modification of the Inverse-Square Law", "Section"],
Cell[TextData[{
"The inverse-square dependence of the gravitational force (and the Coulomb \
force between two charged objects) is related to the ",
Cell[BoxData[
FormBox[
SuperscriptBox["r", "2"], TraditionalForm]]],
" radial dependence of the area of a surface enclosing one object. The force \
of gravity emanates from an object, and the force is directed towards that \
object. The total gravitational \"field lines\" or \"lines of force\" are \
conserved as the enclosing surface moves to different radii. Since the \
surface area scales as ",
Cell[BoxData[
FormBox[
RowBox[{
SuperscriptBox["r", "2"], ","}], TraditionalForm]]],
" the gravitational force must vary as ",
Cell[BoxData[
FormBox[
RowBox[{"1", "/",
SuperscriptBox["r", "2"]}], TraditionalForm]]],
". In this way, the product of the inward gravitational force and the \
surface area of a sphere at radius ",
StyleBox["r",
FontSlant->"Italic"],
" is independent of ",
StyleBox["r",
FontSlant->"Italic"],
". (At least this is true in the classical picture of a Euclidean universe.)\
\n\nAn easy way to explore the consequences of the inverse square dependence \
of gravity is to examine what happens when gravity behaves differently. For \
example, consider the orbit of a planet when the force of gravity is given by \
the equation\n\t",
StyleBox["F",
FontSlant->"Italic"],
"= ",
Cell[BoxData[
FormBox[
RowBox[{"4",
SuperscriptBox["\[Pi]", "2"]}], TraditionalForm]]],
StyleBox["m / ",
FontSlant->"Italic"],
Cell[BoxData[
FormBox[
SuperscriptBox["r", "\[Beta]"], TraditionalForm]]],
"\nwhere \[Beta] is not precisely equal to 2. (In this equation, time is \
measured in years and radius in AU, the average distance between the Sun and \
the Earth.)\n\nExplore what happens to the planet's orbit as \[Beta] changes. \
When \[Beta] = 3, planets are eventually ejected from the solar system. For \
\[Beta] = 2.1, elliptic orbits precess. (What about circular ones?) Is there \
a simple relationship between orbit precession and the value of \[Beta]? \n\n\
[",
StyleBox["For the mathematically inclined,",
FontSlant->"Italic"],
" the precession rate can be estimated easily for nearly circular orbits. \
First, express the equations of motion into polar coordinates, (",
StyleBox["r",
FontSlant->"Italic"],
", \[Theta]). The radial motion depends upon the inward gravitational force \
and the outward centrifugal force. For circular orbits, these two \"forces\" \
balance exactly and ",
Cell[BoxData[
FormBox[
RowBox[{
SubscriptBox["\[PartialD]", "t"],
RowBox[{"r", " ", "~", " ", "0."}]}], TraditionalForm]],
FormatType->"TraditionalForm"],
" For elliptical orbits, ",
Cell[BoxData[
FormBox[
RowBox[{
RowBox[{
SubscriptBox["\[PartialD]", "t"], "r"}], " ", "=", " ", "0"}],
TraditionalForm]],
FormatType->"TraditionalForm"],
" only at the minimum (apogee) and maximum (perigee) distance from the sun. \
Then, using the conservation of angular momentum (since the force is always \
central), the rate of precession can be shown to be equal to difference \
between the ",
StyleBox["angular frequency at perigee",
FontSlant->"Italic"],
" and the ",
StyleBox["frequency of oscillation of the radius",
FontSlant->"Italic"],
" about the average (circular) radius of the orbit. For nearly circular \
orbits, the precession rate is (1- ",
Cell[BoxData[
FormBox[
SqrtBox[
RowBox[{"3", " ", "-", " ", "\[Beta]"}]], TraditionalForm]],
FormatType->"TraditionalForm"],
") times the orbit rate. This requires an understanding of how to ",
StyleBox["linearize",
FontSlant->"Italic"],
" the equations of motion for small deviations from a purely circular orbit.]"
}], "Text",
CellChangeTimes->{{3.474907745937557*^9, 3.474908034398013*^9}, {
3.474908077395359*^9, 3.474908450022564*^9}, {3.474908487134033*^9,
3.474908503740888*^9}}]
}, Open ]],
Cell[CellGroupData[{
Cell["The Precession of Mercury", "Section"],
Cell[TextData[{
"Mercury (the inner most planet) and Pluto (the outermost) have orbits which \
deviate most from circular. Since Mercury is close to the Sun, the variations \
to the usual gravitational force caused by Einstein's general theory of \
relativity cause the elliptical orbit of Mercury to precess. The precession \
of Pluto's elliptic orbit is imperceptible since Pluto is much further from \
the Sun than Mercury.\n\nIn this problem, you are to investigate the \
precession of Mercury due to the first (and most significant) correction to \
Newton's law of gravity due to Einstein. Let the force of gravity be\n\t",
StyleBox["F",
FontSlant->"Italic"],
"= ",
Cell[BoxData[
FormBox[
RowBox[{"4",
SuperscriptBox["\[Pi]", "2"]}], TraditionalForm]]],
StyleBox["m",
FontSlant->"Italic"],
" (1",
StyleBox[" / ",
FontSlant->"Italic"],
Cell[BoxData[
FormBox[
SuperscriptBox["r", "2"], TraditionalForm]]],
" + ",
Cell[BoxData[
FormBox[
RowBox[{"\[Alpha]", "/",
SuperscriptBox["r", "4"]}], TraditionalForm]]],
")\nwhere \[Alpha] \[TildeTilde] 1.1 \[Times] ",
Cell[BoxData[
FormBox[
SuperscriptBox["10",
RowBox[{"-", "8"}]], TraditionalForm]]],
" when radius is measured in units of AU and time in units of years. The \
correction, ",
Cell[BoxData[
FormBox[
RowBox[{"\[Alpha]", "/",
SuperscriptBox["r", "4"]}], TraditionalForm]]],
", is very, very small. Mercury's orbit precesses once every 230,000 years! \
Therefore, in this problem you will need to examine larger values of \
\[Alpha]. The precessional rate will depend on the value of \[Alpha]. When \
\[Alpha] is large, the orbit precession rate will be larger than when \
\[Alpha] is small. (See ",
StyleBox["note",
FontColor->RGBColor[1, 0, 0]],
" below...)\n\nIn order to calculate Mercury's orbit, you will have to \
initialize the orbit with realistic values for its position and velocity. \
Mercury's semimajor radius (the longer radius of its ellipse, ",
StyleBox["a",
FontSlant->"Italic"],
") is about 0.39 AU, and the eccentricity is ",
StyleBox["e ",
FontSlant->"Italic"],
"= 0.206. (Remember, the eccentricity relates the semiminor radius, ",
StyleBox["b, ",
FontSlant->"Italic"],
" to the semimajor radius, ",
StyleBox["a",
FontSlant->"Italic"],
", by ",
StyleBox["b = a ",
FontSlant->"Italic"],
Cell[BoxData[
FormBox[
SqrtBox[
RowBox[{"1", " ", "-", " ",
SuperscriptBox["e", "2"]}]], TraditionalForm]]],
". The distance from the Sun to Mercury is (1 + ",
StyleBox["e) a = ",
FontSlant->"Italic"],
"0.47 AU.) The velocity of Mercury at its greatest distance from the Sun is \
about 8.2 AU/year. This value can actually be determined by consideration of \
conservation of total energy (kinetic plus potential) and conservation of \
angular momentum. \n\nOnce you've initialized Mercury's orbit, examine its \
orbit for various values of \[Alpha]. Carefully consider how to let ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" determine the changing angle of the elliptic axes of the orbit over time. \
See if you can develop a relationship between the value of \[Alpha] and the \
rate of precession of the ellipse. This relation is linear for small values \
of \[Alpha]. By calculating the precession rate for values of \[Alpha] \
between 0.001 and 0.01, you can perform a linear extrapolation to \[Alpha] \
\[TildeTilde] 1.1 \[Times] ",
Cell[BoxData[
FormBox[
SuperscriptBox["10",
RowBox[{"-", "8"}]], TraditionalForm]]],
"and compare your estimate of Mercury's precession to observations. "
}], "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell["The Precession of Moon", "Section",
CellChangeTimes->{{3.4749085869333267`*^9, 3.474908587651841*^9}}],
Cell[TextData[{
"Because the Moon is attracted to both the Earth and the Sun, lunar motion \
is fundamentally a three-body problem. The mathematics of three-body dynamics \
is complex, and the computation of the precession of the Moon's elliptical \
orbit was one of the most challenging of 18th-century physics. Newton \
struggled to explain the lunar orbit, but failed. He calculated a precession \
period of 18 years, which is twice as long as observation. \n\nUsing \
computational methods, the rate of precession of the lunar orbit can be \
understood. The key to the correct solution is to allow both the Earth and \
the Moon to execute orbital motion about each other as the combined Earth and \
Moon system orbits about the Sun. The Sun can be assummed immobile, but the \
Earth and Moon are fully dynamic. You can see Newton's error by neglecting \
the small displacement of the Earth due to the Moon's orbit.\n\nThe relative \
parameters of the the three body system are: Earth: (mass = 5.97 \[Cross] ",
Cell[BoxData[
FormBox[
SuperscriptBox["10", "24"], TraditionalForm]],
FormatType->"TraditionalForm"],
" kg, eccentricity = 0.0167), Moon: (mass = 7.348 \[Cross] ",
Cell[BoxData[
FormBox[
SuperscriptBox["10", "22"], TraditionalForm]],
FormatType->"TraditionalForm"],
" kg, eccentricity = 0.0549). There is a 6.68 degree tilt of the lunar \
orbital plane to that of the Earth-Sun orbital plane. I recommend that you \
ignore this relatively small inclination. "
}], "Text",
CellChangeTimes->{{3.474908827665332*^9, 3.474909292632875*^9}, {
3.474909347990932*^9, 3.4749095591124697`*^9}}]
}, Open ]],
Cell[CellGroupData[{
Cell["Other Ideas for your Project", "Section",
CellChangeTimes->{{3.4749085869333267`*^9, 3.474908587651841*^9}, {
3.474909574098674*^9, 3.474909583039028*^9}}],
Cell[TextData[{
"You are encouraged to investigate any research topic related to the \
solution of coupled nonlinear ordinary differential equations. Please contact \
Prof. Mauel (",
ButtonBox["mauel@columbia.edu",
BaseStyle->"Hyperlink",
ButtonData->{
URL["mailto:mauel@columbia.edu"], None},
ButtonNote->"mailto:mauel@columbia.edu"],
") prior to starting your research project"
}], "Text",
CellChangeTimes->{{3.474908827665332*^9, 3.474909292632875*^9}, {
3.474909347990932*^9, 3.474909682776926*^9}}]
}, Open ]]
}, Open ]]
},
WindowSize->{746, 724},
WindowMargins->{{85, Automatic}, {Automatic, 109}},
PrintingCopies->1,
PrintingPageRange->{1, Automatic},
ShowSelection->True,
FrontEndVersion->"7.0 for Mac OS X x86 (32-bit) (February 18, 2009)",
StyleDefinitions->"Classic.nb"
]
(* End of Notebook Content *)
(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[CellGroupData[{
Cell[567, 22, 47, 0, 58, "Title"],
Cell[617, 24, 83, 1, 54, "Subsubtitle"],
Cell[703, 27, 107, 4, 84, "Subsubtitle"],
Cell[CellGroupData[{
Cell[835, 35, 31, 0, 48, "Section"],
Cell[869, 37, 2707, 47, 345, "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell[3613, 89, 57, 0, 48, "Section"],
Cell[3673, 91, 3923, 97, 393, "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell[7633, 193, 44, 0, 48, "Section"],
Cell[7680, 195, 3601, 90, 411, "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell[11318, 290, 109, 1, 48, "Section"],
Cell[11430, 293, 1621, 28, 210, "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell[13088, 326, 164, 2, 48, "Section"],
Cell[13255, 330, 519, 12, 39, "Text"]
}, Open ]]
}, Open ]]
}
]
*)
(* End of internal cache information *)