UC SCI 111, vM+. 
Notebook 1

Needs["Graphics`Colors`"];
Off[General::spell1];

Mathematica: first steps.

Getting started.

This is a Mathematica notebook. Mathematica notebooks are divided into cells. There are a few things you have to know about cells.

Cells are marked by the blue brackets on the right.

A cell that has a filled blue triangle in its bracket can be opened by double clicking.

A cell(bracket) that contains other cell(bracket)s can be closed by double clicking in the cell bracket.

Exercise 7.
a) Close and open the group of cells you're reading now.
b) Open the next cell and continue.

Text and input.

There are two types of cells:

	text cells
and
	input cells
		
In this course you can recognise these cells from the fonts used. Now you are reading a text cell.

Input cells contain commands that can be send to the kernel of Mathematica. The kernel is the part of Mathematica that does the computations. 

You can evaluate input cells by putting your cursor somewhere in the input cell and pressing

	<shift> <enter>
		
Try to evaluate the contents of the next cell.

4 5

Exercise 8.
a) What is the meaning of the blank?
b) Change the command into 

	4+5
	
and evaluate it.
c) Find out how to create new cells. Compute

	 (4^5)
	 
in a new cell.

Defining  functions.

Note that the kernel of Mathematica knows only about things you have told  her first. So, in the following,  evaluate input cells whenever you encounter them.

A function is defined in the following manner.

f[x_]  :=  E^x

Pay attention to the _  in the left hand side of

	f[x_] := E^x.
             
It is needed  in the definition, but after that you should never use it! In the sequel one just types  f[x],  f[t] or f[34] or whatever.

Exercise 9.
Compute f(a),  f(0),  f(1),  f(-∞),  and f(∞).
Hint: all Mathematica commands start with a capital letter. So guess what the command for ∞ is.

Plotting  functions and Help.

Plotting functions is as easy as defining them.

Plot[{Sin[x], Cos[x]},  {x, 0, 2 Pi}];

Plot is the command that does the job for us. It needs at least two arguments, a function and the domain of the function. In our example the domain is [0,2π]. Instead of giving Plot one function, we gave it a list of functions.

When you want to know more about a command, you can read Help. A complete manual is available via the menu. A short explanation is given by the command 

	?
	
Try it!

?Plot

Exercise 10.
a) Plot 

	f(x) = x (e^(-x))
	
on the domain [0,10]
b) Plot
 
	g(x) = (x^6)-5(x^3)+sin x
	
c) g(x) has three local extrema. Create a picture that shows them by adapting the domain.

Options for  functions.

Mathematica has its own idea of how a nice picture should look. Sometimes this has unwanted and unexpected side effects. Consider

Plot[3  x, {x, -1, 1}];

Note that the scale on the x-axis and y-axis differs. We can repair this by adding an option to the command. It has the mystereous name AspectRatio.

Plot[3 x,  {x, -1, 1},  
        AspectRatio  ->  Automatic];

Automatic is the value of the option. Adding this option forces Mathematica to use the same scale on the x- and the y-axis. Another useful option is PlotRange. It tells Mathematica what range we want to see.

Plot[3  x,  {x, -1, 1}, PlotRange  ->  {0, 1}];

When you want a list of all possible options, type ??. For example,

?? Plot

Exercise 11.
Plot the function f(x)= sin x + x  on the intervals  [0,24π] and [0,1000π], first in the standard way, then with option PlotPoints -> 40. What is your conclusion?

Variables.

As  every other programming language, Mathematica allows the use of variables. The normal way of assigning a value to a variable is like this.

a = 2

