(***********************************************************************
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[ 32950, 968]*)
(*NotebookOutlinePosition[ 34216, 1008]*)
(* CellTagsIndexPosition[ 34172, 1004]*)
(*WindowFrame->Normal*)
Notebook[{
Cell["Differential Equations Lab II", "Title"],
Cell[TextData[{
"This lab ",
StyleBox["explores more 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, shown in \
bolder type. To execute a command, click the mouse in the cell containing \
the command then do a Shift-Enter (or Shift-Return), that is, hold down the \
shift key and hit the Enter key (or Return key).",
Evaluatable->False,
AspectRatioFixed->True]
}], "Text"],
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["Explicit Solutions for Second Order Equations.",
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["Explicit Solutions for Second Order Equations", "Section"],
Cell[CellGroupData[{
Cell[TextData[{
StyleBox[
"To find the exact formula for the solution of simple scalar 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. 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]\),
\(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]\),
\(y2\ = Simplify[\ ComplexExpand[y[t] /. rules[\([1]\)]\ ]]\),
\(\(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]\)]\),
\(ypart = ygen /. {C[1] -> 1, C[2] -> 0}\),
\(\(Plot[{Re[ypart], Im[ypart]}, {t, \(-5\), 5}]; \)\)}], "Input"]
}, Open ]],
Cell["And how about some Legendre functions ?", "Text"],
Cell[BoxData[
\(n = 2; \n
rules = DSolve[\
\((1 - t^2)\)*\(\(y'\)'\)[t]\ - 2*t*\(y'\)[t]\ + \
n*\((n + 1)\) y[t] == 0, \ y[t], \ t]\)], "Input"],
Cell["\<\
This time both solutions are real. The first solution is a \
polynomial and the second is a special function.\
\>", "Text"],
Cell[BoxData[{
\(ygen\ = \ \ y[t] /. rules[\([1]\)]\),
\(ypart1 = ygen /. {C[1] -> 1, C[2] -> 0}\),
\(ypart2 = ygen /. {C[1] -> 0, C[2] -> 1}\),
\(\(Plot[{ypart1, ypart2}, {t, \(-1\), 1}]; \)\)}], "Input"],
Cell["\<\
Try going back and changing the value of n for the Legendre \
equation to some other integer. If you try non-integers, there won't be a \
polynomial solution.\
\>", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Eigenvalues, Eigenvectors -- Coupled Linear Oscillators", "Section"],
Cell[TextData[{
StyleBox["Mathematica",
FontSlant->"Italic"],
" has lots of built-in linear algebra functions which make it easy to work \
with linear systems. A matrix is represented as a list of lists. Each list \
represents one row of the matrix. For example"
}], "Text"],
Cell[BoxData[{
\(A\ = \ {{1, 1}, {1, 0}}\),
\(MatrixForm[A]\)}], "Input"],
Cell[TextData[{
"As you can see, the function ",
StyleBox["MatrixForm", "Input"],
"will print the matrix in the familiar way. Once you have entered the \
matrix you can find its determinant, characteristic polynomial, eigenvalues \
and eigenvectors."
}], "Text"],
Cell[BoxData[{
\(Det[A]\),
\(CharacteristicPolynomial[A, z]\),
\(Eigenvalues[A]\),
\(Eigenvectors[A]\)}], "Input"],
Cell["\<\
Here is a more complicated example in four-dimensions. Two coupled \
oscillators are modelled by coupled second-order equations: x1'' + omega1^2 \
x1 = eps(x2-x1) and x2'' + omega2^2 x2= eps(x1-x2). Here x1 and x2 represent \
the positions of the oscillators. If we introduce the velocities v1 and v2 \
as new variables then x1' = v1 and x2' = v2. The second derivatives of the \
positions become the first derivatives of the velocities. In the end we have \
a four-dimensional first-order linear system with the matrix below.\
\>",
"Text"],
Cell[BoxData[
\(A\ =
\ {{0, 1, 0, 0}, \n\t\t{\(-omega1^2\) - eps, 0, eps, 0}, \n
\t\t{0, 0, 0, 1}, \n\t\t{eps, 0, \(-omega2^2\) - eps, 0}}; \n
MatrixForm[A]\)], "Input"],
Cell["\<\
You can compute the eigenvalues even with the variable entries in \
the matrix.\
\>", "Text"],
Cell[BoxData[{
\(CharacteristicPolynomial[A, z]\),
\(Eigenvalues[A] // Simplify\)}], "Input"],
Cell["\<\
This is fairly complicated so we will pick some particular values \
for the frequencies of the two oscillators and for the coupling constant. \
First we will take the two frequencies equal.\
\>", "Text"],
Cell[BoxData[{
\(omega1 = 1; \nomega2 = 1; \neps\ = \ 0.1; \nMatrixForm[A]\),
\(lambda\ = \ Eigenvalues[A]\),
\(V\ = \ Eigenvectors[A]\)}], "Input"],
Cell["\<\
All of the eigenvalues are imaginary numbers. Each pair of \
imaginary eigenvalues determines a simple periodic motion of the coupled \
oscillator system. These are the so-called normal modes of oscillation. You \
can find these periodic solutions as follows. First generate the complex \
solution corresponding to one of the complex eigenvectors, then extract its \
real part to get a real solution representing a normal mode of \
oscillation.\
\>", "Text"],
Cell[BoxData[{
\(Z[t]\ = \ Exp[lambda[\([1]\)]*t]*V[\([1]\)] // ComplexExpand\),
\(X[t]\ = \ Re[Z[t]] // ComplexExpand\)}], "Input"],
Cell["\<\
Now the four variables in our solution vector are \
{x1[t],v1[t],x2[t],v2[t]}. We will plot the first and third components which \
represent the positions of the oscillators so we can visualize the normal \
mode. The first oscillator is shown in red and the second in green.\
\>",
"Text"],
Cell[BoxData[{
\(x1\ = \ \(X[t]\)[\([1]\)]\),
\(x2 = \ \ \(X[t]\)[\([3]\)]\),
\(\(Plot[{x1, x2}, {t, 0, 10},
PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \)\)}], "Input"],
Cell["\<\
As you can see, for this normal mode of oscillation, the two \
oscillators are exactly out of phase with one another. When one moves to the \
left, the other moves to the right. We can find the other normal mode in the \
same way using the other eigenvalue and corresponding eigenvector.\
\>",
"Text"],
Cell[CellGroupData[{
Cell[BoxData[{
\(Z[t]\ = \ Exp[lambda[\([3]\)]*t]*V[\([3]\)] // ComplexExpand\),
\(X[t]\ = \ Re[Z[t]] // ComplexExpand\),
\(x1\ = \ \(X[t]\)[\([1]\)]\),
\(x2 = \ \ \(X[t]\)[\([3]\)]\),
\(\(Plot[{x1, x2}, {t, 0, 10},
PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \)\)}], "Input"],
Cell["\<\
For this normal mode the two oscillators are exactly in synch with \
one another. These normal modes are the ones for the case when the two \
frequencies omega1 and omega2 are equal. Now we are in a position to look at \
other case, too complicated to be worked out by hand. \
\>", "Text"],
Cell[TextData[{
"In the cell below you can choose different values for the frequencies of \
the separate oscillators and then look at the two normal modes. For example \
omega1 = 1 and omega2 = 2, the normal modes are entirely different than the \
ones above. In one mode, the first oscillator moves with large amplitude \
while the second is almost motionless. In the second mode, the situation is \
reveresed. Try several other choices of omega1 and omega2. For example, you \
could keep ",
StyleBox["omega1=1 ", "Input"],
"and take a sequence of values for ",
StyleBox["omega2 ", "Input"],
"decreasing from 2 back down to 1."
}], "Text"],
Cell[BoxData[{
\(omega1 = 1; \nomega2 = \ 2; \neps\ = \ 0.1; \nMatrixForm[A]\),
\(lambda\ = \ Eigenvalues[A]\),
\(V\ = \ Eigenvectors[A]; \n
X1[t]\ = Re[\ Exp[lambda[\([1]\)]*t]*V[\([1]\)]] // ComplexExpand; \n
X2[t]\ = Re[\ Exp[lambda[\([3]\)]*t]*V[\([3]\)]] // ComplexExpand; \n
x1\ = \ \(X1[t]\)[\([1]\)]; x2 = \ \ \(X1[t]\)[\([3]\)]; \n
Plot[{x1, x2}, {t, 0, 10},
PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \n
x1\ = \ \(X2[t]\)[\([1]\)]; x2 = \ \ \(X2[t]\)[\([3]\)]; \n
Plot[{x1, x2}, {t, 0, 10},
PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \)}], "Input"],
Cell["\<\
The normal modes are just the simplest kinds of solutions for the \
coupled oscillators. A typical solution looks like a linear combination of \
the two normal modes and is quasiperiodic rather than periodic. To plot one \
of these solutions you just have to take a combination of the two complex \
solutions.\
\>", "Text"],
Cell[BoxData[{
\(omega1 = 1; \nomega2 = 1; \neps\ = \ 0.1; \nMatrixForm[A]\),
\(lambda\ = \ Eigenvalues[A]\),
\(V\ = \ Eigenvectors[A]; \n
Z[t]\ = \
Exp[lambda[\([1]\)]*t]*V[\([1]\)]\ +
Exp[lambda[\([3]\)]*t]*V[\([3]\)]; \n
X[t]\ = Re[Z[t]] // ComplexExpand; \nx1\ = \ \(X[t]\)[\([1]\)];
x2 = \ \ \(X[t]\)[\([3]\)]; \n
Plot[{x1, x2}, {t, 0, 60},
PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \)}], "Input"],
Cell["\<\
For the case of equal frequencies shown above, the resulting motion \
is very interesting. First one oscillator has large amplitude while the \
other is relatively quiet, then later the situation is reversed. \
\>",
"Text"]
}, Open ]]
}, Closed]],
Cell[CellGroupData[{
Cell["Linear Flows -- Mapping the Cat", "Section"],
Cell["\<\
For a linear system X' = AX we have seen that the solution is a \
linear function of its ititial condition, namely,
X(t) = exp(At) X(0) where exp(At) is the matrix exponential. Thinking of all \
initial vectors X(0) moving simultaneously in this way we have the concept of \
the flow of the system. In this section we will compute exp(tA) for some \
simple linear systems in the plane and use it to animate the flow.\
\>",
"Text"],
Cell[TextData[{
"There is a built-in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" function called ",
StyleBox["MatrixExp ", "Input"],
"which can actually compute the flow matrix in simple cases.\nFor \
example,"
}], "Text"],
Cell[BoxData[{
\(A\ = \ {{1, 1}, {1, 0}}; \nMatrixForm[A]\),
\(etA\ = \ MatrixExp[t*A] // Simplify; \nMatrixForm[etA]\)}], "Input"],
Cell[TextData[{
"Next we will use the flow matrix to transform a simple shape. The \
simplest shapes are curves represented by ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" Line objects. A Line is just a list of points to be connected with line \
segments. For example, a square could be represented"
}], "Text"],
Cell[BoxData[{
\(square\ = \ \ Line[{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 0}}]\),
\(\(Show[Graphics[square], Axes -> True, AspectRatio -> 1]; \)\)}],
"Input"],
Cell["\<\
To transform the square we just need to multiply each of the points \
in the Line by the flow matrix. We will define a function to do this for us.\
\
\>", "Text"],
Cell[BoxData[
\(Transform[l_Line, A_]\ := \ Line[Map[A . #&, l[\([1]\)]]]\)], "Input"],
Cell["\<\
To compute the image of the square under the time 1 map of the \
flow, set\
\>", "Text"],
Cell[BoxData[
\(t = 1; \nsq1\ = \ Transform[square, etA]; \n
Show[Graphics[{square, sq1}], Axes -> True]; \)], "Input"],
Cell["Of course an animation is even better.", "Text"],
Cell[BoxData[
\(\(Do[
Show[Graphics[Transform[square, etA]], Axes -> True,
PlotRange -> {{0, 6}, {0, 5}}], {t, 0, 1, 0.1}]; \)\)], "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 ",
StyleBox["Cell", "Input"],
" menu and choose the item called ",
StyleBox["Animate Selected Graphics", "Input"],
". You should see a movie of your plots showing how the square is \
transformed. 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", "Input"],
" item in the ",
StyleBox["Edit", "Input"],
" menu. ",
StyleBox["Select GraphicsOptions", "Input"],
" and then ",
StyleBox["Animation", "Input"],
". You will find an item called ",
StyleBox["AnimationDisplayTime", "Input"],
". Click on the old value, enter a new one and hit Enter (or Return). \
Increasing it will slow the movie down."
}], "Text"],
Cell["\<\
It would be much more interesting to animate the images of some \
more complicated shape. The shape below consists of several Line objects \
which together form an approximation to a cat face.\
\>", "Text"],
Cell[BoxData[{
\(c1\ = \ Line[Table[{2*Cos[t], 1.5*Sin[t]}, {t, 0, 2*Pi, Pi/20}]]; \n
c2\ = \ Line[
Table[{1 + 0.3*Cos[t], 0.5 + 0.2*Sin[t]}, {t, 0, 2*Pi, Pi/20}]]; \n
c3\ = \ Line[
Table[{\(-1\) + 0.3*Cos[t], 0.5 + 0.2*Sin[t]}, {t, 0, 2*Pi, Pi/20}]];
\nc4\ = \
Line[{{2*Cos[Pi/6], 1.5*Sin[Pi/6]}, {1.5, 2}, {2*Cos[Pi/3],
1.5*Sin[Pi/3]}}]; \n
c5\ = \ Line[{{2*Cos[5*Pi/6], 1.5*Sin[5*Pi/6]}, {\(-1.5\), 2}, {
2*Cos[2*Pi/3], 1.5*Sin[2*Pi/3]}}]; \n
c6\ = \ Line[{{2.2, \(-0.8\)}, {1, \(-0.6\)}, {2.5, \(-0.5\)}}]\),
\(c7\ = \
Line[{{\(-2.2\), \(-0.8\)}, {\(-1\), \(-0.6\)}, {\(-2.5\), \(-0.5\)}}]\n
\),
\(cat\ = \ {c1, c2, c3, c4, c5, c6, c7}; \n\n
Show[Graphics[cat], PlotRange -> {{\(-3\), 3}, {\(-3\), 3}},
Axes -> True, AspectRatio -> 1]; \)}], "Input"],
Cell["\<\
Since the cat is made of a list of lines, we need to define a \
function for tranforming lists of lines by transforming each line in the \
list.\
\>", "Text"],
Cell[BoxData[{
\(Transform[l_List, A_]\ := \ Map[Transform[#, A]&, l]\),
\(t = 1; \ncat1 = \ Transform[cat, etA]; \n
Show[Graphics[cat1], Axes -> True]; \)}], "Input"],
Cell["And finally, CAT -- The Movie.", "Text"],
Cell[BoxData[{
\(t =. ; \nA\ = \ {{1, 1}, {1, 0}}; \nMatrixForm[A]\),
\(etA\ = \ MatrixExp[t*A] // Simplify; \n\n
Do[Show[Graphics[Transform[cat, etA]], Axes -> True, AspectRatio -> 1,
PlotRange -> {{\(-10\), 10}, {\(-10\), 10}}], {t, 0, 1, 0.1}]; \)}],
"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 ",
StyleBox["Cell", "Input"],
" menu and choose the item called ",
StyleBox["Animate Selected Graphics", "Input"],
". You should see a movie of your plots showing how the cat is \
transformed. To stop the cat 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", "Input"],
" item in the ",
StyleBox["Edit", "Input"],
" menu. ",
StyleBox["Select GraphicsOptions", "Input"],
" and then ",
StyleBox["Animation", "Input"],
". You will find an item called ",
StyleBox["AnimationDisplayTime", "Input"],
". Click on the old value, enter a new one and hit Enter (or Return). \
Increasing it will slow the movie down."
}], "Text"],
Cell[TextData[{
"You might want to go back and try mapping the cat with some other matrices \
A which have different kinds of flows. For example ",
StyleBox["A = {{1,-5},{5,1}}", "Input"],
" has a unstable spiral at the origin while ",
StyleBox["A = {{2,0},{0,1}}", "Input"],
" has an unstable node."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[" 2D Phase Portraits", "Section"],
Cell[TextData[{
"For nonlinear systems, it is not generally possible to find explicit \
formulas for the solutions. However, ",
StyleBox["NDSolve ", "Input"],
"can be used to generate an approximate solution corresponding t o a given \
initial condition. You have to specify the ODE and the initial conditions as \
indicated below. This example is a simple linear system in the plane."
}], "Text"],
Cell[BoxData[
\(rules\ = \
NDSolve[{\(x'\)[t] == \(-x[t]\)/2 - y[t], \(y'\)[t] == \ x[t] - y[t]/2,
x[0] == 1, y[0] == 0}, \n\t\t{x[t], y[t]}, {t, 0, 10}]\)], "Input"],
Cell["\<\
To extract the solution from these substitution rules use the \
following \
\>", "Text"],
Cell[BoxData[
\(X\ = \ {x[t], y[t]} /. rules[\([1]\)]\)], "Input"],
Cell[TextData[{
"The answer is in the form of a vector ",
StyleBox["X", "Input"],
" of two functions, representing the x and y 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_, {x0_, y0_}, time_]\ := \n
\t{x[t], y[t]} /.
\(NDSolve[
ode~Join~{x[0] == x0, y[0] == y0}, {x[t], y[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\ =
\ {\(x'\)[t] == \(-x[t]\)/2 - y[t], \(y'\)[t] == \ x[t] - y[t]/2}; \n
X1\ = \ Orbit[ode, {1, 1}, 10]; \nX2\ = \ Orbit[ode, {1, \(-1\)}, 10];
\nX3\ = \ 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\ = \ {\(x'\)[t] == y[t], \(y'\)[t] == \(-4\)*x[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"],
Cell[TextData[{
"Now for some nonlinear systems. Here is a competing species system. By \
judiciously choosing the intial conditions, you can get a pretty \
comprehensive phase portrait. The choice of a and b below leads to \
coexistence of the two species. Try some other values such as: ",
StyleBox["a=0.8,b=0.6", "Input"],
" and ",
StyleBox["a=1.5,b=1.2", "Input"]
}], "Text",
Evaluatable->False],
Cell["\<\
a=1.2;
b=1.5;
ode = {x'[t]==x[t]*(1-x[t]-y[t]),y'[t]==y[t]*(a-x[t]-b*y[t])}
iclist = {{1.5,0.1},{1.5,0.5},{1.5,1},{1.5,1.5},
\t\t {1,1.5},{0.5,1.5},{0.1,1.5},
\t\t {0.01,0.1},{0.05,0.1},{0.1,0.1},
\t\t {0.1,0.05},{0.1,0.01}};
orbs = OrbitList[ode,iclist,20];
ParametricPlot[Evaluate[orbs],{t,0,20},
\tAxes->True,AspectRatio->1,
\tPlotRange->{{0,1.5},{0,1.5}}];\
\>", "Input"],
Cell[TextData[{
"Now for the famous Van der Pol equations. There is a single attracting \
periodic orbit to which most other orbits converge. Try a few other values \
of ",
StyleBox["mu", "Input"],
", such as ",
StyleBox["mu = 0.1, 1, 10", "Input"]
}], "Text",
Evaluatable->False],
Cell["\<\
mu = 0.6;
ode := {x'[t]==y[t],y'[t]==-x[t]+mu*(1-x[t]^2)*y[t]}
iclist = {{4,4},{4,-4},{-4,4},{-4,-4},{0.1,0.1}};
orbs = OrbitList[ode,iclist,20];
ParametricPlot[Evaluate[orbs],{t,0,20},
\tAxes->True,AspectRatio->1,PlotRange->All];
\t\
\>", "Input"],
Cell[TextData[{
"The nonlinear oscillation of the Van der Pol equation is far from the \
usual sinusoidal oscillations produced by linear systems. To illustrate this \
it is better to look at the graphs of the solutions rather than their orbits. \
We can use our ",
StyleBox[" Orbit ", "Input"],
"function to compute a solution and then plot the x and y components of \
this solution as functions of t."
}], "Text"],
Cell[BoxData[
\(mu\ = \ 10; \n{x, y}\ = \ Orbit[ode, {2, 0}, 20]; \n
Plot[{x, y}, {t, 0, 20},
PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 1, 0]}]; \)], "Input"],
Cell[TextData[{
"Try several different ",
StyleBox["mu ", "Input"],
"values in the cell above. For ",
StyleBox["mu=10 ", "Input"],
"the oscillation looks more like a step function than a sine function."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["3D Phase Portraits -- Lorenz and Rossler Attractors", "Section"],
Cell[TextData[{
"To handle 3D nonlinear systems, only a small changeto the functions \
defined in the last section is required. We need a version of the ",
StyleBox["Orbit ", "Input"],
"function which takes a 3D initial vector."
}], "Text"],
Cell[BoxData[
\(Orbit[ode_, {x0_, y0_, z0_}, time_, opts___]\ := \n
\t{x[t], y[t], z[t]} /.
\(NDSolve[
ode~Join~{x[0] == x0, y[0] == y0, z[0] == z0}, {x[t], y[t],
z[t]}, {t, 0, time}, opts]\)[\([1]\)]\)], "Input"],
Cell[TextData[{
"Also, when plotting the results you have to use ",
StyleBox["ParametricPlot3D. ", "Input"],
"As a simple first example, we will try a linear equation with two \
spiralling eigenvalue and one real one."
}], "Text",
Evaluatable->False],
Cell["\<\
ode = {\tx'[t]==-0.2*x[t]-y[t],
\t\ty'[t]==x[t]-0.2*y[t],
\t\tz'[t]==-0.1*z[t]};
iclist = { {1,1,0},{1,1,1},{-1,-1,-1} };
orbs = OrbitList[ode,iclist,30];
pp = ParametricPlot3D[Evaluate[orbs],{t,0,30},
\tBoxRatios->{1,1,1},PlotRange->All];\
\>", "Input"],
Cell["\<\
This shows three orbits spiralling in toward the origin, one in the \
(x,y) plane, one above it and one below it. You can look at this from a \
different point of view in 3D by using the viewpoint option. Since the \
pictureis named pp, you can just use the Show command to redisplay it. \
\>",
"Text",
Evaluatable->False],
Cell["\<\
Show[pp,ViewPoint->{1,2,0}];
Show[pp,ViewPoint->{0,0,1}];\
\>", "Input"],
Cell["\<\
In 3D it is possible for ODE's to have attractors which are more \
complicated than restpoints or limit cycles. Probably the most famous ODE \
with a \"strange\" attractor is the Lorenz equation. One long orbit of this \
is sufficient since all orbits end up on the attractor anyway.\
\>", "Text",
Evaluatable->False],
Cell["\<\
Unprotect[Out];Clear[Out]; (*this is to free up memory*)
ode = {\tx'[t]== 10*(y[t]-x[t]),
\t\ty'[t]== -x[t]*z[t]+28*x[t]-y[t],
\t\tz'[t]== x[t]*y[t]-(8/3)*z[t]};
orb = Orbit[ode,{1,0,10},50,MaxSteps->5000];
ParametricPlot3D[orb,{t,0,50},PlotPoints->5000];\
\>", "Input"],
Cell["\<\
Next we will make a movie by generating a sequence of pictures of \
longer and longer segments of the orbit. This requires lots of memory so to \
prevent crashes there is some code in the loop to remove old output \
variables.\
\>", "Text",
Evaluatable->False],
Cell["\<\
range = {{-30,30},{-30,30},{0,50}};
Do[\tParametricPlot3D[orb,{t,0,5*i},
AspectRatio->1,PlotRange->range,PlotPoints->3000,Compiled->False];Clear[Out];,\
{i,1,10}];\
\>", "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 ",
StyleBox["Cell", "Input"],
" menu and choose the item called ",
StyleBox["Animate Selected Graphics", "Input"],
". You should see a movie of your plots showing how the orbit gradually \
fills in the attractor. 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", "Input"],
" item in the ",
StyleBox["Edit", "Input"],
" menu. ",
StyleBox["Select GraphicsOptions", "Input"],
" and then ",
StyleBox["Animation", "Input"],
". You will find an item called ",
StyleBox["AnimationDisplayTime", "Input"],
". Click on the old value, enter a new one and hit Enter (or Return). \
Increasing it will slow the movie down."
}], "Text",
Evaluatable->False],
Cell["\<\
Another three-dimensional system with a strange attractor is the \
Rossler system. This attractor looks like a fractal Mobius band.\
\>", "Text",\
Evaluatable->False],
Cell["\<\
ode = {\tx'[t]== -y[t]-z[t],
\t\ty'[t]== x[t]+ 0.398*y[t],
\t\tz'[t]== 2 + z[t]*(x[t]-4)};
orb = Orbit[ode,{1,0,0},200,MaxSteps->5000];
ParametricPlot3D[orb,{t,0,200},PlotPoints->2000,PlotRange->All];\
\>", "Input"]
}, Closed]]
},
FrontEndVersion->"Macintosh 3.0",
ScreenRectangle->{{0, 832}, {0, 604}},
WindowToolbars->{},
CellGrouping->Manual,
WindowSize->{721, 527},
WindowMargins->{{-2, Automatic}, {Automatic, -1}},
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<"
]
(***********************************************************************
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, 46, 0, 98, "Title"],
Cell[1758, 51, 876, 25, 82, "Text"],
Cell[2637, 78, 628, 10, 82, "Text"],
Cell[CellGroupData[{
Cell[3290, 92, 64, 0, 50, "Section"],
Cell[CellGroupData[{
Cell[3379, 96, 890, 30, 65, "Text",
Evaluatable->False],
Cell[4272, 128, 93, 1, 27, "Input"],
Cell[4368, 131, 429, 13, 47, "Text"],
Cell[4800, 146, 27, 0, 27, "Input"],
Cell[4830, 148, 1084, 38, 64, "Text",
Evaluatable->False],
Cell[5917, 188, 66, 1, 27, "Input"],
Cell[5986, 191, 1070, 32, 81, "Text",
Evaluatable->False],
Cell[7059, 225, 203, 4, 43, "Input"],
Cell[7265, 231, 451, 11, 65, "Text"],
Cell[7719, 244, 71, 1, 27, "Input"],
Cell[7793, 247, 31, 0, 30, "Text"],
Cell[7827, 249, 59, 1, 27, "Input"],
Cell[7889, 252, 111, 3, 30, "Text"],
Cell[8003, 257, 275, 5, 59, "Input"],
Cell[8281, 264, 718, 15, 98, "Text"],
Cell[9002, 281, 101, 2, 27, "Input"],
Cell[9106, 285, 166, 4, 46, "Text"],
Cell[9275, 291, 180, 3, 59, "Input"]
}, Open ]],
Cell[9470, 297, 57, 0, 30, "Text"],
Cell[9530, 299, 176, 4, 43, "Input"],
Cell[9709, 305, 133, 3, 30, "Text"],
Cell[9845, 310, 225, 4, 75, "Input"],
Cell[10073, 316, 183, 4, 46, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[10293, 325, 74, 0, 30, "Section"],
Cell[10370, 327, 285, 6, 48, "Text"],
Cell[10658, 335, 83, 2, 43, "Input"],
Cell[10744, 339, 270, 6, 47, "Text"],
Cell[11017, 347, 133, 4, 75, "Input"],
Cell[11153, 353, 561, 9, 94, "Text"],
Cell[11717, 364, 193, 4, 91, "Input"],
Cell[11913, 370, 104, 3, 30, "Text"],
Cell[12020, 375, 102, 2, 43, "Input"],
Cell[12125, 379, 214, 4, 46, "Text"],
Cell[12342, 385, 165, 3, 107, "Input"],
Cell[12510, 390, 471, 8, 78, "Text"],
Cell[12984, 400, 143, 2, 43, "Input"],
Cell[13130, 404, 302, 6, 62, "Text"],
Cell[13435, 412, 201, 4, 75, "Input"],
Cell[13639, 418, 316, 6, 62, "Text"],
Cell[CellGroupData[{
Cell[13980, 428, 320, 6, 107, "Input"],
Cell[14303, 436, 302, 5, 62, "Text"],
Cell[14608, 443, 656, 12, 96, "Text"],
Cell[15267, 457, 639, 11, 235, "Input"],
Cell[15909, 470, 334, 6, 62, "Text"],
Cell[16246, 478, 475, 10, 187, "Input"],
Cell[16724, 490, 238, 5, 46, "Text"]
}, Open ]]
}, Closed]],
Cell[CellGroupData[{
Cell[17011, 501, 50, 0, 30, "Section"],
Cell[17064, 503, 444, 8, 78, "Text"],
Cell[17511, 513, 246, 8, 48, "Text"],
Cell[17760, 523, 142, 2, 75, "Input"],
Cell[17905, 527, 332, 7, 48, "Text"],
Cell[18240, 536, 170, 3, 43, "Input"],
Cell[18413, 541, 173, 4, 46, "Text"],
Cell[18589, 547, 91, 1, 27, "Input"],
Cell[18683, 550, 98, 3, 30, "Text"],
Cell[18784, 555, 128, 2, 59, "Input"],
Cell[18915, 559, 54, 0, 30, "Text"],
Cell[18972, 561, 158, 3, 59, "Input"],
Cell[19133, 566, 914, 21, 129, "Text"],
Cell[20050, 589, 217, 4, 46, "Text"],
Cell[20270, 595, 867, 17, 219, "Input"],
Cell[21140, 614, 168, 4, 30, "Text"],
Cell[21311, 620, 181, 3, 75, "Input"],
Cell[21495, 625, 47, 0, 30, "Text"],
Cell[21545, 627, 289, 5, 139, "Input"],
Cell[21837, 634, 909, 21, 129, "Text"],
Cell[22749, 657, 324, 7, 48, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[23110, 669, 38, 0, 30, "Section"],
Cell[23151, 671, 406, 7, 63, "Text"],
Cell[23560, 680, 189, 3, 43, "Input"],
Cell[23752, 685, 98, 3, 30, "Text"],
Cell[23853, 690, 71, 1, 27, "Input"],
Cell[23927, 693, 256, 6, 48, "Text"],
Cell[24186, 701, 86, 1, 27, "Input"],
Cell[24275, 704, 389, 7, 63, "Text"],
Cell[24667, 713, 212, 5, 59, "Input"],
Cell[24882, 720, 134, 3, 30, "Text"],
Cell[25019, 725, 388, 7, 123, "Input"],
Cell[25410, 734, 541, 10, 80, "Text"],
Cell[25954, 746, 112, 2, 27, "Input"],
Cell[26069, 750, 297, 6, 47, "Text"],
Cell[26369, 758, 315, 5, 107, "Input"],
Cell[26687, 765, 418, 9, 64, "Text",
Evaluatable->False],
Cell[27108, 776, 391, 15, 222, "Input"],
Cell[27502, 793, 293, 8, 47, "Text",
Evaluatable->False],
Cell[27798, 803, 261, 11, 162, "Input"],
Cell[28062, 816, 424, 8, 63, "Text"],
Cell[28489, 826, 178, 3, 59, "Input"],
Cell[28670, 831, 228, 6, 31, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[28935, 842, 70, 0, 30, "Section"],
Cell[29008, 844, 248, 5, 48, "Text"],
Cell[29259, 851, 260, 5, 59, "Input"],
Cell[29522, 858, 259, 6, 47, "Text",
Evaluatable->False],
Cell[29784, 866, 266, 10, 147, "Input"],
Cell[30053, 878, 340, 7, 62, "Text",
Evaluatable->False],
Cell[30396, 887, 82, 3, 42, "Input"],
Cell[30481, 892, 332, 6, 46, "Text",
Evaluatable->False],
Cell[30816, 900, 283, 9, 132, "Input"],
Cell[31102, 911, 272, 6, 46, "Text",
Evaluatable->False],
Cell[31377, 919, 189, 6, 72, "Input"],
Cell[31569, 927, 953, 22, 129, "Text",
Evaluatable->False],
Cell[32525, 951, 180, 5, 30, "Text",
Evaluatable->False],
Cell[32708, 958, 226, 7, 102, "Input"]
}, Closed]]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)