(* 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[ 29386, 1006] NotebookOptionsPosition[ 26283, 905] NotebookOutlinePosition[ 26670, 922] CellTagsIndexPosition[ 26627, 919] WindowFrame->Normal ContainsDynamic->False*) (* Beginning of Notebook Content *) Notebook[{ Cell[CellGroupData[{ Cell["Solutions to Execise 3", "Title"], Cell[TextData[{ "Assigned: January 31, 2006\nDue: (", StyleBox["Thursday", FontColor->RGBColor[1, 0, 0]], ") February 9, 2006" }], "Subsubtitle"], Cell[CellGroupData[{ Cell["Setup your Solution Notebook (in-class)", "Section"], Cell["\<\ In order to solve these execises, you need need to define the equations of \ motion and a particular solution. Copy the equations of motion for an \ orbiting planet and demonstrate that your equations are correct by \ recomputing and plotting Halley's orbit.\ \>", "Text"], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell[BoxData[ RowBox[{ RowBox[{"p", " ", "=", " ", RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", " ", RowBox[{"y", "[", "t", "]"}], ",", " ", RowBox[{"vx", "[", "t", "]"}], ",", " ", RowBox[{"vy", "[", "t", "]"}]}], "}"}]}], ";"}]], "Input"], Cell[BoxData[ RowBox[{"dpdt", " ", "=", " ", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", "t"}], "]"}], " ", "\[Equal]", " ", RowBox[{"vx", "[", "t", "]"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"y", "[", "t", "]"}], ",", "t"}], "]"}], " ", "\[Equal]", " ", RowBox[{"vy", "[", "t", "]"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"vx", "[", "t", "]"}], ",", "t"}], "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"-", "4"}], " ", SuperscriptBox["\[Pi]", "2"], " ", RowBox[{ RowBox[{"x", "[", "t", "]"}], "/", RowBox[{"r", "^", "3"}]}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"D", "[", RowBox[{ RowBox[{"vy", "[", "t", "]"}], ",", "t"}], "]"}], " ", "\[Equal]", " ", RowBox[{ RowBox[{"-", "4"}], " ", SuperscriptBox["\[Pi]", "2"], " ", RowBox[{ RowBox[{"y", "[", "t", "]"}], "/", RowBox[{"r", "^", "3"}]}]}]}]}], "}"}], " ", "/.", " ", RowBox[{"r", " ", "\[Rule]", " ", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"x", "[", "t", "]"}], "^", "2"}], " ", "+", " ", RowBox[{ RowBox[{"y", "[", "t", "]"}], "^", "2"}]}], "]"}]}]}]}]], "Input"], Cell[BoxData[ RowBox[{"p0", " ", "=", " ", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"x", "[", "0", "]"}], " ", "\[Equal]", " ", "x0"}], ",", " ", RowBox[{ RowBox[{"y", "[", "0", "]"}], " ", "\[Equal]", " ", "y0"}], ",", RowBox[{ RowBox[{"vx", "[", "0", "]"}], " ", "\[Equal]", " ", "vx0"}], ",", " ", RowBox[{ RowBox[{"vy", "[", "0", "]"}], " ", "\[Equal]", " ", "vy0"}]}], "}"}]}]], "Input"], Cell[BoxData[ RowBox[{"eqs", " ", "=", " ", RowBox[{"dpdt", " ", "~", " ", "Join", " ", "~", "p0"}]}]], "Input"], Cell[BoxData[ RowBox[{"halley", " ", "=", " ", RowBox[{ RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", " ", "0.59"}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", " ", RowBox[{"vy0", "->", "11.46"}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "76"}], "}"}], ",", RowBox[{"MaxSteps", " ", "->", " ", "2000"}]}], "]"}], " ", "//", " ", "Flatten"}]}]], "Input"], Cell[BoxData[ RowBox[{"ParametricPlot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"x", "[", "t", "]"}], ",", RowBox[{"y", "[", "t", "]"}]}], "}"}], "/.", "\[InvisibleSpace]", "halley"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "76"}], "}"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}]], "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Conservation of Energy (in-class)", "Section"], Cell[TextData[{ "Planetary motion conserves the total energy as well as the angular \ momentum. The total energy is the sum of the kinetic energy and the \ gravitational potential energy.\n\t", StyleBox["Energy ", FontSlant->"Italic"], " = ", Cell[BoxData[ FormBox[ FractionBox["1", "2"], TraditionalForm]]], " ", Cell[BoxData[ FormBox[ SuperscriptBox["v", "2"], TraditionalForm]]], " - 4", Cell[BoxData[ FormBox[ SuperscriptBox["\[Pi]", "2"], TraditionalForm]]], "/ ", StyleBox["r", FontSlant->"Italic"], "\nShow that the total energy of an orbiting body is constant." }], "Text"], Cell[CellGroupData[{ Cell["Hint", "Subsubsection"], Cell["\<\ For this problem, you need to define a function that returns the total energy \ for Halley's comet as a function of time.\ \>", "Text"], Cell[TextData[{ "Start by defining (and testing!) a function of the form:\n\t", StyleBox["energyHalley[t_] := (...) /. halley", FontWeight->"Bold"] }], "Text"], Cell[TextData[{ "This function uses the solution to the differential equations (labeled by \ \"", StyleBox["halley", FontWeight->"Bold"], "\") to substitute for your expressions within the parentheses." }], "Text"], Cell[TextData[{ "Once your function has been defined, you can plot the error as a function \ of time by writing...\n\t", StyleBox["Plot[100 (energyHalley[t]-energyHalley[0])/energyHalley[0], {t, 0, \ 76}, \n \t\tPlotLabel->\"Percent Error of Total Energy\"]", FontWeight->"Bold"] }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell[BoxData[ RowBox[{ RowBox[{"energyHalley", "[", "t_", "]"}], " ", "=", " ", RowBox[{ RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"vx", "[", "t", "]"}], "*", RowBox[{"vx", "[", "t", "]"}]}], " ", "+", " ", RowBox[{ RowBox[{"vy", "[", "t", "]"}], "*", RowBox[{"vy", "[", "t", "]"}]}]}], ")"}], "/", "2"}], " ", "-", " ", RowBox[{"4", RowBox[{ SuperscriptBox["\[Pi]", "2"], "/", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"x", "[", "t", "]"}], "*", RowBox[{"x", "[", "t", "]"}]}], " ", "+", " ", RowBox[{ RowBox[{"y", "[", "t", "]"}], "*", RowBox[{"y", "[", "t", "]"}]}]}], "]"}]}]}]}], " ", "/.", " ", "halley"}]}]], "Input"], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ FractionBox[ RowBox[{"100", " ", RowBox[{"(", RowBox[{ RowBox[{"energyHalley", "[", "t", "]"}], "-", RowBox[{"energyHalley", "[", "0", "]"}]}], ")"}]}], RowBox[{"energyHalley", "[", "0", "]"}]], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "76"}], "}"}], ",", RowBox[{ "PlotLabel", "\[Rule]", "\"\\""}]}], "]"}]], "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["The Planets (in-class)", "Section"], Cell["\<\ The orbital parameters of the planets of our solar system are\ \>", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"pname", " ", "=", " ", RowBox[{"{", RowBox[{ "\"\\"", ",", " ", "\"\\"", ",", " ", "\"\\"", ",", " ", "\"\\"", ",", " ", "\"\\"", ",", " ", "\"\\"", ",", " ", "\"\\"", ",", " ", "\"\\"", ",", " ", "\"\\""}], "}"}]}], ";"}], "\n", RowBox[{ RowBox[{"pradius", " ", "=", " ", RowBox[{"{", RowBox[{ "0.39", ",", " ", "0.72", ",", " ", "1.00", ",", " ", "1.52", ",", " ", "5.20", ",", " ", "9.54", ",", " ", "19.19", ",", " ", "30.06", ",", " ", "39.53"}], "}"}]}], ";"}], "\n", RowBox[{ RowBox[{"peccentricity", " ", "=", " ", RowBox[{"{", RowBox[{ "0.206", ",", " ", "0.007", ",", " ", "0.017", ",", " ", "0.093", ",", " ", "0.048", ",", " ", "0.056", ",", " ", "0.046", ",", " ", "0.010", ",", " ", "0.248"}], "}"}]}], ";"}]}], "Input"], Cell[BoxData[ RowBox[{"TableForm", "[", RowBox[{ RowBox[{"{", RowBox[{"pname", ",", "pradius", ",", " ", "peccentricity"}], "}"}], ",", RowBox[{"TableDirections", "->", RowBox[{"{", RowBox[{"Row", ",", "Column"}], "}"}]}], ",", RowBox[{"TableHeadings", "->", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{ "\"\< \>\"", ",", " ", "\"\\"", ",", "\"\\""}], "}"}], ",", " ", "Automatic"}], "}"}]}]}], "]"}]], "Input"], Cell[TextData[{ "Find the harmonic constant, ", Cell[BoxData[ FormBox[ SuperscriptBox["T", "2"], TraditionalForm]]], "/", Cell[BoxData[ FormBox[ SuperscriptBox["a", "3"], TraditionalForm]]], ", for each of the nine planets. " }], "Text"], Cell[TextData[{ "In order to solve this problem, you need to know the initial conditions for \ each planet. As Kepler showed, the initial conditions determine the geometry \ of the orbit and ", StyleBox["vice-versa", FontSlant->"Italic"], ". " }], "Text"], Cell[TextData[{ "The function that determines the velocity of an object orbiting the sun at \ the ", StyleBox["closest distance to the sun", FontSlant->"Italic"], " is\n\t", StyleBox["vc[a_, e_] := ", FontWeight->"Bold"], Cell[BoxData[ RowBox[{"2", "\[Pi]", " ", SqrtBox[ StyleBox[ FractionBox[ RowBox[{"1", "+", "e"}], RowBox[{"a", " ", "-", " ", RowBox[{"a", " ", "e"}]}]], FontSize->14]], " "}]], FontWeight->"Bold"] }], "Text"], Cell[TextData[{ "where \"", StyleBox["a", FontWeight->"Bold"], "\" is the semimajor radius and \"", StyleBox["e", FontWeight->"Bold"], "\" is the eccentricity." }], "Text"], Cell[CellGroupData[{ Cell["Hint", "Subsubsection"], Cell[TextData[{ "For simplifity, take the planet's mean radius to be equal to the planet's \ semimajor radius, ", StyleBox["a. ", FontSlant->"Italic"], "Then, find the planet's orbital period (including the planet's \ eccentricity). Remember, in our astronomical units, ", Cell[BoxData[ FormBox[ SuperscriptBox["T", "2"], TraditionalForm]]], "/", Cell[BoxData[ FormBox[ SuperscriptBox["a", "3"], TraditionalForm]]], " = 1." }], "Text"], Cell[TextData[{ "You may solve this problem in many ways, but one that I like begins by \ finding the period to orbit from one side of the sun to the other. This is \ one-half of the period, ", StyleBox["T", FontSlant->"Italic"], ". " }], "Text"], Cell[TextData[{ "Let's use a ", StyleBox["Module[...]", FontWeight->"Bold"], " and define a local variable, sol, that is the solution to the equations of \ motion. Then, use FindRoot[...], to find the time when the orbit crosses the \ x-axis. For example,\n\t", StyleBox["findT[a_,e_] :=Module[{sol},\n \t\tsol = NDSolve[...] // \ Flatten;\n \t\t2t /. FindRoot[Evaluate[vx[t]/.sol],{t, a^(3/2)/2}]]", FontWeight->"Bold"] }], "Text"], Cell[TextData[{ "With this function, the orbital periods of the planets in years is obtained \ with...\n\t", StyleBox["pT = Table[findT[pradius[[i]], peccentricity[[i]]], {i, 1, \ 9}]", FontWeight->"Bold"] }], "Text"], Cell[TextData[{ "Further proof of Kepler's Third Law is found with...\n\t", Cell[BoxData[ RowBox[{"pHarmonic", " ", "=", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ SuperscriptBox[ RowBox[{"pT", "[", RowBox[{"[", "i", "]"}], "]"}], "2"], "/", SuperscriptBox[ RowBox[{"pradius", "[", RowBox[{"[", "i", "]"}], "]"}], "3"]}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", "9"}], "}"}]}], "]"}]}]], FontWeight->"Bold"] }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell["\<\ First, we define the initial condition from the geometry of the orbit.\ \>", "Text"], Cell[BoxData[ RowBox[{ StyleBox[ RowBox[{"vc", "[", RowBox[{"a_", ",", " ", "e_"}], "]"}], FontWeight->"Bold"], StyleBox[" ", FontWeight->"Bold"], StyleBox[":=", FontWeight->"Bold"], StyleBox[" ", FontWeight->"Bold"], RowBox[{"2", "\[Pi]", " ", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"(", RowBox[{"1", "+", "e"}], ")"}], "/", RowBox[{"(", RowBox[{"a", " ", "-", " ", RowBox[{"a", " ", "e"}]}], ")"}]}], "]"}]}]}]], "Input"], Cell[TextData[{ "Next, we'll take the planet's mean radius to be equal to the planet's \ semimajor radius, ", StyleBox["a. ", FontSlant->"Italic"], "Then, we must find the planet's orbital period (including the planet's \ eccentricity). We note in our astronomical units, ", Cell[BoxData[ FormBox[ SuperscriptBox["T", "2"], TraditionalForm]]], "/", Cell[BoxData[ FormBox[ SuperscriptBox["a", "3"], TraditionalForm]]], " = 1." }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"findT", "[", RowBox[{"a_", ",", "e_"}], "]"}], " ", ":=", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", RowBox[{"sol", ",", " ", "tGuess"}], "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"tGuess", " ", "=", " ", RowBox[{"Sqrt", "[", RowBox[{"a", "^", "3"}], "]"}]}], ";", RowBox[{"sol", " ", "=", RowBox[{ RowBox[{"NDSolve", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"eqs", " ", "/.", RowBox[{"x0", " ", "->", RowBox[{"a", " ", RowBox[{"(", RowBox[{"1", " ", "-", " ", "e"}], ")"}]}]}]}], " ", "/.", " ", RowBox[{"y0", " ", "->", " ", "0"}]}], " ", "/.", RowBox[{"vy0", "->", RowBox[{"vc", "[", RowBox[{"a", ",", "e"}], "]"}]}]}], "/.", " ", RowBox[{"vx0", " ", "->", " ", "0"}]}], ",", " ", "p", ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "tGuess"}], "}"}], ",", RowBox[{"MaxSteps", " ", "->", " ", "2000"}]}], "]"}], " ", "//", " ", "Flatten"}]}], ";", RowBox[{ RowBox[{"2", "t"}], " ", "/.", " ", RowBox[{"FindRoot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"vx", "[", "t", "]"}], "/.", "sol"}], "]"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", RowBox[{"tGuess", "/", "2"}]}], "}"}]}], "]"}]}]}]}], "]"}]}]], "Input"], Cell["The orbital periods of the planets in years...", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"(", RowBox[{"pT", "=", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{"findT", "[", RowBox[{ RowBox[{"pradius", "[", RowBox[{"[", "i", "]"}], "]"}], ",", RowBox[{"peccentricity", "[", RowBox[{"[", "i", "]"}], "]"}]}], "]"}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "9"}], "}"}]}], "]"}]}], " ", ")"}], "//", " ", "TableForm"}]], "Input"], Cell["Further proof of Kepler's Third Law:", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"(", RowBox[{"pHarmonic", " ", "=", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ SuperscriptBox[ RowBox[{"pT", "[", RowBox[{"[", "i", "]"}], "]"}], "2"], "/", SuperscriptBox[ RowBox[{"pradius", "[", RowBox[{"[", "i", "]"}], "]"}], "3"]}], ",", RowBox[{"{", RowBox[{"i", ",", "1", ",", "9"}], "}"}]}], "]"}]}], " ", ")"}], "//", " ", "TableForm"}]], "Input"], Cell[BoxData[ RowBox[{"TableForm", "[", RowBox[{ RowBox[{"{", RowBox[{ "pname", ",", "pradius", ",", " ", "peccentricity", ",", " ", "pT", ",", " ", "pHarmonic"}], "}"}], ",", RowBox[{"TableDirections", "\[Rule]", " ", RowBox[{"{", RowBox[{"Row", ",", "Column"}], "}"}]}], ",", RowBox[{"TableHeadings", "\[Rule]", " ", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{ "\"\< \>\"", ",", " ", "\"\\"", ",", "\"\\"", ",", " ", "\"\\"", ",", " ", "\"\\""}], "}"}], ",", " ", "Automatic"}], "}"}]}]}], "]"}]], "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Using Euler-Cromer", "Section"], Cell[TextData[{ "Write a procedure to integrate the equations of motion for an orbiting \ object using the Euler-Cromer method. Use your procedure to integrate the \ orbits for the planet Earth and the comet Halley. Examine your solution for a \ few values of the time step, ", StyleBox["dt", FontSlant->"Italic"], ", and discuss the accuracy of this method (with a fixed and constant time \ step) to the adaptive time-step method used by ", StyleBox["NDSolve[\[Ellipsis]]", FontWeight->"Bold"], ". " }], "Text"], Cell[CellGroupData[{ Cell["Hint", "Subsubsection"], Cell[TextData[{ "Define your solution on the (", StyleBox["x, y", FontSlant->"Italic"], ") plane. This implies that you have four unknown functions of time: ", StyleBox["x", FontSlant->"Italic"], "(", StyleBox["t", FontSlant->"Italic"], "), ", StyleBox["y", FontSlant->"Italic"], "(", StyleBox["t", FontSlant->"Italic"], "), ", StyleBox["vx", FontSlant->"Italic"], "(", StyleBox["t", FontSlant->"Italic"], "), and ", StyleBox["vy", FontSlant->"Italic"], "(", StyleBox["t", FontSlant->"Italic"], "). Your procedure needs to map a list of four numbers, ", StyleBox["{x_, y_, vx_, vy_}", FontWeight->"Bold"], ", to four new numbers corresponding to the new orbital position of the of \ the object. Your procedure will look something like:\n\n\t", StyleBox["orbitStep[{x_, y_, vx_, vy_}] := Module[{vxNext, vyNext, \ \[Ellipsis]},\n\t\t\[Ellipsis];\n\t\t{list of four nunbers}]\n\n", FontWeight->"Bold"], "Of course, you also need to remember the form of the acceleration due to \ gravity. You must define the radius between the star and the obit, ", StyleBox["r = ", FontSlant->"Italic"], Cell[BoxData[ FormBox[ SqrtBox[ RowBox[{ SuperscriptBox["x", "2"], " ", "+", " ", SuperscriptBox["y", "2"]}]], TraditionalForm]]], ", and use the equations of motion discussed in class. \n\nFinally, remember \ the initial conditions for Earth are (", StyleBox["x", FontSlant->"Italic"], "(0), ", StyleBox["y", FontSlant->"Italic"], "(0), ", StyleBox["vx", FontSlant->"Italic"], "(0), ", StyleBox["vy", FontSlant->"Italic"], "(0)) = (1.0, 0.0, 0.0, 2\[Pi]), and the initial conditions for Halley are \ ", "(", StyleBox["x", FontSlant->"Italic"], "(0), ", StyleBox["y", FontSlant->"Italic"], "(0), ", StyleBox["vx", FontSlant->"Italic"], "(0), ", StyleBox["vy", FontSlant->"Italic"], "(0)) = (0.59, 0.0, 0.0, 11.46).\n\nGood values of ", StyleBox["dt", FontSlant->"Italic"], " to try are ", StyleBox["dt", FontSlant->"Italic"], " = 0.001 and ", StyleBox["dt", FontSlant->"Italic"], " = 0.01, but be careful: when ", StyleBox["dt", FontSlant->"Italic"], " is small, it takes a ", StyleBox["very long time", FontSlant->"Italic"], " (about 1 minute) to compute the 100,000 steps needed to represent Halley's \ orbit!" }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Solution", "Subsubsection"], Cell[TextData[{ "The global variable \"", StyleBox["dt", FontWeight->"Bold", FontColor->RGBColor[1, 0, 1]], "\" must be set to a numerical value in order for ", StyleBox["orbitStep[...]", FontWeight->"Bold"], " to work." }], "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"fourPi2", " ", "=", " ", RowBox[{"N", "[", RowBox[{"4", " ", RowBox[{"\[Pi]", "^", "2"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"orbitStep", "[", RowBox[{"{", RowBox[{"x_", ",", " ", "y_", ",", " ", "vx_", ",", " ", "vy_"}], "}"}], "]"}], " ", ":=", " ", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", RowBox[{"vxNext", ",", "vyNext", ",", " ", "r3"}], "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"r3", " ", "=", " ", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"x", "^", "2"}], " ", "+", " ", RowBox[{"y", "^", "2"}]}], ")"}], "^", RowBox[{"(", RowBox[{"3", "/", "2"}], ")"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"vxNext", " ", "=", " ", RowBox[{"N", "[", RowBox[{"vx", " ", "-", " ", RowBox[{ StyleBox["dt", FontColor->RGBColor[1, 0, 1]], " ", "fourPi2", " ", RowBox[{"x", "/", "r3"}]}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"vyNext", " ", "=", " ", RowBox[{"N", "[", RowBox[{"vy", " ", "-", " ", RowBox[{ StyleBox["dt", FontColor->RGBColor[1, 0, 1]], " ", "fourPi2", " ", RowBox[{"y", "/", "r3"}]}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"{", RowBox[{ RowBox[{"N", "[", RowBox[{"x", " ", "+", " ", RowBox[{ StyleBox["dt", FontColor->RGBColor[1, 0, 1]], " ", "vxNext"}]}], "]"}], ",", " ", RowBox[{"N", "[", RowBox[{"y", " ", "+", " ", RowBox[{ StyleBox["dt", FontColor->RGBColor[1, 0, 1]], " ", "vyNext"}]}], "]"}], ",", " ", "vxNext", ",", " ", "vyNext"}], "}"}]}]}], "]"}]}]}], "Input"], Cell[BoxData[{ RowBox[{ RowBox[{"initialEarth", " ", "=", " ", RowBox[{"N", "[", RowBox[{"{", RowBox[{"1.0", ",", " ", "0.0", ",", " ", "0.0", ",", " ", RowBox[{"2", " ", "\[Pi]"}]}], "}"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"initialHalley", " ", "=", " ", RowBox[{"{", RowBox[{"0.59", ",", " ", "0.0", ",", " ", "0.0", ",", " ", "11.46"}], "}"}]}], ";"}]}], "Input"], Cell[BoxData[{ RowBox[{ RowBox[{"dt", "=", "0.01`"}], ";", RowBox[{"stepsEarth", "=", RowBox[{"NestList", "[", RowBox[{"orbitStep", ",", "initialEarth", ",", "500"}], "]"}]}], ";"}], "\n", RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"Take", "[", RowBox[{"stepsEarth", ",", "500", ",", "2"}], "]"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}], "\n", RowBox[{"ListPlot", "[", RowBox[{"First", "[", RowBox[{"Transpose", "[", "stepsEarth", "]"}], "]"}], "]"}]}], "Input"], Cell[BoxData[{ RowBox[{ RowBox[{"dt", "=", "0.01`"}], ";", RowBox[{"stepsHalley", "=", RowBox[{"NestList", "[", RowBox[{"orbitStep", ",", "initialHalley", ",", "10000"}], "]"}]}], ";"}], "\n", RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"stepsHalley", "\[LeftDoubleBracket]", RowBox[{"All", ",", RowBox[{"{", RowBox[{"1", ",", "2"}], "}"}]}], "\[RightDoubleBracket]"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}], "\n", RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"stepsHalley", "\[LeftDoubleBracket]", RowBox[{"All", ",", "1"}], "\[RightDoubleBracket]"}], ",", RowBox[{"AxesLabel", "\[Rule]", RowBox[{"{", RowBox[{"\"\\"", ",", "\"\\""}], "}"}]}]}], "]"}]}], "Input"], Cell["Clearly, this is wrong! Halley's orbit can not be tilted!", "Text"], Cell[BoxData[{ RowBox[{ RowBox[{"dt", "=", "0.001`"}], ";", RowBox[{"stepsHalley", "=", RowBox[{"NestList", "[", RowBox[{"orbitStep", ",", "initialHalley", ",", "100000"}], "]"}]}], ";"}], "\n", RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"stepsHalley", "\[LeftDoubleBracket]", RowBox[{"All", ",", RowBox[{"{", RowBox[{"1", ",", "2"}], "}"}]}], "\[RightDoubleBracket]"}], ",", RowBox[{"AspectRatio", "\[Rule]", "Automatic"}]}], "]"}], "\n", RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"stepsHalley", "\[LeftDoubleBracket]", RowBox[{"All", ",", "1"}], "\[RightDoubleBracket]"}], ",", RowBox[{"AxesLabel", "\[Rule]", RowBox[{"{", RowBox[{"\"\\"", ",", "\"\\""}], "}"}]}]}], "]"}]}], "Input"], Cell["This looks better. We can compare below...", "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Comparing Solutions", "Subsubsection"], Cell[BoxData[{ RowBox[{ RowBox[{"dt", " ", "=", " ", "0.001"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"nSkip", " ", "=", " ", "100"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"ecStepsX", " ", "=", " ", RowBox[{"Flatten", "[", " ", RowBox[{"Take", "[", RowBox[{"stepsHalley", ",", " ", RowBox[{"{", RowBox[{"1", ",", "100000", ",", "nSkip"}], "}"}], ",", " ", "1"}], "]"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ StyleBox["difference", FontColor->RGBColor[1, 0, 1]], " ", "=", " ", RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"dt", " ", "i", " ", "nSkip"}], ",", RowBox[{ RowBox[{ "ecStepsX", "\[LeftDoubleBracket]", "i", "\[RightDoubleBracket]"}], " ", "-", " ", RowBox[{"(", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"x", "[", "t", "]"}], " ", "/.", " ", "halley"}], "]"}], " ", "/.", " ", RowBox[{"t", "\[Rule]", " ", RowBox[{"dt", " ", "i", " ", "nSkip"}]}]}], ")"}]}]}], "}"}], ",", " ", RowBox[{"{", RowBox[{"i", ",", " ", "1", ",", " ", "700"}], "}"}]}], "]"}]}], ";"}]}], "Input"], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{"difference", ",", RowBox[{ "PlotLabel", "\[Rule]", "\"\\""}], ",", RowBox[{"AxesLabel", "\[Rule]", RowBox[{"{", RowBox[{"\"\\"", ",", "\"\\""}], "}"}]}]}], "]"}]], "Input"] }, Closed]] }, Closed]] }, Open ]] }, WindowSize->{562, 744}, WindowMargins->{{1, Automatic}, {Automatic, 1}}, PrintingCopies->1, PrintingPageRange->{1, Automatic}, FrontEndVersion->"6.0 for Mac OS X x86 (32-bit) (May 21, 2008)", StyleDefinitions->"Default.nb" ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[590, 23, 39, 0, 73, "Title"], Cell[632, 25, 150, 5, 45, "Subsubtitle"], Cell[CellGroupData[{ Cell[807, 34, 58, 0, 69, "Section"], Cell[868, 36, 282, 5, 59, "Text"], Cell[CellGroupData[{ Cell[1175, 45, 33, 0, 24, "Subsubsection"], Cell[1211, 47, 287, 8, 28, "Input"], Cell[1501, 57, 1475, 45, 104, "Input"], Cell[2979, 104, 441, 12, 28, "Input"], Cell[3423, 118, 116, 2, 28, "Input"], Cell[3542, 122, 626, 16, 63, "Input"], Cell[4171, 140, 421, 12, 46, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[4641, 158, 52, 0, 39, "Section"], Cell[4696, 160, 611, 22, 84, "Text"], Cell[CellGroupData[{ Cell[5332, 186, 29, 0, 24, "Subsubsection"], Cell[5364, 188, 145, 3, 46, "Text"], Cell[5512, 193, 163, 4, 48, "Text"], Cell[5678, 199, 219, 6, 48, "Text"], Cell[5900, 207, 296, 6, 100, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[6233, 218, 33, 0, 18, "Subsubsection"], Cell[6269, 220, 809, 25, 66, "Input"], Cell[7081, 247, 462, 14, 73, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[7592, 267, 41, 0, 39, "Section"], Cell[7636, 269, 85, 2, 27, "Text"], Cell[7724, 273, 938, 22, 114, "Input"], Cell[8665, 297, 509, 15, 80, "Input"], Cell[9177, 314, 252, 10, 30, "Text"], Cell[9432, 326, 260, 7, 43, "Text"], Cell[9695, 335, 483, 18, 85, "Text"], Cell[10181, 355, 181, 8, 27, "Text"], Cell[CellGroupData[{ Cell[10387, 367, 29, 0, 24, "Subsubsection"], Cell[10419, 369, 455, 15, 64, "Text"], Cell[10877, 386, 251, 7, 48, "Text"], Cell[11131, 395, 448, 10, 118, "Text"], Cell[11582, 407, 225, 6, 48, "Text"], Cell[11810, 415, 505, 16, 60, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[12352, 436, 33, 0, 18, "Subsubsection"], Cell[12388, 438, 94, 2, 27, "Text"], Cell[12485, 442, 499, 19, 28, "Input"], Cell[12987, 463, 457, 15, 62, "Text"], Cell[13447, 480, 1587, 44, 131, "Input"], Cell[15037, 526, 62, 0, 27, "Text"], Cell[15102, 528, 475, 14, 46, "Input"], Cell[15580, 544, 52, 0, 27, "Text"], Cell[15635, 546, 481, 15, 50, "Input"], Cell[16119, 563, 640, 17, 97, "Input"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[16808, 586, 37, 0, 39, "Section"], Cell[16848, 588, 522, 12, 91, "Text"], Cell[CellGroupData[{ Cell[17395, 604, 29, 0, 24, "Subsubsection"], Cell[17427, 606, 2333, 90, 343, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[19797, 701, 33, 0, 18, "Subsubsection"], Cell[19833, 703, 243, 9, 27, "Text"], Cell[20079, 714, 1836, 52, 131, "Input"], Cell[21918, 768, 435, 12, 46, "Input"], Cell[22356, 782, 528, 14, 80, "Input"], Cell[22887, 798, 777, 21, 80, "Input"], Cell[23667, 821, 73, 0, 27, "Text"], Cell[23743, 823, 779, 21, 80, "Input"], Cell[24525, 846, 58, 0, 27, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[24620, 851, 44, 0, 24, "Subsubsection"], Cell[24667, 853, 1275, 37, 165, "Input"], Cell[25945, 892, 298, 8, 63, "Input"] }, Closed]] }, Closed]] }, Open ]] } ] *) (* End of internal cache information *)