(************** Content-type: application/mathematica ************** Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 34646, 1079]*) (*NotebookOutlinePosition[ 35709, 1114]*) (* CellTagsIndexPosition[ 35665, 1110]*) (*WindowFrame->Normal*) Notebook[{ Cell["Differential Equations Lab II", "Title"], Cell[TextData[{ "This lab ", StyleBox["explores some of the ODE solving capabilities of ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Mathematica", Evaluatable->False, AspectRatioFixed->True, FontSlant->"Italic"], StyleBox[ ". The notebook is organized into several sections. Each section is just a \ group of ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Mathematica", Evaluatable->False, AspectRatioFixed->True, FontSlant->"Italic"], StyleBox[ " cells containing explanatory text like this and actual commands like the \ one below. To execute a command, click the mouse in the cell containing the \ command then do a Shift-Return, that is, hold down the shift key and hit the \ Return key. Try it on the cell below, which evaluates a complicated \ integral.", Evaluatable->False, AspectRatioFixed->True] }], "Text"], Cell[BoxData[ \(Integrate[\((1 - 2*x^2)\)^\((\(-3\)/2)\), x]\)], "Input"], Cell[TextData[{ "To begin the actual lab work you should open the first group of cells \ below by double-clicking on the bracket to the far right of the title ", StyleBox["Finding Exact Solutions with DSolve.", FontWeight->"Bold"], " The bracket will open up to reveal lots of interesting text and commands. \ Just follow through the section, reading, executing commands and (hopefully) \ thinking about the results. When you get done you can close up the section \ again by double clicking on the expanded bracket at the right edge of the \ cells and move on to the next section. Have fun !" }], "Text"], Cell[CellGroupData[{ Cell["Finding Exact Solutions with DSolve", "Section"], Cell[CellGroupData[{ Cell[TextData[{ StyleBox["To find the exact formula for the solution of simple one-variable \ ODE's, use the ", Evaluatable->False, AspectRatioFixed->True], StyleBox["DSolve", "Input", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[" ", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox["command. Of course, not all equations can be solved. When using \ this command, you have to be careful to use the double equal sign and to \ write all the y's as ", Evaluatable->False, AspectRatioFixed->True], StyleBox["y[t]", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox[" in the equation.", Evaluatable->False, AspectRatioFixed->True] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell["rules = DSolve[ y'[t] == (-1/2)*y[t] + Sin[t], y[t], t]", "Input", AspectRatioFixed->True], Cell[TextData[{ StyleBox["The answer is always a list of one or more subtitution rules. A \ list in ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Mathematica", Evaluatable->False, AspectRatioFixed->True, FontSlant->"Italic"], StyleBox[" is always enclosed in set brackets like this { first thing, \ second thing ,... }. In this case the list contains just one thing, a \ substitution rule. A substitution rule is indicated by an arrow, -> and it \ is also enclosed in set brackets. Hence the double set brackets in the \ output. The first (and in this case the only) element of the list of \ solution rules is called ", Evaluatable->False, AspectRatioFixed->True], StyleBox["rules[[1]]", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox[".", Evaluatable->False, AspectRatioFixed->True] }], "Text"], Cell["rules[[1]]", "Input"], Cell[TextData[{ StyleBox["To extract the solution as a function for plotting or further \ analysis, you have to apply one of the substitution rules to the expression ", Evaluatable->False, AspectRatioFixed->True], StyleBox["y[t]", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[". You apply a rule in ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Mathematica", Evaluatable->False, AspectRatioFixed->True, FontSlant->"Italic"], StyleBox[" using the cryptic", Evaluatable->False, AspectRatioFixed->True], StyleBox[" ", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox["/.", "Input", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[" ", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox[" notation. Give the result another name, ygen, so you can use \ it later.", Evaluatable->False, AspectRatioFixed->True] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell["ygen = y[t]/.rules[[1]]", "Input", AspectRatioFixed->True], Cell[TextData[{ StyleBox["Here ygen is the general solution with the arbitrary constant, ", Evaluatable->False, AspectRatioFixed->True], StyleBox["C[1]", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[". Note that this is a linear ODE and so the general solution \ consists of a particular solution plus and exponential term.To pick out a \ particular solution, just apply another substitution rule, substituting a \ particular constant for ", Evaluatable->False, AspectRatioFixed->True], StyleBox["C[1]", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[".", Evaluatable->False, AspectRatioFixed->True] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell["\<\ ypart = ygen/.C[1]->0 Plot[ypart,{t,-10.0,10}];\ \>", "Input", AspectRatioFixed->True], Cell[TextData[{ StyleBox["To see lots of integral curves , you can plot a list of several \ solutions with different values of ", Evaluatable->False, AspectRatioFixed->True], StyleBox["C[1]", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[". One quick way to do this is to use the ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Table ", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox["command to generate the list of particular solutions. For some \ reason you need the ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Evaluate", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[" command inside the Plot command.", Evaluatable->False, AspectRatioFixed->True] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell["\<\ sollist = Table[ ygen/.C[1]->i/5,{i,-5,5}] Plot[Evaluate[sollist],{t,-10,10}, PlotRange->{-2,2} ];\ \>", "Input", AspectRatioFixed->True], Cell[TextData[{ "You can use the mouse to make the plot bigger. Note that all solutions \ converge to the particular solution we plotted earlier. This is due to the \ exponential decay of the term involving ", StyleBox["C[1]", "Input"], ". In a case like this you would refer to the particular solution as the \ ", StyleBox["steady state", FontWeight->"Bold"], " solution. " }], "Text", Evaluatable->False, AspectRatioFixed->True] }, Open ]], Cell[TextData[{ "Go back to the top of this section and try some other solvable ODE's. For \ example, \n\t\t", StyleBox["y'[t] == y[t]^2", "Input"], "\t or \t", StyleBox["y'[t] == -t*y[t]", "Input"], "\nJust replace the differential equation inside ", StyleBox["DSolve", FontWeight->"Bold"], " with the new equation. You may also want to change the plotting ranges." }], "Text"], Cell[TextData[{ "Some other ODE's are too difficult to solve explicitly. Try the equation \ ", StyleBox["y'[t] == (t^2-y[t]^2)/(1+y[t]^2)", "Input"], ". The program just spits the problem back at you." }], "Text"], Cell[CellGroupData[{ Cell[TextData[StyleBox["You can also use DSolve for second-order ODEs. The \ first example is a typical damped, unforced oscillator.", Evaluatable->False, AspectRatioFixed->True]], "Text", Evaluatable->False, AspectRatioFixed->True], Cell["rules = DSolve[ y''[t] + y'[t] + y[t]==0, y[t], t]", "Input", AspectRatioFixed->True], Cell[TextData[{ StyleBox["As for first-order equations, the answer is a list of one or more \ subtitution rules. The first (and in this case the only) element of the list \ of solution rules is called ", Evaluatable->False, AspectRatioFixed->True], StyleBox["rules[[1]]", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox[".", Evaluatable->False, AspectRatioFixed->True] }], "Text"], Cell["rules[[1]]", "Input"], Cell[TextData[{ StyleBox["To extract the solution as a function for plotting or further \ analysis, you have to apply one of the substitution rules to the expression ", Evaluatable->False, AspectRatioFixed->True], StyleBox["y[t]", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[". You apply a rule in ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Mathematica", Evaluatable->False, AspectRatioFixed->True, FontSlant->"Italic"], StyleBox[" using the cryptic", Evaluatable->False, AspectRatioFixed->True], StyleBox[" ", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox["/.", "Input", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[" ", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox[" notation. Give the result another name, ygen, so you can use \ it later.", Evaluatable->False, AspectRatioFixed->True] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell["ygen = y[t]/.rules[[1]]", "Input", AspectRatioFixed->True], Cell[TextData[{ StyleBox["Here ygen is the general solution with two arbitrary constants, ", Evaluatable->False, AspectRatioFixed->True], StyleBox["C[1] ", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox["and", Evaluatable->False, AspectRatioFixed->True], StyleBox[" C[2]", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[". Note that the solution is given in terms of complex \ exponentials. The exponents are cube roots of -1. To get a real solution \ out of this you have to select the constants to be complex conjugates of one \ another. Since this is rather inconvenient, we will adopt a different \ strategy, namely, specify some real-valued initial conditions directly in the \ ", Evaluatable->False, AspectRatioFixed->True], StyleBox["DSolve ", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox["command.", Evaluatable->False, AspectRatioFixed->True] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[{ \(rules\ = \ DSolve[\ {\(\(y'\)'\)[t]\ + \ \(y'\)[t]\ + \ y[t] == 0, y[0] == 1, \(y'\)[0] == 0}, \ y[t], \ t]\), "\n", \(y1\ = y[t] /. rules[\([1]\)]\)}], "Input"], Cell[TextData[{ "This still doesn't look real but it actually is. You can force ", StyleBox["Mathematica", FontSlant->"Italic"], " to put it into a nicer form by using the commands ", StyleBox["ComplexExpand", "Input"], " and ", StyleBox["Simplify", "Input"], ". The former tries to expand out the expression into its real and \ imaginary parts under the assumption that all the variables in the expression \ are real." }], "Text"], Cell[BoxData[ \(y1\ = \ Simplify[\ ComplexExpand[y1]]\)], "Input"], Cell["Here is a plot.", "Text"], Cell[BoxData[ \(\(Plot[y1, {t, 0, 10}];\)\)], "Input"], Cell["\<\ Another independent real solution can be obtained by changing the \ initial conditions.\ \>", "Text"], Cell[BoxData[{ \(rules\ = \ DSolve[\ {\(\(y'\)'\)[t]\ + \ \(y'\)[t]\ + \ y[t] == 0, y[0] == 0, \(y'\)[0] == 1}, \ y[t], \ t]\), "\n", \(y2\ = Simplify[\ ComplexExpand[y[t] /. rules[\([1]\)]\ ]]\), "\n", \(\(Plot[{y1, y2}, {t, 0, 10}];\)\)}], "Input"], Cell[TextData[{ "Perhaps the best use of ", StyleBox["DSolve", "Input"], " is to find explicit solutions of certain special second-order linear \ equations with non-constant coefficients. In the 18th and 19th centuries, \ many such equations were solve via power series and their solutions were \ names after their discovers. These are the so-called \"special functions\". \ ", StyleBox["Mathematica", FontSlant->"Italic"], " has been trained to recognize many of these equations and can give you \ explicit solutions to them in terms of these special functions. Here is \ Airy's equation ", StyleBox["y''[t]+t*y[t]==0 ", "Input"], "and its solution in terms of \"Airy functions\"." }], "Text"], Cell[BoxData[ \(rules\ = \ DSolve[\ \(\(y'\)'\)[t]\ + t*y[t] == 0, \ y[t], \ t]\)], "Input"], Cell["\<\ Here is what some Airy functions look like. Again, the solution is \ complex so we first extract real and imaginary parts before plotting.\ \>", \ "Text"], Cell[BoxData[{ \(ygen\ = \ \ y[t] /. rules[\([1]\)]\), "\n", \(ypart = ygen /. {C[1] -> 1, C[2] -> 0}\), "\n", \(\(Plot[{Re[ypart], Im[ypart]}, {t, \(-5\), 5}];\)\)}], "Input"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell["Finding Approximate Solutions with NDSolve", "Section"], Cell[TextData[{ "For problems which cannot be solved explicitly you can still approximate \ the solutions numerically on any desired interval of time. For this you can \ use the ", StyleBox["Mathematica", FontSlant->"Italic"], " function ", StyleBox["NDSolve. ", FontWeight->"Bold"], "The following command brings up information about this command." }], "Text"], Cell["?NDSolve", "Input"], Cell[TextData[{ "We will apply it to the ODE which ", StyleBox["DSolve", "Input"], " failed to solve in the last section." }], "Text"], Cell["\<\ rules = NDSolve[ \t\t\t{y'[t] == (t^2-y[t]^2)/(1+y[t]^2),y[0]==1}, \t\t\ty[t],{t,-5,5} \t\t\t]\ \>", "Input"], Cell[TextData[{ "As before, you get a list of substitution rules, this time not given as a \ formula but as a so-called ", StyleBox["Interpolating Function", FontWeight->"Bold"], ". This is some function (whose formulas are hidden from you) which \ approximates the solution over the given time interval. You plot it just as \ you did the particular solutions above." }], "Text"], Cell["\<\ ypart = y[t]/.rules[[1]]; Plot[ypart,{t,-5,5},PlotRange->All];\ \>", "Input"], Cell[TextData[{ "To draw lots of curves, make a list of solutions with different initial \ conditions. This is a bit more complicated than just choosing different \ values for the arbitrary constant as we did above. We will start with an \ empty list and append solutions to it using the ", StyleBox["Append", FontWeight->"Bold"], " command The solutions are obtained by repeating the process we just went \ through. We will make the machine loop through these steps automatically \ using a ", StyleBox["Do", FontWeight->"Bold"], " loop which runs through the initial conditions ", StyleBox["y[0]=-3,-2.5,-2,...,0,...,2.5,3", "Input"], ", that is, ", StyleBox["y[0]=i/2", "Input"], " for ", StyleBox["i=-6 ", "Input"], "to", StyleBox[" 6", "Input"], "." }], "Text"], Cell["\<\ sollist = {}; Do[rules = NDSolve[{y'[t]==(t^2-y[t]^2)/(1+y[t]^2),y[0]==i/2},y,{t,-5,5}]; AppendTo[sollist, y[t]/.rules[[1]]], {i,-6,6}] sollist\ \>", "Input"], Cell["Finally we plot the solutions on the list.", "Text"], Cell["Plot[ Evaluate[sollist], {t,-5,5}];", "Input"], Cell[TextData[{ "Now see if you can make a picture of the integral curves of the equation \ ", StyleBox["y'[t]==Sin[Pi*y[t]]", "Input"], ". Just go back two cells and make some changes. Can you explain the \ resulting picture by thinking about the phase line ?" }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell[" 2D Phase Portraits", "Section"], Cell[TextData[{ StyleBox["NDSolve ", "Input"], "can be used to produce phase portraits for second-order ODEs or \ first-order systems. Here is the second-order damped oscillator y '' + y ' + \ y = 0 written as a first-order system" }], "Text"], Cell[BoxData[ \(rules\ = \ NDSolve[{\(y'\)[t] \[Equal] v[t], \(v'\)[t] == \ \(-v[t]\) - y[t], y[0] == 1, v[0] == 0}, \n\t\t{y[t], v[t]}, {t, 0, 10}]\)], "Input"], Cell["\<\ To extract the solution from these substitution rules use the \ following \ \>", "Text"], Cell[BoxData[ \(X\ = \ {y[t], v[t]} /. rules[\([1]\)]\)], "Input"], Cell[TextData[{ "The answer is in the form of a vector ", StyleBox["X", "Input"], " of two functions, representing the y and v coordinates of the solution. \ Now you can plot it using the function ", StyleBox["ParametricPlot", "Input"] }], "Text"], Cell[BoxData[ \(\(ParametricPlot[X, {t, 0, 10}, AspectRatio -> 1];\)\)], "Input"], Cell[TextData[{ "To get a phase portrait you just have to compute several orbits and plot \ them together. Before doing this it is convenient to automate the process of \ computing an orbit a bit. The following function ", StyleBox["Orbit ", "Input"], "takes an ODE, an initial vector and a time as input and returns the \ correspoding solution ready for plotting." }], "Text"], Cell[BoxData[ \(Orbit[ode_, {y0_, v0_}, time_]\ := \n\t{y[t], v[t]} /. \(NDSolve[ ode~Join~{y[0] \[Equal] y0, v[0] \[Equal] v0}, {y[t], v[t]}, {t, 0, time}]\)[\([1]\)]\)], "Input"], Cell["\<\ With this function defined, we can compute lots of orbits and plot \ them together to form the phase portrait.\ \>", "Text"], Cell[BoxData[{ \(\(ode\ = \ {\(y'\)[t] \[Equal] v[t], \(v'\)[t] == \ \(-v[t]\) - y[t]};\)\), "\n", \(\(X1\ = \ Orbit[ode, {1, 1}, 10];\)\), "\n", \(\(X2\ = \ Orbit[ode, {1, \(-1\)}, 10];\)\), "\n", \(\(X3\ = \ Orbit[ode, {\(-1\), \(-1\)}, 10];\)\), "\n", \(\(X4\ = \ Orbit[ode, {\(-1\), 1}, 10];\)\n\), "\n", \(\(ParametricPlot[{X1, X2, X3, X4}, {t, 0, 10}, AspectRatio -> 1, PlotRange -> All];\)\)}], "Input"], Cell[TextData[{ "In the never ending quest for automation, we could go further and define a \ function which takes a list of initial conditions and returns a list of \ orbits. This is pretty easy to do. First, the expression ", StyleBox["Orbit[ode,#,time]&", "Input"], "denotes a function of one argument (the initial condition) with the ode \ and the time fixed. Then one can apply this new function to each entry in a \ list of initial conditions using the ", StyleBox["Map ", "Input"], "command as shown below." }], "Text"], Cell[BoxData[ \(OrbitList[ode_, iclist_, time_]\ := \ Map[Orbit[ode, #, time] &, iclist]\)], "Input"], Cell[TextData[{ "We will try this out on a different linear system. The iclist is a list \ of vectors representing the initial conditions of the orbits we want to \ compute. For some reason you need the ", StyleBox["Evaluate ", "Input"], "command inside the plotting function." }], "Text"], Cell[BoxData[{ \(\(ode\ = \ {\(y'\)[t] \[Equal] v[t], \(v'\)[t] == \(-4\)*y[t]};\)\), "\n", \(\(iclist\ = \ {{1, 0}, {2, 0}, {3, 0}, {0.5, 0}};\)\), "\n", \(\(orbs\ = \ OrbitList[ode, iclist, 10];\)\n\), "\n", \(\(ParametricPlot[Evaluate[orbs], {t, 0, 10}, \n\tAspectRatio -> 1, PlotRange -> {{\(-6\), 6}, {\(-6\), 6}}];\)\)}], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["The Euler and Runge-Kutta Methods", "Section"], Cell[TextData[{ StyleBox["Although the NDSolve command implements a complicated numerical \ method for you, it is easy and instructive to implement some of the simple \ methods discussed in class. Below we define a function which carries out one \ step of Euler's method for any ODE. As input it takes the function ", Evaluatable->False, AspectRatioFixed->True], StyleBox["f[t,y] ", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox["defining the ODE, the initial point ", Evaluatable->False, AspectRatioFixed->True], StyleBox["{t,y}", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox[" and the step size ", Evaluatable->False, AspectRatioFixed->True], StyleBox["h", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox[". As output it returns the point ", Evaluatable->False, AspectRatioFixed->True], StyleBox["{t+h,y+ h*f[t,y]}", "Input", Evaluatable->False, AspectRatioFixed->True], StyleBox[".", Evaluatable->False, AspectRatioFixed->True] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell["EulerStep[f_,{t_,y_},h_] := {t+h,y+h*f[t,y]}", "Input", AspectRatioFixed->True], Cell[TextData[{ "We can apply this method over and over again to step forward in time as \ far as we want. At each step we want to start at the point obtained from the \ previous step. The ", StyleBox["Mathematica", FontSlant->"Italic"], " function ", StyleBox["NestList", "Input", FontWeight->"Bold"], " automates this in a simple way. It is designed to take a function and \ feed its outputs back into its inputs n times producing a list of all values \ obtained. Throughout the process we want to keep the ODE f and the stepsize \ h fixed. Only the point {t , y} is changing. The expression ", StyleBox["EulerStep[f,#,h]&", "Input", FontWeight->"Bold"], " inside the ", StyleBox["NestList", "Input", FontWeight->"Bold"], " command represents a function with one argument (denoted by the # symbol) \ into which we repeatedly feed back values. " }], "Text"], Cell["\<\ EulerList[f_,{t_,y_},h_,n_] := \tNestList[ EulerStep[f,#,h]&,{t,y},n]\ \>", "Input", AspectRatioFixed->True], Cell["\<\ We will try it out on the simple ODE y ' = y with initial condition \ y[0]=1. Of course we know the exact solution is the exponential function. \ We will try to approximate it over the time interval from 0 to 1 by taking 4 \ steps of size h = 0.25. The y value at the end of the list should be e = \ 2.71828 ...\ \>", "Text"], Cell["\<\ f[t_,y_] := y h = 0.25; n = 4; elist = EulerList[f,{0,1},h,n]\ \>", "Input"], Cell[TextData[{ "To improve the approximation you need to make h smaller and n larger. Go \ back to the last command and try taking ", StyleBox["h=0.1", "Input"], " and ", StyleBox["n=10", "Input"], " and also ", StyleBox["h=0.01", "Input"], " and ", StyleBox["n=100", "Input"], ". \n\nThe lists are rather long and to continue, it is better not to \ print them all out. You can simply print out the last element of the list \ using the ", StyleBox["Last", FontWeight->"Bold"], " command. The command below shows the result of using 1000 steps. \ Experiment to see how many steps it takes to get an approximation for e \ which is correct in the first four decimal places." }], "Text"], Cell["\<\ h = 0.001; n = 1000; Last[ EulerList[f,{0,1},h,n ] ]\ \>", "Input"], Cell[TextData[{ "Next we will compute an approximate solution to a more difficult problem \ and plot it. Plotting just involves applying the ", StyleBox["Line ", "Input", FontWeight->"Bold"], "command to the list produced by ", StyleBox["EulerList", "Input", FontWeight->"Bold"], " . We will try to compute one of the solutions to the ODE ", StyleBox["y'[t]==(t^2-y[t]^2)/(1+y[t])^2 ", "Input"], "that we got in the previous section of this notebook using ", StyleBox["NDSolve", FontWeight->"Bold"], " . We will work on the time interval from -4 to 4." }], "Text"], Cell["\<\ f[t_,y_] := (t^2-y^2)/(1+y^2) n = 4; h = 8.0/n; t0 = -4; y0 = -2; elist = EulerList[f,{t0,y0},h,n]; Show[ Graphics[Line[elist]], Axes->True ];\ \>", "Input"], Cell["\<\ You should go back to the last cell and try a sequence of larger n \ values so you can see the method converging to the actual solution.\ \>", "Text"], Cell[TextData[{ "For fun, you can make a movie showing the convergence of Euler's method by \ making a table of pictures and then using ", StyleBox["Animate Selected Graphics", FontWeight->"Bold"], " command in the ", StyleBox["Cell", FontWeight->"Bold"], " menu. First we produce a table of 20 pictures with n values ranging from \ 4 to 80 and appropriately chosen h values." }], "Text"], Cell[BoxData[ \(Table[\n Show[Graphics[Line[EulerList[f, {t0, y0}, 8.0/\((4*n)\), 4*n]]], \n\t Axes -> True, \ PlotRange -> {{\(-4\), 4}, {\(-2\), 4}}], \n{n, 1, 20}]\)], "Input"], Cell[TextData[{ "To animate them you have to select them by clicking with the mouse on the \ blue bracket on the right of the cell. Once they are selected, go into the \ Cell menu and choose the item called ", StyleBox["Animate Selected Graphics", FontWeight->"Bold"], ". You should see a movie of your 20 plots showing how Euler's method \ converges. To stop the movie just click the mouse. \n\nIf your movie goes \ too fast or too slow you can fix it as follows. Go to the ", StyleBox["Preferences", FontWeight->"Bold"], " item (depending on your version of ", StyleBox["Mathematica", FontSlant->"Italic"], ", it may be in the ", StyleBox["Edit", FontWeight->"Bold"], " menu or somewhere else). Select", StyleBox[" GraphicsOptions", FontWeight->"Bold"], " and then ", StyleBox["Animation", FontWeight->"Bold"], ". You will find an item called ", StyleBox["AnimationDisplayTime", FontWeight->"Bold"], ". Click on the old value, enter a new one and hit ", StyleBox["Return", FontWeight->"Bold"], ". Increasing it will slow the movie down." }], "Text"], Cell[TextData[{ "We will now write our own version of the Runge-Kutta method. The only \ difference is that the method for taking one step is more complicated. To \ carry it out we will use a ", StyleBox["Module ", FontWeight->"Bold"], "command. This allows us to use temporary variable to hold the four \ necessary slope values." }], "Text"], Cell["\<\ RKStep[f_,{t_,y_},h_] := Module[{m1,m2,m3,m4}, \t\t\tm1 = f[t,y]; \t\t\tm2 = f[t+h/2,y+h/2*m1]; \t\t\tm3 = f[t+h/2,y+h/2*m2]; \t\t\tm4 = f[t+h,y+h*m3]; \t\t\t{t+h, y+(m1 + 2*m2 + 2*m3 + m4)*h/6.0 } \t\t\t] \t\t\t RKList[f_,{t_,y_},h_,n_] := \tNestList[ RKStep[f,#,h]&,{t,y},n]\ \>", "Input"], Cell[TextData[{ "This can be used in the same way as ", StyleBox["EulerList", FontWeight->"Bold"], " to produce an approximate solution, The difference is that it does a \ better job for a given step size, h. Here is a comparison of the two methods \ for the ODE we have been looking at with Euler in red and Runge-Kutta in \ green. After you see the plot, try increasing n and watching them converge \ (the formula for h is adjusted automatically when you change n)." }], "Text"], Cell["\<\ n = 4; h = 8.0/n; elist = EulerList[f,{t0,y0},h,n]; rklist = RKList[f,{t0,y0},h,n]; Show[ Graphics[{RGBColor[1,0,0],Line[elist], \t\tRGBColor[0,1,0],Line[rklist]}],Axes->True ];\ \>", "Input"], Cell[TextData[{ "Try increasing the number, ", StyleBox["n,", FontWeight->"Bold"], " of steps in the last cell from 4 to higher values. For large ", StyleBox["n ", FontWeight->"Bold"], StyleBox["the two methods both give a good approximations to the actual \ solution.", FontVariations->{"CompatibilityType"->0}] }], "Text"], Cell[TextData[{ StyleBox[ "Finally, we will compare the errors made by the two methods in \ approximating ", Evaluatable->False, AspectRatioFixed->True], StyleBox["e ", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[ "by solving the differential equation y ' = y with y[0] = 1 for the time \ interval from 0 to 1. We will plot the differences between the values of \ y[1] produced by the two methods and the actual value of ", Evaluatable->False, AspectRatioFixed->True], StyleBox["e ", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[ "for various values of h. To read off the value of y[1] we apply the ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Last", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[" method to get the pair {1 , y [1] } and then apply ", Evaluatable->False, AspectRatioFixed->True], StyleBox["Last", Evaluatable->False, AspectRatioFixed->True, FontWeight->"Bold"], StyleBox[" again.", Evaluatable->False, AspectRatioFixed->True] }], "Text", Evaluatable->False, AspectRatioFixed->True], Cell["\<\ e = N[Exp[1]]; f[t_,y_] := y errortable = Table[ \t\t{1.0/2^n,\t \t\te - Last[Last[ EulerList[f,{0,1},1.0/2^n,2^n]]], \t\te - Last[Last[ RKList[f,{0,1},1.0/2^n,2^n]]]}, \t\t{n,0,10}]; PrependTo[errortable,{\"Stepsize\",\"Euler\",\"Runge-Kutta\"}]; TableForm[errortable]\t\ \>", "Input", AspectRatioFixed->True], Cell["\<\ Recall that the three methods have errors of orders h and h^4 \ respectively. This is reflected in the table where dividing h by 2 causes \ the errors to be divided by factors near 2 and 16 respectively. To explore \ this we will make a table of the ratios of the errors as the step size is \ repeatedly cut in half. We need to extract the second and third columns of \ the table above and make new list of ratios.\ \>", "Text", Evaluatable->False, AspectRatioFixed->True], Cell[BoxData[ \(\(\(\[IndentingNewLine]\)\(\(errortable\ = \ Drop[errortable, 1];\)\[IndentingNewLine] \(eulererrors\ = \ \(Transpose[ errortable]\)[\([2]\)];\)\[IndentingNewLine] \(rkerrors\ = \ \(Transpose[errortable]\)[\([3]\)];\)\[IndentingNewLine] eulererrorratios\ = \ Table[eulererrors[\([i]\)]/eulererrors[\([i + 1]\)], {i, 1, Length[eulererrors] - 1}]\[IndentingNewLine] rkerrorratios\ = \ Table[rkerrors[\([i]\)]/rkerrors[\([i + 1]\)], {i, 1, Length[rkerrors] - 1}]\)\)\)], "Input"], Cell["\<\ As you can see, the error for the Euler method changes by a factor \ of about 2 when the stepsize does whereas the Runge-Kutta errors change by a \ factor of 16.\ \>", "Text"], Cell[TextData[{ "We can also use the Euler and Runge-Kutta methods for first-order systems \ of differential equations. The only change is that the function ", StyleBox["f[t , y]", FontWeight->"Bold"], " needs to be replaced by a vector valued function of a vector argument ", StyleBox[" f[t_, {y_,v_}]", FontWeight->"Bold"], " as shown below. Here is the harmonic oscillator ", StyleBox["y'' + 2 y = 0", FontWeight->"Bold"], StyleBox[" as a first order system. Note that the initial condition is \ also a vector.", FontVariations->{"CompatibilityType"->0}] }], "Text"], Cell[BoxData[{ \(f[t_, {y_, v_}]\ := \ {v, \ \(-2\)*y}\), "\[IndentingNewLine]", \(elist = EulerList[f, {0, {1, 0}}, 0.1, \ 100]\)}], "Input"], Cell["\<\ To plot the phase portrait we need to extract the list of points \ without the time values. The tricky way to do this is to take the second row \ of the transposed array.\ \>", "Text"], Cell[BoxData[ \(ListPlot[\(Transpose[elist]\)[\([2]\)], PlotJoined \[Rule] True, Axes \[Rule] True, AspectRatio \[Rule] 1, PlotRange \[Rule] {{\(-3\), 3}, {\(-3\), 3}}]\)], "Input"], Cell["\<\ The real solution moves on an ellipse, Here we get a curve which \ appears to spiral out. Try increasing the number of steps in the cell below \ to get a more accurate result. You will find that it is rather hard to get \ something which looks right.\ \>", "Text"], Cell[BoxData[{ \(\(n = 100;\)\), "\[IndentingNewLine]", \(\(h\ = \ 10.0/n;\)\), "\[IndentingNewLine]", \(\(elist = EulerList[f, {0, {1, 0}}, h, \ n];\)\), "\[IndentingNewLine]", \(ListPlot[\(Transpose[elist]\)[\([2]\)], PlotJoined \[Rule] True, Axes \[Rule] True, AspectRatio \[Rule] 1, PlotRange \[Rule] {{\(-3\), 3}, {\(-3\), 3}}]\)}], "Input"], Cell["\<\ Here is a similar test of the Runge-Kutta method for systems. Try \ gradually increasing n until you get a reasonable picture.\ \>", "Text"], Cell[BoxData[{ \(\(n = 10;\)\), "\[IndentingNewLine]", \(\(h\ = \ 10.0/n;\)\), "\[IndentingNewLine]", \(\(rklist = RKList[f, {0, {1, 0}}, h, \ n];\)\), "\[IndentingNewLine]", \(ListPlot[\(Transpose[rklist]\)[\([2]\)], PlotJoined \[Rule] True, Axes \[Rule] True, AspectRatio \[Rule] 1, PlotRange \[Rule] {{\(-3\), 3}, {\(-3\), 3}}]\)}], "Input"], Cell["\<\ If you get to this point you could try replacing the harmonic \ oscillator by some other second-order ODE. Or, for a challenge, try making a \ movie showing the convergence of the Euler method phase plots as n gets \ larger.\ \>", "Text"] }, Closed]] }, FrontEndVersion->"4.1 for Macintosh", ScreenRectangle->{{0, 1024}, {0, 746}}, WindowToolbars->{}, CellGrouping->Manual, WindowSize->{699, 530}, WindowMargins->{{37, Automatic}, {Automatic, 48}}, PrintingCopies->1, PrintingPageRange->{1, Automatic}, PrintingOptions->{"PrintingMargins"->{{54, 54}, {72, 72}}, "PrintCellBrackets"->True, "PrintRegistrationMarks"->False, "PrintMultipleHorizontalPages"->False}, PrivateNotebookOptions->{"ColorPalette"->{RGBColor, -1}}, ShowCellLabel->True, ShowCellTags->False, RenderingOptions->{"ObjectDithering"->True, "RasterDithering"->False} ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[1705, 50, 46, 0, 108, "Title"], Cell[1754, 52, 907, 26, 86, "Text"], Cell[2664, 80, 77, 1, 27, "Input"], Cell[2744, 83, 617, 10, 86, "Text"], Cell[CellGroupData[{ Cell[3386, 97, 54, 0, 56, "Section"], Cell[CellGroupData[{ Cell[3465, 101, 824, 26, 68, "Text", Evaluatable->False], Cell[4292, 129, 98, 1, 27, "Input"], Cell[4393, 132, 883, 23, 86, "Text"], Cell[5279, 157, 27, 0, 27, "Input"], Cell[5309, 159, 1083, 37, 50, "Text", Evaluatable->False], Cell[6395, 198, 66, 1, 27, "Input"], Cell[6464, 201, 774, 23, 68, "Text", Evaluatable->False], Cell[7241, 226, 98, 4, 42, "Input"], Cell[7342, 232, 902, 29, 68, "Text", Evaluatable->False], Cell[8247, 263, 150, 5, 57, "Input"], Cell[8400, 270, 452, 12, 68, "Text", Evaluatable->False] }, Open ]], Cell[8867, 285, 398, 10, 68, "Text"], Cell[9268, 297, 222, 5, 50, "Text"], Cell[CellGroupData[{ Cell[9515, 306, 241, 5, 32, "Text", Evaluatable->False], Cell[9759, 313, 93, 1, 27, "Input"], Cell[9855, 316, 426, 12, 50, "Text"], Cell[10284, 330, 27, 0, 27, "Input"], Cell[10314, 332, 1083, 37, 50, "Text", Evaluatable->False], Cell[11400, 371, 66, 1, 27, "Input"], Cell[11469, 374, 1071, 32, 86, "Text", Evaluatable->False], Cell[12543, 408, 209, 4, 43, "Input"], Cell[12755, 414, 451, 11, 68, "Text"], Cell[13209, 427, 71, 1, 27, "Input"], Cell[13283, 430, 31, 0, 32, "Text"], Cell[13317, 432, 58, 1, 27, "Input"], Cell[13378, 435, 111, 3, 32, "Text"], Cell[13492, 440, 286, 5, 59, "Input"], Cell[13781, 447, 718, 15, 104, "Text"], Cell[14502, 464, 105, 2, 27, "Input"], Cell[14610, 468, 165, 4, 32, "Text"], Cell[14778, 474, 191, 3, 59, "Input"] }, Open ]] }, Closed]], Cell[CellGroupData[{ Cell[15018, 483, 61, 0, 36, "Section"], Cell[15082, 485, 381, 10, 50, "Text"], Cell[15466, 497, 25, 0, 27, "Input"], Cell[15494, 499, 140, 4, 32, "Text"], Cell[15637, 505, 119, 5, 72, "Input"], Cell[15759, 512, 390, 8, 68, "Text"], Cell[16152, 522, 88, 4, 57, "Input"], Cell[16243, 528, 803, 21, 104, "Text"], Cell[17049, 551, 170, 7, 102, "Input"], Cell[17222, 560, 58, 0, 32, "Text"], Cell[17283, 562, 52, 0, 27, "Input"], Cell[17338, 564, 281, 6, 50, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[17656, 575, 38, 0, 36, "Section"], Cell[17697, 577, 248, 5, 50, "Text"], Cell[17948, 584, 184, 3, 43, "Input"], Cell[18135, 589, 98, 3, 32, "Text"], Cell[18236, 594, 71, 1, 27, "Input"], Cell[18310, 597, 256, 6, 50, "Text"], Cell[18569, 605, 85, 1, 27, "Input"], Cell[18657, 608, 389, 7, 68, "Text"], Cell[19049, 617, 232, 5, 59, "Input"], Cell[19284, 624, 134, 3, 32, "Text"], Cell[19421, 629, 466, 8, 123, "Input"], Cell[19890, 639, 541, 10, 86, "Text"], Cell[20434, 651, 113, 2, 27, "Input"], Cell[20550, 655, 297, 6, 50, "Text"], Cell[20850, 663, 380, 6, 107, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[21267, 674, 52, 0, 36, "Section"], Cell[21322, 676, 1133, 33, 86, "Text", Evaluatable->False], Cell[22458, 711, 87, 1, 27, "Input"], Cell[22548, 714, 899, 20, 104, "Text"], Cell[23450, 736, 121, 4, 42, "Input"], Cell[23574, 742, 338, 6, 68, "Text"], Cell[23915, 750, 87, 5, 72, "Input"], Cell[24005, 757, 718, 18, 122, "Text"], Cell[24726, 777, 77, 4, 57, "Input"], Cell[24806, 783, 596, 14, 86, "Text"], Cell[25405, 799, 168, 8, 117, "Input"], Cell[25576, 809, 163, 4, 32, "Text"], Cell[25742, 815, 407, 10, 68, "Text"], Cell[26152, 827, 206, 4, 75, "Input"], Cell[26361, 833, 1125, 30, 140, "Text"], Cell[27489, 865, 354, 8, 50, "Text"], Cell[27846, 875, 302, 11, 162, "Input"], Cell[28151, 888, 493, 9, 68, "Text"], Cell[28647, 899, 204, 9, 132, "Input"], Cell[28854, 910, 347, 10, 50, "Text"], Cell[29204, 922, 1219, 40, 68, "Text", Evaluatable->False], Cell[30426, 964, 326, 15, 207, "Input"], Cell[30755, 981, 491, 9, 68, "Text", Evaluatable->False], Cell[31249, 992, 574, 11, 139, "Input"], Cell[31826, 1005, 185, 4, 50, "Text"], Cell[32014, 1011, 602, 14, 68, "Text"], Cell[32619, 1027, 152, 2, 43, "Input"], Cell[32774, 1031, 195, 4, 50, "Text"], Cell[32972, 1037, 199, 3, 43, "Input"], Cell[33174, 1042, 277, 5, 50, "Text"], Cell[33454, 1049, 389, 7, 91, "Input"], Cell[33846, 1058, 151, 3, 32, "Text"], Cell[34000, 1063, 378, 6, 91, "Input"], Cell[34381, 1071, 249, 5, 50, "Text"] }, Closed]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************)