(* Content-type: application/vnd.wolfram.mathematica *) (*** Wolfram Notebook File ***) (* http://www.wolfram.com/nb *) (* CreatedBy='Mathematica 9.0' *) (*CacheID: 234*) (* Internal cache information: NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 157, 7] NotebookDataLength[ 25794, 583] NotebookOptionsPosition[ 25227, 561] NotebookOutlinePosition[ 25573, 576] CellTagsIndexPosition[ 25530, 573] WindowFrame->Normal*) (* Beginning of Notebook Content *) Notebook[{ Cell[BoxData[ RowBox[{ RowBox[{"RK23", "[", RowBox[{ "f_", ",", "u0_", ",", "t0_", ",", "T_", ",", " ", "h0_", ",", " ", "tol_"}], "]"}], ":=", RowBox[{"(", "\[IndentingNewLine]", RowBox[{ RowBox[{"h", "=", "h0"}], ";", "\[IndentingNewLine]", RowBox[{"t", "=", "t0"}], ";", "\[IndentingNewLine]", RowBox[{"y", "=", RowBox[{"{", "0", "}"}]}], ";", "\[IndentingNewLine]", RowBox[{ RowBox[{"y", "[", RowBox[{"[", "1", "]"}], "]"}], "=", "u0"}], ";", "\[IndentingNewLine]", RowBox[{"i", "=", "1"}], ";", "\[IndentingNewLine]", RowBox[{"t", "=", RowBox[{"{", "0", "}"}]}], ";", "\[IndentingNewLine]", RowBox[{"While", "[", RowBox[{ RowBox[{ RowBox[{"t", "[", RowBox[{"[", "i", "]"}], "]"}], "<", "T"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"k1", "=", RowBox[{"h", " ", RowBox[{"f", "[", RowBox[{ RowBox[{"t", "[", RowBox[{"[", "i", "]"}], "]"}], ",", RowBox[{"y", "[", RowBox[{"[", "i", "]"}], "]"}]}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"k2", "=", RowBox[{"h", " ", RowBox[{"f", "[", RowBox[{ RowBox[{ RowBox[{"t", "[", RowBox[{"[", "i", "]"}], "]"}], "+", "h"}], ",", RowBox[{ RowBox[{"y", "[", RowBox[{"[", "i", "]"}], "]"}], "+", "k1"}]}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"k3", "=", RowBox[{"h", " ", RowBox[{"f", "[", RowBox[{ RowBox[{ RowBox[{"t", "[", RowBox[{"[", "i", "]"}], "]"}], "+", FractionBox["h", "2"]}], ",", RowBox[{ RowBox[{"y", "[", RowBox[{"[", "i", "]"}], "]"}], "+", RowBox[{ FractionBox["1", "4"], "k1"}], "+", RowBox[{ FractionBox["1", "4"], "k2"}]}]}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"While", "[", RowBox[{ RowBox[{ RowBox[{"Abs", "[", RowBox[{ FractionBox["1", "3"], RowBox[{"(", RowBox[{"k1", "+", "k2", "-", RowBox[{"2", "k3"}]}], ")"}]}], "]"}], ">", "tol"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"h", "=", RowBox[{"h", SuperscriptBox[ RowBox[{"(", FractionBox["tol", RowBox[{"Abs", "[", RowBox[{ FractionBox["1", "3"], RowBox[{"(", RowBox[{"k1", "+", "k2", "-", RowBox[{"2", "k3"}]}], ")"}]}], "]"}]], ")"}], RowBox[{"1", "/", "3"}]]}]}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{"h", "<", SuperscriptBox["10", RowBox[{"-", "5"}]]}], ",", RowBox[{ RowBox[{ "Print", "[", "\"\\"", "]"}], ";", " ", RowBox[{"Return", "[", "]"}]}]}], "]"}], ";", "\[IndentingNewLine]", RowBox[{"k1", "=", RowBox[{"h", " ", RowBox[{"f", "[", RowBox[{ RowBox[{"t", "[", RowBox[{"[", "i", "]"}], "]"}], ",", RowBox[{"y", "[", RowBox[{"[", "i", "]"}], "]"}]}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"k2", "=", RowBox[{"h", " ", RowBox[{"f", "[", RowBox[{ RowBox[{ RowBox[{"t", "[", RowBox[{"[", "i", "]"}], "]"}], "+", "h"}], ",", RowBox[{ RowBox[{"y", "[", RowBox[{"[", "i", "]"}], "]"}], "+", "k1"}]}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"k3", "=", RowBox[{"h", " ", RowBox[{"f", "[", RowBox[{ RowBox[{ RowBox[{"t", "[", RowBox[{"[", "i", "]"}], "]"}], "+", FractionBox["h", "2"]}], ",", RowBox[{ RowBox[{"y", "[", RowBox[{"[", "i", "]"}], "]"}], "+", RowBox[{ FractionBox["1", "4"], "k1"}], "+", RowBox[{ FractionBox["1", "4"], "k2"}]}]}], "]"}]}]}], ";"}]}], "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", RowBox[{"y", "=", RowBox[{"Append", "[", RowBox[{"y", ",", RowBox[{ RowBox[{"y", "[", RowBox[{"[", "i", "]"}], "]"}], "+", RowBox[{ FractionBox["1", "6"], "k1"}], "+", RowBox[{ FractionBox["1", "6"], "k2"}], "+", RowBox[{ FractionBox["2", "3"], "k3"}]}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"t", "=", RowBox[{"Append", "[", RowBox[{"t", ",", RowBox[{ RowBox[{"t", "[", RowBox[{"[", "i", "]"}], "]"}], "+", "h"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"h", "=", RowBox[{"h", SuperscriptBox[ RowBox[{"(", FractionBox["tol", RowBox[{ FractionBox["1", "3"], RowBox[{"Abs", "[", RowBox[{"(", RowBox[{"k1", "+", "k2", "-", RowBox[{"2", "k3"}]}], ")"}], "]"}]}]], ")"}], RowBox[{"1", "/", "3"}]]}]}], ";", "\[IndentingNewLine]", RowBox[{"i", "++"}]}]}], "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", RowBox[{"{", RowBox[{"t", ",", "y"}], "}"}]}], "\[IndentingNewLine]", ")"}]}]], "Input", CellChangeTimes->{{3.664871610509589*^9, 3.664871614656519*^9}, { 3.664871696727519*^9, 3.664871831755891*^9}, {3.6648718619652395`*^9, 3.664871954563657*^9}, {3.6648720176121984`*^9, 3.6648720224646025`*^9}, { 3.6648720694027576`*^9, 3.6648720698631115`*^9}, {3.6648721058325005`*^9, 3.664872134340661*^9}, {3.6648721934073696`*^9, 3.664872217575444*^9}, { 3.6648722598593082`*^9, 3.664872259970392*^9}, {3.6648722916597743`*^9, 3.664872305987896*^9}, {3.6648859919958224`*^9, 3.664886012330188*^9}, { 3.6685006275952377`*^9, 3.6685006998925457`*^9}, {3.6685007378007164`*^9, 3.668500927249287*^9}, {3.6685009586664133`*^9, 3.6685009625164323`*^9}, { 3.6685010818597136`*^9, 3.6685013564686213`*^9}, {3.668501545838112*^9, 3.6685015539951315`*^9}, {3.66850278262576*^9, 3.6685028976983414`*^9}, 3.668502998349098*^9, {3.6685030341751475`*^9, 3.6685030387584095`*^9}, { 3.6685031475376315`*^9, 3.668503149301732*^9}, {3.6685031834636865`*^9, 3.6685032900097804`*^9}, {3.6685033278939476`*^9, 3.6685034179800997`*^9}, {3.6685039380828485`*^9, 3.6685039437231708`*^9}, {3.6685048254836044`*^9, 3.668504833079039*^9}, 3.668608102483925*^9, {3.668608278807457*^9, 3.6686083806469235`*^9}, { 3.6686085820739594`*^9, 3.6686085849690523`*^9}, {3.6691016955336423`*^9, 3.669101721614134*^9}, {3.669101797545477*^9, 3.6691018039028406`*^9}, 3.669101871676717*^9, {3.669101919200435*^9, 3.669101919777468*^9}, { 3.6691027582564263`*^9, 3.669102787679109*^9}, {3.6691028185438747`*^9, 3.669102875848152*^9}, {3.6691047896316147`*^9, 3.6691048618467445`*^9}, 3.669104941499301*^9, 3.669105166041144*^9, 3.6691073882942495`*^9, 3.6691074399552045`*^9, {3.6691075138724318`*^9, 3.6691075359486947`*^9}, { 3.6691076052846603`*^9, 3.6691076431918287`*^9}, 3.669107808274271*^9, { 3.6691079528905425`*^9, 3.6691079622910805`*^9}, 3.669107992394802*^9}], Cell[CellGroupData[{ Cell[BoxData[{ RowBox[{ RowBox[{"\[Lambda]", "=", "1000000"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"f", "[", RowBox[{"t_", ",", "y_"}], "]"}], ":=", RowBox[{"\[Lambda]", RowBox[{"(", RowBox[{ RowBox[{"-", "y"}], "+", RowBox[{"Sin", "[", "t", "]"}]}], ")"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{"t", "=", RowBox[{ RowBox[{"RK23", "[", RowBox[{"f", ",", "0", ",", "0", ",", "10", ",", "0.1", ",", "0.1"}], "]"}], "[", RowBox[{"[", "1", "]"}], "]"}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{"y", "=", RowBox[{ RowBox[{"RK23", "[", RowBox[{"f", ",", "0", ",", "0", ",", "10", ",", "0.1", ",", "0.1"}], "]"}], "[", RowBox[{"[", "2", "]"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"plotAppr", "=", RowBox[{"ListPlot", "[", RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"t", "[", RowBox[{"[", RowBox[{"i", "+", "1"}], "]"}], "]"}], ",", " ", RowBox[{"y", "[", RowBox[{"[", RowBox[{"i", "+", "1"}], "]"}], "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"i", ",", "0", ",", RowBox[{ RowBox[{"Length", "[", "y", "]"}], "-", "1"}]}], "}"}]}], "]"}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"c", "=", FractionBox["\[Lambda]", RowBox[{"1", "+", SuperscriptBox["\[Lambda]", "2"]}]]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"plotExact", "=", RowBox[{"Plot", "[", RowBox[{ RowBox[{ RowBox[{"c", " ", SuperscriptBox["E", RowBox[{ RowBox[{"-", "\[Lambda]"}], " ", "t"}]]}], "+", RowBox[{ FractionBox[ SuperscriptBox["\[Lambda]", "2"], RowBox[{"1", "+", SuperscriptBox["\[Lambda]", "2"]}]], RowBox[{"Sin", "[", "t", "]"}]}], "-", RowBox[{ FractionBox["\[Lambda]", RowBox[{"1", "+", SuperscriptBox["\[Lambda]", "2"]}]], RowBox[{"Cos", "[", "t", "]"}]}]}], ",", RowBox[{"{", RowBox[{"t", ",", "0", ",", "10"}], "}"}], ",", RowBox[{"PlotStyle", "\[Rule]", "Red"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{"Show", "[", RowBox[{"plotAppr", ",", "plotExact"}], "]"}], "\[IndentingNewLine]"}], "Input", CellChangeTimes->{{3.669107455475092*^9, 3.66910749687746*^9}}], Cell[BoxData["\<\"The step size is too small. Stiff equation suspected\"\>"], \ "Print", CellChangeTimes->{3.669107996827056*^9}], Cell[BoxData[ RowBox[{ StyleBox[ RowBox[{"Part", "::", "partd"}], "MessageName"], RowBox[{ ":", " "}], "\<\"Part specification \[NoBreak]\\!\\(Null \ \[LeftDoubleBracket] 1 \[RightDoubleBracket]\\)\[NoBreak] is longer than \ depth of object. \\!\\(\\*ButtonBox[\\\"\[RightSkeleton]\\\", ButtonStyle->\\\ \"Link\\\", ButtonFrame->None, \ ButtonData:>\\\"paclet:ref/message/General/partd\\\", ButtonNote -> \ \\\"Part::partd\\\"]\\)\"\>"}]], "Message", "MSG", CellChangeTimes->{3.6691079968800583`*^9}], Cell[BoxData[ RowBox[{"Null", "\[LeftDoubleBracket]", "1", "\[RightDoubleBracket]"}]], "Output", CellChangeTimes->{ 3.6691078845636344`*^9, {3.669107939816795*^9, 3.669107996882059*^9}}], Cell[BoxData["\<\"The step size is too small. Stiff equation suspected\"\>"], \ "Print", CellChangeTimes->{3.6691079968830585`*^9}], Cell[BoxData[ RowBox[{ StyleBox[ RowBox[{"Part", "::", "partd"}], "MessageName"], RowBox[{ ":", " "}], "\<\"Part specification \[NoBreak]\\!\\(Null \ \[LeftDoubleBracket] 2 \[RightDoubleBracket]\\)\[NoBreak] is longer than \ depth of object. \\!\\(\\*ButtonBox[\\\"\[RightSkeleton]\\\", ButtonStyle->\\\ \"Link\\\", ButtonFrame->None, \ ButtonData:>\\\"paclet:ref/message/General/partd\\\", ButtonNote -> \ \\\"Part::partd\\\"]\\)\"\>"}]], "Message", "MSG", CellChangeTimes->{3.6691079969160604`*^9}], Cell[BoxData[ GraphicsBox[{{{}, {{}, {RGBColor[0.368417, 0.506779, 0.709798], PointSize[ 0.012833333333333334`], AbsoluteThickness[1.6], PointBox[{{0.00008434326656410108, 2.}}]}, {}}, {}}, {{}, {}, {RGBColor[1, 0, 0], AbsoluteThickness[1.6], Opacity[1.], LineBox[CompressedData[" 1:eJwVmHk8VO8Xxy2lyDqDss2MpY3iS1SWOicqpY2kRMkSKbJFqVRkKUuhbKmE lKVIklC5ktBCpUXWZN9mrp2x/u7vn7mv9+uZ+zznnM/nOffeR9Hefb8jHw8P jw718//r/cv2AprHwzfzqOQ29/Ue2JxkNaB3iHUSaJnplfqMk1C/J8yLxToP 8KgyV5B5Hs4ZHxZdyAoF36sXwi2YobC9HI9MMO/Auq154mLUNfuqduAgMxMs g6I26Slkwj+7HGk2swh+DH04ECZVBBfNXGz7mR/hsZevVtZ4FQhy7350VW8A Nkvrr4BPPXiWXYkY8GoH7b1xjvPb2mBSIia7MbwdHi95/NfTqg0u22XUfEpr B7NsOZ16tzaI4Pkunvm7HTR25GFIfBtkbFaKddTrAIFsmwDLnjZoLS6/38Lb CQc2jo3VXmsH0xeCz75FdYGiSWbY57cdoJ52+3t+Ti+cDfvy3pC/GyQT/pso /dgL7sM5PxZKd8N0eLV8dUcvLEw95Fi4qhuqzixy7pTtg+26d/707+kGe6OL c1LX+kB6pCw8LL4bYtvs1c4d6wf36VUv7q7sgRmWVrCuOBtkXp+ZeqHXC22S X7O2qbGBwLJde3f1wsfFrt/MtrNhYcds2G/rXogbfCR3yo8N01JL0nL9euG/ Upm8xB42vPzCnQx92wsOx3j/TpdyINhzUr5Svw869zrL/904CPV3x/wUVPtB ssJV5fK2QRALv/3EdUM/GG7yXKOwfxC217zty9raD8lqFwysXQYhOvHz+VGb fjiyOOJI3f1B+PQka7z4Vj/Uvnt27yvfEBydMT8tOtEPJdoT8qVfhiB4TvbW nbwBYD+ZVjlWPwQdXVuyXr8dAHllnrVznUMgOCzd+rlqAM5LCG4ymKfuN9wT WtAyANoc2aOFmsNQe7uHni/Ihqz0zfdz44YhsO/2ZpMjbIiXDVFIsR0Bi4Wn /hWMsGG8vSZ9r9sIRLxL59szy4YD2Uu1Zi6OQJCuhW61AAdomLn9UPwIfI7c +TlchgMRjl/cRWpGwNlK6T3PZg7459LKzuuPQmbXJ6PyqxxoOW+9e+XOUfgV IJw0GcaBTUZpv38eHAXTIJ3H0rc5MPVLp1/DaxQaYm5niD3kgPe0pWRXxihc kHL9MfqOA87bHzjtlx4DK+mSrpBpDlSI9QzOK4/BLXXPi+Z8JCyv/+9ituYY 2GyodRIUJKHNpSxy8Z4x2Ngsrq4pTcKR6M7CksAxeFy2leuuQUKxlbqRa/QY 0GpfvQ3RIUFG5Vy1zIMxuCfs9iJEn4TfBYvbvIvHoL1SVdHQmATTJrUlakNj UOWvH/n4KAk5j7xj/8yNQeiFTCcBBxJE3N8yQ4THwSpEdu8+ZxI+8e7T/rdy HAK/PNx0z4uErSu9jsbbjIPgZ9MNSiEkpA4Wd291HYdTVcdWfAojgbeY32v4 /DhsuWy/2DqShJLdsSF7YschUXQ6fkMCCbpehc/4v4zD1GLNiMYMEr7N69t4 1o+DmDQnsvApCU43COG/XePwJ33O6EouCbcyKk4W806A86GQvI+vSOhr+ans uWECAotj7Pd/IMHf1bK2ZesEjH7sSPGrIkF6qtF/9/4JGAjlXor6TIKRVHvz ytMTsEHkF7/PdxLu7RqKb0mdgJ2XAkU0m0jQqvfevjt3Amppxdz0FhKqnCZH i95OwJaCpWcW/yNhNGDeLPbPBGRHPRII6SQhTCyQh79rAjIWP1V42E0C677A M4+RCfhAC3md2UvCnkIR4d2ik8DxFQs5wyahY1t0cZHcJNgX0XPXkyRc+CF5 cuXqSUg7GGvdPkiCuF3C0tj1k/D0qNO188MkPObIVfBtnYQr+c5a3BESDPwe eHuYTcKk7E1r+zESagWVlVtsJiFI7Q/vy3ES5lRU/YvOT0LqrEaZNJeEmLxs 9ZXXJkEIohqWT5GgiprNMTGT4KBW7MqYJqG0Oj+cL3USUp7MneOZIeGg9UY9 j2eTUBdSMv2Z4oGe1z3NbyYh3f7UhP8sCVfPQvyuT1R876+eYs6RsGzB+21F dZPwRSThQDrFz6K3j67onISIIJOXy+ZJ2M78lBozPAmmu0wizlDc9HSPGR8P F2IdD9W/othL7/u8uwgXFF9tuddJ8eKqAznNslxgLuj9OUNxksWfI7tWcUF/ yCBwlmLtduslRTpcSOfflN5N8SePv0UrjLignDKJrym2nbN3jjHlwqBYpKkv xePhXdJ8NlzIP6vxm0FxhMypD+4uXLBtWlabQ8WrlM4+0+zLBf5ViVtXUVyk 7am0K4QL4SyaejiVr2nZ6LfC21y44UrerKPq0bXP98qKFC6Q0bknRCj2a55e G5PDhZix9sI1VD1pLleaeN9wYahhZaAOVe+MSb5w949UvA6qlSspPTaHhOg2 /+bC3M3MYP5JEn7ShXpMOrhwQvcMUUXpdyrlRlzhEBf65Rx9zlH68mhIbFsx z4Xgz64ZYqMkxL2JGbktPAURhdcOR1P+eF9319R95RScNTXQM6P8FGNvMTq2 bgrc6/Kscyi/ObFFE/xwCjQ8HvRP9JMgxH/1b+jhKbDNkbE27aH0uKEfKO40 Bd2FHQbHukjIXja2It5rChYkJUdbdVD5qzu7pYVPgVhcb4J0K1WvIkWaWsIU rOmj761rpvaDUePL52lTkHLn5MWgRhLiD++dLXk7BfabLbam/SahJUg7op6c ApXJuSIutf9yxTgatjNTsD2H15rxkfJTYvqPrsXT8Eay3V+jgup3ubJyo4rT cD/lsv3SUhJcG3myRA9Mw0za0bzZPKp+TsV7Ym2nIV/V+FnIM2q/DJ0Zkjs9 DVM1JTjzhIQXAt0bV4dMw7UAmmp2GglTmtWVRoXT0Pd32mJFLAlf3oS4fCqf hh2mL9ZLRlP+MkYxs+/TcGE+pYQTQcKWoy8O2vRNw2tN65qTQSRcv36n01d+ Bjya9m7wpfqd9F9H/hz/GUg4Gf940oSEnpPMdO0bMzAwFrctbxvVf0f/mLy+ MwPl63iLLZGEo0K7b1XlzcAML/znSvXjNB0txfaOGRAIlXiezSBBM2Ju87Kd s3BfYhEric0BrZMzx7MPzkLHAwGdxC4OaG+fCjM8PgsjOboL/f9yYD3P+G/X K7Ng7/u+f/F3Dhh4s93evZyFma7Hx8VfcGDHkaakU0pz8FCe8XPKiwMmug3l 8xpzoLXoyvLYUxzYJf2nL2bTHOz2PxMuY8+Bvd9+rCcs50CqKjG71ZQDB7Z+ rqFHzoHFDIfrvJYDtmuK595Mz8H0gQSbE21ssBcsVNkvOA9DT1W5j/6wwaHr pUm39Dy03boo8bWGDU7Jz+MktOaBEfahrrWYDa6SmepOzvOQ4mMezH+LDedm 4m3Efs2D5Gxa7y4DNhh36f397xAP5l8tU2IFDIDTj58D147yoM+hoOwonwEI KnWbanHgQf8AH0vy5AC8S0yVuuHBg4Vtq9ac2z8A+vuEdveG8uCLZrnE50oD oFFYX5jyhgd5HLnqf4h+WBrme4umxIvtg0t8lNl90L22YNvoAC8+J/zDIqV6 QenSaUe1EV5c8bpUa5lALxz5ohJsz+XFGpndpdHjPfDt1O3ybwv58Ojkx9xD dT1QnO5hlM3gw2Zx3kGNhB6IUFTb4mTKh4prm3bKLu0BTalkgz8v+DDN9Zip nFA3XJwJXffmAj+qCEcMSFR2wCKJnpMHVBfiJaboG/BugYXVC42XvBXA/O/J kNZSC9M12a63PBdjkL6cwTqBAug8lCPDYAjhgoXrN/W5E8SW626nq7KX4P1r qaZfBX4SsXOFgkbHRLA8Nkb9w5pWQnjlB7dYRxEc1yv/GW3YSgTt+/6j20UE /WajzQ9YthJeyb33I3xF8FjHi7TiwFZir6Gs1u9oEUzZsUZ2SUMrIRBy0epk uQiGaa2Iehb0j/AV2fwkcrUoavH7TGV+byOs5d/vbhoRRfUht5w4406i7/TQ 8MspUVRKe1Wy8Egn4Usw70TyiqFEWu1yR49OIs7er9NQTAwbNn1+M5rQSdRm 6FzOVBXDoXBbN4veTsJEJ+PZWTsxXDS9b4n4tS5Cb+8NmsRXMbSfb+IuLOgm Pj54Xdj3SwyFaAV/HKq6CcuhXpvyJjFcbBBinNvQTfjEbH9yrk8M9WhuDyTn qXH2vCLJJ46xeY3vBY17iJh1RcUWDHGUEa5yTPzRQwiXqg0oWYjjcSOHyemW XqLtjlnVFUtxLA0+xlfF6SVenTmX1mQtjqvcTxYGzPUS9ivfH4mzF8el4Ztn 3ir0Ea9uWlULeojjL56hHEvrPsL+aGjOYJg4ijXE7BL50UcUTnV5lJSKI5eR 7a2a30/c/Cm8R66cWu9H4rn4d/2EQ47Wat9KcQw/8KRpsKafELW//E+zRhyf ZgenO/dS4x+l9j9qFMc7d+z0XRQGCNGEresixsUxROOFQ37AAOGokzp2eI0E FlmEN8vos4mwm1JilzQksHh2xybjbWwip/v6qmQtCRyTnL1rv49NjN9xt+7a KIGXv/hscnBgE9fnDN55bZPA3q/mjuFhbCKrsi4i3EYCl80+oXn9YhM1rF2P c+wkcOPOtVF+LWxi+HwJ8f24BI46OOd4d7MJ/TWPhpe6SGCamMCrtVw28SXK yzLtnAQ++q1hVybHITiHRVTeRklgFp9YgOdRDkF74b+p9bYEHteZizB35BDr hUcP8sdLYMlQQZfSaQ5xpaQhdOd9CexwOewR6schxJUzyF+ZEpjoUmDpk8gh tAYMX3PKJNCs9mhIeS2HOLit4KdEhQQav8szPl/PIS4kreZof5TAbDNbMblW DvHeVFzx4lcJDPt7JWU5m0MceNkcsqhJAs0VztsdFCCJc5d997PGJPDqGpEZ JR2SWDXgJ3xtUgJjdjB3yeuTxJ/DARXsaQncJOS4cB5JYqNOmN5rPhouZMVl +e8micn+u0oHxWl4JcN1q7c9SWQcTm56Q6fhw+pr5wtPkIRlZVqc8lIaTpjH cjtdSaIoNVtoSIGGybZJn2fOksSFw8RwuBoN55el7JcJIwnVyvdPh9VpKCQk dSH9Jkk0aFc5HdaioV31gW+M2yShL/69YYUuDdVijs5U3SWJmYq2snfbabgv SHdG6ClJPNXu9ltlQvHMJqO2ZyRxJLV/feQeGjYELkpJfUESby+NZh05QMPV n1hxLcUkcVlb4PaEHQ1rR28wjKpIQj1VaI+NIw2Plzr/sP9MEi1iYos+ONPw zszycJcakoD+pRdvudOQh+/W1IafJDFoKa/D9aKhjp9YzvRvkkiuYJHHztJw 6RfhIxn1JMGTutph7SUaHrVpfPCmhSRyxdQVYvxp6HWRR0vlH0nYXtKqmwqk oeGKPa+820ninaX+rqow6qtr+ZuYb90k4VUBCzVu0pCR+aq3sZcklLS3ErHR NDSXvaH2rZ8krortWeeQQEOnVA2PMyRJaF0yY3+8S8NyW2sPpSGSaOuzSP/v AQ2/HmVYFg+TxC1LK7v4VBp+yjJboT9KEkYVNnJzj2hYLPy3Pm2MJEbXOfw6 nklDjdlSz4lxkkhLORH5+SkNp7L42FqTJHFAzHWnVi4Ne6cIk0NcklhwyYP/ zgtq/evLw45PkUR+n/fb+QIaij/2yLKaJonjlufPORXTcKsNmaU7QxKSFZc0 q9/SMP0MTzjvLEmUr7vav+4dDePs6CZ5FHunhDxKLKdh6JGM7l1zJKEiFn6M t4qGZ/ebOtRQ/NMvUsb5Mw2V+5++0Z0nieC+2z9qaii9rnuM3KA4qK7z0/R3 Gq44HPW88v9cvqFs1S8aruV6bGVTHPg8tMjiD5Vv0dSTOYqvJjXmXm2k/Juq 1DNNcUD42oxnLTRM9eyb7KDY3/fKg6Z/NCTu2rYUUXzF8XucYCcNm3clJ1yg +PJ+5Zvre2ho5fB95UqKL4FPsEM/DZvOr4ogqHj91lT6RXFouLy564MRxRdl ZLzfDtFwIMrq1wsq3wsCLi59ozTM3dtSLEbx+ZE39ksnaXg5uP7sIapevq2i Vlunafj0cKlgOFXPc9W2Zp5zNLyVRvfOoup9tjhvRxIvHfVGbQteUnr4pC/A zwvoKBBu+P0ppZd3zMENk4vouPvYsQ+RE5Sf3KaW7xel47TsCgMZSm9P690K VyTouNj9aCkxQhIeO5Ikn0rSseVPH2s/5Rd3JUN+ATk6HnfNuaBH+ctNLGZK i0HHXZF1PlFsimc6h44p0lG9yMCslvKj6+/Q1qKVdLR7sCZbpocknMO+vz2t TceDyoFyr/6SxIlzyi8TN9CxXk5Sz7WZJJyO+zyt1KOj6kMRI5FGyi+bZe4q bqHjxl1bZOnUfrIbtvX9tYeO5ZMPhNSp/Wj7N8+dz4yOJ4olMg2p/Wr7ZcEJ jQN0jC8wBKMPJGHzOMMi1IqO9KpLVosJkrCyGtTadIKO919bBaU+J4nDxoaq p07RsdXMY9nSHKo/accoxp+m4/4/7/PPZZHEQdGN1Cs0HX85m4vNppKEedkV dpo/He89kotcRvWf3apiGcKJdKycL56Kcafmq3LPWH2fjsJGI1fbT1H5OH3L 2J5Mx+plfXxMJ8o/qdGZ/o/pyOP1dOzIEWp/yEo+Gc2jo+/ZDVOjxiQhv0Qm p+kzHW+aFxkmyJHE6szzOdwaOppejC1+JUUS640bcqRr6Sh7p+pouRhJmAbe fWb6h44Np7O0H/FT/p5WeF7eQUfDgzIRwX0cYqBfKf/pLB3Hn9duLc/nENzQ wPyPPJIo3hqlX5jNIQRWdeR38Utilvzje3cecwjW8UcvWUKSmHpW8ZlKAvV8 aFr5KkZaEtcPt0UducAh3nxZU+ynIYn9djd3HNLnEDey15futpNEU4OF9YnP 2MQK21V22ccl8a1jcuH5R2yihCbLJ+osiRMGZj5Gd9nE4LlZo69ukphnrfQw PoRNHDAsrzTzk0StMvFkoyNsQr7O7OvBeOr/azK/5PKziac8p1tsq6l4/K6f WGY0QGx7YXPl3TdJLF98b7n9+gGixdGUpfRTEsMUL1XErx4gxD+vs29vkERX 5f9+VYgNEN6xU52OPZIYXK61Wqexn9BXvc4+xS+FZlG0ISu3fuKT+cMZb10p HDJMEBYP6SPygtsa7hpI4e4873sD5/qIxFeKRWUghSZ/9r4tPNlHnJJL9hHf LoWBqrt/Mff0EYLt9zhP9kvhc9qPQ6r0PmKHV1xrm4sUDga73954r5f4EBVa bpokhW+qLfZfTO0hSqvdI9bwS+PIVVL8Q0AXoS2JS+IFpHF7EO3DPrcuIsNK PIxXSBrz7AwVPll1EZHdudd+i0tjr0DV8I11XYQN7/BVf4Y0SvTNvjfv6CRm tL3P/9SVRnKIkddh2Eno3vN18vOQxrqStc8CR9uJvFMB+KVZGq8tHImVWfuP YOyCjW75S1GmvaTC8ehvwmj+SOADv2V4/UDGq7CYcuL4b2K40EIGX5547dDU mAmbA66+tGPJYoVT8nfj81/hlvvQa+FuWRRT2p/f5NQIbf6Puh6UyOE15o0k Qc92OPvtvkf5Ozn0X2PZNB3UDkKsuKmecjmU1QkSaUtoBy0iRFTrsxxaq3cs Cytth8DZE+vL6+SQvjvrkrFYB6hcUA3pGZTDchHXbtesDnA682y5prI8Sk6k SorVdcJUWfozixXyeKv4fohzbyfcpCXrXlgtj8F72J75051Q8Dxqz3sNeRT9 5HhNmdUFiwY9fSwM5FFKOvLMUecuyHDV/nDeQh7bz62/u368C/odi46XXZPH 30bft2gt6gHhT4aRquHymGlcmLZapgfWqH8punVTHsfSdkeIqfWA63iLqEOs PNanJ758vLcHBkIWFPE/lMeTb+wX88ZSnL5PZFsJxZ5xsVaMXuD0dr6sHJXH spQzOruU+0B0r3urxqQ8Sthc5J3W6gP1vEmhhGl5rEgX/BJv2AduF5bYOvMp 4EWJCxoJdn3AEdIUEhRXQD31kbG8pD4g1fxsTNQU8POUxb/tUv0weJq2qNpO ARvdfi9LGumHY2a9btPHFfBa/zZV7nw/1GiX/l7trIBl74hgQ+EBeDp9+nGI mwI6RVdOP1AZAOfQj9vwogIG15mKBB8YgNY0/6AXMQoo9SRdeyh3APZdPzTw L14BT278YbX2zQCUuKgfEL9LzZ83a36ocgDuaTUpn06h4huQe+jdPACWpRve r8hRQJe8xKEfi9nwtZHDm1ipgL7ZkTbMI2zYTHw4WfVJAT2uxl3WcWLD09R7 38erFXDvCTEBXQ82hJ00STH/qYAshojW4iA2bJ94hCL/FDC5pXxfRRYbSmhH r/hPKeCclpPFm1E25Oz8POO0loE3OGTsiQscUGvJPR38HwOrX+UVTQZwIMMr ruXhOgaKdGianQvlQOpdu9JWXQaua17/0SCBA7HsiUCr7Qz8tfOmbU8+B/yi VZbsPcZA7x26LiV9HOAuF/JztWdgyzOiVGWIA2eLyYEwRwZqBLhs953ggEd7 cU2lCwM9TJtqu/hJcNAxu7XFl4HxGYbX/smRYFJ/SUYnmuJT3wxMdpBQedoh zDyGmm9IvzVgDwlb+XZOe8YzsB0fXXy8n4TNapLNOfcZeJiculFwhATNS1nJ q7IYqHd8/vOQOwk5ElESxtkM9BKxaIjyJkH1sc9Vx1wGvrbprZI/T4LKVzye WsDAbUs+608HkLBMsW6l/HsGqrb0rq29RULsyzcJehUMXPjOaqQujgSaSarg 4Y8MfGy4LuhTIgnCZ073x35lIL/ujTmfVBLmy/mfiTUxsJh+S6orlwS/w71M 9b8MNLe+O2aWT8IkuyZqdxsD08+zIrNekTAsnegV2sPA2Y8hXYolJHQ5/7ee f4yBNp65g6s+kUBKpBs7TjKw7pte0vpq6v5ihcOV0wycOLsiYu03EgRFlviF 8zEx8lABq+4XFV9BwA32QiYOP9rzPvQPCbLHJpP2CTLxleOapOWNJKx53vmO Ls7Eb3GHx5a1kqBjdeSHD52Jbr4Kbh5tVD35f3TUSTPxWg592fP/n9dZlC66 p8DEgyd0uOweEg7PrZeZYzHR4K/l2r4+EuzTs1VtVZi4IiEp4usACWe4iXtU 1Jg4e1C70mKQyj9V4liIOhOfc+cejA+RELTrukePJhMrl/1L8h8h4cboXICJ DhNPt5SWjfz/PPK+z+2nG5moWRMjZDpOwoPtA2miBkwkPjj4xE6QkEHaF3gA E2Nnt/KVT5LwPKG+staQiRmCtrmNXBKKt5jWa2+n1q8fuNw0RcL7voq+uJ1M 3Px63emKaRK+3N40M7mbiRsP3bhwZ4aEnwb5ItamTDzw0v6xxSwJzZ2qzLfm TOTXXzs4RXHXzZT/mIeYqPT07KHQORI4G5YZBlgxcatafSPPPAkTrTfN248y UVhv3s+OYp6whY7b7Jh4qOih3hOKBdf5nU0/zsRLvdGSLRTTmoavCToz0f3U I6EpimWDT95xcaH00P8kP0+xsnprVrUbExe9nTGh3rdhTd3BNxpeTDTM3hlb /v/zav/q6mgfJrqo1UwHUbxp9da/I75MFC94dVn9/+fjtcWDFn5MNFHbynhH xbvvoiZf4RUmhle9azCg2FIlgy4byERn1cj8FCo/u2rGcr8QJprvlntKUvU4 dTZ2fUsoE82ud75bQfEZpvAOvMHExLhL49up+vlVXT2cGsVErRuBJqZUfYM8 uacWxDDx4t/0EkOq/jdkPfyc4pmYd/SJOYPSJ+59142qRCY6aR8Taqf0THI9 +kA1iYk9/WntUWMk5JaYlHHSqHg63/JmDJNQdOLdD9MMJtbMPNohQfmlTHxj Z94Tyn9JOwscSEov++WLz+UxMYGbLFTeT+m15J5M/UsmBu0cGPneS+mVT1PT L2Li0r1Kgh+7KX0EePbOE0yk7dR869VO1ftJw+1r1Uws+Jbsm19PgtePpLZn 35joGKMnoFdHwsdpe80/P5iYdKfhzZOfJPju7q9Z1cBE28sJT8y/klDHnhb8 1EX5VSG75OV7EtSlSy2Heyl/V71/f6+UhJDNQemybCZaj03Mu78lYX2k8DaX ESbqN+Usbyyg+luCKArwsNCQLLMtzyDh7FOPE4m8LEw51mgp8oiE1NLam+r8 LIzw33THMIWEqZ645oMCFN+2mfG8Q0KmLvNiujALz8Vtie64TsLiRo2CHTIs /HZ7rbHcccofZHRzkywLhb/oTa0+RoIt/+gCT3kWOrxQUFxuRUKhWqH5HSYL 60VGfrfuI+GEHw71Lqfmv39qk5seCRUKZmvCtVj4SLFuqH0J1Z80X5gztVnI +Cvt/XohCYztUhdf6LBwprAz5co8B3zc6j82bmRh0ZoKsZphDqgQdifWAAv1 T+5MC/zDgUBbr9TqXSxcFexmYZHMgdE5+cTAPSzk8vTzWVPPE6f7ldG6+1j4 4arFHpMoDpg0yAc82s/CvqXJ3zqvcIB2oPLYpcMs3Hlmn56EDQeSjeUV1p5g 4YPILBmQpsa7KiTbnVkYHfDv0j9hDgQFeQrfOcVCl1rJztP8HDjxrmJ6gRsL My+o/tk6yAZ1fc+GJm8WpopmSQ9UseGtekV8xFUWLtJSP7jiLDVe7RFpGMTC ER3JXa4ubEh2kbs2GcxC52+yRfds2RCU4XH2eCgLY35uYr01YcNuJTkLgygW Oo2tn5lRYEODlAdt4D4LdcSZq/1LB2BiRubGriIWDps0e38fo94/RnaFcotZ 2HZbN/xlTz/09F0KTn/DwuslAcPBjf3Q8OffJf5SFrY8lq7vL+2HkvxMt9cV LKz6yRStDuuHIFddU7WfLFRaqx3Ilu0H8SZL+hKShVdMA0sG1vTBqtfxdz6p KKLC/q+CW2R7oMrTI996hSIW98T990awB5xX7fw6sFIRDRalx63idkNm7NQC MTVF9Imw3vmrrhtUPY56mGsqYlndKucrMd2gtlzZuGmTIjZc3pDCXNIN6pE5 o+yDirju0jO/BnYn6Dh82CcRqojXNn05nRneDqX3pJZ9GVNEWWHa4PbrjbCq VCpsr7USNrIenBbQ+wp7Fg4o1JYoIS3nRYtrehYMNna8/6iqjKpb5Gc8N5UT 7vdu0EMjlbHupvHpGOPfxBdeFdIkWhnNZsq/G/r9JlRPFH8Svq2Mt/uvNPzN /U10aHYHRMUpY7vK1vWdMnWEZSUOxt9XxofzspnL+usIw+GRL4+fKKP0r1d6 86H1hPQOq5DySmXkORO7JL+wiTiTPWgX8lEZX4z9OdLb10R8o13btOOzMs5d X2W9WKGZCG9+Mfq5RhnFqzaU8wU0EzxnRBx+/FJGLbkZ8U7jFqLvfim0tSvj sqCv0pu//yWM+Q/JpXUq4yeuPmOOp5VIc2aPO3Yro1GVwW9TmVbi2DrZnN4+ Zfww8bzv285W4mfVGfmhIWXUfbUoUTerldBUF5rMG1HG39Z3P60rayVu3k7+ 4T2mjHzf1i1RbGgldtpUh01OKqPo60Svn4L/iMfvHZyKp5SxbHmcd6LiP2LB 6qktfjPUescHDcx1/xF2N6MUNs8po47V4aI5039EycgK7vw8Fb9GeN0953/E /wC2as3U "]]}}}, AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948], Axes->{True, True}, AxesLabel->{None, None}, AxesOrigin->{0, 0}, DisplayFunction->Identity, Frame->{{False, False}, {False, False}}, FrameLabel->{{None, None}, {None, None}}, FrameTicks->{{Automatic, Automatic}, {Automatic, Automatic}}, GridLines->{None, None}, GridLinesStyle->Directive[ GrayLevel[0.5, 0.4]], Method->{}, PlotRange->{{0, 0.00016868653312820215`}, {0, 4.}}, PlotRangeClipping->True, PlotRangePadding->{{ Scaled[0.02], Scaled[0.02]}, { Scaled[0.02], Scaled[0.05]}}, Ticks->{Automatic, Automatic}]], "Output", CellChangeTimes->{ 3.6691078845636344`*^9, {3.669107939816795*^9, 3.6691079969520626`*^9}}] }, Open ]] }, WindowSize->{1259, 951}, WindowMargins->{{202, Automatic}, {-69, Automatic}}, FrontEndVersion->"10.1 for Microsoft Windows (64-bit) (March 23, 2015)", StyleDefinitions->"Default.nb" ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[557, 20, 7904, 194, 739, "Input"], Cell[CellGroupData[{ Cell[8486, 218, 2445, 76, 229, "Input"], Cell[10934, 296, 130, 2, 23, "Print"], Cell[11067, 300, 512, 11, 21, "Message"], Cell[11582, 313, 193, 4, 31, "Output"], Cell[11778, 319, 132, 2, 23, "Print"], Cell[11913, 323, 512, 11, 21, "Message"], Cell[12428, 336, 12783, 222, 275, "Output"] }, Open ]] } ] *) (* End of internal cache information *)