(***********************************************************************
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[ 96279, 2547]*)
(*NotebookOutlinePosition[ 97206, 2577]*)
(* CellTagsIndexPosition[ 97162, 2573]*)
(*WindowFrame->Normal*)
Notebook[{
Cell[CellGroupData[{
Cell[TextData[StyleBox["Dynamical Systems and Chaos", "Subtitle"]], "Title"],
Cell[CellGroupData[{
Cell["Initialization", "Section",
InitializationCell->True],
Cell[TextData[{
"This section contains the definitions of a package of ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" functions for computing and plotting orbits of dynamical systems. These \
functions will be used throughout the lab. To evaluate them go up to the \
kernel menu and in the ",
StyleBox["Evaluation",
FontWeight->"Bold"],
" submenu choose the command ",
StyleBox["Evaluate Initialization",
FontWeight->"Bold"],
". It is not necessary to look at the detailed definitions in this section \
but each one is preceded by some words of explanation."
}], "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[TextData[StyleBox[
"Dynamical Systems, Orbits and Iteration", "Subsection"]], "Subsubsection",
InitializationCell->True],
Cell[TextData[{
"One of the interesting aspects of ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" is that the various mathematical objects which arise (numbers, functions, \
etc.) are stored with a kind of label in front. For example the number 5 is \
really ",
StyleBox["Integer[5]", "Input"],
" internally. The label ",
StyleBox["Integer", "Input"],
" is what makes it possible to handle operations differently for integers \
than, say, for real numbers. To make the same kind of distinctions possible \
in this package of functions, we are going to introduce new kinds of objects \
appropriate to dynamical systems theory. Specifically we are going to use \
dynamical system objects, of the form ",
StyleBox["DS[", "Input"],
"function here",
StyleBox["]", "Input"],
" and orbit objects, of the form ",
StyleBox["Orbit[", "Input"],
"list here",
StyleBox["]", "Input"],
". This extra layer of notation will usually not be visible while using \
this package."
}], "Text",
InitializationCell->True],
Cell[TextData[{
"A discrete dynamical system is just a mapping from some set to itself. We \
will represent it in the form ",
StyleBox["f = DS[f1]", "Input"],
" where ",
StyleBox["f1", "Input"],
" is the mapping. For example we could have a mapping of the real line \
like ",
StyleBox["f1[x_] := 4*x*(1-x)", "Input"],
" or a mapping of the plane to itself, such as ",
StyleBox["f1[{x_,y_}]:=", "Input"],
StyleBox["f[{x_,y_}]:={1.28-x^2+0.3*x,y}", "Input",
FontWeight->"Bold"],
". Then the command ",
StyleBox["f = DS[f1]", "Input"],
" turns the mapping into a \"dynamical system object\". You can think of ",
StyleBox["DS", "Input"],
" as an operator which does nothing, except label the mapping as a \
dynamical system."
}], "Text",
InitializationCell->True],
Cell[TextData[{
"We can use this labelling to teach ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" how to iterate. If ",
StyleBox["f=DS[f1]", "Input"],
" is a dynamical system object, we want to be able to use the notation ",
StyleBox["f[x]", "Input"],
" to evaluate the mapping and the notation ",
StyleBox["f^n[x]", "Input"],
" to evaluate the n-th iterate of the mapping. The following commands make \
this possible. The mapping ",
StyleBox["f1", "Input"],
" is recovered from ",
StyleBox["f", "Input"],
" as ",
StyleBox["f[[1]]", "Input"],
", that is, the first element inside ",
StyleBox["f", "Input"],
". The notation ",
StyleBox["DS/:", "Input"],
" in front of the definition of ",
StyleBox["f^n", "Input"],
" assures that the definition is associated with the \"operator\" ",
StyleBox["DS", "Input"],
" rather than with the exponentiation operator. Iteration is basically the \
built-in function ",
StyleBox["Nest", "Input"],
"."
}], "Text",
InitializationCell->True],
Cell[BoxData[
\(\(Off[General::spell1]\ (*
first\ turn\ off\ annoying\ spelling\ messages*) \)\)], "Input",
InitializationCell->True],
Cell[BoxData[{
\(f_DS[x_]\ := \ \(f[\([1]\)]\)[x]\ \ \n\),
\(DS /: \ f_DS^n_Integer[x_]\ := \ Nest[f[\([1]\)], x, n]\)}], "Input",
InitializationCell->True],
Cell[TextData[{
"An orbit will be a list of states enclosed in an ",
StyleBox["Orbit[]", "Input"],
" label. We will use the notation ",
StyleBox["F^{n}[x0]", "Input"],
" to represent the orbit ",
StyleBox["Orbit[{x0,x1,...,xn}]", "Input"],
" of some initial point and ",
StyleBox["F^{m,n}", "Input"],
" to represent ",
StyleBox["Orbit[{xm,x1,...,xn}]", "Input"],
". Orbits are computed using the built-in function ",
StyleBox["NestList", "Input"],
"."
}], "Text",
InitializationCell->True],
Cell[BoxData[{
\(DS /: \ f_DS^{n_Integer}[x_]\ := \
NestList[f[\([1]\)], x, n] // Orbit\n\),
\(DS /: \ f_DS^{m_Integer, n_Integer}[x_]\ := \
NestList[f[\([1]\)], Nest[f[\([1]\)], x, m], n - m] // Orbit\)}],
"Input",
InitializationCell->True],
Cell[TextData[{
"Suppose that ",
StyleBox["orb", "Input"],
" is an orbit object ",
StyleBox["Orbit[{x0,x1,...,xn}]", "Input"],
". We will introduce the notations ",
StyleBox["orb[i]", "Input"],
" for ",
StyleBox["xi, orb[]", "Input"],
" for the list ",
StyleBox["{x0,x1,...,xn}", "Input"],
" without the enclosing ",
StyleBox["Orbit[ ] ", "Input"],
"bracket, ",
StyleBox["orb[{m}]", "Input"],
" for the partial list ",
StyleBox["{x0,x1,...,xm},", "Input"],
" and ",
StyleBox["orb[{m,n}]", "Input"],
" for ",
StyleBox["{xm,x1,...,xn}", "Input"],
". Since orbits are indexed starting with ",
StyleBox["x0", "Input"],
" and ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" lists are indexed from index 1, an ",
StyleBox["n+1", "Input"],
" is needed."
}], "Text",
InitializationCell->True],
Cell[BoxData[{
\(l_Orbit[]\ := \ l[\([1]\)]\n\t\t\),
\(l_Orbit[n_Integer]\ := \ l[\([1, n + 1]\)]\n\),
\(l_Orbit[{n_Integer}]\ := \ Take[l[\([1]\)], n + 1]\n\),
\(l_Orbit[{m_Integer, n_Integer}]\ := \
Take[l[\([1]\)], {m + 1, n + 1}]\)}], "Input",
InitializationCell->True],
Cell[TextData[{
"Sometimes we want to find the iterates or orbits of a list of initial \
states. The following definitions extend the built-in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" function Map to handle this situation."
}], "Text",
InitializationCell->True]
}, Closed]],
Cell[CellGroupData[{
Cell["Periodic Orbits", "Subsection",
InitializationCell->True],
Cell[TextData[{
"This group of functions is concerned with finding periodic orbits and \
computing their multipliers. For the purposes of these functions, an orbit \
of period k will be represented by an ",
StyleBox["Orbit", "Input"],
" of length k+1 of the form ",
StyleBox["Orbit[{p0,p1,p2,....,pk}]", "Input"],
" with ",
StyleBox["pk=p0", "Input"],
"."
}], "Text",
InitializationCell->True],
Cell["\<\
The first function just tests an orbit to see if it seems to be \
periodic. The test is just to see when it comes back with a certain small \
distance of itself. For vectors the distance is measured using the maximum \
absolute value of the differences of the coordinates instead of with the \
Euclidean distance.\
\>", "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
StyleBox[
\(MinimalPeriod[orb_Orbit, tol_: 10^\((\(-6\))\)]\ := \n\t
With[{diff\ = \
Map[Max[Abs[#]]&,
Drop[orb[] - Table[orb[0], {Length[orb[]]}], 1]]}, \n\t\t
Min[Position[diff, _?\((Function[d, d <= tol])\)]]]\n\t\t\ \),
InitializationCell->True,
AspectRatioFixed->True],
\(MinimalPeriod::usage\ = \
"\"\)}], "Input",\
InitializationCell->True],
Cell[BoxData[
\("MinimalPeriod[orb,tol] tries to find the minimal period of orb = \
Orbit[{x0,x1,...}] by searching for the an item xk in the sequence which \
agrees with x0 to the given tolerance. The xi may be vectors."\)], "Output"]
}, Open ]],
Cell[TextData[{
"If you have a periodic orbit you can compute the multipliers by \
multiplying the derivatives at all points of the orbit. ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" can compute derivatives if you supply the formula and the following \
functions make use of this. The first function just computes the derivative \
of a given function at a given point. In the vector case, the derivative is \
the Jacobian matrix. There are two versions of the function for that \
reason."
}], "Text",
InitializationCell->True],
Cell[BoxData[{
\(\nScalarQ[x_]\ := \ \(! \((Head[x] === List)\)\) // LogicalExpand\n\t
\),
\(DerAtPoint[f_DS, x0_?ScalarQ]\ := Module[{x}, D[f[x], x] /. x -> x0]\n
\),
\(DerAtPoint[f_DS, x0_List]\ :=
Module[{x, xx}, x\ = \ Array[xx, Length[x0]];
Transpose[
Map[\ D[f[x], #]&, \ x] /.
Table[xx[i] -> x0[\([i]\)], {i, 1, Length[x0]}]]\n\t\t]\)}],
"Input",
InitializationCell->True],
Cell[TextData[{
"Next we multiply the derivatives at {p0,p1,...,p(k-1)} for a given \
periodic orbit {p0,p1,...,pk}. For the vector case we need to use matrix \
multiplication, called ",
StyleBox["Dot", "Input",
FontWeight->"Bold"],
"."
}], "Text",
InitializationCell->True],
Cell[BoxData[{
\(ScalarOrbitQ[orb_]\ := \
\((Head[orb] === Orbit)\) && \((Length[orb[]] >= 1)\) &&
ScalarQ[orb[1]] // LogicalExpand\n\),
\(VectorOrbitQ[orb_]\ := \
\((Head[orb] === Orbit)\) && \((Length[orb[]] >= 1)\) &&
VectorQ[orb[1]] // LogicalExpand\n\),
\(DerAlongOrbit[f_DS, orb_?ScalarOrbitQ] := \ \n\t
If[Length[orb[]] === 1, 2,
Apply[Times, Map[DerAtPoint[f, #]&, Drop[Reverse[orb[]], 1]]]\n\t\t]\n
\),
\(DerAlongOrbit[f_DS, orb_?VectorOrbitQ]\ := \n\t
If[Length[orb[]] === 1, IdentityMatrix[Length[orb[1]]], \ \n\t\t
Apply[Dot, Map[DerAtPoint[f, #]&, Drop[Reverse[orb[]], 1]]]\n\t\t]
\)}], "Input",
InitializationCell->True],
Cell["\<\
Finally, the function below computes the multiplier of a given \
periodic orbit. In the vector case, it returns the eigenvalues of the \
product of the derivative matrices along the orbit.\
\>", "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
\(Multiplier[f_DS, orb_?ScalarOrbitQ]\ := \ DerAlongOrbit[f, orb]\n\),
\(Multiplier[f_DS, orb_?VectorOrbitQ]\ := \ \t\n\t\t
With[{A = DerAlongOrbit[f, orb]}, Eigenvalues[A]]\n\),
\(Multiplier::usage\ = \
"\"\)}], "Input",
InitializationCell->True],
Cell[BoxData[
\("Multiplier[f,orb] computes the multiplier or multipliers of a periodic \
orbit, orb = Orbit[{x0,...,xq}] of a dynamical system f by multiplying the \
derivatives of f at x0,...,x(q-1). For maps in dimension greater than 1 the \
derivatives will be matrices and Multiplier returns the list of eigenvalues."\
\)], "Output"]
}, Open ]],
Cell[TextData[{
"Next we have a fairly complicated function called ",
StyleBox["SeekPeriodicPoint ", "Input"],
"which searches for periodic points. It tries to solve the equation f^k(x) \
= x using ",
StyleBox["Mathematica",
FontSlant->"Italic"],
"'s built-in ",
StyleBox["FindRoot", "Input"],
" command which is based on Newton's method. The search is conducted \
starting at a given point, but if that fails nearby points are chosen at \
random. So we need to define a ",
StyleBox["RandomPoint", "Input"],
" function first, which produced a number or vector chosen at random from a \
given range. In the scalar case, the range should be an interval given as a \
vector ",
StyleBox["{xmin,xmax}", "Input"],
". In the vector case the range will be specified by a list of two vectors \
-- a vector of minima and a vector of maxima. For example ",
StyleBox["{{xmin,ymin,zmin},{xmax,ymax,zmax}}", "Input"],
" represents a cube in three dimensions."
}], "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
StyleBox[\(RandomPoint[r_?VectorQ]\ := \tRandom[Real, r]; \n\n
RandomPoint[r : {{___}, {___}}]\ := \
Map[Random[Real, #]&, Transpose[r]]\n\),
InitializationCell->True,
AspectRatioFixed->True],
\(RandomPoint::usage\ = \
"\"\)}], "Input",
InitializationCell->True],
Cell[BoxData[
\("RandomPoint[range] returns a random vector in the given range. For an \
n-dimensional point, the range should be a pair of the form \
{{x1min,x2min,...},{x1max,x2max,...}}."\)], "Output"]
}, Open ]],
Cell[TextData[{
"There are several options for ",
StyleBox["SeekPeriodicPoint", "Input"],
" which constrain the search procedure."
}], "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
RowBox[{
StyleBox[
\(Options[SeekPeriodicOrbit]\ =
\ {\n\tRadius -> Automatic, \n\tMaxTries -> 10, \n\t
Tolerance -> 0.000001, \n\tMaxIterations -> 100\n\t}; \n
spooptpattern\ = \
\((Radius -> _)\) | \((MaxTries -> _)\) | \((Tolerance -> _)\); \n\n
\nSeekPeriodicOrbit[f_DS, k_Integer, xguess_, opts___]\ := \
Module[{fk, guess, sol, x, xx, x0, po, ic, mt, mi, tol, ra, xrange,
fropts}, \n\t\t
mt\ = \ \(MaxTries /. {opts}\) /. Options[SeekPeriodicOrbit]; \n
\t\ttol\ = \
\(Tolerance /. {opts}\) /. Options[SeekPeriodicOrbit]; \n\t\t
ra\ = \ \(Radius /. {opts}\) /. Options[SeekPeriodicOrbit]; \n
\t\tfropts\ = \ DeleteCases[{opts}, spooptpattern]; \n
\t\t (*\ Try\ the\ given\ initial\ point\ first\ *) \n\t\t
guess\ = \ xguess; \n\t\t
If[VectorQ[guess], \n\t\t\tx = Array[xx, Length[guess]];
ic = Transpose[{x, guess}], \n\t\t\tic = {{x, guess}}]; \n\t\t
fk\ = \ f^k[x] - x // Simplify; \n\t\t
sol\ = \ \(FindRoot[fk, ##]&\)\ @@\ Flatten[{ic, fropts}, 1]; \n
\t\tx0\ = \ x /. sol; \n\t\tpo\ = \ f^{k}[x0]; \n\t\t
If[MinimalPeriod[po] === k, Return[po]]; \n
\t\t (*\ Failing\ this, \ try\ random\ guesses\ nearby\ *) \n\t\t
If[ra === Automatic, \n\t\t\t
ra = 2*Max[Abs[xguess - f^k[xguess]]]; \n
\t\t\t (*ra\ = \ {0.1, 0.1}; *) \n\t\t\t
xrange\ = \ {xguess - ra, xguess + ra}]; \n\t\t
While[mt > 1, \n\t\t\tguess\ = \ RandomPoint[xrange]; \n\t\t\t
If[VectorQ[guess], \n\t\t\t\tic = Transpose[{x, guess}], \n
\t\t\t\tic = {{x, guess}}\n\t\t\t\t]; \n\t\t\t
x0\ = \ x /.
\((\(FindRoot[fk, ##]&\)\ @@Flatten[{ic, fropts}, 1])\); \n
\t\t\tpo\ = \ f^{k}[x0]; Print[po]; \n\t\t\t
If[MinimalPeriod[po] === k, Return[po]]; \n\t\t\t\(mt--\)]; \n
\t\t (*\ Failing\ this, \ return\ the\ empty\ list\ *) \n\t\t
Return[Orbit[{}]]; \n\t\t]\),
InitializationCell->True,
AspectRatioFixed->True], "\n"}],
RowBox[{
StyleBox[
\(SeekPeriodicOrbit::usage\ = \
"\"\),
InitializationCell->True,
AspectRatioFixed->True], "\n"}]}], "Input",
InitializationCell->True],
Cell[BoxData[
\("SeekPeriodicOrbit[f,k,xguess,opts] tries to find a periodic orbit of \
period k for the dynamical system f near the point xguess using the \
Mathematica function FindRoot. Works well only for small k and low \
dimensional systems. If it succeeds,it returns Orbit[{x0,x1,...,xk}] of \
length k+1. If it fails,it returns Orbit[{}]."\)], "Output"]
}, Open ]]
}, Closed]],
Cell[CellGroupData[{
Cell["Unstable Manifolds", "Subsection",
InitializationCell->True],
Cell["\<\
Here is a function for computing unstable manifolds of 2D maps. \
Supposing that we have already located a period q point x0, we first test it \
to make sure it is a saddle point. Then we locate the unstable eigenvector, \
construct a short line segment in that direction, subdivide it into many \
smaller segments and then extend these by iteration of the dynamical \
system.\
\>", "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
\(IsSaddle[f_DS, x0_, q_]\ := \
With[\n\t{mu\ = \ Multiplier[f, f^{0, q}[x0]]}, \n\t
If[Im[mu] != {0, 0}, Return[False], ]; \n\t
If[\tAbs[mu[\([1]\)]] > 1 && Abs[mu[\([2]\)]] < 1 || \n\t\ \ \t
Abs[mu[\([2]\)]] > 1 && Abs[mu[\([1]\)]] < 1, True, False]\n\t]\n
\t\),
\(Subdivide[l_, i_, n_: 2]\ := \
Take[l, i - 1]~Join~\n\t\t
Table[N[l[\([i]\)]*\((1 - t/n)\) + l[\([i + 1]\)]*t/n], {t, 0, n}]~
Join~\n\t\t\t\tDrop[l, i + 1]\n\),
\(ListOrbit[f_, l_List, n_]\ := \
Flatten[Transpose[\n\t\t\tMap[\ NestList[f, #, n]&, l]], 1]\n\),
\(Options[UnstableManifold]\ =
\ {\n\tPeriod -> 1, \n\tSubdivisions -> 20, \n\tRadius -> 0.01}; \n\n
UnstableManifold[f_DS, x0_, n_, opts___]\ := \ \n\t
Module[{q, su, ra, A, mu, ev, vu, l, l1, l2, g}, \n\t
q\ = \ \(Period /. {opts}\) /. Options[UnstableManifold]; \n\t
su\ = \ \(Subdivisions /. {opts}\) /. Options[UnstableManifold]; \n\t
ra\ = \ \(Radius /. {opts}\) /. Options[UnstableManifold]; \n\t
If[\(! IsSaddle[f, x0, q]\), Print["\"]; Return[{}]];
\n\tA\ = \ DerAlongOrbit[f, f^{0, q}[x0]]; \n\t
mu\ = \ Eigenvalues[A]; \n\tev\ = \ Eigenvectors[A]; \n\t
vu\ = \ If[Abs[mu[\([1]\)]] > 1, ev[\([1]\)],
mu[\([1]\)] = mu[\([2]\)]; ev[\([2]\)]]; \n\t
g\ = \ If[mu[\([1]\)] > 0, f^q[#]&, f^\((2*q)\)[#]&]; \n\t
l1\ = \ ListOrbit[g, Subdivide[{x0 + ra*vu, g[x0 + ra*vu]}, 1, su],
n]; \n\tl2\ = \
ListOrbit[g, Subdivide[{x0 - ra*vu, g[x0 - ra*vu]}, 1, su], n]; \n\t
Orbit[Reverse[l2]~Join~l1]\n\t]\n\t\),
\(UnstableManifold::usage\ = \
"\q. A \
number of points are chosen along the unstable eigenvector close to x0. The \
number of points is determined by the option Subdivisions which defaults to \
20. The distance from x0 is controlled by the option Radius which defaults to \
0.01. The argument n is the number of iterations used to extend the manifold \
outside this radius.\>"\)}], "Input",
InitializationCell->True],
Cell[BoxData[
\("UnstableManifold[f,x0,n,opts] returns a list of points on the unstable \
manifold of the periodic point x0 for the map f. By default, x0 is a fixed \
point; period q can be handled by setting the option Period->q. A number of \
points are chosen along the unstable eigenvector close to x0. The number of \
points is determined by the option Subdivisions which defaults to 20. The \
distance from x0 is controlled by the option Radius which defaults to 0.01. \
The argument n is the number of iterations used to extend the manifold \
outside this radius."\)], "Output"]
}, Open ]]
}, Closed]],
Cell[CellGroupData[{
Cell["Plotting", "Subsection",
InitializationCell->True],
Cell[TextData[{
"This section contains some functions for plotting the orbits of dynamical \
systems in interesting ways. There are many different types of plots, all of \
which take similar kinds of options, such as the ranges for the variables and \
the sizes and colors of the dots and lines. So all the plotting functions \
below are based on a ",
StyleBox["GenericPlot ", "Input"],
"function which handles all of these details. Then the only difference \
between the various plots is some preliminary transformations which they \
perform on the orbits. And you can even make up new types of plots by \
defining new transformations."
}], "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell["Generic Plotting Functions", "Subsubsection",
InitializationCell->True],
Cell[TextData[{
"First there are some functions for turning a list of 2D or 3D vectors into \
a ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" graphics object -- either a list of ",
StyleBox["Point", "Input"],
" 's, a single ",
StyleBox["Line", "Input"],
" or a ",
StyleBox["Polygon", "Input"],
". There are optional arguments to control the color and size. There are \
two versions of each function, one for 2D and one for 3D."
}], "Text",
InitializationCell->True],
Cell[BoxData[{
RowBox[{
\(ScalarListQ[l_]\ := \
\((Head[l] === List)\) && \((Length[l] >= 1)\) &&
ScalarQ[l[\([1]\)]] // LogicalExpand\), "\n"}],
RowBox[{
\(TwoDListQ[l_]\ := \
\((Head[l] === List)\) && \((Length[l] >= 1)\) &&
\((Head[l[\([1]\)]] === List)\) &&
\((Length[l[\([1]\)]] === 2)\) // LogicalExpand\), "\n"}],
RowBox[{
\(ThreeDListQ[l_]\ := \
\((Head[l] === List)\) && \((Length[l] >= 1)\) &&
\((Head[l[\([1]\)]] === List)\) &&
\((Length[l[\([1]\)]] === 3)\) // LogicalExpand\), "\n", "\n"}],
StyleBox[
\(MakePoints[l_?TwoDListQ, color_: {0, 0, 0}, size_: 0.01]\ := \n\t
Graphics[{Apply[RGBColor, color], PointSize[size], Map[Point, l]}]\n
\),
InitializationCell->True,
AspectRatioFixed->True],
StyleBox[
\(MakeLine[l_?TwoDListQ, color_: {0, 0, 0}, thickness_: 0.005]\ := \n\t
Graphics[{Apply[RGBColor, color], Thickness[thickness], Line[l]}]\n
\),
InitializationCell->True,
AspectRatioFixed->True],
StyleBox[
\(MakePolygon[l_?TwoDListQ, color_: {0, 0, 0}, size_: 0.01]\ := \n\t
Graphics[{Apply[RGBColor, color], PointSize[size], Polygon[l]}]\n\t
\n (*\ 3 D\ versions\ *) \t\t\t\t\t\t\t\t\t\),
InitializationCell->True,
AspectRatioFixed->True],
StyleBox[
\(MakePoints[l_?ThreeDListQ, color_: {0, 0, 0}, size_: 0.01]\ := \n\t
Graphics3D[{Apply[RGBColor, color], PointSize[size], Map[Point, l]}]\n
\),
InitializationCell->True,
AspectRatioFixed->True],
StyleBox[
\(MakeLine[l_?ThreeDListQ, color_: {0, 0, 0}, thickness_: 0.005]\ := \n
\tGraphics3D[{Apply[RGBColor, color], Thickness[thickness], Line[l]}]
\n\),
InitializationCell->True,
AspectRatioFixed->True],
StyleBox[
\(MakePolygon[l_?ThreeDListQ, color_: {0, 0, 0}, size_: 0.01]\ := \n\t
Graphics3D[{Apply[RGBColor, color], PointSize[size], Polygon[l]}]\),
InitializationCell->True,
AspectRatioFixed->True]}], "Input",
InitializationCell->True],
Cell["\<\
Here are a few definitions of popular RGB colors so you don't have \
to supply explicit R,G,and B values.\
\>", "Text",
InitializationCell->True],
Cell[BoxData[
\(\(ColorTable\ =
\ {Black -> {0, 0, 0}, White -> {1, 1, 1}, Red -> {1, 0, 0}, \n\t
Green -> {0, 1, 0}, Blue -> {0, 0, 1}, Yellow -> {1, 1, 0},
Purple -> {0.5, 0, 1}, \n\tBlueGreen -> {0, 1, 1},
LightGray -> {0.75, 0.75, 0.75}, Gray -> {0.5, 0.5, 0.5}, \n\t
DarkGray -> {0.25, 0.25, 0.25}}; \)\)], "Input",
InitializationCell->True],
Cell[TextData[{
"Some of the plots included decorations like a copy of a graph of a \
function, a diagonal line, or a circle. These also need to be converted to ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" graphics objects."
}], "Text",
InitializationCell->True],
Cell[BoxData[{
RowBox[{
\(MakeDiagonal[x_, color_: {0, 0, 0}, thickness_: 0.005]\ := \ \n\t
Graphics[{Apply[RGBColor, color], Thickness[thickness], \n
\t\t\t\t\t\t\t\t\t\t
Line[{{x[\([1]\)], x[\([1]\)]}, {x[\([2]\)], x[\([2]\)]}}]}]\),
"\n"}],
RowBox[{
StyleBox[
\(MakeCircle[rad_: 1.0, pts_: 100, color_: {0, 0, 0},
thickness_: 0.005]\ := \n\t
With[{d\ = \ N[2 Pi/pts]},
Graphics[{Apply[RGBColor, color], Thickness[thickness],
Line[Table[{Cos[d*i], Sin[d*i]}, {i, 0, pts}]]}]]\),
InitializationCell->True,
AspectRatioFixed->True], "\n", "\t\t\t\t\t\t\t\t\t\t"}],
\(MakeFunctionGraph[f_, x_, color_: {0, 0, 0}, thickness_: 0.005]\ := \ \n
\tWith[{r\ = \
Range[x[\([1]\)], x[\([2]\)],
0.005*\((x[\([2]\)] - x[\([1]\)])\)]},
Graphics[{Apply[RGBColor, color], Thickness[thickness],
Line[Transpose[{r, Map[f, r]}]]}]]\)}], "Input",
InitializationCell->True],
Cell["\<\
Here is the generic plotting function. The first version is for \
plotting a single orbit. There are numerous options\
\>", "Text",
InitializationCell->True],
Cell[BoxData[
\(\(Options[GenericPlot]\ =
\ {\n\t\tPoints -> True, \n\t\tLines -> True, \n\t\t
PointColor -> Black, \n\t\tLineColor -> Black, \n\t\t
PointSize -> 0.01, \n\t\tLineSize -> 0.002, \n\t\t
Transformation -> IdentityTransformation, \n\t\tMoreGraphics -> {}\n
\t\t}; \ngenoptpattern\ = \
\((Points -> _)\) | \((Lines -> _)\) | \((PointColor -> _)\) | \n\t\t
\((LineColor -> _)\) | \((PointSize -> _)\) | \((LineSize -> _)\) | \n
\t\t\((Transformation -> _)\) | \((MoreGraphics -> _)\); \n\n
GenericPlot[orb_, opts___]\ := \ \n\t\t
Module[{showopts, pf, lf, pc, lc, ps, ls, pr, tr, mg, \n\t\t\t\txmin,
xmax, trorb, gr}, \n\t\t\t
pf\ = \ \(Points /. {opts}\) /. Options[GenericPlot]; \n\t\t\t
pc\ = \ \(\(PointColor /. {opts}\) /. Options[GenericPlot]\) /.
ColorTable; \n\t\t\t
ps\ = \ \(PointSize //. {opts}\) /. Options[GenericPlot]; \n\t\t\t
lf\ = \ \(Lines /. {opts}\) /. Options[GenericPlot]; \n\t\t\t
lc\ = \ \(\(LineColor /. {opts}\) /. Options[GenericPlot]\) /.
ColorTable; \n\t\t\t
ls\ = \ \(LineSize /. {opts}\) /. Options[GenericPlot]; \n\t\t\t
pr\ = \ PlotRange /. {opts}; \n\t\t\t
tr\ = \ \(Transformation /. {opts}\) /. Options[GenericPlot]; \n
\t\t\tmg\ = \ \(MoreGraphics /. {opts}\) /. Options[GenericPlot]; \n
\t\t\tshowopts\ = \
DeleteCases[{opts}, genoptpattern]~Join~\n
\t\t\t\t\t\t{Axes -> True, AspectRatio -> 1}; \n\t\t\t
xmin\ = \ Min[orb[]]; xmax\ = \ Max[orb[]]; \n\t\t\t
xmin\ -= \ 0.1*\((xmax - xmin)\); xmax\ += \ 0.1*\((xmax - xmin)\);
\n\t\t\tIf[Head[pr] == List, xmin = pr[\([1, 1]\)];
xmax = pr[\([1, 2]\)], ]; \n\t\t\ttrorb\ = \ tr[orb]; \n\t\t\t
gr\ = \ {\n\t\t\t\tIf[lf, MakeLine[trorb, lc, ls], {}], \n\t\t\t\t
If[pf, MakePoints[trorb, pc, ps], {}]\n\t\t\t\t}~Join~mg; \n
\t\t\tShow[gr, showopts]\n\t\t\t]\n\)\)], "Input",
InitializationCell->True],
Cell["The default identity transformation is defined here.", "Text",
InitializationCell->True],
Cell[BoxData[
\(IdentityTransformation[orb_]\ := \ orb[]\)], "Input",
InitializationCell->True],
Cell["Here is the version for plotting a list of orbits.", "Text",
InitializationCell->True],
Cell[BoxData[{
RowBox[{
\(OrbitListQ[l_]\ :=
\((Head[l] === List)\) && \((Length[l] >= 1)\) &&
\((Head[l[\([1]\)]] === Orbit)\) // LogicalExpand\), "\n"}],
StyleBox[
\(GenericPlot[orbs_?OrbitListQ, opts___]\ := \ \n\t\t
Module[{num, showopts, pf, lf, pc, lc, ps, ls, pr, tr, mg, \n\t\t\t\t
xmin, xmax, trorbs, gr}, \n\t\t\t
pf\ = \ \(Points /. {opts}\) /. Options[GenericPlot]; \n\t\t\t
pc\ = \ \(\(PointColor /. {opts}\) /. Options[GenericPlot]\) /.
ColorTable; \n\t\t\t
ps\ = \ \(PointSize //. {opts}\) /. Options[GenericPlot]; \n\t\t\t
lf\ = \ \(Lines /. {opts}\) /. Options[GenericPlot]; \n\t\t\t
lc\ = \ \(\(LineColor /. {opts}\) /. Options[GenericPlot]\) /.
ColorTable; \n\t\t\t
ls\ = \ \(LineSize /. {opts}\) /. Options[GenericPlot]; \n\t\t\t
pr\ = \ PlotRange /. {opts}; \n\t\t\t
tr\ = \ \(Transformation /. {opts}\) /. Options[GenericPlot]; \n
\t\t\tmg\ = \ \(MoreGraphics /. {opts}\) /. Options[GenericPlot];
\n\t\t\t\t (*\ make\ the\ options\ into\ lists\ if\ necessary\ *) \n
\t\t\tnum\ = \ Length[orbs]; \n\t\t\t
If[Head[pf] =!= List, pf = Table[pf, {num}]]; \n\t\t\t
If[VectorQ[pc], pc\ = \ Table[pc, {num}]]; \n\t\t\t
If[Head[ps] =!= List, ps = Table[ps, {num}]]; \n\t\t\t
If[Head[lf] =!= List, lf = Table[lf, {num}]]; \n\t\t\t
If[VectorQ[lc], lc = Table[lc, {num}]]; \n\t\t\t
If[Head[ls] =!= List, ls = Table[ls, {num}]]; \n\t\t\t
showopts\ = \
DeleteCases[{opts}, genoptpattern]~Join~\n
\t\t\t\t\t\t{Axes -> True, AspectRatio -> 1}; \n\t\t\t
xmin\ = \ Min[Union[orbs]]; xmax\ = \ Max[Union[orbs]]; \n\t\t\t
xmin\ -= \ 0.2*\((xmax - xmin)\);
xmax\ += \ 0.2*\((xmax - xmin)\); \n\t\t\t
If[Head[pr] == List, xmin = pr[\([1, 1]\)];
xmax = pr[\([1, 2]\)], ]; \n\t\t\ttrorbs\ = \ Map[tr, orbs]; \n
\t\t\tgr\ = \
Flatten[{\n\t\t\t\t
MapThread[If[#2, MakeLine[#1, #3, #4], {}]&, \n
\t\t\t\t\t\t\t\t{trorbs, lf, lc, ls}], \n\t\t\t\t
MapThread[If[#2, MakePoints[#1, #3, #4], {}]&, \n
\t\t\t\t\t\t\t\t{trorbs, pf, pc, ps}]\n\t\t\t\t}]~Join~
mg; \n\t\t\tShow[gr, showopts]\n\t\t\t]\n\),
FontWeight->"Bold"]}], "Input",
InitializationCell->True]
}, Closed]],
Cell[CellGroupData[{
Cell["Plots for One-dimensional Orbits", "Subsubsection",
InitializationCell->True],
Cell[TextData[{
"Each plot is mainly determined by the preliminary transformation it does \
to its orbits. The ",
StyleBox["CobwebPlot", "Input"],
" transforms a scalar orbit ",
StyleBox["Orbit[{x0,x1,...}] ", "Input"],
"to a list of pairs ",
StyleBox["{{x0,x0},{x0,x1},{x1,x1},{x1,x2},...} ", "Input"],
"suitable for plotting in the plane. It also contains a copy of the \
diagonal and a copy of the graph of the dynamical system, added with the ",
StyleBox["MoreGraphics", "Input"],
" option of ",
StyleBox["GenericPlot", "Input"],
"."
}], "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
RowBox[{
\(CobwebTransformation[orb_?ScalarOrbitQ]\ := \
With[{dup = Flatten[Transpose[{orb[], orb[]}]]}, \n\t\t\t\t\t
Transpose[{Drop[dup, \(-1\)], Drop[dup, 1]}]]\), "\n"}],
StyleBox[
\(CobwebPlot[f_DS, orb_, xrange_, opts___]\ := \ \n\t
GenericPlot[orb, opts, Transformation -> CobwebTransformation,
Points -> False, \n\t\tPlotRange -> {xrange, xrange}, \n\t\t
PlotLabel -> StringForm["\", Evaluate[f["\"]]], \n\t\t
MoreGraphics -> {MakeDiagonal[xrange], \n\t\t\t\t\t\t
MakeFunctionGraph[f, xrange]}\n\t\t\t\t]\n\t\t\t\t\),
InitializationCell->True,
AspectRatioFixed->True],
StyleBox[
\(CobwebPlot::usage\ = \
"\{Black,Red}] will make a plot \
showing the two orbits in the given colors.\>"\),
InitializationCell->True,
AspectRatioFixed->True]}], "Input",
InitializationCell->True],
Cell[BoxData[
\("CobwebPlot[f,orb,xrange,opts] produces a graphical analysis or cobweb \
plot showing the graph of the function f, the diagonal and the orbit orb. \
xrange is a list giving the left and right boundaries of the plot. opts is \
any number of options. A list of orbits can be given instead of a single \
orbit. For example if f is a dynamical system and orb1 and orb2 are orbits \
which have already been computed, then \
CobwebPlot[f,{orb1,orb2},{0.0,1.0},LineColor->{Black,Red}] will make a plot \
showing the two orbits in the given colors."\)], "Output"]
}, Open ]],
Cell[TextData[{
"Another useful plot for scalar dynamical systems is the time series plot. \
This time the scalar orbit ",
StyleBox["Orbit[{x0,x1,...}] ", "Input"],
"is transformed into the list of pairs ",
StyleBox["{{0,x0},{1,x1},{2,x2},...}", "Input"]
}], "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
RowBox[{
\(TimeSeriesTransformation[orb_?ScalarOrbitQ]\ := \
Map[Flatten, Transpose[{Range[0, Length[orb[]] - 1], orb[]}]]\),
"\n", "\n"}],
StyleBox[
\(TimeSeriesPlot[orb_, opts___]\ := \ \n\t
GenericPlot[orb, opts, Transformation -> TimeSeriesTransformation,
PointSize -> 0.02, AspectRatio -> 0.5,
AxesLabel -> {"\", "\"}]\n\t\t\t\t\),
InitializationCell->True,
AspectRatioFixed->True],
StyleBox[
\(TimeSeriesPlot::usage\ = \
"\{Black,Red}] will make a plot showing \
the two orbits in the given colors.\>"\),
InitializationCell->True,
AspectRatioFixed->True]}], "Input",
InitializationCell->True],
Cell[BoxData[
\("TimeSeriesPlot[orb,opts] produces a time series plot showing the graph \
of the orbit orb as a function of time index n. opts is any number of \
options. A list of orbits can be given instead of a single orbit. For \
example if orb1 and orb2 have already been computed using Orbit[...], then \
TimeSeriesPlot[{orb1,orb2},LineColor->{Black,Red}] will make a plot showing \
the two orbits in the given colors."\)], "Output"]
}, Open ]],
Cell[TextData[{
"The final plots for scalar data are two and three dimensional delay plots. \
In the two-dimensional case, the scalar orbit ",
StyleBox["Orbit[{x0,x1,...}] ", "Input"],
"is transformed into the list of pairs ",
StyleBox["{{x0,xd},{x1,x(d+1)},{x2,x(d+2)},...}", "Input"],
" where ",
StyleBox["d", "Input"],
" is a delay which defaults to ",
StyleBox["d=1", "Input"],
"."
}], "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
RowBox[{
\(DelayTransformation[orb_?ScalarOrbitQ, dim_: 2, delay_: 1]\ := \ \n
\tWith[{ln = Length[orb[]] - \((dim - 1)\)*delay},
Map[Flatten,
Transpose[
Table[Take[orb[], {1 + i*delay, ln + i*delay}], {i, 0,
dim - 1}]]]\n\t]\), "\n", "\t\t"}],
StyleBox[
\(DelayPlot[orb_, dim_Integer: 2, delay_Integer: 1, opts___]\ := \n\t
Switch[dim, \ \n\t\t2,
GenericPlot[orb, opts, Lines -> False, \n\t\t\t
Transformation -> \((DelayTransformation[#, dim, delay]&)\), \n
\t\t\tAxesLabel -> {"\", "\"}], \n\t\t3,
GenericPlot[orb, opts, Lines -> False, \n\t\t\t
Transformation -> \((DelayTransformation[#, dim, delay]&)\), \n
\t\t\tAxesLabel -> {"\", "\", "\"}]\n
\t\t\t]\n\t\t\t\),
InitializationCell->True,
AspectRatioFixed->True],
RowBox[{
StyleBox[
\(DelayPlot::usage\ = \
"\"\),
InitializationCell->True,
AspectRatioFixed->True], "\n"}]}], "Input",
InitializationCell->True],
Cell[BoxData[
\("DelayPlot[orb,dim,delay,opts] produces a delay plot showing the orbit \
orb embedded in dim dimensions, using a given delay. opts is any number of \
options. dim defaults to 2 and delay to 1 so, for example, the command \
DelayPlot[orb] will make a plot showing {orb_n,orb_n+1} in the plane. A list \
of orbits can be given instead of a single orbit."\)], "Output"]
}, Open ]],
Cell[BoxData[
\(CircleOrbit[orb_]\ := \
With[{tp\ = \ N[2 Pi]}, \n\t\t\t\t\t\t
Map[{Cos[tp*#], Sin[tp*#]}&, orb]]\)], "Input",
InitializationCell->True]
}, Open ]],
Cell[CellGroupData[{
Cell["Plots for Vector Orbits", "Subsubsection",
InitializationCell->True],
Cell["\<\
To plot a vector orbit, just project out two of the coordinates \
(the first two by default) and plot some dots.\
\>", "Text",
InitializationCell->True],
Cell[BoxData[{
RowBox[{
\(ProjectionTransformation[orb_?VectorOrbitQ, indices_: {1, 2}]\ := \
\nMap[Part[#, indices]&, orb[]]\), "\n"}],
StyleBox[
\(PhasePlot[orb_, opts___]\ := \n\t
GenericPlot[orb, opts, Transformation -> ProjectionTransformation, \n
\t\t\t\tPointSize -> 0.007, AspectRatio -> 1, \n\t\t\t\t
Lines -> False]\),
InitializationCell->True,
AspectRatioFixed->True]}], "Input",
InitializationCell->True]
}, Closed]],
Cell[CellGroupData[{
Cell["Bifurcation Plot", "Subsubsection",
InitializationCell->True],
Cell[CellGroupData[{
Cell[BoxData[{
\( (*\
A\ function\ for\ plotting\ attractor\ bifurcation\ \n\tdiagrams\ for\
one - parameter\ families\ of\ maps\ *) \n\t\n
Options[BifurcationPlot]\ =
\ {\n\tParameterPoints -> 100, \n\tOrbitPoints -> 100, \n\t
Delay -> 100, \n\tClearOut -> True, \n\tRoundDigits -> 8, \n\t
PlotRange -> All, \n\tTransformation -> IdentityTransformation}; \n\t
\nbpoptpattern\ = \
\((ParameterPoints -> _)\) | \((OrbitPoints -> _)\) | \n\t
\((Delay -> _)\) | \((ClearOut -> _)\) | \((RoundDigits -> _)\) |
\((Transformation -> _)\); \n\t\n
BifurcationPlot[f_DS, ic_, pr_, opts___]\ := \n\t
Module[{showopts, pp, op, de, co, rd, plr, tr, fc, orb, Rnd, Clp,
ptlist = {}}, \n\t\t
pp\ = \ \(ParameterPoints /. {opts}\) /. Options[BifurcationPlot]; \n
\t\top\ = \ \(OrbitPoints /. {opts}\) /. Options[BifurcationPlot]; \n
\t\tde\ = \ \(Delay /. {opts}\) /. Options[BifurcationPlot]; \n\t\t
co\ = \ \(ClearOut /. {opts}\) /. Options[BifurcationPlot]; \n\t\t
rd\ = \ \((\(RoundDigits /. {opts}\) /. Options[BifurcationPlot])
\)*1.0; \n\t\t
plr\ = \ \(PlotRange /. {opts}\) /. Options[BifurcationPlot]; \n\t\t
tr\ = \ \(Transformation /. {opts}\) /. Options[BifurcationPlot]; \n
\t\tshowopts\ = \
DeleteCases[{opts}, bpoptpattern]~Join~\n
\t\t\t\t\t\t{Axes -> True, AspectRatio -> 1}; \n\t\t
If[co, Unprotect[Out]; Clear[Out]; Protect[Out]]; \n
\t\t (*fc\ = \ Compile[{x, c}, Evaluate[f[x, c]]]; *) \n\t\t
Rnd\ = \ Function[x, Round[10^rd*x]*10^\((\(-rd\))\)]; \n\t\t
Clp\ = \
If[ListQ[plr], \n\t\t\t
Function[l,
Select[l, \((plr[\([1]\)] <= #\ && \ # <= plr[\([2]\)])\)&]],
\n\t\t\tFunction[l, l]]; \n\t\t
For[c = pr[\([1]\)], c <= pr[\([2]\)],
c = c + \((pr[\([2]\)] - pr[\([1]\)])\)/pp, \n\t\t\t
orb\ = \ tr[f^{de, de + op}[ic[c]]]; \n\t\t\t
orb\ = \ Clp[Union[Rnd[orb]]]; \n\t\t\t
ptlist\ = \
ptlist~Join~\n\t\t\t\tTranspose[{Table[c, {Length[orb]}], orb}];
\n\t\t\t]; \n\t\t (*Print[Length[ptlist]]; *) \n\t\t
Show[Graphics[{PointSize[0.001], \n\t\t\t\tMap[Point, ptlist]}],
showopts]; \n\t\t]\n\t\t\t\),
\(BifurcationPlot::usage\ = \
"\300,OrbitPoints->200,Delay->100, \
PlotRange->{0,1}] will produce an orbit diagram with the parameter c ranging \
over the interval {0,4} in 300 steps. For each step an orbit with \
initialcondition ic[c] is computed. The first 100 points are not plotted and \
the next 200 are. Only points lying in the range 0<=x<=1 are plotted. This \
function is memory intensive.\>"\n\t\t\)}], "Input",
InitializationCell->True],
Cell[BoxData[
\("BifurcationPlot[f,ic,pr,opts] produces an orbit diagram for a one \
parameter family of maps depending on a parameter, c. Here f should be a \
dynamical system object whose underlying formula contains the variable c. \
ic[c] is a function of one variable specifying the initial point to use in \
iterating. pr is the range of c values to use. There are several options. \
For example, if f has been defined as the logistic map f= DS[f1] where f1[x_] \
:= c*x*(1-x) and ic[c_] := 0.5 the command: \
BifurcationPlot[f,ic,{0,4},ParameterPoints->300,OrbitPoints->200,Delay->100, \
PlotRange->{0,1}] will produce an orbit diagram with the parameter c ranging \
over the interval {0,4} in 300 steps. For each step an orbit with \
initialcondition ic[c] is computed. The first 100 points are not plotted and \
the next 200 are. Only points lying in the range 0<=x<=1 are plotted. This \
function is memory intensive."\)], "Output"]
}, Open ]]
}, Open ]]
}, Closed]],
Cell[CellGroupData[{
Cell["Iterated Function Systems", "Subsection"],
Cell["\<\
Here we define an IFS object representing an iterated function \
system. To apply and IFS to a point we simply apply each function to the \
point and return the list of answers. To iterate this we need to define what \
it means to apply and IFS to a list of vectors. Finally ,we define how to \
use an IFS to transform a line or a list of lines\
\>", "Text"],
Cell[BoxData[{
\(SwitchHead[f_, h_]\ := \ f /. f[\([0]\)] -> h\n\),
\(f_IFS[x : {__Line}]\ := Flatten[Map[f, x], 1]\),
\(f_IFS[x : {__List}]\ := \ Flatten[Map[f, x], 1]\),
\(f_IFS[x_List /; VectorQ[x]]\ := \ Map[#[x]&, SwitchHead[f, List]]\),
\(f_IFS[Line[l_]]\ := \ Map[Line, Transpose[Map[f, l]]]\n\),
\(IFS /: \ f_IFS^n_Integer[x_]\ := \ Nest[f, x, n]\n\),
\(IFS /: \ f_IFS^n_Integer[x : {__Line}]\ := \ Nest[f, x, n]\)}], "Input",\
InitializationCell->True]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell["Lab 1 --- Dynamics of One-Dimensional Maps", "Section"],
Cell[TextData[{
"In this lab you will use ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" to compute and plot orbits of one-dimensional dynamical systems. The \
goal will be to illustrate some of the features described in class, such as \
attracting fixed points, periodic points and chaos. To make things easier, \
you will use a package of ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" programs that I have written which are contained in the ",
StyleBox["Initialization",
FontWeight->"Bold"],
" section of this notebook. To get ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" started and to enter all these programs into memory, go to the ",
StyleBox["Kernel",
FontWeight->"Bold"],
" menu and inside the ",
StyleBox["Evaluation",
FontWeight->"Bold"],
" submenu you should select the item called ",
StyleBox["Evaluate Initialization",
FontWeight->"Bold"],
". It is not necessary to look in detail at these programs now. You will \
see how to use them as we go along."
}], "Text"],
Cell[TextData[{
"This lab is organized into several sections whose titles are visible \
below. To begin the first section just double-click on the bracket to the \
far right of the title ",
StyleBox["Dynamical systems and orbits. ",
FontWeight->"Bold"],
" The bracket will expand to reveal lots of ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" cells containing explanatory text and commands. Just work your way \
through the section reading the explanations and executing the commands. \
When you finish a section you can close it up again by double clicking on the \
expanded bracket at the extreme right. Then go on to the next section. \
Enjoy !"
}], "Text"],
Cell[CellGroupData[{
Cell["Dynamical systems and orbits", "Subsection"],
Cell[TextData[{
"To use the dynamical systems package, you need to enter a dynamical \
system. The cell below defines the logistic function ",
StyleBox["f1[x]", "Input"],
" and then uses it to define a \"dynamical system object\" called ",
StyleBox["f", "Input"],
". This extra step lets the system know that the package I wrote should be \
applied to compute iterates, orbits, etc. To execute a ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" command, click the mouse somewhere in the \"cell\" containing the command \
and do a ",
StyleBox["Shift-Return",
FontWeight->"Bold"],
", that is, hold down the ",
StyleBox["Shift",
FontWeight->"Bold"],
" key and then hit the ",
StyleBox["Return",
FontWeight->"Bold"],
" key."
}], "Text"],
Cell[BoxData[{
\(f1[x_]\ := \ a*x*\((1 - x)\)\),
\(f\ = \ DS[f1]\)}], "Input"],
Cell[TextData[{
"To compute, say, the 12-th iterate of ",
StyleBox["f", "Input"],
" applied to the initial state ",
StyleBox["x0 = 0.1", "Input"],
", we can use the standard exponential notation below. First we will set \
the value of ",
StyleBox["a", "Input"],
" (the semicolon prevents the a value from being printed)."
}], "Text"],
Cell[BoxData[
\(a\ = \ 2.9; \nf^12[0.1]\)], "Input"],
Cell["\<\
You can edit previous commands and reexecute them. Go back to the \
last cell and compute some larger iterate of the same point, say the 1200-th \
one or the 12000-th one. \
\>", "Text"],
Cell["\<\
To compute an entire orbit, instead of just one point on it, just \
enclose the exponent in list brackets as shown below.\
\>", "Text"],
Cell[BoxData[
\(f^{100}[0.1]\)], "Input"],
Cell[TextData[{
"The return value is an \"orbit object\", which is just a list of points \
contained inside an ",
StyleBox["Orbit[ ]", "Input"],
" symbol. The points represent the orbit ",
StyleBox["x0,x1,...", "Input"],
". Once you have computed an orbit object you can use a simple bracket \
notation to print out entries of it. But you need to give the orbit a name \
first."
}], "Text"],
Cell[BoxData[
\(\(orb\ = \ f^{100}[0.1]; \)\)], "Input"],
Cell[TextData[{
"Now you can use, for example, orb[27] to read off ",
StyleBox["x27", "Input"],
". You can also put a range of indices inside to get a sublist or use an \
empty bracket to print the orbit as a plain list."
}], "Text"],
Cell[BoxData[{
\(orb[27]\),
\(orb[{0, 3}]\),
\(orb[]\)}], "Input"],
Cell[TextData[{
"It seems that this orbit is attracted to a fixed point of ",
StyleBox["f ", "Input"],
"located at ",
StyleBox["p=0.65517. ", "Input"],
"This is even clearer when we plot the orbit. There are many built-in \
plotting functions in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" and I have included several more in the package. For example, the cell \
below produces a time-series plot."
}], "Text"],
Cell[BoxData[
\(\(TimeSeriesPlot[orb, PlotRange -> All]; \)\)], "Input"],
Cell["\<\
You can use the mouse to resize graphic images for clearer viewing.\
\
\>", "Text"],
Cell["Another interesting plot is shown below.", "Text"],
Cell[BoxData[
\(\(CobwebPlot[f, orb, {0, 1}]; \)\)], "Input"]
}, Open ]],
Cell[CellGroupData[{
Cell["Periodic Attractors", "Subsection"],
Cell[TextData[{
"It is easy to use simple iteration and plotting to search for attracting \
periodic orbits. In the cell below, we plot the time series for an orbit of \
the logistic map with ",
StyleBox["a=2.5", "Input"],
"."
}], "Text"],
Cell[BoxData[
\(a\ = \ 2.5; \norb\ = \ f^{100}[0.1]; \n
TimeSeriesPlot[orb, PlotRange -> All]; \)], "Input"],
Cell[TextData[{
"The commands in the previous cell can be used to explore the changing \
dynamics of the logistic map as the parameter ",
StyleBox["a", "Input"],
" is changed. Go back and change the value of ",
StyleBox["a", "Input"],
" and see what happens. Try the values ",
StyleBox["a=2.9,3.1,3.4,3.5,3.55,3.83,3.845.", "Input"],
"In each case you should be able to identify the period of the attractor. \
"
}], "Text"],
Cell[TextData[{
"Cobweb plots are also interesting. Try several ",
StyleBox["a ", "Input"],
"values in the command below."
}], "Text"],
Cell[BoxData[
\(a\ = \ 2.5; \norb\ = \ f^{100}[0.1]; \nCobwebPlot[f, orb, {0, 1}];
\)], "Input"],
Cell[TextData[{
"For high period attractors it is sometimes hard to see the attracting \
periodic orbit because the \"transient\" part of the orbit messes up the \
picture. A nice way to avoid this is to ignore the first part of the orbit \
to give it time to be attracted to the periodic orbit. Another variation on \
the iteration commands makes this easy. The command below plots iterates ",
StyleBox["x1000,..x1050", "Input"],
" along the orbit, ignoring the first 1000 iterates."
}], "Text"],
Cell[BoxData[{
\(a\ = \ 3.1; \norb\ = \ f^{1000, 1050}[0.1]\),
\(\(CobwebPlot[f, orb, {0, 1}]; \)\)}], "Input"],
Cell[TextData[{
"Have a look at some of the other ",
StyleBox["a", "Input"],
" values using this technique."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Periodic Orbits via the equation f^k[x] = x", "Subsection"],
Cell[TextData[{
"Of course iteration will only find attracting periodic orbits. To find \
repellers you have to solve the equation ",
StyleBox["f^k[x]=x", "Input"],
". You can find a formula for ",
StyleBox["f^k[x]", "Input"],
" by just using the variable ",
StyleBox["x", "Input"],
" as your initial condition. The first line below unsets ",
StyleBox["a", "Input"],
" so it is a variable again."
}], "Text"],
Cell[BoxData[{
\(a =. \),
\(f[x]\),
\(f^2[x]\),
\(f^3[x]\)}], "Input"],
Cell[TextData[{
"If ",
StyleBox["f[x]", "Input"],
" is a polynomial, as it is for the logistic map, you can try to solve ",
StyleBox["f^k[x]=x", "Input"],
"by factorization , at least for small values of ",
StyleBox["k", "Input"],
"."
}], "Text"],
Cell[BoxData[{
\(Factor[f^2[x] - x]\),
\(Factor[f^3[x] - x]\),
\(Factor[f^4[x] - x]\)}], "Input"],
Cell[TextData[{
"Unfortunately, the polynomials are still pretty big. For specific values \
of ",
StyleBox["a", "Input"],
" you can use ",
StyleBox["Mathematica",
FontSlant->"Italic"],
"'s NSolve function to look for solutions."
}], "Text"],
Cell[BoxData[{
\(a\ = \ 3.5; \nNSolve[f^2[x] - x == 0, x]\),
\(NSolve[f^3[x] - x == 0, x]\),
\(NSolve[f^4[x] - x == 0, x]\)}], "Input"],
Cell[TextData[{
"These results show where the period 1, 2 and 4 orbits are for this value \
of ",
StyleBox["a", "Input"],
". Apparently there are no period 3 orbits since the roots of ",
StyleBox["f^3[x]-x", "Input"],
" other than the fixed points are all complex numbers. Try ",
StyleBox["a=3.83", "Input"],
", which seemed to have a period 3 attractor in the plots from the \
previous section. You will actually find two period-3 orbits -- one is a \
repeller which was not previously detected. Try ",
StyleBox["a=4", "Input"],
" as well."
}], "Text"],
Cell[TextData[{
"Another way to study periodic points is to plot the graphs of the \
iterates and look for places where they cross the diagonal. For this the \
built-in plotting function ",
StyleBox["Plot", "Input"],
" can be used. You can plot several iterates on the same picture. \
Plotting the 0-th iterate is one way to draw the \"diagonal\"."
}], "Text"],
Cell[BoxData[
\(a = 3.83; \nPlot[{f^0[x], f[x], f^2[x]}, {x, 0, 1}, AspectRatio -> 1];
\nPlot[{f^0[x], f[x], f^3[x]}, {x, 0, 1}, AspectRatio -> 1]; \)], "Input"],
Cell[TextData[{
"The graph of ",
StyleBox["f^2[x]", "Input"],
" intersects the diagonal four times, at the fixed points and at two other \
points. The latter form an orbit of period 2 (a repeller). On the other \
hand, the graph of ",
StyleBox["f^3[x]", "Input"],
" crosses the diagonal only at the fixed points. Now compare ",
StyleBox["a=3.83.", "Input"]
}], "Text"],
Cell[TextData[{
"The abundance of periodic points for ",
StyleBox["a=4", "Input"],
" is visible from the graphs of its iterates. For example, the fifth \
iterate has a graph consisting of 32 monotonic pieces, each crossing the \
diagonal."
}], "Text"],
Cell[BoxData[
\(a = 4.0; \nPlot[{f^0[x], f[x], f^5[x]}, {x, 0, 1}, AspectRatio -> 1];
\)], "Input"],
Cell["Try looking at some higher iterates of this map.", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Finding Periodic Orbits and Multipliers Numerically", "Subsection"],
Cell[TextData[{
"Another way to look for periodic points is to use a special function \
included with this package called ",
StyleBox["SeekPeriodicOrbit", "Input"],
". Instead of trying to find all possible solutions to the ",
StyleBox["f^k[x]=x", "Input"],
", it looks for one solution close to a given point which you supply. \
This is sometimes feasible when a more general solution is not, say when ",
StyleBox["f[x]", "Input"],
" is not a polynomial or when ",
StyleBox["k", "Input"],
" is a bit larger. "
}], "Text"],
Cell[TextData[{
"We will try to find both of the period three orbits for the logistic map \
with ",
StyleBox["a=3.83",
FontWeight->"Bold"],
". Here is the graph of the third iterate and a blowup of that graph near \
the center"
}], "Text"],
Cell[BoxData[
\(a = 3.83; \nPlot[{f^0[x], f^3[x]}, {x, 0, 1}, AspectRatio -> 1]; \n
Plot[{f^0[x], f^3[x]}, {x, 0.4, 0.6}, AspectRatio -> 1]; \)], "Input"],
Cell["\<\
From these graphs you can see some good guesses for where to search \
for the period 3 points. In the blown up picture, the point on the left is \
on the attracting period 3 orbit and the point on the right is on the \
repeller. We will use 0.5 as a starting point for a search.\
\>", "Text"],
Cell[BoxData[
\(orb1 = \ SeekPeriodicOrbit[f, 3, 0.5]\)], "Input"],
Cell["\<\
This is apparently the attractor. To find the repeller we search a \
bit to the right starting at 0.55.\
\>", "Text"],
Cell[BoxData[
\(orb2 = \ SeekPeriodicOrbit[f, 3, 0.55]\)], "Input"],
Cell[TextData[{
"Thus, using two different initial guesses leads to two different period 3 \
orbits which must be the ones we are after. You can test them for stability \
using a package function called ",
StyleBox["Multiplier.", "Input"]
}], "Text"],
Cell[BoxData[{
\(Multiplier[f, orb1]\),
\(Multiplier[f, orb2]\)}], "Input"],
Cell[TextData[{
"This confirms that ",
StyleBox["orb1", "Input"],
" is an attractor and ",
StyleBox["orb2", "Input"],
" is a repeller. Finally, we can display them both in a cobweb plot using \
red for the attractor and blue for the repeller."
}], "Text"],
Cell[BoxData[
\(\(CobwebPlot[f, {orb1, orb2}, {0, 1}, LineColor -> {Red, Blue}]; \)\)],
"Input"],
Cell["\<\
This demonstrates another feature of the plotting functions in this \
package. You can plot several orbits just by supplying a list of orbits as \
the argument instead of a single orbit. Also there are various options which \
you can manipulate. To get more information about a command you can \
type:\
\>", "Text"],
Cell[BoxData[
\(\(?CobwebPlot\)\)], "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell["Sensitive Dependence on Initial Conditions", "Subsection"],
Cell[TextData[{
"In this section we will have a look at the dynamics of some dynamical \
systems which exhibit chaos. A first sign of chaos is just the absence of a \
periodic attractor. Here is a cobweb plot of the logistic maps with ",
StyleBox["a=4.0", "Input"],
"."
}], "Text"],
Cell[BoxData[
\(a = 4.0; \norb\ = \ f^{1000, 1200}[0.1]; \nCobwebPlot[f, orb, {0, 1}];
\)], "Input"],
Cell[TextData[{
"In spite of the fact that we have skipped the first 1000 points on the \
orbit to give it a chance to converge to an attractor, there is no evidence \
of convergence -- the orbit seems to wander all over the interval [0,1]. Go \
back and try a similar plot for ",
StyleBox["a=3.6785", "Input"],
" and ",
StyleBox["a=3.927737", "Input"],
". You will find similar behavior but on a slightly smaller interval. \
Finally, have a look at ",
StyleBox["a=3.57345", "Input"],
" which seems to hop back and forth between two intervals of chaos. Just \
from the pictures, it is not clear that these maps are really chaotic. \
Perhaps there is an attracting periodic orbit with a very large period so the \
resulting cobweb plot looks complicated. For special values like ",
StyleBox["a=4", "Input"],
" we can prove that no attractors exist. Also there is a theorem which \
says that there are many other such ",
StyleBox["a", "Input"],
" values (a so-called \"positive measure set\" of them) without specifying \
exactly which ones."
}], "Text"],
Cell["\<\
It is easy to demonstrate sensitive dependence on initial \
conditions by simply comparing orbits with nearby initial conditions. We \
will compute two such orbits and plot them together on the same graph.\
\>",
"Text"],
Cell[BoxData[
\(a = 4.0; \norb1\ = \ f^{50}[0.1]; \norb2\t = \ f^{50}[0.10001]; \n
TimeSeriesPlot[{orb1, orb2}, LineColor -> {Red, Green}]; \)], "Input"],
Cell[TextData[{
"You can see that after only about 12 or 13 iterations, the two orbits \
separate. After that they are more or less indepedent, sometimes close and \
sometimes far apart. Here is a cobweb plot demonstrating the same thing. \
Perhaps we just need to take the initial conditions closer. Go back and try \
changing the initial condition of ",
StyleBox["orb2", "Input"],
" to ",
StyleBox["0.1000001.", "Input"],
" Keep adding zeroes and see what happens (machine precision is 15 places \
or so). "
}], "Text"],
Cell["\<\
If you try the same experiment for a parameter value with an \
attractor, the results are quite different.\
\>", "Text"],
Cell[BoxData[
\(a = 3.845; \norb1 = \ f^{70}[0.1]; \norb2 = \ f^{70}[0.10001]; \n
TimeSeriesPlot[{orb1, orb2}, LineColor -> {Red, Green}]; \)], "Input"],
Cell[TextData[{
"After a while the orbits get attracted to the period 3 attractor. \
However, at any given iteration they are far apart because they are \"out of \
phase\". In other words, one orbit has every third iterate ",
StyleBox["x0,x3,x6,... ", "Input"],
"converging to one point, say ",
StyleBox["p0", "Input"],
", on the period-3 orbit while the other orbit has ",
StyleBox["x0,x3,x6,... ", "Input"],
"converging to ",
StyleBox["p1", "Input"],
". However, this can be remedied by taking the two initial points closer. \
Go back up to the previous cell and try the initial condion ",
StyleBox["0.100001", "Input"],
". Now the long term behaviors are the same. "
}], "Text"],
Cell[TextData[{
"We can demonstrate exponential divergence and the idea of Lyapunov \
exponents by looking quantitatively at the differences bewteen two nearby \
orbits. The plot below looks at the logaritbms of the absolute values of \
the differences bewteen two nearby orbits for the chaotic logistic map with \
",
StyleBox["a=4", "Input"],
"."
}], "Text"],
Cell[BoxData[{
\(a = 4.0; \norb1 = \ f^{50}[0.1]; \norb2 = \ f^{50}[0.10001]; \n
diff\ = \ Abs[orb1[] - orb2[]]\),
\(\(ListPlot[Log[diff]]; \)\)}], "Input"],
Cell["\<\
Note that the first part of this plot looks roughly like a straight \
line. We can estimate its slope as follows:\
\>", "Text"],
Cell[BoxData[
\(slope\ = \ \((Log[diff[\([11]\)]] - Log[diff[\([1]\)]])\)/10\)],
"Input"],
Cell[TextData[{
"Compare this with the Lyapunov exponent. This is approximated by the \
logarithm of the absolute value of the product of the derivatives along the \
orbit divided by the length of the orbit. The product of the derivatives \
along the orbit is just what the function ",
StyleBox["Multiplier", "Input"],
" is designed to compute. It can also be applied to nonperiodic orbits."
}], "Text"],
Cell[BoxData[{
\(mu\ = \ Multiplier[f, orb1]\),
\(lyap\ = \ Log[\ Abs[mu]]/50\)}], "Input"],
Cell[TextData[{
"Note that the lyapunov exponent is almost equal to the slope of the graph \
we found above. This shows how the Lyapunov exponent along the orbit ",
StyleBox["orb1", "Input"],
" measures the exponential rate at which nearby orbits (like ",
StyleBox["orb2", "Input"],
") diverge from it."
}], "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Some More Things to Try", "Subsection"],
Cell["\<\
If you have more time you could try some of the following \
experiments. For these you will have to type the commands yourself.\
\>",
"Text"],
Cell[TextData[{
"1. Define a new dynamical system ",
StyleBox["g = DS[Cos]", "Input"],
" based on the cosine function. Plot graphs some of its iterates and \
explain what you find."
}], "Text"],
Cell[TextData[{
"2. Using ",
StyleBox["SeekPeriodicOrbit",
FontWeight->"Bold"],
", find several periodic points of the map 4x(1-x). Compute their \
multipliers using ",
StyleBox["Multiplier. ",
FontWeight->"Bold"],
" The results are surprisingly simple."
}], "Text"],
Cell[TextData[{
"3. Using iteration and plotting, convince yourself that the parameter \
value where the logistic map ax(1-x) first has an attracting period 2 orbit \
is ",
StyleBox["a=3", "Input"],
". Then by increasing ",
StyleBox["a", "Input"],
", find a slightly larger value where a period 4 attractor appears. \
Continue to periods 8 and 16, if you can. After finding these values ( call \
them a2, a4, a8, a16 ) have a look at the ratios (a4-a2)/(a8-a4) and \
(a8-a4)/(a16-a8). These approximate Feigenbaum's number. The pattern \
continues."
}], "Text"]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell["Lab 2 --- Two-Dimensional Dynamics and Bifurcations", "Section"],
Cell["\<\
In this lab we will look at some dynamical phenomena peculiar to \
two dimensions and also plot bifurcations for one and two dimensional \
systems.\
\>", "Text"],
Cell[CellGroupData[{
Cell["Orbits of Linear Maps", "Subsection"],
Cell[TextData[{
"We will plot some orbits of simple linear maps of the plane. In ",
StyleBox["Mathematica",
FontSlant->"Italic"],
", vectors are represented as lists of numbers enclosed in set brackets. \
For example, {1,2} is a vector in the plane. A matrix is given as a list of \
vectors which constitute the rows of the matrix. The most general possible \
2x2 matrix would be {{a,b},{c,d}}. The built-in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" function ",
StyleBox["MatrixForm", "Input"],
" just prints the matrix in a nice way."
}], "Text"],
Cell[BoxData[
\(A\ := \ {{a, b}, {c, d}}; \nMatrixForm[A]\)], "Input"],
Cell[TextData[{
"There are many ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" functions to analyze matrices. After setting specific values for \
a,b,c,d, a few useful commands are demonstated below. The chosen values \
represent the Fibonacci matrix."
}], "Text"],
Cell[BoxData[{
\(a\ = \ 1; \ b = 1; \nc = 1; \ d = 0; \n"\" MatrixForm[A]\),
\("\"\ MatrixForm[Inverse[A]]\),
\(Eigenvalues[A]\),
\(Eigenvectors[A]\)}], "Input"],
Cell["\<\
To view a matrix as a dynamical system, we define a linear function \
which takes as input a vector {x,y} and produces as output, the vector \
obtained by multiplying by the matrix A. Then we define a dynamical system \
object using this function. As before, this allows the program to recognize \
the commands for iteration and orbit computation.\
\>", "Text"],
Cell[BoxData[{
\(l1[{x_, y_}]\ := \ A . {x, y}\),
\(L\ = \ DS[l1]\n\),
\(L[{x, y}]\),
\(L[{1, 5}]\)}], "Input"],
Cell["\<\
Just as for one-dimensional maps, you can use the standard exponent \
notation to find an iterate of this map applied to an initial vector x0. To \
compute an orbit just put a range of values in the exponent. The result is \
then an Orbit object containing a list of vectors which are the individual \
points of the orbit. In this case, the coordinates of these points are the \
Fibonacci numbers.\
\>", "Text"],
Cell[BoxData[{
\(x0\ = \ {1, 0}; \nL^100[x0]\),
\(orb\ = \ L^{0, 50}[x0]\)}], "Input"],
Cell[TextData[{
"Next we plot the orbit using a vector plotting function called ",
StyleBox["PhasePlot", "Input"],
". Since the origin is a saddle point, the orbit just heads to infinity in \
the direction of the unstable subspace."
}], "Text"],
Cell[BoxData[
\(\(PhasePlot[orb, PlotRange -> All, PointSize -> 0.02, Lines -> True];
\)\)], "Input"],
Cell["\<\
To make an interesting plot of a saddle point, you need to start \
close to the stable subspace. The eigenvectors were computed above and the \
stable eigenvector is {(1-Sqrt[5])/2,1}. The point x0 below is \
approximately, but not exactly, equal to this. The orbit hops back and forth \
near the stable subspace (it's a flip saddle) and then takes off along the \
unstable direction.\
\>", "Text"],
Cell[BoxData[{
\(x0\ = \ {\(-6.18\), 10}; \norb\ = \ L^{0, 25}[x0]\),
\(\(PhasePlot[orb, PlotRange -> All, PointSize -> 0.02, Lines -> True];
\)\)}], "Input"],
Cell[TextData[{
"It is more interesting to look at a case where the origin is an attractor. \
Here is a spiral sink, i.e., the eigenvalues are complex with norm less than \
1. By changing the numbers ",
StyleBox["r ", "Input"],
"and ",
StyleBox["theta ", "Input"],
"below you can produce some really nice pictures. "
}], "Text"],
Cell[BoxData[{
\(\nr\ = \ 0.95; \ theta\ = 2*Pi/7; \na\ = \ r*Cos[theta];
b\ = \ \(-r\)*Sin[theta]; \nc\ = \ r*Sin[theta]; d\ = \ r*Cos[theta];
\nEigenvalues[A]\),
\(Abs[Eigenvalues[A]]\),
\(x0\ = \ {1, 0}; \norb\ = \ L^{0, 50}[x0]; \n
PhasePlot[orb, PlotRange -> All, PointSize -> 0.02, Lines -> True]; \)}],
"Input"],
Cell["\<\
Here is a sink with real eigenvalues. Several different orbits \
are shown to illustrate how they all approach the origin tangent to the \
weaker eigenspace (in this case, the x-axis). Try switching the values of a \
and d so that the other eigenvector becomes the weaker one. \
\>", "Text"],
Cell[BoxData[{
\(a = 0.8; b = 0; \nc = 0; d = 0.7; \nEigenvalues[A]\),
\(x1\ = \ {1, \(-1\)}; orb1 = L^{0, 20}[x1]; \nx2\ = \ {1, \(-0.1\)};
orb2 = L^{0, 20}[x2]; \nx3\ = \ {1, 1}; orb3 = L^{0, 20}[x3]; \n
x4\ = \ {\(-1\), 1}; orb4 = L^{0, 20}[x4]; \nx5\ = \ {\(-1\), 0.1};
orb5 = L^{0, 20}[x5]; \nx6\ = \ {\(-1\), \(-1\)}; orb6 = L^{0, 20}[x6];
\n\nPhasePlot[{orb1, orb2, orb3, orb4, orb5, orb6},
PlotRange -> {{\(-1\), 1}, {\(-1\), 1}}, PointSize -> 0.02,
Lines -> True]; \n\)}], "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell["Orbits the Henon Map", "Subsection"],
Cell["\<\
We begin by entering the formulas for the Henon map and declaring \
it to be a dynamical system object. A mapping of the plane such as the \
Henon map takes a vector as input and produces another vector as output. It \
can be entered as shown below.\
\>", "Text"],
Cell[BoxData[{
\(h1[{x_, y_}]\ := {a - x^2 + b*y, x}\),
\(henon\ = \ DS[h1]\)}], "Input"],
Cell["\<\
Next we set the values of the parameters, define an initial point, \
compute an orbit and plot it. Since we will be plotting long orbits we will \
plot them with small dots not connected by lines. In this case we get the \
Henon attractor.\
\>", "Text"],
Cell[BoxData[
\(a = 1.4; \ b\ = \ 0.3; \nx0\ = \ {1, 1}; \n
orb\ = \ henon^{0, 5000}[x0]; \n
PhasePlot[orb, PlotRange -> {{\(-1.6\), 2}, {\(-1.6\), 2}}]; \)], "Input"],
Cell["\<\
To help understand this attractor we will plot a sequence of \
pictures where the value of the parameter a is gradually increased. At each \
stage, we plot the omega limit set of the orbit with initial conditions \
{1,1}. This is done by skipping the first 100 points along the orbit to \
allow time to converge to whatever attractor may be present.
\tTo make the simpler attractors visible, we can use larger dots in the plot. \
Starting from the value a = 0.3 for which there is an attracting fixed \
point, we progress through a sequence of period-doubling bifurcations \
arriving at a strange attractor with several parts. There are interludes \
where a periodic attractor returns eventually leading to the attractor found \
above. Try stepping through the a values suggested below.\
\>", "Text"],
Cell[BoxData[
\(a = 0.3;
\ \ (*\ try\ a\ = \ 0.4, \ 0.9, \ 0.95, \ 1.05, \ 1.055, \ 1.07, \ 1.1,
\ 1.25, \ 1.27, \ 1.3, \ 1.4*) \ \ \nb\ = \ 0.3; \nx0\ = \ {1, 1}; \n
orb\ = \ henon^{100, 2000}[x0]; \n
PhasePlot[orb, PlotRange -> {{\(-1.6\), 2}, {\(-1.6\), 2}},
PointSize -> 0.02]; \)], "Input"],
Cell["\<\
For a close-up look at the attractor, you can plot a very long \
orbit and plot blow-ups of it by changing the PlotRange option.\
\>", "Text"],
Cell[BoxData[
\(a = 1.4; \ b\ = \ 0.3; \nx0\ = \ {1, 1}; \n
orb\ = \ henon^{0, 20000}[x0]; \n
PhasePlot[orb, PlotRange -> {{0.3, 0.9}, {0.6, 1.2}}]; \n
PhasePlot[orb, PlotRange -> {{0.4, 0.6}, {0.9, 1.1}}]; \)], "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell["Bifurcations and Orbit Diagrams", "Subsection"],
Cell["\<\
In this section we will have a look at the amazing details of the \
orbit diagram for the logistic map and also construct a similar diagram for \
the two-dimensional Henon map. First we define the logistic dynamical \
system. The bifurcation plotting software is set up to vary a parameter \
called c so we will use that letter in the formula.\
\>", "Text"],
Cell[BoxData[{
\(f1[x_]\ := \ c*x*\((1 - x)\)\),
\(f\ = \ DS[f1]\)}], "Input"],
Cell["\<\
To construct an orbit diagram we need to start with some initial \
condition. In general this could vary with c so we define an initial \
condition function ic[c] which gives the desired starting point. In this \
case it is a constant function at the critical point x0 = 0.5.\
\>", "Text"],
Cell[BoxData[
\(ic[c_]\ := \ 0.5\)], "Input"],
Cell["\<\
Now the function bifurcation plot can be used to draw the orbit \
diagram. Here is some information about this function, whose actual \
definition is contained in the initialization section of this notebook.\
\>",
"Text"],
Cell[BoxData[
\(\(?BifurcationPlot\)\)], "Input"],
Cell[BoxData[
\(BifurcationPlot[f, ic, {0, 4}, ParameterPoints -> 300,
OrbitPoints -> 200, Delay -> 100, PlotRange -> {0, 1}]\)], "Input"],
Cell["\<\
Of course the most interesting part of the diagram is from c=2 to \
c=4. You can zoom in on that part by changing the range for c.\
\>", "Text"],
Cell[BoxData[
\(BifurcationPlot[f, ic, {2, 4}, ParameterPoints -> 300,
OrbitPoints -> 200, Delay -> 100, PlotRange -> {0, 1}]\)], "Input"],
Cell["\<\
It is endlessly amusing to blow up small regions of the diagram and \
to ponder the incredible complexity contained in it. For example here is a \
tiny window where a period 5 orbit is born in a saddle-node bifurcation and \
then proceeds to period double producing attractors of period 10, 20 , 40 , \
etc.\
\>", "Text"],
Cell[BoxData[
\(BifurcationPlot[f, ic, {3.73, 3.75}, ParameterPoints -> 300,
OrbitPoints -> 200, Delay -> 100, PlotRange -> {0, 1}]\)], "Input"],
Cell[TextData[{
"You can blow up just one of the five period-doubling regions by changing \
the ",
StyleBox["PlotRange", "Input"],
" as well as the range for c. Here a blow-up of the part of the last \
diagram near the center revealing further smaller windows. For deep blow-ups \
it is sometimes good to increase the Delay to give more time to converge to \
the attractor and the number of OrbitPoints, since not all of them will \
appear in the chosen Plot Range."
}], "Text"],
Cell[BoxData[
\(BifurcationPlot[f, ic, {3.742, 3.745}, ParameterPoints -> 300,
OrbitPoints -> 400, Delay -> 200, PlotRange -> {0.45, 0.55}]\)], "Input"],
Cell["\<\
I suggest that you go back to the previous cell and explore other \
regions of the orbit diagram by making changes to the c range and the Plot \
Range.\
\>", "Text"],
Cell["\<\
Next we will make a similar plot for a one-parameter family of \
Henon maps. Here is the family which was studied using simple forward \
iteration in the last section. Once again, we have called the parameter c.\
\
\>", "Text"],
Cell[BoxData[{
\(h1[{x_, y_}]\ := \ {c + 0.2*y\ - x^2, x}\),
\(h\ = \ DS[h1]\)}], "Input"],
Cell["\<\
Now the orbits for the Henon map consist of vectors. To get a 2D \
plot we have to plot c versus one or the other of the x, y coordinates (or \
some combination of them). Here is a function which projects a 2D vector \
onto its x coordinate and another which takes a vector orbit and produces its \
list of x coordinates\
\>", "Text"],
Cell[BoxData[{
\(xproj[{x_, y_}]\ := \ x\),
\(XProj[orb_Orbit]\ := \ Map[xproj, orb[]]\)}], "Input"],
Cell["\<\
This can be supplied as an option to the BifurcationPlot function \
to construct an orbit diagram for the Henon family. We also need to supply \
an initial conction somewhere in the basin of the attractor. The point {1,1} \
will do.\
\>", "Text"],
Cell[BoxData[
\(ic[c_]\ := \ {1.0, 1.0}\)], "Input"],
Cell[BoxData[
\(BifurcationPlot[h, ic, {0, 1.6}, ParameterPoints -> 300,
OrbitPoints -> 200, Delay -> 100, Transformation -> XProj]\)], "Input"],
Cell["\<\
By blowing-up parts of this picture, you can find subattractors of \
the Henon attractor with various periods.\
\>", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Stable and Unstable Manifolds", "Subsection"],
Cell["\<\
Next we will plot some stable and unstable manifolds for 2D \
diffeomorphisms. The first step is to find a saddle fixed or periodic point \
for the map. Of course we begin with the Henon map.\
\>", "Text"],
Cell[BoxData[{
\(h1[{x_, y_}]\ := {a - x^2 + b*y, x}\),
\(henon\ = \ DS[h1]\)}], "Input"],
Cell[TextData[{
"After setting some familiar values for the parameters we search for a \
fixed point using the function ",
StyleBox["SeekPeriodicOrbit.", "Input"],
" This function looks for a solution of the fixed point equations near a \
given initial point. If it finds one it returns an ",
StyleBox["Orbit", "Input"],
" object containing the point and its first iterate (so you can see it \
really is fixed)."
}], "Text"],
Cell[BoxData[
\(a = 0.8; \nb\ = \ 0.2; \n
orb\ = \ SeekPeriodicOrbit[henon, 1, {1, 1}]\)], "Input"],
Cell["\<\
You can test to see if it is a saddle point by computing its \
multipliers.\
\>", "Text"],
Cell[BoxData[
\(Multiplier[henon, orb]\)], "Input"],
Cell["\<\
Since one of them has absolute value bigger than one and one has \
absolute value less than one, it is indeed a saddle point. The first entry \
of the orbit is the actual fixed point which we will call X0.\
\>", "Text"],
Cell[BoxData[
\(X0\ = \ orb[0]\)], "Input"],
Cell[TextData[{
"Now we can use the magic function ",
StyleBox["UnstableManifold ", "Input"],
"to compute a piece of the unstable manifold near the fixed point."
}], "Text"],
Cell[BoxData[
\(\(?UnstableManifold\)\)], "Input"],
Cell[BoxData[
\(Wu\ = \ UnstableManifold[henon, X0, 1]\)], "Input"],
Cell["\<\
The result is a list of 40 points appoximately on the unstable \
manifold of the saddle, 20 on each side of the fixed point. We can plot \
this to see what is happening.\
\>", "Text"],
Cell[BoxData[
\(\(PhasePlot[Wu, Lines -> True, Points -> False, LineColor -> Blue,
PlotRange -> {{\(-1\), 2}, {\(-1\), 2}}]; \)\)], "Input"],
Cell[TextData[{
"This is only a small piece of the unstable manifold. To extend it just \
change the third argument of the ",
StyleBox["UnstableManifold ", "Input"],
"function to increase the number of iterations it uses to extend the \
manifolds. "
}], "Text"],
Cell[BoxData[
\(Wu\ = \ UnstableManifold[henon, X0, 3]; \n
PhasePlot[Wu, Lines -> True, Points -> False, LineColor -> Blue, \n
PlotRange -> {{\(-1\), 2}, {\(-1\), 2}}]; \)], "Input"],
Cell["\<\
Go back and gradually increase the number of iterations used up to \
about 10 to see the manifold grow.\
\>", "Text"],
Cell["\<\
To get a stable manifold, we just compute the unstable manifold for \
the inverse of the map. So first we have to input the formulas for the \
inverse of the Henon map.\
\>", "Text"],
Cell[BoxData[{
\(hinv1[{x_, y_}]\ := \ {y, \((x + y^2 - a)\)/b}\),
\(henoninv\ = \ DS[hinv1]\)}], "Input"],
Cell["\<\
The next command will compute pieces of both the stable and \
unstable manifolds and plot them on the same plot. The unstable manifold is \
blue and the stable manifold is green.\
\>", "Text"],
Cell[BoxData[
\(a\ = 0.8; \nb = \ 0.2; \n
orb\ = \ SeekPeriodicOrbit[henon, 1, {1, 1}]; \nX0\ = \ orb[0]; \n
Wu\ = \ UnstableManifold[henon, X0, 6, Subdivisions -> 200]; \n
Ws\ = \ UnstableManifold[henoninv, X0, 4]; \n
PhasePlot[{Ws, Wu}, Lines -> True, Points -> False,
LineColor -> {Green, Blue}, \nPlotRange -> {{\(-2\), 2}, {\(-2\), 2}}];
\)], "Input"],
Cell["\<\
As you can see, the stable manifold grows a lot faster (under \
backward iteration) than the unstable manifold does. Now go back and \
gradually increase the value of the parameter a over the range [0.8,1.6] to \
watch how the manifolds change. If you get ragged picture try increasing the \
subdivisions used to compute the unstable manifold or decreasing the number \
of iterations used to extend it.\
\>", "Text"],
Cell["\<\
Now the classic Henon attractor occurs for the parameters given \
below. Here is a picture showing an orbit on the attractor as well as the \
stable and unstable manifolds of the saddle fixed point inside the attractor. \
\
\>", "Text"],
Cell[BoxData[
\(a\ = \ 1.4; \nb = \ 0.3; \n
orb\ = \ SeekPeriodicOrbit[henon, 1, {1, 1}]; \nX0\ = \ orb[0]; \n
Wu\ = \ UnstableManifold[henon, X0, 6, Subdivisions -> 100]; \n
Ws\ = \ UnstableManifold[henoninv, X0, 4]; \n
orb\ = \ henon^{100, 1000}[{1.0, 1.0}]; \n
PhasePlot[{Ws, Wu, orb}, Lines -> {True, True, False},
Points -> {False, False, True}, LineColor -> {Green, Blue, Black},
PointSize -> 0.01, \nPlotRange -> {{\(-2\), 2}, {\(-2\), 2}}]; \)],
"Input"],
Cell["\<\
As you can see, the attractor is virtually indistinguishable from \
the unstable manifold of the saddle point.\
\>", "Text"],
Cell["\<\
Another interesting 2D diffeomorphism is the so-called standard \
mapping of the annulus. It is a kind of discrete dynamical system analogue \
of the forced, undamped pendulum differential equation. The x-variable is an \
angle (analogous to the position of the pendulum) so the state space is \
really a cylinder or annulus instead of a plane. First we enter the formulas \
for the map and its inverse. The Mod command in the first component forces \
the angular coordinate x to always lie in the interval from 0 to 2*Pi\
\>",
"Text"],
Cell[BoxData[{
\(st1[{x_, y_}]\ := \ {Mod[x + y + a*Sin[x], 2*Pi], y + a*Sin[x]}\),
\(stmap\ = \ DS[st1]\),
\(stinv1[{x_, y_}]\ := \ {Mod[x - y, 2*Pi], y - a*Sin[x - y]}\),
\(stinv\ = \ DS[stinv1]\)}], "Input"],
Cell["First we just plot some orbits.", "Text"],
Cell[BoxData[
\(a\ = \ 0.1; \norb1\ = \ stmap^{0, 1000}[{0.1, 0}]; \n
orb2\ = \ stmap^{0, 1000}[{0, 0.1}]; \n
orb3\ = \ stmap^{0, 1000}[{1, 0}]; \norb4\ = \ stmap^{0, 1000}[{2, 0}];
\norb5\ = \ stmap^{0, 1000}[{0, 1}]; \n
orb6\ = \ stmap^{0, 1000}[{0, 2}]; \n
orb7\ = \ stmap^{0, 1000}[{0, \(-1\)}]; \n
orb8\ = \ stmap^{0, 1000}[{0, \(-2\)}]; \n
PhasePlot[{orb1, orb2, orb3, orb4, orb5, orb6, orb7, orb8},
PointSize -> 0.002]; \)], "Input"],
Cell["\<\
As you can see the plot resembles the familiar picture for the \
pendulum. Now go back to the last cell and try some other values for the \
parameter a. For example, try a = 0.5. You will nottice a small chaotic \
band associated with the orbits which started out near the origin. Later we \
will see that this is an example of homoclinic chaos produced by a \
transversal homoclinic point to the saddle fixed point at the origin. Try \
some larger values of a such as a = 0.8, 1.0, 1.2, 1.5. You will find that \
the chaotic zone gradually takes over thestate space.\
\>", "Text"],
Cell["\<\
To understand the origin of the chaos we will plot the stable and \
unstable manifolds of the fixed point at the origin. It turns out that the \
best pictures are obtained by thinking of the points
(0,0), (2*Pi,0), (4*Pi,0) etc as distinct. For this reason we will first \
redefine the dynamical system so that it does not force the x-coordinate to \
lie in the interval [0, 2*Pi].\
\>", "Text"],
Cell[BoxData[{
\(st2[{x_, y_}]\ := \ {x + y + a*Sin[x], y + a*Sin[x]}\),
\(stmap2\ = \ DS[st2]\),
\(sti2[{x_, y_}]\ := \ {x - y, y - a*Sin[x - y]}\),
\(stinv2\ = \ DS[sti2]\)}], "Input"],
Cell["First we check that the origin is a saddle point.", "Text"],
Cell[BoxData[{
\(a = 0.8; \norb\ = \ SeekPeriodicOrbit[stmap2, 1, {0, 0}]\),
\(Multiplier[stmap2, orb]\)}], "Input"],
Cell["\<\
Now we will plot the stable and unstable manifolds of the origin \
and also those of its translate at (2*Pi,0). You can easily see the \
transverse intersections of the stable and unstable manifolds and also how \
these manifolds accumulate on themselves near the fixed point.\
\>", "Text"],
Cell[BoxData[
\(a = 0.8; \niters\ = \ 3; \n
wu\ = \ UnstableManifold[stmap2, {0, 0}, 7, Radius -> 0.5]; \n
ws\ = \ UnstableManifold[stinv2, {0, 0}, 7, Radius -> 0.5]; \n
wu1\ = \ UnstableManifold[stmap2, {2*Pi, 0}, 7, Radius -> 0.5]; \n
ws1\ = \ UnstableManifold[stinv2, {2*Pi, 0}, 7, Radius -> 0.5]; \n\n
PhasePlot[{wu, ws, wu1, ws1}, Lines -> True, \ Points -> False, \n
LineColor -> {Blue, Green, Blue, Green},
PlotRange -> {{\(-2\), 8}, {\(-3\), 3}}]; \n\n
PhasePlot[{wu, ws, wu1, ws1}, Lines -> True, \ Points -> False, \n
LineColor -> {Blue, Green, Blue, Green},
PlotRange -> {{\(-1\), 1}, {\(-1\), 1}}]; \)], "Input"],
Cell["Go back and try increasing a to a = 1.0 and a = 1.2.", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Iterated Function Systems and Fractals", "Subsection"],
Cell["\<\
In this section we will use iterated function systems to construct \
some of the famous fractals. First we have to enter a collection of \
similarity transformations and collect them together into an IFS object. \
Here are the three maps used to make the Sierpinski triangle.\
\>", "Text"],
Cell[BoxData[
\(sf1[{x_, y_}]\ := \ {x/2, y/2}; \n
sf2[{x_, y_}]\ := \ {0.5 + x/2, y/2}; \n
sf3[{x_, y_}]\ := \ {0.25 + x/2, 0.25*Sqrt[3] + y/2}; \n
sierpinski\ = IFS[\ sf1, sf2, sf3]\)], "Input"],
Cell["\<\
The program know how to apply some iterate of the IFS to a given \
initial point. The result is a list of points obtained by applying all the \
maps to the given point. In this case the nth iterate would produce a list \
of 3^n points.\
\>", "Text"],
Cell[BoxData[{
\(sierpinski[{0, 0}]\),
\(sierpinski^2[{0, 0}]\),
\(sierpinski^3[{0, 0}]\)}], "Input"],
Cell["\<\
You can plot such an iterate to get an approximate picture of the \
attractor of the iterated function system.\
\>", "Text"],
Cell[BoxData[
\(\(ListPlot[sierpinski^7[{0, 0}], AspectRatio -> 1,
PlotRange -> {{0, 1}, {0, 1}}]; \)\)], "Input"],
Cell[TextData[{
"It is interesting to use an initial shape other than a point. Simple \
shapes can be represented as a ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" Line object. Here are some examples"
}], "Text"],
Cell[BoxData[
\(point\ = \ Line[{{0, 0}, {0, 0}}]; \n
segment\ = \ Line[{{0, 0}, {1, 0}}]; \n
triangle\ = \ Line[{{0, 0}, {1, 0}, {0.5, N[Sqrt[3]/2]}, {0, 0}}]; \n
square\ = \ Line[{{0, 0}, {0, 1}, {1, 1}, {1, 0}, {0, 0}}]; \n
diamond\ = \
Line[{{0.5, 0}, {0, 0.5}, {\(-0.5\), 0}, {0, \(-0.5\)}, {0.5, 0}}]; \n\n
Show[Graphics[{triangle, square, diamond}], AspectRatio -> 1,
Axes -> True]; \)], "Input"],
Cell["Now we can apply the IFS to one of these shapes.", "Text"],
Cell[BoxData[{
\(frac\ = \ sierpinski[triangle]\),
\(\(Show[Graphics[frac], AspectRatio -> 1,
PlotRange -> {{0, 1}, {0, 1}}]; \)\)}], "Input"],
Cell["\<\
Here is a movie showing the convergence of the iterates to the \
attractor.\
\>", "Text"],
Cell[BoxData[
\(Do[frac\ = \ sierpinski^i[triangle]; \n
\(Show[Graphics[frac], AspectRatio -> 1, PlotRange -> {{0, 1}, {0, 1}}];
\), {i, 1, 7}]\)], "Input"],
Cell[TextData[{
"To animate them you have to select them by clicking with the mouse on the \
second 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["\<\
Next we can demonstrate the fact that the attractor is independent \
of the initial shape by plotting a similar picture starting from a square or \
a diamond.\
\>", "Text"],
Cell[BoxData[
\(iterations\ \ = \ 3; \n\nfrac\ = \ sierpinski^iterations[square]; \n
Show[Graphics[frac], AspectRatio -> 1, PlotRange -> {{0, 1}, {0, 1}}]; \n
\nfrac\ = \ sierpinski^iterations[diamond]; \n
Show[Graphics[frac], AspectRatio -> 1, PlotRange -> {{0, 1}, {0, 1}}];
\)], "Input"],
Cell["\<\
Go back and increase the number of iterations to see that, in the \
end, it all comes out the same.\
\>", "Text"],
Cell["\<\
Now we will have a look at a couple other fractals. Here are some \
Sierpinski carpets using 5 and 8 squares. First we define 9 maps and then \
take subsets of them to construct the IFS.\
\>", "Text"],
Cell[BoxData[
\(s5f1[{x_, y_}]\ := \ {x/3, y/3}; \n
s5f2[{x_, y_}]\ := \ {x/3 + 1/3, y/3}; \n
s5f3[{x_, y_}]\ := \ {x/3 + 2/3, y/3}; \n
s5f4[{x_, y_}]\ := \ {x/3, y/3 + 1/3}; \n
s5f5[{x_, y_}]\ := \ {x/3 + 1/3, y/3 + 1/3}; \n
s5f6[{x_, y_}] := \ {x/3 + 2/3, y/3 + 1/3}; \n
s5f7[{x_, y_}]\ := \ {x/3, y/3 + 2/3}; \n
s5f8[{x_, y_}]\ := \ {x/3 + 1/3, y/3 + 2/3}; \n
s5f9[{x_, y_}]\ := \ {x/3 + 2/3, y/3 + 2/3}; \n
square5\ = \ IFS[s5f1, s5f3, s5f5, s5f7, s5f9]; \n
square8\ = \ IFS[s5f1, s5f2, s5f3, s5f4, s5f6, s5f7, s5f8, s5f9]; \)],
"Input"],
Cell["\<\
And now the fractals. With an IFS like square8 which has 8 \
functions in it, you need to be careful about how many iterations you use. \
Even n = 4 iterations gives you 8^4 = 4096 separate squares in your \
picture.\
\>", "Text"],
Cell[BoxData[
\(frac\ = square5^4[square]; \n
Show[Graphics[frac], AspectRatio -> 1, PlotRange -> {{0, 1}, {0, 1}}]; \n
\nfrac\ = \ square8^4[square]; \n
Show[Graphics[frac], AspectRatio -> 1, PlotRange -> {{0, 1}, {0, 1}}];
\)], "Input"],
Cell["\<\
Finally, we will plot the Koch curve. This time the functions \
involve some rotation.\
\>", "Text"],
Cell[BoxData[
\(kf1[{x_, y_}]\ := \ {x/3, y/3}; \n
kf4[{x_, y_}]\ := \ {x/3 + 2/3, y/3}; \n
kf2[{x_, y_}]\ :=
\ {x/6 - Sqrt[3]*y/6 + 1/3, \n\t\t\tSqrt[3]*x/6 + y/6}; \n
kf3[{x_, y_}]\ :=
\ {x/6 + Sqrt[3]*y/6 + 1/2, \n\t\t\t
\(-Sqrt[3]\)*x/6 + y/6 + 1/\((2*Sqrt[3])\)}; \n
koch\ = \ IFS[kf1, kf2, kf3, kf4]; \)], "Input"],
Cell["And here's the movie:", "Text"],
Cell[BoxData[
\(Do[frac\ = koch^i[segment]; \n
\(Show[Graphics[frac], AspectRatio -> 1,
PlotRange -> {{0, 1}, {\(-1\)/2, 1/2}}];\), {i, 1, 5}]\)], "Input"],
Cell["You can animate as before.", "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell["Delay Embedding and Attractor Reconstruction", "Subsection"],
Cell["\<\
In this section we will compute some orbits of vector dynamical \
systems, and then try to reconstruct the dynamics using a scalar \
\"measurement function\". For example here is a simple rotation map of the \
plane producing an orbit which is dense in a circle.\
\>", "Text"],
Cell[BoxData[{
\(f1[{x_, y_}]\ :=
\ {Cos[2*Pi*omega]*x - Sin[2*Pi*omega]*y,
Sin[2*Pi*omega]*x + Cos[2*Pi*omega]*y}\),
\(\(f\ = \ DS[f1]; \)\)}], "Input"],
Cell["\<\
Here is a rotation by one radian, an irrational fraction of a whole \
rotation.\
\>", "Text"],
Cell[BoxData[
\(omega\ = \ 1/\((2*Pi)\); \norb\ = \ f^{0, 1000}[{1.0, 0.0}]; \n
PhasePlot[orb]; \)], "Input"],
Cell["\<\
Next we construct a scalar time series by applying some scalar \
measurement function to each point of this orbit. We will try the \
measurement function h = x+y^2\
\>", "Text"],
Cell[BoxData[
\(h[{x_, y_}]\ := \ x + y^2; \ntimeseries\ = \ Orbit[Map[h, orb[]]];
\)], "Input"],
Cell["\<\
Now that we have a scalar time series we replot it as a series of \
delay vectors and compare it with the original plot.\
\>", "Text"],
Cell[BoxData[
\(\(DelayPlot[timeseries, 2, PointSize -> 0.004]; \)\)], "Input"],
Cell["\<\
This is not homeomorphic to the original circle. Moving up to \
three dimensions solves the problem however.\
\>", "Text"],
Cell[BoxData[
\(\(DelayPlot[timeseries, 3, PointSize -> 0.004]; \)\)], "Input"],
Cell["\<\
Now for a more interesting orbit. Here is an orbit from the Henon \
attractor.\
\>", "Text"],
Cell[BoxData[{
\(h1[{x_, y_}]\ := {a - x^2 + b*y, x}\),
\(henon\ = \ DS[h1]; \na = 1.4; b = 0.3; \n
orb\ = \ henon^{100, 10000}[{1, 1}]; \n
PhasePlot[orb, PointSize -> 0.002]; \)}], "Input"],
Cell["As a measurement function we will use x+y", "Text"],
Cell[BoxData[{
\(h[{x_, y_}]\ := \ x + y\),
\(timeseries\ = \ Orbit[Map[h, orb[]]]; \n
DelayPlot[timeseries, 2, PointSize -> 0.002]; \)}], "Input"],
Cell["\<\
Again going to 3D eliminates the self-intersections producing an \
object homeomorphic to the original attractor.\
\>", "Text"],
Cell[BoxData[
\(\(DelayPlot[timeseries, 3, PointSize -> 0.002]; \)\)], "Input"],
Cell["\<\
Finally we consider a higher dimensional example. You can \
construct a 4D system analogous to the rotation of the circle, which involves \
two separate rotations. The resulting orbit will generally be dense in the \
surface of a torus (the product space of two circles).\
\>", "Text"],
Cell[BoxData[{
\(f1[{x1_, y1_, x2_, y2_}]\ :=
\ {Cos[2*Pi*omega1]*x1 - Sin[2*Pi*omega1]*y1,
Sin[2*Pi*omega1]*x1 + Cos[2*Pi*omega1]*y1, \n\t\t
Cos[2*Pi*omega2]*x2 - Sin[2*Pi*omega2]*y2,
Sin[2*Pi*omega2]*x2 + Cos[2*Pi*omega2]*y2}\),
\(\(f\ = \ DS[f1]; \)\)}], "Input"],
Cell[BoxData[
\(omega1\ = \ 1/\((2*Pi)\); \nomega2\ = \ 1.618; \n
orb\ = \ f^{0, 5000}[{1.0, 0.0, 1.0, 0.0}]; \norb[{0, 20}]\)], "Input"],
Cell["\<\
The orbit is a list of 4D points so we can't plot it directly at \
all. We will use a measurement function to reduce it to a scalar time \
series.\
\>", "Text"],
Cell[BoxData[{
\(h[{x1_, y1_, x2_, y2_}]\ := \ x1 + x2\),
\(\(timeseries\ = \ Orbit[Map[h, orb[]]]; \)\)}], "Input"],
Cell["\<\
Here are the delay plots in 2 and in 3 dimensions. You can almost \
see the torus but it seems that even in 3D there are still \
self-intersections. According to the theory we should really be using an \
embedding dimension of 5 !\
\>", "Text"],
Cell[BoxData[
\(DelayPlot[timeseries, 2, PointSize -> 0.002]; \n
DelayPlot[timeseries, 3, PointSize -> 0.002]; \)], "Input"]
}, Closed]]
}, Closed]]
}, Open ]]
},
FrontEndVersion->"Macintosh 3.0",
ScreenRectangle->{{0, 832}, {0, 604}},
AutoGeneratedPackage->Automatic,
WindowSize->{529, 400},
WindowMargins->{{Automatic, 124}, {Automatic, 7}},
PrintingCopies->1,
PrintingPageRange->{1, Automatic},
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[CellGroupData[{
Cell[1731, 51, 76, 0, 91, "Title"],
Cell[CellGroupData[{
Cell[1832, 55, 61, 1, 50, "Section",
InitializationCell->True],
Cell[1896, 58, 631, 15, 100, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[2552, 77, 127, 2, 44, "Subsubsection",
InitializationCell->True],
Cell[2682, 81, 1042, 25, 146, "Text",
InitializationCell->True],
Cell[3727, 108, 800, 21, 114, "Text",
InitializationCell->True],
Cell[4530, 131, 1039, 30, 134, "Text",
InitializationCell->True],
Cell[5572, 163, 146, 3, 43, "Input",
InitializationCell->True],
Cell[5721, 168, 169, 3, 59, "Input",
InitializationCell->True],
Cell[5893, 173, 519, 15, 82, "Text",
InitializationCell->True],
Cell[6415, 190, 271, 6, 91, "Input",
InitializationCell->True],
Cell[6689, 198, 851, 30, 101, "Text",
InitializationCell->True],
Cell[7543, 230, 305, 6, 139, "Input",
InitializationCell->True],
Cell[7851, 238, 284, 7, 48, "Text",
InitializationCell->True]
}, Closed]],
Cell[CellGroupData[{
Cell[8172, 250, 65, 1, 30, "Subsection",
InitializationCell->True],
Cell[8240, 253, 410, 11, 64, "Text",
InitializationCell->True],
Cell[8653, 266, 367, 7, 78, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[9045, 277, 658, 14, 187, "Input",
InitializationCell->True],
Cell[9706, 293, 241, 3, 71, "Output"]
}, Open ]],
Cell[9962, 299, 549, 11, 96, "Text",
InitializationCell->True],
Cell[10514, 312, 449, 11, 187, "Input",
InitializationCell->True],
Cell[10966, 325, 287, 8, 47, "Text",
InitializationCell->True],
Cell[11256, 335, 735, 15, 315, "Input",
InitializationCell->True],
Cell[11994, 352, 241, 5, 46, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[12260, 361, 603, 10, 235, "Input",
InitializationCell->True],
Cell[12866, 373, 342, 5, 101, "Output"]
}, Open ]],
Cell[13223, 381, 1021, 22, 184, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[14269, 407, 499, 10, 171, "Input",
InitializationCell->True],
Cell[14771, 419, 209, 3, 71, "Output"]
}, Open ]],
Cell[14995, 425, 176, 5, 31, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[15196, 434, 2832, 51, 971, "Input",
InitializationCell->True],
Cell[18031, 487, 367, 5, 116, "Output"]
}, Open ]]
}, Closed]],
Cell[CellGroupData[{
Cell[18447, 498, 68, 1, 30, "Subsection",
InitializationCell->True],
Cell[18518, 501, 430, 8, 78, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[18973, 513, 2349, 39, 875, "Input",
InitializationCell->True],
Cell[21325, 554, 589, 8, 70, "Output"]
}, Open ]]
}, Closed]],
Cell[CellGroupData[{
Cell[21963, 568, 58, 1, 30, "Subsection",
InitializationCell->True],
Cell[22024, 571, 687, 12, 127, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[22736, 587, 79, 1, 42, "Subsubsection",
InitializationCell->True],
Cell[22818, 590, 496, 14, 81, "Text",
InitializationCell->True],
Cell[23317, 606, 2163, 50, 683, "Input",
InitializationCell->True],
Cell[25483, 658, 157, 4, 46, "Text",
InitializationCell->True],
Cell[25643, 664, 394, 7, 123, "Input",
InitializationCell->True],
Cell[26040, 673, 284, 8, 48, "Text",
InitializationCell->True],
Cell[26327, 683, 1058, 22, 315, "Input",
InitializationCell->True],
Cell[27388, 707, 171, 4, 46, "Text",
InitializationCell->True],
Cell[27562, 713, 2096, 34, 731, "Input",
InitializationCell->True],
Cell[29661, 749, 96, 1, 30, "Text",
InitializationCell->True],
Cell[29760, 752, 102, 2, 27, "Input",
InitializationCell->True],
Cell[29865, 756, 94, 1, 30, "Text",
InitializationCell->True],
Cell[29962, 759, 2545, 44, 731, "Input",
InitializationCell->True]
}, Closed]],
Cell[CellGroupData[{
Cell[32544, 808, 85, 1, 28, "Subsubsection",
InitializationCell->True],
Cell[32632, 811, 602, 15, 98, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[33259, 830, 1427, 26, 443, "Input",
InitializationCell->True],
Cell[34689, 858, 583, 8, 176, "Output"]
}, Open ]],
Cell[35287, 869, 303, 7, 64, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[35615, 880, 1075, 22, 331, "Input",
InitializationCell->True],
Cell[36693, 904, 448, 6, 131, "Output"]
}, Open ]],
Cell[37156, 913, 444, 12, 81, "Text",
InitializationCell->True],
Cell[CellGroupData[{
Cell[37625, 929, 1515, 30, 491, "Input",
InitializationCell->True],
Cell[39143, 961, 387, 5, 116, "Output"]
}, Open ]],
Cell[39545, 969, 176, 4, 43, "Input",
InitializationCell->True]
}, Open ]],
Cell[CellGroupData[{
Cell[39758, 978, 76, 1, 28, "Subsubsection",
InitializationCell->True],
Cell[39837, 981, 164, 4, 46, "Text",
InitializationCell->True],
Cell[40004, 987, 490, 11, 155, "Input",
InitializationCell->True]
}, Closed]],
Cell[CellGroupData[{
Cell[40531, 1003, 69, 1, 28, "Subsubsection",
InitializationCell->True],
Cell[CellGroupData[{
Cell[40625, 1008, 3403, 57, 1211, "Input",
InitializationCell->True],
Cell[44031, 1067, 948, 13, 266, "Output"]
}, Open ]]
}, Open ]]
}, Closed]],
Cell[CellGroupData[{
Cell[45040, 1087, 47, 0, 30, "Subsection"],
Cell[45090, 1089, 371, 6, 78, "Text"],
Cell[45464, 1097, 505, 9, 187, "Input",
InitializationCell->True]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[46018, 1112, 61, 0, 30, "Section"],
Cell[46082, 1114, 1054, 27, 136, "Text"],
Cell[47139, 1143, 689, 14, 114, "Text"],
Cell[CellGroupData[{
Cell[47853, 1161, 50, 0, 46, "Subsection"],
Cell[47906, 1163, 779, 21, 118, "Text"],
Cell[48688, 1186, 89, 2, 43, "Input"],
Cell[48780, 1190, 348, 9, 64, "Text"],
Cell[49131, 1201, 57, 1, 43, "Input"],
Cell[49191, 1204, 198, 4, 46, "Text"],
Cell[49392, 1210, 145, 3, 46, "Text"],
Cell[49540, 1215, 45, 1, 27, "Input"],
Cell[49588, 1218, 402, 9, 80, "Text"],
Cell[49993, 1229, 61, 1, 27, "Input"],
Cell[50057, 1232, 240, 5, 47, "Text"],
Cell[50300, 1239, 80, 3, 59, "Input"],
Cell[50383, 1244, 438, 11, 65, "Text"],
Cell[50824, 1257, 76, 1, 27, "Input"],
Cell[50903, 1260, 93, 3, 30, "Text"],
Cell[50999, 1265, 56, 0, 30, "Text"],
Cell[51058, 1267, 65, 1, 27, "Input"]
}, Open ]],
Cell[CellGroupData[{
Cell[51160, 1273, 41, 0, 46, "Subsection"],
Cell[51204, 1275, 245, 6, 47, "Text"],
Cell[51452, 1283, 118, 2, 59, "Input"],
Cell[51573, 1287, 442, 10, 80, "Text"],
Cell[52018, 1299, 143, 4, 31, "Text"],
Cell[52164, 1305, 107, 2, 59, "Input"],
Cell[52274, 1309, 510, 9, 95, "Text"],
Cell[52787, 1320, 122, 2, 59, "Input"],
Cell[52912, 1324, 127, 4, 31, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[53076, 1333, 65, 0, 30, "Subsection"],
Cell[53144, 1335, 430, 11, 64, "Text"],
Cell[53577, 1348, 89, 4, 75, "Input"],
Cell[53669, 1354, 263, 8, 48, "Text"],
Cell[53935, 1364, 111, 3, 59, "Input"],
Cell[54049, 1369, 257, 8, 50, "Text"],
Cell[54309, 1379, 150, 3, 75, "Input"],
Cell[54462, 1384, 579, 13, 98, "Text"],
Cell[55044, 1399, 372, 7, 63, "Text"],
Cell[55419, 1408, 170, 2, 59, "Input"],
Cell[55592, 1412, 386, 9, 65, "Text"],
Cell[55981, 1423, 259, 6, 47, "Text"],
Cell[56243, 1431, 108, 2, 43, "Input"],
Cell[56354, 1435, 64, 0, 30, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[56455, 1440, 73, 0, 30, "Subsection"],
Cell[56531, 1442, 547, 12, 98, "Text"],
Cell[57081, 1456, 253, 7, 48, "Text"],
Cell[57337, 1465, 162, 2, 59, "Input"],
Cell[57502, 1469, 304, 5, 62, "Text"],
Cell[57809, 1476, 70, 1, 27, "Input"],
Cell[57882, 1479, 128, 3, 30, "Text"],
Cell[58013, 1484, 71, 1, 27, "Input"],
Cell[58087, 1487, 256, 5, 64, "Text"],
Cell[58346, 1494, 84, 2, 43, "Input"],
Cell[58433, 1498, 267, 7, 47, "Text"],
Cell[58703, 1507, 103, 2, 43, "Input"],
Cell[58809, 1511, 328, 6, 78, "Text"],
Cell[59140, 1519, 48, 1, 27, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[59225, 1525, 64, 0, 30, "Subsection"],
Cell[59292, 1527, 290, 6, 63, "Text"],
Cell[59585, 1535, 110, 2, 59, "Input"],
Cell[59698, 1539, 1080, 21, 162, "Text"],
Cell[60781, 1562, 233, 5, 62, "Text"],
Cell[61017, 1569, 162, 2, 75, "Input"],
Cell[61182, 1573, 537, 11, 95, "Text"],
Cell[61722, 1586, 130, 3, 30, "Text"],
Cell[61855, 1591, 160, 2, 75, "Input"],
Cell[62018, 1595, 710, 15, 114, "Text"],
Cell[62731, 1612, 370, 8, 79, "Text"],
Cell[63104, 1622, 170, 3, 91, "Input"],
Cell[63277, 1627, 138, 3, 46, "Text"],
Cell[63418, 1632, 97, 2, 27, "Input"],
Cell[63518, 1636, 412, 7, 79, "Text"],
Cell[63933, 1645, 102, 2, 43, "Input"],
Cell[64038, 1649, 325, 7, 64, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[64400, 1661, 45, 0, 30, "Subsection"],
Cell[64448, 1663, 155, 4, 46, "Text"],
Cell[64606, 1669, 200, 5, 47, "Text"],
Cell[64809, 1676, 285, 9, 50, "Text"],
Cell[65097, 1687, 576, 12, 95, "Text"]
}, Closed]]
}, Closed]],
Cell[CellGroupData[{
Cell[65722, 1705, 70, 0, 30, "Section"],
Cell[65795, 1707, 171, 4, 46, "Text"],
Cell[CellGroupData[{
Cell[65991, 1715, 43, 0, 46, "Subsection"],
Cell[66037, 1717, 582, 13, 99, "Text"],
Cell[66622, 1732, 75, 1, 43, "Input"],
Cell[66700, 1735, 280, 7, 64, "Text"],
Cell[66983, 1744, 199, 4, 107, "Input"],
Cell[67185, 1750, 373, 6, 78, "Text"],
Cell[67561, 1758, 132, 4, 91, "Input"],
Cell[67696, 1764, 424, 7, 78, "Text"],
Cell[68123, 1773, 97, 2, 59, "Input"],
Cell[68223, 1777, 251, 5, 47, "Text"],
Cell[68477, 1784, 110, 2, 43, "Input"],
Cell[68590, 1788, 411, 7, 78, "Text"],
Cell[69004, 1797, 174, 3, 75, "Input"],
Cell[69181, 1802, 340, 8, 63, "Text"],
Cell[69524, 1812, 357, 7, 171, "Input"],
Cell[69884, 1821, 306, 5, 62, "Text"],
Cell[70193, 1828, 538, 8, 235, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[70768, 1841, 42, 0, 30, "Subsection"],
Cell[70813, 1843, 277, 5, 62, "Text"],
Cell[71093, 1850, 100, 2, 43, "Input"],
Cell[71196, 1854, 264, 5, 62, "Text"],
Cell[71463, 1861, 182, 3, 75, "Input"],
Cell[71648, 1866, 817, 12, 158, "Text"],
Cell[72468, 1880, 331, 6, 123, "Input"],
Cell[72802, 1888, 152, 3, 46, "Text"],
Cell[72957, 1893, 239, 4, 91, "Input"]
}, Closed]],
Cell[CellGroupData[{
Cell[73233, 1902, 53, 0, 30, "Subsection"],
Cell[73289, 1904, 369, 6, 78, "Text"],
Cell[73661, 1912, 89, 2, 43, "Input"],
Cell[73753, 1916, 301, 5, 62, "Text"],
Cell[74057, 1923, 50, 1, 27, "Input"],
Cell[74110, 1926, 234, 5, 62, "Text"],
Cell[74347, 1933, 53, 1, 27, "Input"],
Cell[74403, 1936, 149, 2, 43, "Input"],
Cell[74555, 1940, 155, 3, 46, "Text"],
Cell[74713, 1945, 149, 2, 43, "Input"],
Cell[74865, 1949, 332, 6, 78, "Text"],
Cell[75200, 1957, 155, 2, 59, "Input"],
Cell[75358, 1961, 489, 9, 95, "Text"],
Cell[75850, 1972, 163, 2, 59, "Input"],
Cell[76016, 1976, 176, 4, 46, "Text"],
Cell[76195, 1982, 239, 5, 62, "Text"],
Cell[76437, 1989, 102, 2, 43, "Input"],
Cell[76542, 1993, 347, 6, 78, "Text"],
Cell[76892, 2001, 111, 2, 43, "Input"],
Cell[77006, 2005, 257, 5, 62, "Text"],
Cell[77266, 2012, 57, 1, 27, "Input"],
Cell[77326, 2015, 155, 2, 59, "Input"],
Cell[77484, 2019, 134, 3, 46, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[77655, 2027, 51, 0, 30, "Subsection"],
Cell[77709, 2029, 219, 4, 46, "Text"],
Cell[77931, 2035, 100, 2, 43, "Input"],
Cell[78034, 2039, 436, 9, 80, "Text"],
Cell[78473, 2050, 109, 2, 59, "Input"],
Cell[78585, 2054, 99, 3, 30, "Text"],
Cell[78687, 2059, 55, 1, 27, "Input"],
Cell[78745, 2062, 230, 4, 62, "Text"],
Cell[78978, 2068, 48, 1, 27, "Input"],
Cell[79029, 2071, 179, 4, 47, "Text"],
Cell[79211, 2077, 54, 1, 27, "Input"],
Cell[79268, 2080, 72, 1, 27, "Input"],
Cell[79343, 2083, 196, 4, 46, "Text"],
Cell[79542, 2089, 151, 2, 43, "Input"],
Cell[79696, 2093, 269, 6, 63, "Text"],
Cell[79968, 2101, 197, 3, 75, "Input"],
Cell[80168, 2106, 127, 3, 46, "Text"],
Cell[80298, 2111, 193, 4, 46, "Text"],
Cell[80494, 2117, 117, 2, 43, "Input"],
Cell[80614, 2121, 203, 4, 46, "Text"],
Cell[80820, 2127, 393, 7, 171, "Input"],
Cell[81216, 2136, 428, 7, 94, "Text"],
Cell[81647, 2145, 250, 5, 62, "Text"],
Cell[81900, 2152, 509, 9, 203, "Input"],
Cell[82412, 2163, 134, 3, 46, "Text"],
Cell[82549, 2168, 552, 9, 110, "Text"],
Cell[83104, 2179, 233, 4, 91, "Input"],
Cell[83340, 2185, 47, 0, 30, "Text"],
Cell[83390, 2187, 489, 9, 203, "Input"],
Cell[83882, 2198, 600, 9, 110, "Text"],
Cell[84485, 2209, 407, 7, 78, "Text"],
Cell[84895, 2218, 209, 4, 75, "Input"],
Cell[85107, 2224, 65, 0, 30, "Text"],
Cell[85175, 2226, 126, 2, 59, "Input"],
Cell[85304, 2230, 302, 5, 62, "Text"],
Cell[85609, 2237, 683, 11, 299, "Input"],
Cell[86295, 2250, 68, 0, 30, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[86400, 2255, 60, 0, 30, "Subsection"],
Cell[86463, 2257, 301, 5, 62, "Text"],
Cell[86767, 2264, 216, 4, 75, "Input"],
Cell[86986, 2270, 261, 5, 62, "Text"],
Cell[87250, 2277, 115, 3, 59, "Input"],
Cell[87368, 2282, 134, 3, 46, "Text"],
Cell[87505, 2287, 125, 2, 43, "Input"],
Cell[87633, 2291, 230, 6, 50, "Text"],
Cell[87866, 2299, 446, 8, 171, "Input"],
Cell[88315, 2309, 64, 0, 30, "Text"],
Cell[88382, 2311, 159, 3, 59, "Input"],
Cell[88544, 2316, 99, 3, 30, "Text"],
Cell[88646, 2321, 175, 3, 75, "Input"],
Cell[88824, 2326, 920, 21, 163, "Text"],
Cell[89747, 2349, 182, 4, 46, "Text"],
Cell[89932, 2355, 314, 5, 155, "Input"],
Cell[90249, 2362, 123, 3, 30, "Text"],
Cell[90375, 2367, 211, 4, 46, "Text"],
Cell[90589, 2373, 598, 12, 203, "Input"],
Cell[91190, 2387, 242, 5, 62, "Text"],
Cell[91435, 2394, 261, 5, 123, "Input"],
Cell[91699, 2401, 111, 3, 30, "Text"],
Cell[91813, 2406, 367, 8, 123, "Input"],
Cell[92183, 2416, 37, 0, 30, "Text"],
Cell[92223, 2418, 177, 3, 75, "Input"],
Cell[92403, 2423, 42, 0, 30, "Text"]
}, Closed]],
Cell[CellGroupData[{
Cell[92482, 2428, 66, 0, 30, "Subsection"],
Cell[92551, 2430, 287, 5, 62, "Text"],
Cell[92841, 2437, 180, 4, 75, "Input"],
Cell[93024, 2443, 103, 3, 30, "Text"],
Cell[93130, 2448, 119, 2, 59, "Input"],
Cell[93252, 2452, 188, 4, 46, "Text"],
Cell[93443, 2458, 107, 2, 43, "Input"],
Cell[93553, 2462, 145, 3, 46, "Text"],
Cell[93701, 2467, 83, 1, 27, "Input"],
Cell[93787, 2470, 133, 3, 46, "Text"],
Cell[93923, 2475, 83, 1, 27, "Input"],
Cell[94009, 2478, 103, 3, 30, "Text"],
Cell[94115, 2483, 210, 4, 91, "Input"],
Cell[94328, 2489, 58, 0, 30, "Text"],
Cell[94389, 2491, 162, 3, 59, "Input"],
Cell[94554, 2496, 137, 3, 46, "Text"],
Cell[94694, 2501, 83, 1, 27, "Input"],
Cell[94780, 2504, 297, 5, 62, "Text"],
Cell[95080, 2511, 310, 6, 107, "Input"],
Cell[95393, 2519, 148, 2, 75, "Input"],
Cell[95544, 2523, 171, 4, 46, "Text"],
Cell[95718, 2529, 127, 2, 43, "Input"],
Cell[95848, 2533, 256, 5, 62, "Text"],
Cell[96107, 2540, 132, 2, 43, "Input"]
}, Closed]]
}, Closed]]
}, Open ]]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)