Once we have evaluated this expression, Mathematica knows about a. Again use `?',  when you want to know what Mathematica knows about a variable or function.

?a

When you want to get rid of a definition you can use the command Clear.

Clear[a]

?a

When you want to get rid of all definitions, use

Clear["@"]

Exercise 12.
a) Plot 

	f(x) =a x (e^(-x))
	
together with

	g(x) = x
	
on the domain [0,10], for a = -1, 1, e and (e^2).

More on Mathematica.

Thoughout this course, we will start every notebook with some comments on how to use Mathematica. When you want more of an overview: at the end of the printout of the notebooks we added an appendix that summarizes Mathematica's most important concepts and comments.

A bifurcation diagram.

States.

We continue our analysis of the behaviour of the recurrence relation

	x(n+1) = a x(n) (e^(- x(n)))		(1.1)
		
We use the computer in order to investigate how the state changes when we increase n. We start by telling Mathematica about our relation. We call the right hand side f and define it as a function.

f[x_]  :=  a x E^(-x)

Exercise 13.
a) Choose a value for the parameter a.
b) Choose an initial value, say x = 10. 
c) Verify that we can execute one step of the recurrence relation (1.1)  by putting 

	x = f(x)
	
Try this in Mathematica.
d) Use Mathematica's commands Do and Print in order to print the first ten states for your choice of the parameter a and initial value x. 
Hint 1: use the Help facilities!
Hint 2: The command N gives a numerical approximation of the number its given.
e) Did you find a steady state? If not, reduce the value of a.

Cobwebbing.

In this paragraph we do exactly the same as in the last one, but now we represent the steps graphically. The tool that we use is cobwebbing. Let's start by recalling the recipe for cobwebbing, as given in the course notes.

Recipe 1.5.1. Draw in one figure the (45^∘) line y = x and the graph of f, i.e. the curve y = f(x), for relevant values of x. Choose (Subscript[x, 0]) on the x-axis. Go vertically to the graph of f, then horizontally to the (45^∘) line, then vertically to the graph of f, then horizontally to the (45^∘) line, ..., etc. ...

Note that we have done the first step of the recipe in exercise 12.

We define f (again?!), choose some value for a and plot f.

f[x_]  :=  a x E^(-x)
a = 20;
plot = 
        Plot[{f[x], x}, {x, 0, 8}, AspectRatio  ->  Automatic];
x0  =  1.5;

Some startvalue, like (Subscript[x, 0]) = 1.5, has to be chosen as well.

We have written a small program that performs one cobweb step at a time. By one step we mean going horizontally and vertically once...

\(nexthook := \ 
      Module[
		{}, 
		y0 = f[x0]; 
		markx0 = 
          Graphics[{Red, PointSize[0.015], Point[{x0, 0}]}]; 
		hook = 
          Graphics[{Blue, 
              Line[{{x0, x0}, {x0, y0}, {y0, y0}}]}]; 
		plot = {plot, 
            markx0, hook}; 
		x0 = y0; 
		Show[plot];
		]\)

Now try it, by evaluating nexthook repeatedly.

nexthook

Exercise 14.
Check your answers to the theoretical exercises by cobwebbing. That is, in each of the following cases, characterized by a range of values to which the parameter a must belong, make a cobweb picture that illustrates the statement.
a) The trivial steady state is stable.
b) The trivial steady state is unstable.
c) The non-trivial steady state is stable.
d) The non-trivial steady state is unstable.
e) An attracting 2-cycle.

Orbits

Cobwebbing is a graphical tool. Graphical operations are time consuming and a picture that is getting too complex doesn't really help. Therefore we continue our analysis using numbers only, for the moment.

For this purpose we recall the definition of orbit. Fix some initial value x0. The sequence of points {x(n)} in the state space is the orbit through x0.

Exercise 15.
a) Show that

	x(n) = f(f(....f(x0)...))
	
where f is successively applied n times, on the right hand side.
b) What happens with the orbit in the case you happen to start in a steady state?

We can use the line of Mathematica commands that we produced in Exercise 13c. for investigating orbits. Do you agree?

c) Start by defining a=7.4 again. You should be able to see the 2-cycle appear. Now gently increase a. For what value of a does the first change in  behaviour appear?

Hint: Mathematica keeps answers exact whenever it is able to do so. This can lead to strange/unwanted results. You can avoid this by feeding it a numerical number, like 10.0 instead of an exact number like 10.

d) How would you characterize the change of attaractors as observed in c)?

Chaotic behaviour.

Exercise 16.
a) Increase a further. How do the orbits look?

We have seen that the behaviour of solutions of the recurrence relation (1.1) varies with the value of a. For small positive a we find steady states. When we increase a, we find 2-cycles and even 4-cycles. For bigger values of a, the behaviour of the recurrence relation becomes chaotic. Therefore it seems a good idea to plot a diagram of the orbits of  (1.1)  as a function of a.  We let  ln(a) increase from 0 to 4.4,  fix  (Subscript[x, o])= 0.1 and plot x(n) for 40 ≤ n ≤ 90 for every value of ln(a) bigger then 2.5.

\(oneorbit[lna_, interval_]\  := \ 
      Module[
		
		{a\  = \ E^lna, 
		\ \ x0\  = \ 
            0.1, 
		\ \ f, \ 
          xn, \ 
		\ \ from\  = \ interval[\([1]\)], 
		\ \ to\  = \ 
            interval[\([2]\)]\ \ \ \ \ \ \ \ }, 
		
		f[x_]\  = \ 
          a\ x\ E^\((\(-x\))\)\ ; 
		xn = 
          Nest[f, x0, from]; 
	\ \ \ \ \ \ \ \ \ \ NestList[f, xn, 
          to - from]
]\)

\(bifurcation[stepx_]\  := \ 
      Module[
		
		{begin\  = \ 0, 
			end\ \ \ \ \  = \ 
            4.4, 
			from\  = \ 40, 
			to\  = \ 
            90, 
			ll, 
			plotdata\  = \ {}, 
			orbit, \ 
          i\ \ \ \ \ \ \ \ \ \ \ \ \ 	}, 
	
	 (*\ 
          The\ first\ piece\ is\ a\ straight\ line\ *) 	
	
        plotdata\  = \ 
          Table[{lna, lna}, {lna, begin, 1.9, 2\ stepx}]; 	
	
	 (*\ 
          Here\ we\ have\ a\ two - cycle\ *) 
	
        For[lna = 1.9, \ lna\  < \ 2.5, \ 
          lna\  = \ lna\  + \ stepx, 
	\ \ orbit\  = \ 
            oneorbit[lna, \ {from, from + 1}]; 
		AppendTo[\ 
            plotdata, {lna, orbit[\([1]\)]}]; 
		AppendTo[\ 
            plotdata, {lna, orbit[\([2]\)]}];
\ \ ]; 
		
	 (*\ 
          This\ is\ the\ chaotic\ part\ *) 
	ll\  = \ to - from + 1; 
	
        For[lna\  = \ 2.5, \ lna\  <= \ end, \ 
          lna\  = \ lna\  + \ stepx, 
		orbit\  = \ 
            oneorbit[lna, \ {from, to}]; 
		plotdata\  = \ 
            Join[\ plotdata, \ 
              Table[{lna, orbit[\([i]\)]}, {i, 1, 
                  ll}]];
\ \ ]; 
		
\ \ \ ListPlot[plotdata\ , 
          Axes\  -> \ True, \ 
          PlotRange\  -> \ {\(-1\), 20}, 
	\ \ Background\  -> \ 
            RGBColor[0, 0, 0], \ 
	\ \ \ DefaultColor\  -> 
            RGBColor[1.0, 0.9, 0], 
	\ \ \ Prolog\  -> \ 
            AbsolutePointSize[1.4]\ ]
]\)

bifurcation[0.01];

This will give you some idea of the enormous orbit complexity that is possible for even such innocent looking recurrence relations.

b) What is the meaning of the parameter stepx in the procedure bifurcation?
c) What will be the increase in the number of computations when you halve the value of this parameter? 
d) Choose a value for this parameter in such a way that your picture looks
nice, but doesn't take ages.
e) Do you see the period doubling route to chaos?

To surprise you even further we mention a remarkable theorem, due to Sharkovskii.

Theorem. Assume that f is  a continous map of an interval to itself. Assume that f has a period p orbit. For any q such that p < q, f has a period q orbit, where the ordening is given by

3 < 5 < 7 < 9 < ... < 2.3 < 2.5 < ...  < (2^2).3 < (2^2).5 < ... < (2^3).3 < (2^3).5 < ... < (2^4) < (2^3) < (2^2) < 2 < 1

In particular the existence of a period 3 orbit implies the existence of a periodic orbit of any period!