In this section, we illustrate how Mathematica
can
be used for research or course work in nonlinear dynamics. Mathematica
is a
complete system for doing symbolic, graphical, and numerical manipulations on
symbols and data, and is commonly available on a wide range of machines from
microcomputers (386s and Macs) to
mainframes and supercomputers.
Mathematica is strong in two- and three-dimensional graphics, and, where appropriate, we would encourage its use in a nonlinear dynamics course. It can serve at least three important functions: (1) a means of generating complex graphical representations of data from experiments or simulations; (2) a method for double-checking complex algebraic manipulations first done by hand; and (3) a general system for writing both numerical and symbolic routines for the analysis of nonlinear equations.
What follows is text from a typical Mathematica session, typeset for legibility, used to produce some of the graphics and to double-check some of the algebraic results presented in this chapter. Of course, a real Mathematica session would not be so heavily commented.
(* This is a Mathematica comment statement.
Mathematica ignores
everything between star parentheses. *)
(* To try out this Mathematica session yourself, type everything
in
bold that is not enclosed in the comment statements.
The
following
line is Mathematica's answer typed in italic.
Mathematica's output
from graphical commands are not
printed here, but are left for the
reader to discover. *)
(* This notebook is written on a Macintosh. 7/22/90 nbt. *)
(* In this notebook we will analytically solve for the
period one and
period two orbits in the quadratic map and
plot their locations as a
function of the control
parameter, lambda. *)
(* First we define the quadratic function with the Mathematica
Function[{x,y, ...}, f(x,y, ...)] command, which takes
two
arguments, the first of which is a list of variables,
and the second
of which is the function.
The basic data type in Mathematica
is a list, and all lists
are enclosed between braces,
{x1, x2, x3, ...}.
Note that all arguments and functions in
Mathematica are
enclosed in square brackets f[], which differs from
the
standard mathematical notation of parentheses f().
This is, in
part, because square brackets are easier to
reach on the keyboard.
Also, note that in defining a
variable one must always put a space
around it, so xy
is equal to a single variable named "xy", while
x y with
a space between is equal to two variables, x and y. *)
(* To evaluate a Mathematica expression tap the enter key, not
the
return key. Now to our first Mathematica command: *)
f = Function[{lambda, x}, lambda x (1 - x)]
Function[{lambda, x}, lambda x (1 - x)]
(* Mathematica should respond by saying Out[1], which tells us
that
Mathematica successfully processed the first
command and has put the
result in the variable Out[1], as
well as the variable we created, f.
To evaluate the quadratic
map, we now can feed f two arguments,
the first
of which is lambda, and the second of which is x. *)
f[4, 1/2]
1
(* Mathematica should respond with the function Out[2],
which
contains the quadratic map evaluated at lambda = 4 and
x = 1/2. To
evaluate a list of values for x we could let
the x variable be a
list, i.e., a series of numbers enclosed
in braces {x1, x2, x3, ...}.
To try this command, type: *)
f[4, {0, 0.25, 0.5, 0.75, 1}]
{0, 0.75, 1., 0.75, 0}
(* To plot the quadratic map we use the Mathematica plot command,
Plot[f, {x, xmin, xmax}]. For instance a plot of the
quadratic map
for lambda = 4 is given by: *)
Plot[f[4,x], {x, 0, 1}]
(* It is as easy as pie to take a composite function in
Mathematica,
just type f[f[x]], so to plot f(f(x)) for
the quadratic map we
simply type: *)
Plot[f[4,f[4,x]], {x, 0, 1}]
(* Now let's find the locations of the period one orbits,
given by
the roots of f(x) = x. To find the roots,
we use the Mathematica
command Solve[eqns, vars]. Notice
that we are going to rename the
parameter "lambda" to
"a" just to save some space when printing the
answer.
The double equals "==" in Mathematica is equivalent to the
single equal "=" of mathematics. *)
Solve[f[a, x1] == x1, x1]
{\x1 --;SPMgt; 
}, {x1 --;SPMgt; 0}\
(* As expected, Mathematica finds two roots. Let's make a
function
out of the first root so that we can plot it later using
Plot. To do
this we need the following sequence of somewhat
cryptic commands: *)
r1 = %[[1]]
{x1 --;SPMgt; }
(* The roots are saved in a list of two items. To pull out the
first item we used the % command, a Mathematica variable
that always
holds the value of the last expression. In this
case it holds the
list of two roots, and the double square
bracket notation %[[1]]
tells Mathematica we want the
first
item in the list of two items.
Now we must pull out the last
part of the expression, 1 - 1/a, with
the replacement
command Replace[expr, rules]: *)
x1 = Replace[x1, r1]
(* To plot the location of the period one orbit we just use
the Plot
command again. *)
Plot[x1, {a, 0.9, 4}]
(* To find the location of the period two orbit,
we solve
for f(f(x)) = x. *)
Solve[f[a,f[a, x2]] == x2, x2]
{\x2 --;SPMgt; 0}, {x2 --;SPMgt; },
{x2 --;SPMgt; 
},
{x2 --;SPMgt;
}\
(* We find four roots, as expected. Before proceeding further, it's
a good idea
to try and
simplify the algebra for the last two new
roots by applying the
Simplify command to the last expression. *)
rt = Simplify[%]
{\x2 --;SPMgt; 0}, {x2 --;SPMgt; },
{x2 --;SPMgt; 
},
{x2 --;SPMgt; 
}\
(* And we can now pull out the positive and negative roots of the
period two orbit. *)
x2plus = Replace[x2, rt[[3]]]
x2minus = Replace[x2, rt[[4]]]
/2
(* As a last step, we can plot the location of the period one
orbit
and both branches of the period two orbit by making
a list of
functions to be plotted. *)
Plot[{x1, x2plus, x2minus}, {a, 1, 4}]
(* This is how we originally plotted the orbit stability diagram
in the
text. Mathematica can go on to find higher-order periodic
orbits
by numerically finding the roots of the nth composite of f,
if no exact solution
exists. *)