(***********************************************************************

                    Mathematica-Compatible Notebook

This notebook can be used on any computer system with Mathematica 3.0,
MathReader 3.0, or any compatible application. The data for the notebook 
starts with the line of 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[     28675,        932]*)
(*NotebookOutlinePosition[     29941,        972]*)
(*  CellTagsIndexPosition[     29897,        968]*)
(*WindowFrame->Normal*)



Notebook[{
Cell["Differential Equations Lab I", "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 Explicit 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 Explicit 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"]
}, 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 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["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 yourself.  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 it's outputs back into it's 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.1;
n = 10; 
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"],
  ".  The 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 in the ",
  StyleBox["Edit",
    FontWeight->"Bold"],
  " menu. 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[{
  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.\
\>", "Text",\

  Evaluatable->False,
  AspectRatioFixed->True]
}, Closed]],

Cell[CellGroupData[{

Cell["Picard Iteration", "Section"],

Cell[TextData[{
  StyleBox[
  "Next we will set up a function which does Picard iteration and create a \
movie illustrating the convergence of the successive approximations.  Except \
for the ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["Evaluate ",
    Evaluatable->False,
    AspectRatioFixed->True,
    FontWeight->"Bold"],
  StyleBox["command, the ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["Picard ",
    Evaluatable->False,
    AspectRatioFixed->True,
    FontWeight->"Bold"],
  StyleBox["function below is just a translation into ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["Mathematica",
    Evaluatable->False,
    AspectRatioFixed->True,
    FontSlant->"Italic"],
  StyleBox[" of Picard's formula.  You feed it the current approximation ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["yn", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[", the ODE ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["f[t,y]", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[", and the initial conditions ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["t0,y0", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[".  It returns the next approximation, ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["yn+1", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[
  ".  This is an example of a function whose return value is a function, a \
possibility disallowed in most programming languages.",
    Evaluatable->False,
    AspectRatioFixed->True]
}], "Text"],

Cell[BoxData[
    \(Picard[yn_, f_, t0_, y0_]\  := \ 
      Function[t, \n\t\t\t\t
        Evaluate[y0\  + \ Integrate[f[u, yn[u]], {u, t0, t}]]]\)], "Input"],

Cell[TextData[{
  "Now we'll try it out on the ODE  y' = y using the initial conditions x0=0 \
and y0 = 1.  We will call the successive approximations ",
  StyleBox["pic[0], pic[1]", "Input"],
  " and so on. To evaluate them, you use expressions like ",
  StyleBox["pic[0][t]", "Input"],
  "."
}], "Text"],

Cell[BoxData[{
    \(t0\  = \ 0; \ny0\  = \ 1; \n\(pic[0]\)[x_]\  := \ y0\ ; \n
    f[t_, y_]\  := \ y; \npic[1]\  = \ Picard[pic[0], f, t0, y0]; \n
    pic[2]\  = \ Picard[pic[1], f, t0, y0]; \n
    pic[3]\  = \ Picard[pic[2], f, t0, y0]; \n\(pic[1]\)[t]\), 
    \(\(pic[2]\)[t]\), 
    \(\(pic[3]\)[t]\)}], "Input"],

Cell["\<\
Now we will generate several graphs of these approximations and \
animate them to make the movie.\
\>", "Text"],

Cell[BoxData[{
    \(Do[\ pic[i]\  = \ Picard[pic[i - 1], f, 0, 1]; \n\t
      Plot[\(pic[i]\)[t], {t, \(-5\), 5}, PlotRange -> {\(-5\), 5}], \n
      \t{i, 1, 12}]\), 
    \(Print["\<Done.\>"]\)}], "Input"],

Cell[TextData[{
  StyleBox[
  "Now try making a movie of the successive approximations to the solution of \
another ODE.  For example,   ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["f[t_,y_] := -Sin[t]*y", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[
  "  makes an amusing show.  For this one, it is better to change the ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["PlotRange", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[" inside the ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["Do",
    Evaluatable->False,
    AspectRatioFixed->True,
    FontWeight->"Bold"],
  StyleBox[" loop to ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["{-1,1}", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[" and the range on ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["t", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[" to ",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox["{t,-10,10}", "Input",
    Evaluatable->False,
    AspectRatioFixed->True],
  StyleBox[".",
    Evaluatable->False,
    AspectRatioFixed->True]
}], "Text"]
}, Open  ]],

Cell[CellGroupData[{

Cell["Blowup in Finite Time", "Section"],

Cell[TextData[{
  "A good test case for numerical methods is the logistic equation   y ' = \
y(1-y).  However, it also gives a few problems due to the fact that some \
solutions don't exist for all time.  For example, here is an attempt to \
compute the solution with initial condition ",
  StyleBox["y[0]==2", "Input"],
  " using NDSolve.\n"
}], "Text"],

Cell[BoxData[{
    \(rules\  = \ 
      NDSolve[\n\t\t\t{\(y'\)[t]\  == \ y[t]*\((1 - y[t])\), y[0] == 2}, \n
        \t\t\ty[t], {t, \(-5\), 5}\n\t\t\t]\n\), 
    \(ypart\  = \ y[t] /. rules[\([1]\)]; \n\n
    Plot[ypart, {t, \(-5\), 5}, PlotRange -> All]; \)}], "Input"],

Cell["\<\
You will get some scary looking output due to the fac t that the \
solution does not exist over the given time interval.  In particular it blows \
up for some negative time  between -5 and 0.  Go back to the previous cell \
and change  the -5 to something less ambitious.  Try to estimate the maximal \
interval of existence for this solution.\
\>", "Text"],

Cell[TextData[{
  "It would be a problem to plot many integral curves for this ODE because \
each choice of initial condition could lead to a different interval of \
existence.  We are going to avoid this problem by introducing a trick -- \
changing the differential equation outside a certain range of ",
  StyleBox["y", "Input"],
  " values. Suppose we agree that solutions outside the range ",
  StyleBox["-10<y<10", "Input"],
  " are no longer interesting.  Then it should do no harm to change the \
equation out there to slow the solutions down.  In fact we will make the new \
ODE zero outside this range."
}], "Text"],

Cell[BoxData[
    \(f[t_, y_]\  := \ If[Abs[y] <= 10, y*\((1 - y)\), 0]\)], "Input"],

Cell[TextData[{
  "In other words, if ",
  StyleBox["|y|<10", "Input"],
  " we use the correct formula ",
  StyleBox["y*(1-y)", "Input"],
  ".  Otherwise we use the formula ",
  StyleBox["0", "Input"],
  ".  Now we can try again."
}], "Text"],

Cell[BoxData[{
    \(rules\  = \ 
      NDSolve[\n\t\t\t{\(y'\)[t]\  == \ f[t, y[t]], y[0] == 2}, \n\t\t\t
        y[t], {t, \(-5\), 5}\n\t\t\t]\n\), 
    \(ypart\  = \ y[t] /. rules[\([1]\)]; \n\n
    Plot[ypart, {t, \(-5\), 5}, PlotRange -> All]; \)}], "Input"],

Cell["\<\
You can see how the solution which  previously exploded to infinity \
is now stopped in its tracks at y=10.  Using the modified ODE we can get a \
nice picture of several integral curves as we did before for less problematic \
ODE's.\
\>", "Text"],

Cell[BoxData[
    \(sollist\  = \ {}; \n
    Do[rules\ \  = \ 
        NDSolve[{\(y'\)[t] == f[t, y[t]], y[0] == i/4}, y, {t, \(-5\), 5}]; \n
      AppendTo[sollist, \ y[t] /. rules[\([1]\)]], \n{i, \(-6\), 6}]\)], 
  "Input"],

Cell["\<\
If you restrict the plot range when you display the results, no one \
will ever know !\
\>", "Text"],

Cell[BoxData[
    \(\(Plot[\ Evaluate[sollist], \ {t, \(-5\), 5}, 
      PlotRange -> {\(-3\), 3}]; \)\)], "Input"]
}, Closed]]
},
FrontEndVersion->"Macintosh 3.0",
ScreenRectangle->{{0, 832}, {0, 604}},
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{544, 515},
WindowMargins->{{56, Automatic}, {Automatic, 15}},
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},
MacintoshSystemPageSetup->"\<\
00<0004/0B`000002mT8o?mooh<<H`Ui0fL5:0?l0080004/0B`000000]P2:001
0000I00000400`<300000BL?00400@<?kE0000000000060801T1T00000000000
00000000000000000000000000000000\>"
]


(***********************************************************************
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[1709, 49, 45, 0, 98, "Title"],
Cell[1757, 51, 907, 26, 98, "Text"],
Cell[2667, 79, 77, 1, 27, "Input"],
Cell[2747, 82, 620, 10, 112, "Text"],

Cell[CellGroupData[{
Cell[3392, 96, 57, 0, 50, "Section"],

Cell[CellGroupData[{
Cell[3474, 100, 830, 28, 65, "Text",
  Evaluatable->False],
Cell[4307, 130, 98, 1, 27, "Input"],
Cell[4408, 133, 885, 23, 81, "Text"],
Cell[5296, 158, 27, 0, 27, "Input"],
Cell[5326, 160, 1084, 38, 64, "Text",
  Evaluatable->False],
Cell[6413, 200, 66, 1, 27, "Input"],
Cell[6482, 203, 777, 24, 66, "Text",
  Evaluatable->False],
Cell[7262, 229, 98, 4, 42, "Input"],
Cell[7363, 235, 908, 31, 66, "Text",
  Evaluatable->False],
Cell[8274, 268, 150, 5, 57, "Input"],
Cell[8427, 275, 453, 12, 65, "Text",
  Evaluatable->False]
}, Open  ]],
Cell[8895, 290, 398, 10, 65, "Text"],
Cell[9296, 302, 222, 5, 48, "Text"]
}, Closed]],

Cell[CellGroupData[{
Cell[9555, 312, 61, 0, 30, "Section"],
Cell[9619, 314, 377, 10, 64, "Text"],
Cell[9999, 326, 25, 0, 27, "Input"],
Cell[10027, 328, 140, 4, 31, "Text"],
Cell[10170, 334, 119, 5, 72, "Input"],
Cell[10292, 341, 390, 8, 64, "Text"],
Cell[10685, 351, 88, 4, 57, "Input"],
Cell[10776, 357, 803, 21, 99, "Text"],
Cell[11582, 380, 170, 7, 117, "Input"],
Cell[11755, 389, 58, 0, 30, "Text"],
Cell[11816, 391, 52, 0, 27, "Input"],
Cell[11871, 393, 282, 6, 47, "Text"]
}, Closed]],

Cell[CellGroupData[{
Cell[12190, 404, 52, 0, 30, "Section"],
Cell[12245, 406, 1145, 34, 80, "Text",
  Evaluatable->False],
Cell[13393, 442, 87, 1, 27, "Input"],
Cell[13483, 445, 901, 20, 114, "Text"],
Cell[14387, 467, 121, 4, 42, "Input"],
Cell[14511, 473, 338, 6, 62, "Text"],
Cell[14852, 481, 87, 5, 72, "Input"],
Cell[14942, 488, 714, 18, 98, "Text"],
Cell[15659, 508, 77, 4, 57, "Input"],
Cell[15739, 514, 596, 14, 82, "Text"],
Cell[16338, 530, 168, 8, 117, "Input"],
Cell[16509, 540, 163, 4, 46, "Text"],
Cell[16675, 546, 407, 10, 66, "Text"],
Cell[17085, 558, 204, 4, 75, "Input"],
Cell[17292, 564, 1006, 27, 134, "Text"],
Cell[18301, 593, 354, 8, 64, "Text"],
Cell[18658, 603, 302, 11, 162, "Input"],
Cell[18963, 616, 493, 9, 80, "Text"],
Cell[19459, 627, 204, 9, 132, "Input"],
Cell[19666, 638, 1219, 40, 84, "Text",
  Evaluatable->False],
Cell[20888, 680, 326, 15, 207, "Input"],
Cell[21217, 697, 280, 7, 46, "Text",
  Evaluatable->False]
}, Closed]],

Cell[CellGroupData[{
Cell[21534, 709, 35, 0, 30, "Section"],
Cell[21572, 711, 1686, 54, 116, "Text"],
Cell[23261, 767, 158, 3, 43, "Input"],
Cell[23422, 772, 305, 7, 64, "Text"],
Cell[23730, 781, 317, 6, 171, "Input"],
Cell[24050, 789, 121, 3, 30, "Text"],
Cell[24174, 794, 207, 4, 75, "Input"],
Cell[24384, 800, 1251, 44, 65, "Text"]
}, Open  ]],

Cell[CellGroupData[{
Cell[25672, 849, 40, 0, 50, "Section"],
Cell[25715, 851, 354, 7, 79, "Text"],
Cell[26072, 860, 272, 5, 139, "Input"],
Cell[26347, 867, 367, 6, 62, "Text"],
Cell[26717, 875, 624, 11, 95, "Text"],
Cell[27344, 888, 84, 1, 27, "Input"],
Cell[27431, 891, 242, 8, 31, "Text"],
Cell[27676, 901, 263, 5, 139, "Input"],
Cell[27942, 908, 257, 5, 46, "Text"],
Cell[28202, 915, 226, 5, 75, "Input"],
Cell[28431, 922, 110, 3, 30, "Text"],
Cell[28544, 927, 115, 2, 27, "Input"]
}, Closed]]
}
]
*)




(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)