(* 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[ 37498, 1198] NotebookOptionsPosition[ 33676, 1080] NotebookOutlinePosition[ 34091, 1098] CellTagsIndexPosition[ 34048, 1095] WindowFrame->Normal ContainsDynamic->False*) (* Beginning of Notebook Content *) Notebook[{ Cell[CellGroupData[{ Cell["Pendulum", "Title"], Cell["\<\ APAM 1601 Columbia University\ \>", "Subsubtitle"], Cell[CellGroupData[{ Cell["Introduction", "Section"], Cell["\<\ The equation of motion for a simple pendulum is a good way to illustrate how \ a computer can solve a differential equation.\ \>", "Text"], Cell[GraphicsData["CompressedBitmap", "\<\ eJztHWuPHEdxbmdn9u7Wj7skkJxxYp8JyZHgEHuFICAUAVECEoG8DCRIgAlZ 4kghwTGY91cbi/cnL+93cokQIB5RzvkJ9xfuL/gvDF3VVTU1M9Uzs+d7Ec9K tz1bXV2vrq6u6e6Ze/j0uWefef70uTNPnz760NnTLz575umXjj74wlkHimei aGbJ/d19NILrLIr4Cz+H4It+0PWTWKTZxsZG9sADD2TUIFteXs4uXLiwprEG 2cmTJ7NTp04hgatXr+I1YJ8/fz56Q6MmmTClj0OPegAt4QEJYM6f9fV1A6+P 7BwbLC9fvowoUQyX3DJ63bd4iigDVVYHuKytrUVrGiVGYlANwr2pq1Ks0gYZ j8cAe6OIxRiA7czFvDQWStqTqwKBvhCAxo4IqwWaslr/rcoM1QsLC5GF/nq1 k0PXnxep3mwLQCn6ALgSAHyC5FQOEB0wkZAuViXhqhSu1lQVUt4frtoHV4Ny FTIfwtWs6oWPC5n5cNV+IX1FVQ2ElzZPqVWlKhG9dNUnDZNVSO8w0sMG0nD7 kD5lIM3vMtKnDaS5XUZ6xECa3WWkRw2kwS4jPWYgpduH9LiBlOwy0hMGUn+X kU4ZSPF0SLdRDwQylsvtadTnM0eojUtIkE/Gw9HBoB3A4deZ9gxrsqLbFDdA cnmND0iU4wAMUiPi9hmDW6/CLZhb3UrNOUVydvMzgINxGgPIz7bn1ipDS+k7 ooQKSpfHQGoGNX1q4rn6zLXXAIOugFQoIdK6G+EPus6pQl3hnKSVTkfIglYO Oaf6ErqLu8wnZEyBMrWcNEvsa3FkoN5kfy//QGQn95LuWF1d9dRVFwHsq+3V Kea7Q4Ip/49KuWpsXJHkLKbTICFCMACJPUmaa8NeyAq7XmHJP0vuUxoaSu5g Ms4DhKscuEl67zTgc/ArJQL0W4kei6hOBYGxv8KvcXvxA3cJs9QUupqN5waQ JXdFAx4KrAEMAx61MdFlrcgPsXlbo5c8YVbzZ9cUy0wrLHsFBDkWVsHKwn5O WZ8+kbJvq5soHl5sEz1wySOU9yMvpzQRlE6DXlL+w/pwgIZGCQ0SHlbkMu1U ab678yzn0DtZF46k8IGBBtJKMHJxHWAuVKCP6PAI7bkv2YRqPHOwAFIxMed4 CKTUAGnS7iNUAVLwfFN2BgitMfUVa8Z2h7oAD7H5/d5tXBcxZRYWgg27GU0X +Oeu9WTQJx8FU4xGYyxBOCqR5UFARHYrxISdKFdpTquJMjLO8ePHORry1PpE PlZSzS0buU6F0nUu2qhPamjbIr95AL3dtDwbZH9JVhh3KQ0ITm84tK6srEgg oWTjceoVAINhlHgOCkHX+yjDevDdM+RFz5oD0M218h5Q8pZTL5rP8RfJ9xj1 LfU1yuf+hihaIpZk64Kt8X4Bb+FvMSS5pOnG7OSoWx+JDlDf8XhSKseYOfYQ p4+Mx+RI6KUZLTEshXk+WuLprQvf88JrMplkFy9erJQkR0yNWO+oEMl74sEv a46p5kh/SK0HpHvAA+SYCHf3y3+z5n7y5YGCyuLVK57NI6qTuDdckx5arI+U lS6e4UAYalUB5oQjM3s9C/1asuzLWoBEIqbSM0IlnYsxec9yHctXX31Vwwk3 9TCiwXZztP+puc0XuKGlJr5khYh6DBfODMyQdPW8Em1jVBdUvT3sRH4lYbai qnagMTONs4kbTigH5g5eENXTvj+wr2AEFDkuN0nRLzgWuqe3NdsV7UwSrLvr K1euIIxkK/hdPrryeQOYHGuSIi3Zwg9YcbQEuYIEUKIEeKvkrpwCZBt2RO4U oKIdADrlnU2C+KqRjBfQKFbeD1xB8wFLIVJJLyGOHx446id59NFD/h9h/Ul3 P7y9hRUn0tvDyB4kVXGcYvNqZ/y9OuDZ6M779fAhdaTLXdmTgbCuPIRs3hdF iR8NNlzN89P3sNDB48IoAz+baHt6rXRwI6sUYgvrRQZ9CIohpWZsTd9GRUrm Br3lx1TivVjFbx3EtAFrnOhjUBzImDdORxOrP9XoKowiz5fnlNZ8cQ7eR6mt B4k7cQjXzihdFusoQnOGHc/e1cB9KLNXtKgHkXYuNC/LoeUpDBQdSKyY1mQH jOr50lYsBhUZnIpeXRUjmJ9r8S9P7KOi1qBELBFiKseg3KPYZ//2pO6HYj9A 55XwzJklLE7VVQ//jyf2YShSL01GycI+GWJ9n5K6K9pp+pAgoRpYtSiw94SN eZ8gzRtINxmwu8LE3i/ErInx3eGG75OGd0zXcCQNVwwk2rU7IUiW7IT0XkGy JrJZA0Y9dbyoNHrQXBgdRUhF1SC6pQ8pvSIk8ApxLAFrSNwhQt8tJCwp7gyT uF1IoAX6ISlqSBwTEveISlb2SANsWfjcK+gWx/sIxvfZWZYV7nRpLSb6mpKi RAIZDELyMAO+A9cM4DffLxGDU7XSfoBgcKNKaxuCB3fItBgTPa+IJaH+sobP QQWjZQxhAJGIF9OIwZMGCUtua67gteXAAjAv4fFaDSwpMIwXUPw85mF8Dw/3 ylsm3bBueRd/q3oSL18ao9Adq96m3HrrJLSWcZPCMi7VRywGGxRu99l4vPLi cFm0pwyWln/fQ7DmJVklg1qmV/tgxPkwmZFXOWjz5oMtbYXNceTzMkrdWiv9 sSHUcqPKGAqC+TtS8FrQAmjD8hD3hl4CSRR/XgKx4mxQi3kdVOoWYUtaMBjW uGq14A90lJNSJK4uPBV3tabSAjOm8NJrnA2HQ3Edy0vLXqO81CuTopPzEIUP 9KNDL7cEg7Ccet/M0sdydlwI2qejc3XRNd+E0s5UDmegeY+swDtuRoDAeMzb bOBHLAnbE4J0kV4qy5j8oTVUNkYleuIaXtsuvRG+ExnMc3p2qS4NFweagzOY xw5FqbJ5YKdMm4C2K9XG0RyZAcC8zA8f8AJHuByCoOfLo6YUn/UunZWIWNaI DRg6Tip2QVv1QnhR87K36yfeLAX9de/SlEy5w0LJiEDQZ+rexzifIROJUNVu s2cJvRkYNlBlk6ovBkCjYH/GOlxVaYSX5EFWkA3GG7udkfmwKfREMyTxqDkH uJ6iUo5ClXkqKu4pWlbAbg0u8h8k/crL+HDN+8+UKoqDGGniYkloHWLAZLR/ q85tJUaLfG6wRj+GGBzfJ6i2ZZK8qKp47OnUjXInQzh74rKEw1qcNl4rMtUB NqZvvSHHt666u/ulttrpW3X3rCGI9h3dNRTFZT9f7eRVRl/TLGHd7QZGYeHU isan+6sbDP9gscn5CvbTW4N9YhkI6npXMCxx63VqllQ7y4JXUTYTq6GSNhCt u/pN8re8RXsahye941hOPiB1vNtgZSUfjSL1pVMin6FyB/BtQKbcrinFs9aB LEP9zUvAXq/zkfxcbSAx1nuDU/L7wnToX2xp4xp0i/qRcJfsDRJfaqn3lCQs KV7bFPrplgJuHbolzK0G7NKWkfhyS6F3noSlyOqm0J9uKWANukX9lenQ3xFW fwtIfKWljlOSsKQ4dO0kaBv9mZZCbx26Jcxfp0OvOXAxbingzpOwFKk5rfJ1 o+ov11J11qj6c7jKOtJzaVPoLxlV1r3P1qH/qbHqnFH1R1/lT3L1JQej5Tmr +wotEknVeV8O8iTLcd4W1oZJqWOlWA/J4aZI5Zkc3BK33ZgjEnggi48tqpPn pihNdKL8XgaSSmvL4Q9FdGWENuhNiv7eoz83jVK6TQsFbgzbgEk0KdWCRJOi NzSQiFsqX0PnvFH1O1/1LaPqt1tdtWjASLRvG1ULu4X+glFlmfo3YXQr8Px6 U+gvthRmSvSDBuxSmIQl4IEwCWtGs6TYGyR+FUa39P6lR7cmVIv69qLvN2A1 M76lz74wCWsytqTYeRKWIvTYmZUmWBx3Dd2S3ZoSSP1vtOS4TSTw9hOXHs9c OzFcrMR5/LkwMcs+d5UbfrOlFJtuuOQbLi8v942rwsNYvMYG65Grq6s9sZhf rDuAk7/ankAG8Fst62Nd6Xm8tEgo+GgeHxGApUJmQ7JEenOO4U4cJV5qblzw cY2m5wFhDRI2b2/cJuv6RXLYjdBb/uoZv1n9jJ+oD6YAOJ+T8+uh3F5nY/m5 P32cMMvTTSbvn/7zPc5nTokVdwBzB2PzURvfZF7IKNaIp/dA2eglIfKnMmqs bI0aC3ZbkaxKTZufJdO7i/CR5xXrjcjN1AN2crAGmDO5cEudw+sVcPJWU2x/ d2tmvZZH1tuq0hkBC/GmgX4csoVxrMdVG5D7Bb+vk25XDKH7V/V7/hBjvXqD ZoZRflIj2gnNy34XYhroWv3EpuURkfGwsrWbz4alJwwLvR8eOLApxTDeUZLz ArbVrFs4y2q4gxqXgjXzbRtTEr1JyRGp7AWxx/enxdUhFHnUWZ0mktgK18XH wsunA3IKc4UzCXwejN2Xj3qyrsaTqdPaboVquDN1okDHB8qCx8aVE8Z4GDXK 8wc3//BJI/V0qmh3+PBhhvmhNYtVxhFOJOFIcVqgn3XlB1StdMKJ11OUax6F ndKCJ8Q6t8gVrzyyEPDrnjBZa9jfWyQbFQ8X8AqeOqjEDKy1Bkvu1IDdWWGq HzN48NoZ4PHgOGeQir9zvgSd4jSiA5zuKt+Fv2uL+JPX8BPWfkW+mLzSZBFN KYbVkYkBo7ui77TU5Be7hW7p83OP/t2W1K3zWJf2HAk8T2c5jHUWTy/qDUIN f+aRcLjPNiC1ooT9E3xC4adtkKwDgxUrTNnwe1BgjX4liv5AnP6Rv/g+FL1N IhfZ/8Qj/cCosjQg9B8aVT9uWVU8/ZTLWV97qURWKegEXVxczI4dO1YoR6MI rp3R3S+HE6jH2h5iJhYdqksFVqpvbmqL1551136r2oPz6PL/Tf6u/S60n/FR ZgHxYgkb8CuD71QhufQc6kajUVPp2c4KqBZtroAWQCVNYgOWhth01V11V70T 1YO6ke6GqbuagUsKJBkAljwUAdieHryDVvDnrnG00x0ygXwEgPe9+N8aXXCQ obyfqVLfVXfVXfUeqK4M+D6CYw2OJETI6OexH/UNWEIwjg+u7GAdbGdhM3A5 46/QexPy18J7+0bjHrksg1WVHzu0uEovGRSGHayDbT0s0TD2RXZocmN+s+FI vZpvlNJmlwITagfrYDsFq7glgslzh3tO3A52vcPmvH+Ox2MJr/6Vbt5/83fW cdOurqvr6vZeXQpXsrblrpcUQN7cql4AmsH3oqqY0LsWhYirWoIvXnrP+LWd /IZM/SrPruzKrnzrlnBWbJtI43t2Jeq465s9YDJx8SuvcnWhd8l3ZVd25Vu6 hP9DsQOs8uzHAZbgS1IoECKTV5V3ZVd25XVVRtFNW00So8oMXEnccddLCuBu ztR/+5D/bMEwgnfgDtyBr0dwUgYjVKKJAyzBVxiwvr6+BF9hwMbGxhJ8hQA9 fPgqk+fI+GmwruzKrrwuSngUdhtZ5PHMAZbgSwDuc4gypZn/AVGbBEA=\ \>"], "Graphics", CellMargins->{{92, Inherited}, {Inherited, Inherited}}, Evaluatable->False, TextAlignment->Left, ImageSize->{281, 299}, ImageMargins->{{0, 0}, {0, 0}}, ImageRegion->{{0, 1}, {0, 1}}], Cell[TextData[{ "The position of a pendulum is given by the angle \[Theta]. Then,\n\t", Cell[BoxData[ FormBox[ RowBox[{ SubscriptBox["\[PartialD]", "t"], " ", "\[Theta]"}], TraditionalForm]]], " = \[Omega](", StyleBox["t", FontSlant->"Italic"], ")\n\t", Cell[BoxData[ FormBox[ RowBox[{ SubscriptBox["\[PartialD]", "t"], " ", "\[Omega]"}], TraditionalForm]]], " = - (", StyleBox["g/L)", FontSlant->"Italic"], " Sin[\[Theta]]\nwhen the angle is small |\[Theta]| << 1, then Sin[\[Theta]] \ ~ \[Theta], and\n\t", Cell[BoxData[ FormBox[ RowBox[{ SubscriptBox["\[PartialD]", "t"], " ", RowBox[{ SubscriptBox["\[PartialD]", "t"], " ", "\[Theta]"}]}], TraditionalForm]]], " = - (", StyleBox["g/L", FontSlant->"Italic"], ") \[Theta]\n\[Theta](", StyleBox["t", FontSlant->"Italic"], ") = ", Cell[BoxData[ FormBox[ SubscriptBox["\[Theta]", "0"], TraditionalForm]]], "Sin[", Cell[BoxData[ FormBox[ SubscriptBox["\[Omega]", "0"], TraditionalForm]]], StyleBox["t", FontSlant->"Italic"], " + ", Cell[BoxData[ FormBox[ SubscriptBox["\[Phi]", "0"], TraditionalForm]]], "] where ", Cell[BoxData[ FormBox[ SuperscriptBox[ SubscriptBox["\[Omega]", "0"], "2"], TraditionalForm]]], " = ", StyleBox["g/L", FontSlant->"Italic"], ", and ", Cell[BoxData[ FormBox[ SubscriptBox["\[Theta]", "0"], TraditionalForm]]], "and ", Cell[BoxData[ FormBox[ SubscriptBox["\[Phi]", "0"], TraditionalForm]]], " are initial conditions." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["Integration with a Computer", "Section"], Cell[TextData[{ "The numerical solution for the pendulum's motion requires the integration \ of two coupled first order ordinary differential equations (ODEs). This is \ the best problem to begin the exploration of computation physics. \nWe will \ first explore Euler's method: integration of the equations of motion by \ approximating the first-derivative with a ", StyleBox["forward difference", FontSlant->"Italic"], ". We'll see that Euler's method fails and becomes unstable. \nNext, we'll \ explore the Euler-Cromer method which significantly improves the stability of \ the integration. It is similar to the ", StyleBox["Leap-Frog", FontSlant->"Italic"], " method which is based on ", StyleBox["central-difference", FontSlant->"Italic"], " approximation to the first-derivative.", "\nFinally, we'll look into the numerical tools built-in to ", StyleBox["Mathematica", FontSlant->"Italic"], "." }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"dt", " ", "=", " ", "0.01"}], ";", " ", RowBox[{"MachineNumberQ", "[", "dt", "]"}]}]], "Input", CellLabel->"In[25]:="], Cell[BoxData[ RowBox[{ RowBox[{"Off", "[", RowBox[{"General", "::", "\"\\""}], "]"}], ";", " ", RowBox[{"(*", " ", RowBox[{"turn", " ", "off", " ", "spell", " ", "checking"}], " ", "*)"}]}]], "Input", CellLabel->"In[26]:="], Cell[CellGroupData[{ Cell["Euler's Method", "Subsection"], Cell[TextData[{ "Using two ", StyleBox["forward-difference", FontSlant->"Italic"], " approximations\[Ellipsis]" }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"stepE", "[", RowBox[{"{", RowBox[{"\[Theta]_", ",", " ", "\[Omega]_"}], "}"}], "]"}], " ", ":=", " ", RowBox[{"N", "[", RowBox[{"{", RowBox[{ RowBox[{"\[Theta]", " ", "+", " ", RowBox[{"\[Omega]", " ", "dt"}]}], ",", " ", RowBox[{"\[Omega]", " ", "-", " ", RowBox[{ RowBox[{"Sin", "[", "\[Theta]", "]"}], " ", "dt"}]}]}], "}"}], "]"}]}]], "Input", CellLabel->"In[27]:="], Cell[BoxData[ RowBox[{ RowBox[{"stepE", "[", RowBox[{"{", RowBox[{"\[Theta]0", ",", " ", "\[Omega]0"}], "}"}], "]"}], " "}]], "Input", CellLabel->"In[28]:="], Cell[BoxData[ RowBox[{"stepE", "[", RowBox[{"{", RowBox[{ RowBox[{"N", "[", RowBox[{"75", " ", "Degree"}], "]"}], ",", "0.0"}], "}"}], "]"}]], "Input", CellLabel->"In[29]:="], Cell[BoxData[ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"steps", " ", "=", " ", RowBox[{"NestList", "[", RowBox[{"stepE", ",", " ", RowBox[{"{", RowBox[{ RowBox[{"N", "[", RowBox[{"75", " ", "Degree"}], "]"}], ",", "0.0"}], "}"}], ",", " ", "5000"}], "]"}]}], ";"}], ")"}], " ", "//", " ", "Timing"}]], "Input",\ CellLabel->"In[30]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", "steps", "]"}]], "Input", CellLabel->"In[31]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"First", "[", RowBox[{"Transpose", "[", "steps", "]"}], "]"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\<\[Theta] vs. time step\>\""}]}], "]"}]], "Input", CellLabel->"In[32]:="], Cell[TextData[{ "If ", StyleBox["dt", FontSlant->"Italic"], " gets too large, then the entire method fails." }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"dt", " ", "=", " ", "0.1"}], ";", RowBox[{"steps", " ", "=", " ", RowBox[{"NestList", "[", RowBox[{"stepE", ",", " ", RowBox[{"{", RowBox[{ RowBox[{"N", "[", RowBox[{"75", " ", "Degree"}], "]"}], ",", "0.0"}], "}"}], ",", " ", "5000"}], "]"}]}], ";"}]], "Input", CellLabel->"In[33]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"First", "[", RowBox[{"Transpose", "[", "steps", "]"}], "]"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\<\[Theta] vs. time step\>\""}]}], "]"}]], "Input", CellLabel->"In[34]:="] }, Closed]], Cell[CellGroupData[{ Cell["Central Difference (Leap-Frog Method)", "Subsection"], Cell[TextData[{ "The ", StyleBox["leap-frog", FontSlant->"Italic"], " method approximates the differential equations with ", StyleBox["central-difference", FontSlant->"Italic"], " approximations.\n\t", Cell[BoxData[ FormBox[ RowBox[{ SubscriptBox["\[PartialD]", "t"], " ", "\[Theta]"}], TraditionalForm]]], " ~ ", Cell[BoxData[ FormBox[ StyleBox[ FractionBox[ RowBox[{ RowBox[{"\[Theta]", "[", RowBox[{"t", "+", RowBox[{"\[Delta]t", "/", "2"}]}], "]"}], " ", "-", " ", RowBox[{"\[Theta]", "[", RowBox[{"t", "-", RowBox[{"\[Delta]t", "/", "2"}]}], "]"}]}], "\[Delta]t"], FontSize->16], TraditionalForm]]], " = \[Omega](", StyleBox["t", FontSlant->"Italic"], ") + O(", Cell[BoxData[ FormBox[ SuperscriptBox["\[Delta]t", "2"], TraditionalForm]]], ")\n\t", Cell[BoxData[ FormBox[ RowBox[{ SubscriptBox["\[PartialD]", "t"], " ", "\[Omega]"}], TraditionalForm]]], " ~ ", Cell[BoxData[ FormBox[ StyleBox[ FractionBox[ RowBox[{ RowBox[{"\[Omega]", "[", RowBox[{"t", "+", RowBox[{"\[Delta]t", "/", "2"}]}], "]"}], " ", "-", " ", RowBox[{"\[Omega]", "[", RowBox[{"t", "-", RowBox[{"\[Delta]t", "/", "2"}]}], "]"}]}], "\[Delta]t"], FontSize->16], TraditionalForm]]], "= - Sin[\[Theta](", StyleBox["t", FontSlant->"Italic"], ")]+ O(", Cell[BoxData[ FormBox[ SuperscriptBox["\[Delta]t", "2"], TraditionalForm]]], ")\nLooking carefully at the difference equations, evaluation of \[Theta] \ steps \"leap\" over \[Omega], and evaluations of \[Omega] steps \"leap\" over \ \[Theta]. This presents two problems: first, the method is not ", StyleBox["self-starting", FontSlant->"Italic"], ", and, second, the values of \[Theta] and \[Omega] are not known at the \ same time step. This method has the advantage of being ", StyleBox["second-order", FontSlant->"Italic"], " with reduced truncation error.\nIn the following, we'll use an Euler \ \"half-step\" to start the integration by computing \[Theta](\[Delta]t/2). \ Then, we'll \"leap\" with \[Omega] followed by \[Theta] using \"full \ steps\"." }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"stepLF", "[", RowBox[{"{", RowBox[{"\[Theta]_", ",", " ", "\[Omega]_"}], "}"}], "]"}], " ", ":=", " ", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", "\[Omega]next", "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Omega]next", " ", "=", " ", RowBox[{"N", "[", RowBox[{"\[Omega]", " ", "-", " ", RowBox[{ RowBox[{"Sin", "[", "\[Theta]", "]"}], " ", "dt"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"{", RowBox[{ RowBox[{"N", "[", RowBox[{"\[Theta]", " ", "+", " ", RowBox[{"dt", " ", "\[Omega]next"}]}], "]"}], ",", " ", "\[Omega]next"}], "}"}]}]}], "]"}]}]], "Input", CellLabel->"In[35]:="], Cell[BoxData[{ RowBox[{ RowBox[{"dt", " ", "=", " ", "0.01"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"w0", " ", "=", " ", "0.0"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Theta]0", " ", "=", " ", RowBox[{"N", "[", RowBox[{"75", " ", "Degree"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Theta]1", " ", "=", " ", RowBox[{"N", "[", RowBox[{"\[Theta]0", " ", "+", " ", RowBox[{"w0", " ", RowBox[{"(", RowBox[{"dt", "/", "2"}], ")"}]}]}], "]"}]}], ";", " ", RowBox[{"(*", " ", RowBox[{"Euler", " ", "Step", " ", "to", " ", "start"}], " ", "*)"}], RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"steps", " ", "=", " ", RowBox[{"NestList", "[", RowBox[{"stepLF", ",", RowBox[{"{", RowBox[{"\[Theta]1", ",", " ", "w0"}], "}"}], ",", " ", "5000"}], "]"}]}], ";"}], ")"}], " ", "//", " ", "Timing"}]}]}], "Input", CellLabel->"In[36]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"First", "[", RowBox[{"Transpose", "[", "steps", "]"}], "]"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\<\[Theta] vs. time step\>\""}]}], "]"}]], "Input", CellLabel->"In[40]:="] }, Closed]], Cell[CellGroupData[{ Cell["Euler-Cromer Method", "Subsection"], Cell[TextData[{ "This method is practically the same as the ", StyleBox["leap-frog", FontSlant->"Italic"], " method. The only difference is to drop the pretense of leaping and shift \ the two variables, {\[Theta]", ", \[Omega]}, so that they \"align\" (instead of \"stagger\") every time \ step. In Euler-Cromer, we actually make \"two steps\": first a step in the \ frequency, \[Omega], and then in the angle, \[Theta]; however, we do not need \ to \"start\" the integration with a \"half-step\"." }], "Text"], Cell[BoxData[ RowBox[{ RowBox[{"stepEC", "[", RowBox[{"{", RowBox[{"\[Theta]_", ",", " ", "\[Omega]_"}], "}"}], "]"}], " ", ":=", " ", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", "\[Omega]next", "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Omega]next", " ", "=", " ", RowBox[{"N", "[", RowBox[{"\[Omega]", " ", "-", " ", RowBox[{ RowBox[{"Sin", "[", "\[Theta]", "]"}], " ", "dt"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"{", RowBox[{ RowBox[{"N", "[", RowBox[{"\[Theta]", " ", "+", " ", RowBox[{"dt", " ", "\[Omega]next"}]}], "]"}], ",", " ", "\[Omega]next"}], "}"}]}]}], "]"}]}]], "Input", CellLabel->"In[41]:="], Cell[BoxData[ RowBox[{ RowBox[{"dt", " ", "=", " ", "0.01"}], ";", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"steps", " ", "=", " ", RowBox[{"NestList", "[", RowBox[{"stepEC", ",", " ", RowBox[{"{", RowBox[{ RowBox[{"N", "[", RowBox[{"75", " ", "Degree"}], "]"}], ",", "0.0"}], "}"}], ",", " ", "5000"}], "]"}]}], ";"}], ")"}], " ", "//", " ", "Timing"}]}]], "Input", CellLabel->"In[42]:="], Cell[TextData[{ StyleBox["(Mathematica", FontSlant->"Italic"], " requires more time when using the Euler-Cromer method, but this can be \ fixed by pre-compiling the module. We'll discuss this later in the course.)" }], "Text"], Cell[BoxData[ RowBox[{"ListPlot", "[", "steps", "]"}]], "Input", CellLabel->"In[43]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"First", "[", RowBox[{"Transpose", "[", "steps", "]"}], "]"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\<\[Theta] vs. time step\>\""}]}], "]"}]], "Input", CellLabel->"In[44]:="], Cell["\<\ The Euler-Cramer method is also more stable for large time steps.\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"dt", " ", "=", " ", "0.1"}], ";", RowBox[{ RowBox[{ RowBox[{"steps", " ", "=", " ", RowBox[{"NestList", "[", RowBox[{"stepEC", ",", " ", RowBox[{"{", RowBox[{ RowBox[{"N", "[", RowBox[{"75", " ", "Degree"}], "]"}], ",", "0.0"}], "}"}], ",", " ", "5000"}], "]"}]}], ";"}], " ", "//", " ", "Timing"}]}]], "Input", CellLabel->"In[45]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"First", "[", RowBox[{"Transpose", "[", "steps", "]"}], "]"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\<\[Theta] vs. time step\>\""}]}], "]"}]], "Input", CellLabel->"In[46]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"First", "[", RowBox[{"Transpose", "[", RowBox[{"steps", "\[LeftDoubleBracket]", RowBox[{"Range", "[", RowBox[{"1", ",", "500"}], "]"}], "\[RightDoubleBracket]"}], "]"}], "]"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\<\[Theta] vs. time step\>\""}]}], "]"}]], "Input", CellLabel->"In[47]:="] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell["Using NDSolve", "Section"], Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], " can numerically integrate differential equations using ", StyleBox["NDSolve[...]", FontWeight->"Bold"], ". ", StyleBox["NDSolve", FontWeight->"Bold"], " is probably the most sophisticated and comprehensive command within ", StyleBox["Mathematica", FontSlant->"Italic"], ".", " Since there are many methods to finding numerical solutions to \ differential equations, ", StyleBox["NDSolve", FontWeight->"Bold"], " has many options." }], "Text"], Cell[TextData[{ "In order to find the solution to the equations of motion for the pendulum, \ we need to define:\n\t1. The unknown functions of time, ", "\[Theta]", "[", StyleBox["t", FontSlant->"Italic"], "] and \[Omega]", "[", StyleBox["t", FontSlant->"Italic"], "].", "\n\t2. The equations of motions, and\n\t3. The equations for the initial \ conditions." }], "Text"], Cell[BoxData[ RowBox[{"eqs", " ", "=", " ", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"\[Theta]", "'"}], "[", "t", "]"}], " ", "==", " ", RowBox[{"w", "[", "t", "]"}]}], ",", " ", RowBox[{ RowBox[{ RowBox[{"w", "'"}], "[", "t", "]"}], " ", "==", " ", RowBox[{"-", " ", RowBox[{"Sin", "[", RowBox[{"\[Theta]", "[", "t", "]"}], "]"}]}]}]}], "}"}]}]], "Input", CellLabel->"In[48]:="], Cell[BoxData[ RowBox[{"initial", " ", "=", " ", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"\[Theta]", "[", "0", "]"}], " ", "\[Equal]", RowBox[{"N", "[", RowBox[{"75", " ", "Degree"}], "]"}]}], ",", " ", RowBox[{ RowBox[{"w", "[", "0", "]"}], " ", "==", " ", "0"}]}], "}"}]}]], "Input",\ CellLabel->"In[49]:="], Cell[BoxData[ RowBox[{"?", "NDSolve"}]], "Input", CellLabel->"In[50]:="], Cell[BoxData[ RowBox[{"Options", "[", "NDSolve", "]"}]], "Input", CellLabel->"In[51]:="], Cell[BoxData[ RowBox[{ RowBox[{"(", RowBox[{"sol", " ", "=", " ", RowBox[{"NDSolve", "[", " ", RowBox[{ RowBox[{"eqs", " ", "~", " ", "Join", " ", "~", " ", "initial"}], ",", " ", RowBox[{"{", RowBox[{ RowBox[{"\[Theta]", "[", "t", "]"}], ",", " ", RowBox[{"w", "[", "t", "]"}]}], "}"}], ",", " ", RowBox[{"{", RowBox[{"t", ",", " ", "0", ",", " ", "60"}], "}"}]}], "]"}]}], " ", ")"}], "//", " ", "Timing"}]], "Input", CellLabel->"In[52]:="], Cell[BoxData[ RowBox[{"Plot", "[", RowBox[{ RowBox[{"Evaluate", "[", RowBox[{ RowBox[{"\[Theta]", "[", "t", "]"}], "/.", "\[InvisibleSpace]", "sol"}], "]"}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "60"}], "}"}]}], "]"}]], "Input", CellLabel->"In[53]:="], Cell["To appreciate the numerical solution, examine", "Text"], Cell[BoxData[ RowBox[{ RowBox[{ RowBox[{ RowBox[{"\[Theta]", "[", "t", "]"}], " ", "/.", " ", "sol"}], "//", " ", "FullForm"}], " ", "//", " ", "Shallow"}]], "Input", CellLabel->"In[54]:="], Cell["\<\ The interpolating function contains all of the numerical solution. For \ example, the time points are found by\ \>", "Text"], Cell[BoxData[ RowBox[{ RowBox[{"Part", "[", " ", RowBox[{ RowBox[{ RowBox[{"\[Theta]", "[", "t", "]"}], " ", "/.", " ", "sol"}], " ", ",", " ", "1", ",", " ", "0", ",", " ", "3", ",", " ", "1"}], "]"}], " ", "//", " ", "Shallow"}]], "Input", CellLabel->"In[55]:="], Cell[BoxData[ RowBox[{"stepsND", " ", "=", " ", RowBox[{"(", RowBox[{ RowBox[{"Part", "[", " ", RowBox[{ RowBox[{ RowBox[{"\[Theta]", "[", "t", "]"}], " ", "/.", " ", "sol"}], " ", ",", " ", "1", ",", " ", "0", ",", " ", "3", ",", " ", "1"}], "]"}], " ", "//", " ", "Length"}], ")"}]}]], "Input", CellLabel->"In[56]:="], Cell[BoxData[ RowBox[{"ListPlot", "[", RowBox[{ RowBox[{"Delete", "[", RowBox[{ RowBox[{ RowBox[{"RotateLeft", "[", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"\[Theta]", "[", "t", "]"}], "/.", "\[InvisibleSpace]", "sol"}], ")"}], "\[LeftDoubleBracket]", RowBox[{"1", ",", "0", ",", "3", ",", "1"}], "\[RightDoubleBracket]"}], "]"}], "-", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"\[Theta]", "[", "t", "]"}], "/.", "\[InvisibleSpace]", "sol"}], ")"}], "\[LeftDoubleBracket]", RowBox[{"1", ",", "0", ",", "3", ",", "1"}], "\[RightDoubleBracket]"}]}], ",", RowBox[{"-", "1"}]}], "]"}], ",", RowBox[{"PlotRange", "\[Rule]", "All"}], ",", RowBox[{"PlotLabel", "\[Rule]", "\"\