(* 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[ 57212, 1883] NotebookOptionsPosition[ 51969, 1722] NotebookOutlinePosition[ 52382, 1740] CellTagsIndexPosition[ 52339, 1737] WindowFrame->Normal ContainsDynamic->False*) (* Beginning of Notebook Content *) Notebook[{ Cell[CellGroupData[{ Cell["The Many Body Problem", "Title"], Cell["\<\ AP1601 Columbia University\ \>", "Subsubtitle"], Cell[CellGroupData[{ Cell["Introduction", "Section"], Cell["\<\ The equations of motion for two bodies interacting with a central force can \ be solved exactly. However, when more than two bodies are interacting, \ solutions are possible only in spe cial cases or with simplifying \ assumptions. When we assumed that the mass of the sun was much larger than \ the planets, we were making one of these simplifying assumptions.\ \>", "Text"], Cell["\<\ In this notebook, we will explore the solution of several objects interacting \ through gravitational forces. \ \>", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Representing Many Objects", "Section"], Cell["How should we represent many objects?", "Text"], Cell[TextData[{ "I find it most convenient to represent the objects as ", StyleBox["arrays", FontSlant->"Italic"], " indexed by an integer. Arrays can be randomly accessed by using a integer \ within the range defined for that array." }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Off", "[", RowBox[{ RowBox[{"General", "::", "\"\\""}], ",", RowBox[{"General", "::", "\"\\""}]}], "]"}], ";"}], ")"}], " ", RowBox[{"(", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Needs", "[", "\"\\"", "]"}], ";", RowBox[{"Needs", "[", "\"\\"", "]"}]}], ")"}], ";"}], ")"}]}]], "Input", CellLabel->"In[1]:="], Cell[BoxData[{ RowBox[{ StyleBox[ RowBox[{"red", " ", "=", " ", RowBox[{"RGBColor", "[", RowBox[{"1", ",", "0", ",", "0"}], "]"}]}], FontColor->RGBColor[1, 0, 0]], ";"}], "\[IndentingNewLine]", RowBox[{ StyleBox[ RowBox[{"blue", " ", "=", " ", RowBox[{"RGBColor", "[", RowBox[{"0", ",", "0", ",", "1"}], "]"}]}], FontColor->RGBColor[0, 0, 1]], ";"}]}], "Input", CellLabel->"In[2]:="], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"a", "[", "1", "]"}], " ", "=", " ", "3"}], ";", " ", RowBox[{ RowBox[{"a", "[", "4", "]"}], " ", "=", " ", "2"}], ";", " ", RowBox[{ RowBox[{ RowBox[{"a", "[", "7", "]"}], "[", "t_", "]"}], " ", "=", " ", RowBox[{"t", "^", "2"}]}], ";"}]], "Input", CellLabel->"In[4]:="], Cell[BoxData[ RowBox[{"?", "a"}]], "Input", CellLabel->"In[5]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{ RowBox[{"a", "[", "7", "]"}], "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5"}], "}"}]}], "]"}]], "Input", CellLabel->"In[6]:="], Cell[BoxData[ RowBox[{"bList", " ", "=", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"b", "[", "i", "]"}], "[", "t_", "]"}], " ", "=", " ", RowBox[{"t", "^", "i"}]}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "10"}], "}"}]}], "]"}]}]], "Input",\ CellLabel->"In[7]:="], Cell[BoxData[ RowBox[{"?", "b"}]], "Input", CellLabel->"In[8]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"b", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", "10"}], "}"}]}], "]"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "2"}], "}"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"Table", "[", RowBox[{ RowBox[{"Hue", "[", FractionBox["i", "10"], "]"}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", "10"}], "}"}]}], "]"}]}]}], "]"}]], "Input",\ CellLabel->"In[9]:="] }, Closed]], Cell[CellGroupData[{ Cell["Equations of Motion", "Section"], Cell[CellGroupData[{ Cell["The Definitions", "Subsection"], Cell["Imagine a collection of objects of equal mass.", "Text"], Cell[BoxData[ RowBox[{ RowBox[{ StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]], " ", "=", " ", "5"}], ";"}]], "Input", CellChangeTimes->{ 3.443438998430773*^9, {3.443439059202115*^9, 3.4434390592721567`*^9}, { 3.443439286746182*^9, 3.4434392888118477`*^9}, 3.443439493588707*^9}, CellLabel->"In[183]:="], Cell[BoxData[ RowBox[{"mass", " ", "=", " ", RowBox[{ RowBox[{"{", StyleBox["10.", FontColor->RGBColor[1, 0, 0]], "}"}], " ", "~", " ", "Join", "~", " ", RowBox[{"Table", "[", RowBox[{ StyleBox["1.", FontColor->RGBColor[0, 0, 1]], ",", " ", RowBox[{"{", RowBox[{"i", ",", "2", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}]}]}]], "Input", CellChangeTimes->{3.443438892470387*^9}, CellLabel->"In[184]:="], Cell[BoxData[ RowBox[{ RowBox[{"p", " ", "=", " ", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"vx", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"vy", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"vz", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}]], "Input", CellChangeTimes->{3.443438896353829*^9}, CellLabel->"In[185]:="], Cell[TextData[{ "\"", StyleBox["p", FontWeight->"Bold"], "\" is a list of all of the unknown dynamical functions of time." }], "Text"], Cell[BoxData["p"], "Input", CellLabel->"In[186]:="], Cell[TextData[{ "The force (per unit mass) on object \"", StyleBox["i", FontWeight->"Bold"], "\" from object \"", StyleBox["j", FontWeight->"Bold"], "\" is found from..." }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"distance", "[", RowBox[{"i_", ",", "j_"}], "]"}], " ", ":=", " ", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], " ", "-", " ", RowBox[{ RowBox[{"x", "[", "j", "]"}], "[", "t", "]"}]}], ")"}], "^", "2"}], " ", "+", " ", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], " ", "-", " ", RowBox[{ RowBox[{"y", "[", "j", "]"}], "[", "t", "]"}]}], ")"}], "^", "2"}], " ", "+", " ", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}], " ", "-", " ", RowBox[{ RowBox[{"z", "[", "j", "]"}], "[", "t", "]"}]}], ")"}], "^", "2"}]}], "]"}]}]], "Input", CellLabel->"In[187]:="], Cell[BoxData[{ RowBox[{ RowBox[{ RowBox[{"forcex", "[", RowBox[{"i_", ",", "j_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"x", "[", "j", "]"}], "[", "t", "]"}], " ", "-", " ", RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}]}], ")"}], RowBox[{ RowBox[{"mass", "[", RowBox[{"[", "j", "]"}], "]"}], "/", RowBox[{ RowBox[{"distance", "[", RowBox[{"i", ",", "j"}], "]"}], "^", "3"}]}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"forcey", "[", RowBox[{"i_", ",", "j_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"y", "[", "j", "]"}], "[", "t", "]"}], " ", "-", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}]}], ")"}], RowBox[{ RowBox[{"mass", "[", RowBox[{"[", "j", "]"}], "]"}], "/", RowBox[{ RowBox[{"distance", "[", RowBox[{"i", ",", "j"}], "]"}], "^", "3"}]}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"forcez", "[", RowBox[{"i_", ",", "j_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"z", "[", "j", "]"}], "[", "t", "]"}], " ", "-", " ", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], ")"}], RowBox[{ RowBox[{"mass", "[", RowBox[{"[", "j", "]"}], "]"}], "/", RowBox[{ RowBox[{"distance", "[", RowBox[{"i", ",", "j"}], "]"}], "^", "3"}]}]}]}], ";"}]}], "Input", CellLabel->"In[188]:="], Cell["\<\ We must sum the forces from all objects, except we must exclude self-forces.\ \>", "Text"], Cell[BoxData[ RowBox[{"?", "Complement"}]], "Input", CellLabel->"In[191]:="], Cell[BoxData[ RowBox[{ RowBox[{"dp", " ", "=", " ", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"vx", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"vy", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"vz", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{"Plus", " ", "@@", " ", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"forcex", "[", RowBox[{"i", ",", "#"}], "]"}], "&"}], " ", "/@", " ", RowBox[{"Complement", "[", RowBox[{ RowBox[{"Range", "[", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]], "]"}], ",", " ", RowBox[{"{", "i", "}"}]}], "]"}]}], ")"}]}], ",", "\[IndentingNewLine]", RowBox[{"Plus", " ", "@@", " ", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"forcey", "[", RowBox[{"i", ",", "#"}], "]"}], "&"}], " ", "/@", " ", RowBox[{"Complement", "[", RowBox[{ RowBox[{"Range", "[", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]], "]"}], ",", " ", RowBox[{"{", "i", "}"}]}], "]"}]}], ")"}]}], ",", "\[IndentingNewLine]", RowBox[{"Plus", " ", "@@", " ", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"forcez", "[", RowBox[{"i", ",", "#"}], "]"}], "&"}], " ", "/@", " ", RowBox[{"Complement", "[", RowBox[{ RowBox[{"Range", "[", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]], "]"}], ",", " ", RowBox[{"{", "i", "}"}]}], "]"}]}], ")"}]}]}], "}"}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}]], "Input", CellChangeTimes->{{3.443438902625885*^9, 3.443438905642428*^9}}, CellLabel->"In[192]:="], Cell["A particular initial condition...", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"p0", "=", RowBox[{"Flatten", "[", RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "0", "]"}], "\[Equal]", RowBox[{"RandomReal", "[", RowBox[{"{", RowBox[{ RowBox[{"-", "1.`"}], ",", "1.`"}], "}"}], "]"}]}], ",", RowBox[{ RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "0", "]"}], "\[Equal]", RowBox[{"RandomReal", "[", RowBox[{"{", RowBox[{ RowBox[{"-", "1.`"}], ",", "1.`"}], "}"}], "]"}]}], ",", RowBox[{ RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "0", "]"}], "\[Equal]", RowBox[{"RandomReal", "[", RowBox[{"{", RowBox[{ RowBox[{"-", "1.`"}], ",", "1.`"}], "}"}], "]"}]}], ",", RowBox[{ RowBox[{ RowBox[{"vx", "[", "i", "]"}], "[", "0", "]"}], "\[Equal]", RowBox[{"RandomReal", "[", RowBox[{"{", RowBox[{ RowBox[{"-", "1.`"}], ",", "1.`"}], "}"}], "]"}]}], ",", RowBox[{ RowBox[{ RowBox[{"vy", "[", "i", "]"}], "[", "0", "]"}], "\[Equal]", RowBox[{"RandomReal", "[", RowBox[{"{", RowBox[{ RowBox[{"-", "1.`"}], ",", "1.`"}], "}"}], "]"}]}], ",", RowBox[{ RowBox[{ RowBox[{"vz", "[", "i", "]"}], "[", "0", "]"}], "\[Equal]", RowBox[{"RandomReal", "[", RowBox[{"{", RowBox[{ RowBox[{"-", "1.`"}], ",", "1.`"}], "}"}], "]"}]}]}], "}"}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}], "]"}]}], ";"}]], "Input", CellChangeTimes->{3.443438907599663*^9}, CellLabel->"In[193]:="], Cell[BoxData[ RowBox[{ RowBox[{"eqs", " ", "=", " ", RowBox[{ RowBox[{"Thread", "[", RowBox[{ RowBox[{"D", "[", RowBox[{"p", ",", "t"}], "]"}], " ", "==", " ", "dp"}], "]"}], "~", " ", "Join", " ", "~", " ", "p0"}]}], ";"}]], "Input", CellLabel->"In[194]:="], Cell[TextData[{ "\"", StyleBox["eqs", FontWeight->"Bold"], "\" are the equations of motion!" }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Example Solution", "Subsection"], Cell[BoxData[ RowBox[{ RowBox[{"Print", "[", RowBox[{"\"\\"", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "]"}], ";", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"sol1", " ", "=", " ", RowBox[{ RowBox[{"NDSolve", "[", RowBox[{"eqs", " ", ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "5"}], "}"}], ",", " ", RowBox[{"MaxSteps", " ", "\[Rule]", " ", "80000"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}], ")"}], " ", "//", " ", "Timing"}]}]], "Input", CellChangeTimes->{3.443438911505159*^9, 3.44343932996006*^9}, CellLabel->"In[195]:="], Cell[TextData[{ "How long does ", StyleBox["Mathematica", FontSlant->"Italic"], " require to solve these equations as the number of planetary objects \ change? Try this with ", StyleBox["nBodies = 5, 10, 15", FontWeight->"Bold"], ". (", StyleBox["Be warned", FontWeight->"Bold", FontColor->RGBColor[1, 0, 0]], ": the solution time will change greatly depending upon the initial \ condition! If two objects are too close, then the solution requires very many \ time steps to maintain accuracy.)" }], "Text", CellChangeTimes->{{3.443439217759988*^9, 3.44343926854702*^9}, { 3.443439416737691*^9, 3.4434394703682117`*^9}}] }, Closed]], Cell[CellGroupData[{ Cell["Visualizing the Solution", "Subsection"], Cell["To view the solution, we need 3D graphics.", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Needs", "[", "\"\\"", "]"}], ";", RowBox[{"Needs", "[", "\"\\"", "]"}]}], ")"}], ";"}]], "Input", CellLabel->"In[166]:="], Cell[BoxData[ RowBox[{ RowBox[{"threeRange", " ", "=", " ", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"-", "3.0"}], ",", " ", "3.0"}], "}"}], ",", " ", RowBox[{"{", RowBox[{ RowBox[{"-", "3.0"}], ",", " ", "3.0"}], "}"}], ",", " ", RowBox[{"{", RowBox[{ RowBox[{"-", "3.0"}], ",", " ", "3.0"}], "}"}]}], "}"}]}], ";"}]], "Input", CellLabel->"In[167]:="], Cell[BoxData[ RowBox[{ RowBox[{"objects", "[", "tt_", "]"}], " ", ":=", " ", RowBox[{"Graphics3D", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"PointSize", "[", "0.04", "]"}], ",", StyleBox["red", FontColor->RGBColor[1, 0, 0]], ",", " ", RowBox[{"Point", "[", StyleBox[ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "1", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "1", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"z", "[", "1", "]"}], "[", "t", "]"}]}], "}"}], FontColor->RGBColor[1, 0, 0]], "]"}], ",", " ", RowBox[{"PointSize", "[", "0.02", "]"}], ",", StyleBox["blue", FontColor->RGBColor[0, 0, 1]]}], "}"}], " ", "~", " ", "Join", " ", "~", RowBox[{"Table", "[", RowBox[{ RowBox[{"Point", "[", StyleBox[ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], FontColor->RGBColor[0, 0, 1]], "]"}], " ", ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "2", ",", " ", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}]}], " ", "/.", " ", "sol1"}], " ", "/.", " ", RowBox[{"t", " ", "\[Rule]", " ", "tt"}]}], ",", " ", RowBox[{"PlotRange", "->", "threeRange"}]}], "]"}]}]], "Input", CellChangeTimes->{3.443438917623988*^9}, CellLabel->"In[168]:="], Cell[BoxData[ RowBox[{"Show", "[", RowBox[{"objects", "[", "1.`", "]"}], "]"}]], "Input", CellLabel->"In[169]:="] }, Closed]], Cell[CellGroupData[{ Cell["Annimation", "Subsection"], Cell[BoxData[ RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"Show", "[", RowBox[{"objects", "[", "t", "]"}], "]"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0.0", ",", " ", "5.0", ",", " ", "0.1"}], "}"}]}], "]"}], ";"}]], "Input", CellLabel->"In[170]:="], Cell[TextData[{ "We made a oversight. In order to follow the moving objects, we must ", StyleBox["change our frame of reference", FontSlant->"Italic"], " so that we are moving and located at the \"center of gravity\" of the many \ bodies." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Averages", "Subsection"], Cell["We can find the average radius and the average position.", "Text"], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"avgRadius", "[", "tt_", "]"}], " ", "=", " ", RowBox[{ RowBox[{"(", RowBox[{"Plus", " ", "@@", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], "^", "2"}], " ", "+", " ", RowBox[{ RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], "^", "2"}], " ", "+", " ", RowBox[{ RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}], "^", "2"}]}], "]"}], " ", "/.", " ", "sol1"}], "/.", " ", RowBox[{"t", " ", "\[Rule]", " ", "tt"}]}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "nBodies"}], "}"}]}], "]"}]}], " ", ")"}], "/", "nBodies"}]}], ";"}]], "Input", CellLabel->"In[171]:="], Cell[TextData[{ "The ", StyleBox["Apply[...]", FontWeight->"Bold"], " function is used in its \"infix\" form..." }], "Text"], Cell[BoxData[ RowBox[{"Plus", " ", "@@", " ", RowBox[{"{", RowBox[{"a1", ",", " ", "a2", ",", " ", "a3", ",", " ", "a4"}], "}"}]}]], "Input", CellLabel->"In[172]:="], Cell[BoxData[ RowBox[{"avgRadius", "[", "0.2", "]"}]], "Input", CellLabel->"In[173]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"avgRadius", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5"}], "}"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input", CellLabel->"In[174]:="], Cell["The average position is a vector...", "Text"], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"avgPosition", "[", "tt_", "]"}], " ", "=", " ", RowBox[{ RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], "/.", " ", "sol1"}], " ", "/.", " ", RowBox[{"t", " ", "\[Rule]", " ", "tt"}]}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "nBodies"}], "}"}]}], "]"}]}], ")"}], "/", "nBodies"}]}], ";"}]], "Input", CellLabel->"In[175]:="], Cell[BoxData[ RowBox[{"Plus", " ", "@@", " ", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{"x1", ",", "y1", ",", "z1"}], "}"}], ",", RowBox[{"{", RowBox[{"x2", ",", "y2", ",", "z2"}], "}"}], ",", RowBox[{"{", RowBox[{"x3", ",", "y3", ",", "z3"}], "}"}]}], "}"}]}]], "Input", CellLabel->"In[176]:="], Cell[BoxData[ RowBox[{"avgPosition", "[", "0.2", "]"}]], "Input", CellLabel->"In[177]:="], Cell[BoxData[ RowBox[{"sp1", "=", RowBox[{"ListPointPlot3D", "[", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"avgPosition", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5", ",", "0.1`"}], "}"}]}], "]"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"{", RowBox[{"blue", ",", RowBox[{"PointSize", "[", "0.02`", "]"}]}], "}"}]}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]}]], "Input", CellLabel->"In[178]:="], Cell["The center of mass (COM) is the density-weighted position...", "Text"], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"avgCOM", "[", "tt_", "]"}], " ", "=", " ", RowBox[{ RowBox[{"(", RowBox[{"Plus", " ", "@@", RowBox[{"(", RowBox[{"mass", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], "/.", " ", "sol1"}], " ", "/.", " ", RowBox[{"t", " ", "\[Rule]", " ", "tt"}]}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "nBodies"}], "}"}]}], "]"}]}], ")"}]}], ")"}], "/", RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", "mass"}], ")"}]}]}], ";"}]], "Input", CellLabel->"In[179]:="], Cell[BoxData[ RowBox[{"sp2", "=", RowBox[{"ListPointPlot3D", "[", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"avgCOM", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5", ",", "0.1`"}], "}"}]}], "]"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"{", RowBox[{"red", ",", RowBox[{"PointSize", "[", "0.02`", "]"}]}], "}"}]}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]}]], "Input", CellLabel->"In[180]:="], Cell[BoxData[ RowBox[{"Show", "[", RowBox[{"sp1", ",", "sp2", ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input", CellLabel->"In[181]:="] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Removing the COM", "Section"], Cell["\<\ Since the mass of each object is identical, we can easily sum (and average) \ the initial velocities and positions of the objects in order to find the \ center of gravity and the center of mass velocity.\ \>", "Text"], Cell[BoxData[ RowBox[{"Print", "[", RowBox[{"\"\\"", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], StyleBox["]", FontColor->RGBColor[1, 0, 1]]}]], "Input", CellChangeTimes->{{3.44343893024223*^9, 3.443438945284687*^9}}, CellLabel->"In[118]:="], Cell[CellGroupData[{ Cell["Defining the COM", "Subsection"], Cell[BoxData[ RowBox[{"v0", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"RandomReal", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"-", "1."}], ",", "1."}], "}"}], ",", "3"}], "]"}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}]}]], "Input", CellChangeTimes->{{3.443438853101368*^9, 3.4434388543793707`*^9}, 3.443438927448036*^9}, CellLabel->"In[119]:="], Cell[BoxData[ RowBox[{"vcom", " ", "=", " ", RowBox[{ RowBox[{"MapThread", "[", RowBox[{"Plus", ",", " ", RowBox[{"v0", " ", "mass"}]}], "]"}], "/", RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", "mass"}], ")"}]}]}]], "Input", CellLabel->"In[120]:="], Cell[BoxData[ RowBox[{ RowBox[{"r0", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"RandomReal", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"-", "1."}], ",", "1."}], "}"}], ",", "3"}], "]"}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}]}], ";"}]], "Input", CellChangeTimes->{{3.4434388557158003`*^9, 3.443438856707183*^9}, 3.4434389551900997`*^9}, CellLabel->"In[121]:="], Cell[BoxData[ RowBox[{"rcom", " ", "=", " ", RowBox[{ RowBox[{"MapThread", "[", RowBox[{"Plus", ",", " ", RowBox[{"r0", " ", "mass"}]}], "]"}], "/", RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", "mass"}], ")"}]}]}]], "Input", CellLabel->"In[122]:="], Cell[BoxData[ RowBox[{ RowBox[{"p0", " ", "=", " ", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "1"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "1", "]"}], "]"}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "2"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "2", "]"}], "]"}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "3"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "3", "]"}], "]"}]}]}], ",", " ", RowBox[{ RowBox[{ RowBox[{"vx", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0", "[", RowBox[{"[", RowBox[{"i", ",", "1"}], "]"}], "]"}], " ", "-", " ", RowBox[{"vcom", "[", RowBox[{"[", "1", "]"}], "]"}]}]}], ",", " ", RowBox[{ RowBox[{ RowBox[{"vy", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0", "[", RowBox[{"[", RowBox[{"i", ",", "2"}], "]"}], "]"}], " ", "-", " ", RowBox[{"vcom", "[", RowBox[{"[", "2", "]"}], "]"}]}]}], ",", RowBox[{ RowBox[{ RowBox[{"vz", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0", "[", RowBox[{"[", RowBox[{"i", ",", "3"}], "]"}], "]"}], "-", " ", RowBox[{"vcom", "[", RowBox[{"[", "3", "]"}], "]"}]}]}]}], "}"}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}]], "Input", CellChangeTimes->{3.4434389580593843`*^9}, CellLabel->"In[123]:="], Cell[BoxData[ RowBox[{ RowBox[{"eqs2", " ", "=", " ", RowBox[{ RowBox[{"Thread", "[", RowBox[{ RowBox[{"D", "[", RowBox[{"p", ",", "t"}], "]"}], " ", "==", " ", "dp"}], "]"}], "~", " ", "Join", " ", "~", " ", "p0"}]}], ";"}]], "Input", CellLabel->"In[124]:="] }, Closed]], Cell[CellGroupData[{ Cell["Example Solution", "Subsection"], Cell[BoxData[ RowBox[{ RowBox[{"Print", "[", RowBox[{"\"\\"", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "]"}], ";", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"sol2", " ", "=", " ", RowBox[{ RowBox[{"NDSolve", "[", RowBox[{"eqs2", " ", ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "5"}], "}"}], ",", " ", RowBox[{"MaxSteps", " ", "\[Rule]", " ", "40000"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}], ")"}], " ", "//", " ", "Timing"}]}]], "Input", CellChangeTimes->{3.443438959832045*^9}, CellLabel->"In[125]:="] }, Closed]], Cell[CellGroupData[{ Cell["Annimation", "Subsection"], Cell[BoxData[ RowBox[{ RowBox[{"objects2", "[", "tt_", "]"}], " ", ":=", " ", RowBox[{"Graphics3D", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"PointSize", "[", "0.04", "]"}], ",", StyleBox["red", FontColor->RGBColor[1, 0, 0]], ",", " ", RowBox[{"Point", "[", StyleBox[ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "1", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "1", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"z", "[", "1", "]"}], "[", "t", "]"}]}], "}"}], FontColor->RGBColor[1, 0, 0]], "]"}], ",", " ", RowBox[{"PointSize", "[", "0.02", "]"}], ",", StyleBox["blue", FontColor->RGBColor[0, 0, 1]]}], "}"}], " ", "~", " ", "Join", " ", "~", RowBox[{"Table", "[", RowBox[{ RowBox[{"Point", "[", StyleBox[ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], FontColor->RGBColor[0, 0, 1]], "]"}], " ", ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "2", ",", " ", "nBodies"}], "}"}]}], "]"}]}], " ", "/.", " ", "sol2"}], " ", "/.", " ", RowBox[{"t", " ", "\[Rule]", " ", "tt"}]}], ",", " ", RowBox[{"PlotRange", "->", "threeRange"}]}], "]"}]}]], "Input", CellLabel->"In[126]:="], Cell[BoxData[ RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"Show", "[", RowBox[{"objects2", "[", "t", "]"}], "]"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0.0", ",", " ", "5.0", ",", " ", "0.1"}], "}"}]}], "]"}], ";"}]], "Input", CellLabel->"In[127]:="], Cell[BoxData[ RowBox[{"ParametricPlot3D", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", RowBox[{"i", "\[Rule]", "1"}]}], "/.", "\[InvisibleSpace]", "sol2"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5"}], "}"}], ",", RowBox[{"PlotPoints", "\[Rule]", "200"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input",\ CellLabel->"In[128]:="], Cell[BoxData[ RowBox[{"ParametricPlot3D", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", RowBox[{"i", "\[Rule]", "3"}]}], "/.", "\[InvisibleSpace]", "sol2"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5"}], "}"}], ",", RowBox[{"PlotPoints", "\[Rule]", "200"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input",\ CellLabel->"In[129]:="], Cell[TextData[{ "Now, what happens when ", StyleBox["nBodies -> 2", FontWeight->"Bold"], "?" }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Averages", "Subsection"], Cell["We can find the average radius and the average position.", "Text"], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"avgRadius", "[", "tt_", "]"}], " ", "=", " ", RowBox[{ RowBox[{"(", RowBox[{"Plus", " ", "@@", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], "^", "2"}], " ", "+", " ", RowBox[{ RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], "^", "2"}], " ", "+", " ", RowBox[{ RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}], "^", "2"}]}], "]"}], " ", "/.", " ", "sol2"}], "/.", " ", RowBox[{"t", " ", "\[Rule]", " ", "tt"}]}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "nBodies"}], "}"}]}], "]"}]}], " ", ")"}], "/", "nBodies"}]}], ";"}]], "Input", CellLabel->"In[130]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"avgRadius", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5"}], "}"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input", CellLabel->"In[131]:="], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"avgPosition", "[", "tt_", "]"}], " ", "=", " ", RowBox[{ RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], "/.", " ", "sol2"}], " ", "/.", " ", RowBox[{"t", " ", "\[Rule]", " ", "tt"}]}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "nBodies"}], "}"}]}], "]"}]}], ")"}], "/", "nBodies"}]}], ";"}]], "Input", CellLabel->"In[132]:="], Cell[BoxData[ RowBox[{"sp1", "=", RowBox[{"ListPointPlot3D", "[", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"avgPosition", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5", ",", "0.1`"}], "}"}]}], "]"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"{", RowBox[{"blue", ",", RowBox[{"PointSize", "[", "0.02`", "]"}]}], "}"}]}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]}]], "Input", CellLabel->"In[133]:="], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{"avgCOM", "[", "tt_", "]"}], " ", "=", " ", RowBox[{ RowBox[{"(", RowBox[{"Plus", " ", "@@", RowBox[{"(", RowBox[{"mass", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], "/.", " ", "sol2"}], " ", "/.", " ", RowBox[{"t", " ", "\[Rule]", " ", "tt"}]}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "nBodies"}], "}"}]}], "]"}]}], ")"}]}], ")"}], "/", RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", "mass"}], ")"}]}]}], ";"}]], "Input", CellLabel->"In[134]:="], Cell[BoxData[ RowBox[{"sp2", "=", RowBox[{"ListPointPlot3D", "[", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"avgCOM", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "5", ",", "0.1"}], "}"}]}], "]"}], ",", RowBox[{"PlotStyle", "\[Rule]", RowBox[{"{", RowBox[{"red", ",", RowBox[{"PointSize", "[", "0.02", "]"}]}], "}"}]}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]}]], "Input", CellChangeTimes->{{3.443438772661255*^9, 3.4434387816125727`*^9}}, CellLabel->"In[135]:="], Cell["The COM remains \"zero\" within a small numerical error.", "Text"], Cell[BoxData[ RowBox[{"Show", "[", RowBox[{"sp1", ",", "sp2", ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input", CellLabel->"In[136]:="] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Nearby Orbits", "Section"], Cell["\<\ Let's choose two nearby initial conditions: the second differing only by a \ slightly different velocity.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"tMax", " ", "=", " ", "5.0"}], ";"}]], "Input", CellLabel->"In[137]:="], Cell[BoxData[ RowBox[{"v0a", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"RandomReal", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"-", "1.`"}], ",", "1.`"}], "}"}], ",", "3"}], "]"}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}]}]], "Input", CellChangeTimes->{3.443438969757483*^9}, CellLabel->"In[138]:="], Cell[BoxData[ RowBox[{"v0b", " ", "=", " ", RowBox[{ RowBox[{"Drop", "[", RowBox[{"v0a", ",", RowBox[{"-", "1"}]}], "]"}], " ", "~", " ", "Join", "~", " ", RowBox[{"{", RowBox[{ RowBox[{"Last", "[", "v0a", "]"}], " ", RowBox[{"{", RowBox[{"1.001", ",", "1", ",", "1"}], "}"}]}], "}"}]}]}]], "Input", CellLabel->"In[139]:="], Cell[BoxData[{ RowBox[{ RowBox[{"vcoma", " ", "=", " ", RowBox[{ RowBox[{"MapThread", "[", RowBox[{"Plus", ",", " ", RowBox[{"v0a", " ", "mass"}]}], "]"}], "/", RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", "mass"}], ")"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"vcomb", " ", "=", " ", RowBox[{ RowBox[{"MapThread", "[", RowBox[{"Plus", ",", " ", RowBox[{"v0b", " ", "mass"}]}], "]"}], "/", RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", "mass"}], ")"}]}]}], ";"}]}], "Input", CellLabel->"In[140]:="], Cell[BoxData[ RowBox[{ RowBox[{"r0", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"RandomReal", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"-", "1.`"}], ",", "1.`"}], "}"}], ",", "3"}], "]"}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", StyleBox["nBodies", FontColor->RGBColor[1, 0, 1]]}], "}"}]}], "]"}]}], ";"}]], "Input", CellChangeTimes->{3.4434389726629953`*^9}, CellLabel->"In[142]:="], Cell[BoxData[ RowBox[{ RowBox[{"rcom", " ", "=", " ", RowBox[{ RowBox[{"MapThread", "[", RowBox[{"Plus", ",", " ", RowBox[{"r0", " ", "mass"}]}], "]"}], "/", RowBox[{"(", RowBox[{"Plus", " ", "@@", " ", "mass"}], ")"}]}]}], ";"}]], "Input", CellLabel->"In[143]:="], Cell[BoxData[{ RowBox[{ RowBox[{"p0a", " ", "=", " ", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "1"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "1", "]"}], "]"}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "2"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "2", "]"}], "]"}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "3"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "3", "]"}], "]"}]}]}], ",", " ", RowBox[{ RowBox[{ RowBox[{"vx", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0a", "[", RowBox[{"[", RowBox[{"i", ",", "1"}], "]"}], "]"}], " ", "-", " ", RowBox[{"vcoma", "[", RowBox[{"[", "1", "]"}], "]"}]}]}], ",", " ", RowBox[{ RowBox[{ RowBox[{"vy", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0a", "[", RowBox[{"[", RowBox[{"i", ",", "2"}], "]"}], "]"}], " ", "-", " ", RowBox[{"vcoma", "[", RowBox[{"[", "2", "]"}], "]"}]}]}], ",", RowBox[{ RowBox[{ RowBox[{"vz", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0a", "[", RowBox[{"[", RowBox[{"i", ",", "3"}], "]"}], "]"}], "-", " ", RowBox[{"vcoma", "[", RowBox[{"[", "3", "]"}], "]"}]}]}]}], "}"}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "nBodies"}], "}"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"p0b", " ", "=", " ", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "1"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "1", "]"}], "]"}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "2"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "2", "]"}], "]"}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"r0", "[", RowBox[{"[", RowBox[{"i", ",", "3"}], "]"}], "]"}], " ", "-", " ", RowBox[{"rcom", "[", RowBox[{"[", "3", "]"}], "]"}]}]}], ",", " ", RowBox[{ RowBox[{ RowBox[{"vx", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0b", "[", RowBox[{"[", RowBox[{"i", ",", "1"}], "]"}], "]"}], " ", "-", " ", RowBox[{"vcomb", "[", RowBox[{"[", "1", "]"}], "]"}]}]}], ",", " ", RowBox[{ RowBox[{ RowBox[{"vy", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0b", "[", RowBox[{"[", RowBox[{"i", ",", "2"}], "]"}], "]"}], " ", "-", " ", RowBox[{"vcomb", "[", RowBox[{"[", "2", "]"}], "]"}]}]}], ",", RowBox[{ RowBox[{ RowBox[{"vz", "[", "i", "]"}], "[", "0", "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"v0b", "[", RowBox[{"[", RowBox[{"i", ",", "3"}], "]"}], "]"}], "-", " ", RowBox[{"vcomb", "[", RowBox[{"[", "3", "]"}], "]"}]}]}]}], "}"}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "nBodies"}], "}"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}]}], "Input", CellLabel->"In[144]:="], Cell[BoxData[{ RowBox[{ RowBox[{"eqsA", " ", "=", " ", RowBox[{ RowBox[{"Thread", "[", RowBox[{ RowBox[{"D", "[", RowBox[{"p", ",", "t"}], "]"}], " ", "==", " ", "dp"}], "]"}], "~", " ", "Join", " ", "~", " ", "p0a"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"eqsB", " ", "=", " ", RowBox[{ RowBox[{"Thread", "[", RowBox[{ RowBox[{"D", "[", RowBox[{"p", ",", "t"}], "]"}], " ", "==", " ", "dp"}], "]"}], "~", " ", "Join", " ", "~", " ", "p0b"}]}], ";"}]}], "Input", CellLabel->"In[146]:="], Cell[BoxData[{ RowBox[{ RowBox[{"solA", " ", "=", " ", RowBox[{ RowBox[{"NDSolve", "[", RowBox[{"eqsA", " ", ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "tMax"}], "}"}], ",", " ", RowBox[{"MaxSteps", " ", "\[Rule]", " ", "40000"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"solB", " ", "=", " ", RowBox[{ RowBox[{"NDSolve", "[", RowBox[{"eqsB", " ", ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "tMax"}], "}"}], ",", " ", RowBox[{"MaxSteps", " ", "\[Rule]", " ", "40000"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";"}]}], "Input", CellLabel->"In[148]:="], Cell[CellGroupData[{ Cell["Orbit Separation", "Subsubsection"], Cell[BoxData[ RowBox[{ RowBox[{"dr", "[", "t_", "]"}], " ", "=", " ", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], " ", "/.", " ", RowBox[{"i", "\[Rule]", "nBodies"}]}], " ", "/.", " ", "solB"}], ")"}], " ", "-", " ", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"y", "[", "i", "]"}], "[", "t", "]"}], ",", " ", RowBox[{ RowBox[{"z", "[", "i", "]"}], "[", "t", "]"}]}], "}"}], " ", "/.", " ", RowBox[{"i", "\[Rule]", "nBodies"}]}], " ", "/.", " ", "solA"}], ")"}]}]}]], "Input", CellLabel->"In[150]:="], Cell[BoxData[ RowBox[{"ParametricPlot3D", "[", RowBox[{ RowBox[{"dr", "[", "t", "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "tMax"}], "}"}], ",", RowBox[{"PlotPoints", "\[Rule]", "200"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input", CellLabel->"In[151]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{"Log", "[", RowBox[{"10", ",", SqrtBox[ RowBox[{ RowBox[{"dr", "[", "t", "]"}], ".", RowBox[{"dr", "[", "t", "]"}]}]]}], "]"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "tMax"}], "}"}], ",", RowBox[{"PlotPoints", "\[Rule]", "200"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\\""}], ",", RowBox[{"Frame", "\[Rule]", "True"}], ",", RowBox[{"PlotRange", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"-", "6"}], ",", "0"}], "}"}]}]}], "]"}]], "Input", CellLabel->"In[152]:="], Cell[TextData[{ "Two objects that begin nearby deviate ", StyleBox["expoentially", FontSlant->"Italic"], " fast!! (This is a measure of chaos.)" }], "Text", CellChangeTimes->{{3.443439127917616*^9, 3.443439177417303*^9}}] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Summary", "Section"], Cell[TextData[{ "We have shown how to use ", StyleBox["Mathematica", FontSlant->"Italic"], " to solve the nonlinear, coupled differential equations for many \ interacting objects in three dimensions. This is a very complicated problem. \ For ", StyleBox["n", FontSlant->"Italic"], " interacting objects, there are 6 \[Cross] ", StyleBox["n", FontSlant->"Italic"], " equations to solve. Additionally, computation of the force acting on any \ object grows as (", StyleBox["n", FontSlant->"Italic"], " - 1). This implies that computation of the interacting forces grows as ", StyleBox["n ", FontSlant->"Italic"], "\[Cross] (", StyleBox["n", FontSlant->"Italic"], " -1) ~ ", Cell[BoxData[ FormBox[ SuperscriptBox["n", "2"], TraditionalForm]]], ". Finally, these equations can not be solved in closed form, and they \ exhibit sensitivity to very small changes in the initial conditions. " }], "Text"], Cell[CellGroupData[{ Cell["Extra notes:", "Subsubsection"], Cell[TextData[{ "Schemes for many body problems with computations scaling like ", StyleBox["N", FontSlant->"Italic"], " log ", StyleBox["N", FontSlant->"Italic"], " are possible using a recursive \"cell averaging\" procedure. \nSay we \ divide ", StyleBox["N", FontSlant->"Italic"], " bodies into ", StyleBox["M", FontSlant->"Italic"], " cells with about ", StyleBox["N/M", FontSlant->"Italic"], " bodies per cell. We need to sort these bodies so that all of the cells are \ separated while the bodies within a cell are relatively close. Then, the \ acceleration of a body within a cell equals to the accelerations due to the \ other (", StyleBox["M", FontSlant->"Italic"], " - 1) cells plus the accelerations due to other bodies within the cell. The \ number of operations is roughly of order ", Cell[BoxData[ FormBox[ SuperscriptBox["M", "2"], TraditionalForm]]], "+ ", Cell[BoxData[ FormBox[ SuperscriptBox["N", "2"], TraditionalForm]]], "/", StyleBox["M", FontSlant->"Italic"], ". When ", StyleBox["M", FontSlant->"Italic"], " ~ ", Cell[BoxData[ FormBox[ SuperscriptBox["N", RowBox[{"2", "/", "3"}]], TraditionalForm]]], ", then the number of operations are minimized, saving considerable \ computational time." }], "Text"] }, Closed]] }, Closed]] }, Open ]] }, WindowSize->{683, 777}, WindowMargins->{{0, Automatic}, {Automatic, 0}}, PrintingCopies->1, PrintingPageRange->{1, Automatic}, ShowSelection->True, FrontEndVersion->"6.0 for Mac OS X x86 (32-bit) (May 21, 2008)", StyleDefinitions->"TutorialBook.nb" ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[590, 23, 38, 0, 51, "Title"], Cell[631, 25, 57, 3, 41, "Subsubtitle"], Cell[CellGroupData[{ Cell[713, 32, 31, 0, 86, "Section"], Cell[747, 34, 385, 6, 77, "Text"], Cell[1135, 42, 134, 3, 23, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[1306, 50, 44, 0, 54, "Section"], Cell[1353, 52, 53, 0, 23, "Text"], Cell[1409, 54, 251, 6, 41, "Text"], Cell[1663, 62, 472, 16, 43, "Input"], Cell[2138, 80, 425, 13, 42, "Input"], Cell[2566, 95, 338, 10, 26, "Input"], Cell[2907, 107, 67, 2, 24, "Input"], Cell[2977, 111, 221, 7, 24, "Input"], Cell[3201, 120, 346, 11, 28, "Input"], Cell[3550, 133, 67, 2, 24, "Input"], Cell[3620, 137, 626, 20, 85, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[4283, 162, 38, 0, 54, "Section"], Cell[CellGroupData[{ Cell[4346, 166, 37, 0, 34, "Subsection"], Cell[4386, 168, 62, 0, 23, "Text"], Cell[4451, 170, 323, 8, 26, "Input"], Cell[4777, 180, 497, 15, 26, "Input"], Cell[5277, 197, 926, 26, 61, "Input"], Cell[6206, 225, 139, 5, 23, "Text"], Cell[6348, 232, 52, 1, 24, "Input"], Cell[6403, 235, 191, 8, 23, "Text"], Cell[6597, 245, 911, 30, 64, "Input"], Cell[7511, 277, 1619, 54, 120, "Input"], Cell[9133, 333, 100, 2, 23, "Text"], Cell[9236, 337, 78, 2, 24, "Input"], Cell[9317, 341, 2171, 58, 166, "Input"], Cell[11491, 401, 49, 0, 23, "Text"], Cell[11543, 403, 1943, 57, 150, "Input"], Cell[13489, 462, 295, 9, 26, "Input"], Cell[13787, 473, 109, 5, 23, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[13933, 483, 38, 0, 34, "Subsection"], Cell[13974, 485, 692, 19, 63, "Input"], Cell[14669, 506, 637, 17, 59, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[15343, 528, 46, 0, 34, "Subsection"], Cell[15392, 530, 58, 0, 23, "Text"], Cell[15453, 532, 222, 7, 26, "Input"], Cell[15678, 541, 438, 15, 24, "Input"], Cell[16119, 558, 1848, 49, 153, "Input"], Cell[17970, 609, 118, 3, 26, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[18125, 617, 32, 0, 34, "Subsection"], Cell[18160, 619, 298, 9, 26, "Input"], Cell[18461, 630, 256, 6, 41, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[18754, 641, 30, 0, 34, "Subsection"], Cell[18787, 643, 72, 0, 23, "Text"], Cell[18862, 645, 997, 29, 81, "Input"], Cell[19862, 676, 129, 5, 23, "Text"], Cell[19994, 683, 178, 5, 24, "Input"], Cell[20175, 690, 89, 2, 26, "Input"], Cell[20267, 694, 273, 8, 26, "Input"], Cell[20543, 704, 51, 0, 23, "Text"], Cell[20597, 706, 844, 24, 64, "Input"], Cell[21444, 732, 337, 10, 24, "Input"], Cell[21784, 744, 91, 2, 26, "Input"], Cell[21878, 748, 538, 15, 63, "Input"], Cell[22419, 765, 76, 0, 23, "Text"], Cell[22498, 767, 990, 28, 79, "Input"], Cell[23491, 797, 515, 14, 45, "Input"], Cell[24009, 813, 172, 4, 24, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[24230, 823, 35, 0, 54, "Section"], Cell[24268, 825, 227, 4, 41, "Text"], Cell[24498, 831, 309, 8, 26, "Input"], Cell[CellGroupData[{ Cell[24832, 843, 38, 0, 34, "Subsection"], Cell[24873, 845, 492, 15, 26, "Input"], Cell[25368, 862, 276, 8, 26, "Input"], Cell[25647, 872, 523, 16, 26, "Input"], Cell[26173, 890, 276, 8, 26, "Input"], Cell[26452, 900, 2586, 74, 135, "Input"], Cell[29041, 976, 296, 9, 26, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[29374, 990, 38, 0, 34, "Subsection"], Cell[29415, 992, 672, 19, 63, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[30124, 1016, 32, 0, 34, "Subsection"], Cell[30159, 1018, 1743, 46, 153, "Input"], Cell[31905, 1066, 299, 9, 26, "Input"], Cell[32207, 1077, 747, 22, 63, "Input"], Cell[32957, 1101, 747, 22, 63, "Input"], Cell[33707, 1125, 109, 5, 23, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[33853, 1135, 30, 0, 34, "Subsection"], Cell[33886, 1137, 72, 0, 23, "Text"], Cell[33961, 1139, 997, 29, 81, "Input"], Cell[34961, 1170, 273, 8, 26, "Input"], Cell[35237, 1180, 844, 24, 64, "Input"], Cell[36084, 1206, 538, 15, 63, "Input"], Cell[36625, 1223, 990, 28, 79, "Input"], Cell[37618, 1253, 581, 15, 45, "Input"], Cell[38202, 1270, 72, 0, 23, "Text"], Cell[38277, 1272, 172, 4, 24, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[38498, 1282, 32, 0, 54, "Section"], Cell[38533, 1284, 129, 3, 23, "Text"], Cell[38665, 1289, 107, 3, 24, "Input"], Cell[38775, 1294, 443, 14, 26, "Input"], Cell[39221, 1310, 369, 11, 26, "Input"], Cell[39593, 1323, 585, 18, 44, "Input"], Cell[40181, 1343, 473, 15, 26, "Input"], Cell[40657, 1360, 300, 9, 26, "Input"], Cell[40960, 1371, 4960, 141, 261, "Input"], Cell[45923, 1514, 572, 17, 45, "Input"], Cell[46498, 1533, 754, 19, 112, "Input"], Cell[CellGroupData[{ Cell[47277, 1556, 41, 0, 33, "Subsubsection"], Cell[47321, 1558, 1012, 32, 62, "Input"], Cell[48336, 1592, 336, 9, 45, "Input"], Cell[48675, 1603, 668, 19, 81, "Input"], Cell[49346, 1624, 227, 6, 23, "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[49622, 1636, 26, 0, 54, "Section"], Cell[49651, 1638, 927, 28, 96, "Text"], Cell[CellGroupData[{ Cell[50603, 1670, 37, 0, 33, "Subsubsection"], Cell[50643, 1672, 1286, 45, 143, "Text"] }, Closed]] }, Closed]] }, Open ]] } ] *) (* End of internal cache information *)