(* 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[ 12249, 452] NotebookOptionsPosition[ 10613, 393] NotebookOutlinePosition[ 11003, 410] CellTagsIndexPosition[ 10960, 407] WindowFrame->Normal ContainsDynamic->False*) (* Beginning of Notebook Content *) Notebook[{ Cell[CellGroupData[{ Cell["Execises (3)", "Title"], Cell[TextData[{ "Assigned: February 12, 2009\nDue: (", StyleBox["Thursday", FontColor->RGBColor[1, 0, 0]], ") February 19, 2009" }], "Subsubtitle", CellChangeTimes->{{3.443437892000667*^9, 3.443437908796715*^9}}], 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"] }, 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]] }, 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]] }, 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]] }, Closed]] }, Open ]] }, WindowSize->{757, 803}, WindowMargins->{{102, Automatic}, {Automatic, 62}}, 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, 29, 0, 73, "Title"], Cell[622, 25, 218, 6, 45, "Subsubtitle"], Cell[CellGroupData[{ Cell[865, 35, 58, 0, 69, "Section"], Cell[926, 37, 282, 5, 43, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[1245, 47, 52, 0, 39, "Section"], Cell[1300, 49, 611, 22, 84, "Text"], Cell[CellGroupData[{ Cell[1936, 75, 29, 0, 24, "Subsubsection"], Cell[1968, 77, 145, 3, 27, "Text"], Cell[2116, 82, 163, 4, 43, "Text"], Cell[2282, 88, 219, 6, 43, "Text"], Cell[2504, 96, 296, 6, 59, "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[2849, 108, 41, 0, 39, "Section"], Cell[2893, 110, 85, 2, 27, "Text"], Cell[2981, 114, 938, 22, 80, "Input"], Cell[3922, 138, 509, 15, 46, "Input"], Cell[4434, 155, 252, 10, 30, "Text"], Cell[4689, 167, 260, 7, 43, "Text"], Cell[4952, 176, 483, 18, 69, "Text"], Cell[5438, 196, 181, 8, 27, "Text"], Cell[CellGroupData[{ Cell[5644, 208, 29, 0, 24, "Subsubsection"], Cell[5676, 210, 455, 15, 46, "Text"], Cell[6134, 227, 251, 7, 43, "Text"], Cell[6388, 236, 448, 10, 91, "Text"], Cell[6839, 248, 225, 6, 43, "Text"], Cell[7067, 256, 505, 16, 48, "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[7621, 278, 37, 0, 39, "Section"], Cell[7661, 280, 522, 12, 59, "Text"], Cell[CellGroupData[{ Cell[8208, 296, 29, 0, 24, "Subsubsection"], Cell[8240, 298, 2333, 90, 278, "Text"] }, Closed]] }, Closed]] }, Open ]] } ] *) (* End of internal cache information *)