UC SCI 111, January 23, 2004
Notebook 1

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.

Cells may be grouped together, with a grouping bracket added on the far right to show the scope of the group. Cell groups may be opened and closed by double clicking on this grouping bracket. Closed cell groups have a filled blue triangle on the top of the grouping bracket.

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

Text and input.

There are many types of cells. Two common cell types are:

	text cells
and
	input cells
		
In this course you can recognise these cells from the fonts used. Now you are reading a text cell. Text cells use the "Times" font, and are used for writing text.

Input cells use the boldface "Courier" font, and contain commands that can be sent to the kernel of Mathematica. The kernel is the part of Mathematica that does the computations. 

In its simplist sense, Mathematica is like a calculator. You can evaluate input cells by putting your cursor somewhere in the input cell and pressing

	<shift> <enter>

You may also evaluate an input cell by pressing the  <enter>  key on the lower right corner of your keyboard.
 
Mathematica displays the output in an output cell, using the plain "Courier" font.

Try to evaluate the contents of the next cell.

4*5

Exercise 2.
a) What is the meaning of  *  in the above input?
b) Change the  *  to a space and re-evaluate the cell. You see that there are two ways to carry out multiplication.
c) Change the command into 

	4+5
	
and evaluate it.
c) Slowly move the cursor above this (red) cell, but below the output cell above until the cursor becomes a horizontal bar. Click once. A horzontal bar should appear across the entire notebook. If you type from the keyboard now, you will create a new input cell at the location of the horizontal bar. Create a new input cell, and use the "Palette" on the right side of your screen to type 

	 (2002^30)
	 
in this cell. Evaluate the cell.
d) Experiment with the so-called "Basic Input" Palette on your own. Note that it contains the mathematical constants π (approximately 3.14), e (approximately 2.72), and i (the square root of -1) on the first row of tiny buttons near the middle of the palette. It also contains ∞ (infinity).  There are many greek letters on the palette, but only π carries a numerical value. You may also type directly "Pi" for π, "E" for e, "I" for i, and "Infinity" for ∞.

Using commands.

Most of the thousands of built-in Mathematica commands are used by typing their names directly into an input cell. As a first example, consider the command whose name is N. It is used for numerically approximating an exact number. To see it in action evauluate these cells:

Pi

N[Pi]

N[Pi, 500]

Here is an example of another built-in command, this one for the sine function:

Sin[Pi/12]

The following conventions are utilized for all built-in commands:

1) The names of built-in commands begin with a capital letter.
2) The argument(s) for commands are contained in square brackets, and if there is more than one argument, they are separated by commas.

In this course you will be exposed to a few dozen commands. They will be introduced as they are needed.

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_]  := Exp[x]

Pay attention to the _  in the left hand side of

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

Exercise 3.
Compute f(a),  f(0),  f(1),  f(-Infinity),  and f(Infinity).

Plotting  functions and Help.

Plotting functions is as easy as defining them.

Plot[Sin[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 the interval [0, 2π]. One may plot several functions at once by typing a list of functions for the first argument:

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

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

	?
	
Try it:

?Plot

Exercise 4.
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-axes. 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

Mathematica has a very complete online help system as well. It is accessed by the Help menu on the top right of your window. Highlight the name of a command on your screen that you wish to learn about (an option for the Plot command in the above output, for instance), and then choose "Find in Help Browser..." in the Help menu. You can use the Help Browser to look up any and all Mathematica related items. A good general strategy is to first push the "Master Index" radio button, then type a keyword into the search field and push the "Go To:" button. Then follow a hyperlink to the explanatory text. Note that a complete edition of The Mathematica Book is contained in the Browser.

Exercise 5.
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

One very common source of errors in using Mathematica is forgetting to clear a variable when you are done with it. You may have forgoten that a=2, but Mathematica has not.

When you want to get rid of all definitions, use

ClearAll

You are also permitted to assign values to entire words, such as

age = 21

You must exercise caution, however! Recall that a space indicates multiplication, and the possibility of assignments such as the one made above shows that adjacent letters (without a space) do NOT indicate multiplication. For instance, enter these cells:

α = 2

α x

αx

This last cell contains a single variable with the name αx, which has no assigned value.

Suppose we want to plot (x^3) - a*(Underscript[x, 1])+1 in the interval [-3,3] for various values of a. First carry out the following command:

Plot[x^3 - a*x + 1, {x, -3, 3}]

We do not get a result because we did not give a a value yet. Now assign a value for a as in

a=2

and make the plot again. Do this for different values of a.

Some sample problems.

To get a feel for what you can do with Mathematica we give a few exercises.

Problem 0.

In this subsection we play with iterations, which are important for the first assignment. We consider the function f(x)=a x(1 - x) for various values of a and study what happens if we apply this function repeatedly to a number between 0 and 1. In other words, we start with a number (Subscript[x, 0]) and compute (Subscript[x, 1]) = f((Subscript[x, 0])), (Subscript[x, 2])=f((Subscript[x, 1])),  (Subscript[x, 3]) = f((Subscript[x, 2])), .... In general,

(Subscript[x, n + 1]) = f((Subscript[x, n]))

First we define the function f with a = 1.5

a =  2; f[x_] = a*x*(1 - x)

We can compute the first 50 values of (Subscript[x, n]) with (Subscript[x, 0]) = 0.5 using the command NestList as follows

xn = NestList[f, 0.5, 50]

We can make a plot of these values using the ListPlot command:

ListPlot[xn, PlotRange -> {0, 1}]

Do the above experiment with some other values of a between 2 and 4. In particular we advise you to take some values between 3 and 4.

Problem 1.

The first problem is to determine the smallest positive number p so that tan(p)=p.
Problems of this sort sometimes occur in so-called boundary value problems, such as quantum mechanical particles in a potential well. It is impossible to give an exact answer to this problem and what we like to do is find a value of p which is correct to, say, 6 decimal places.
Let us first make a combined plot of the functions x and tan(x).

Plot[{Tan[x], x}, {x, 4.4, 4.6}, PlotRange -> {-10, 10}]

Of course we know that tan(0)=0, but we are interested in positive solutions p. From the plot above you get an idea of the approximate value of p. It is somewhere between 4 and 5.

Exercise 7.
By repeatedly decreasing the x-interval in the above plot, make an estimate of the value of p to 3 decimal places.

In Mathematica there is a built there are a number of built in functions to solve equations. One of them is called Solve (naturally...). So let us try that

Solve[Tan[x] == x, x]

We do not get an answer. Unfortunately, the problem is too hard for the function Solve. The reason is that tan(x) is a so-called transcendental function. Moreover, tan(x)=x has an infinite number of solutions (do you understand why?). To get an answer we must give Mathematica an extra hint of the whereabouts of the smallest solution. If we do that the function FindRoot does the trick for us in many cases.

FindRoot[Tan[x] == x, {x, 4.5}]

If you want the solution in 20 decimal places you can add the option WorkingPrecision->20. Do this!

Exercise 8.
Determine the second smallest positive number q such that tan(q)=q to 10 decimal places.

Problem 2.

Mathematica is also a symbolic calculator. That is, it can perform operations on formulas. An important operation is differentiation of a function. Here you see the derivative of ((xlog(x))^2).

D[x*Log[x]^2, x]

Exercise 9.
Compute the derivatives of (x/(1 + (x^2))), sin((x^6)) , (x^x), first by hand, and then using Mathematica.

Derivatives can be very handy if you want to compute local maxima and minima of a function. Here we take a simple example. Define the function f(x) by

f[x_] = ((x^4) - 2*(x^2) + x + 1)/((x^4) + 1)

Exercise 10.
Find the number of local maxima and minima of f(x) by making a number of suitable plots and make an estimate, to 1 decimal place, of their location.

To get a better idea of the location of the local extrema, we must determine the zeros of the derivative of f(x). Here is the derivative which we denote by df(x).

df[x_] = D[f[x], x]

Exercise 11.
Determine the location of the local extrema of f(x) to 10 decimal places by determining the zeros of df(x).

Problem 3.

Let us consider the motion of a valve on a bicycle wheel. This motion consists of two components, the forward speed of the bicycle and the circular motion of the wheel. Suppose at time t=0 the valve starts in its lowest position and the bicycle move to the right with speed v m/s. Suppose also that the radius of the wheel is r meters. We assume that the motion of the bicycle wheel takes place in the xy-plane, the x-axis representing the ground. Suppose also that at t=0 the valve is at the origin.

Convince yourself that the motion of the wheel is given by

x(t)=v t-r sin(vt/r)
y(t)=r-r cos(vt/r)

By rescaling units we can see to it that r=v=1. So let us define

x[t_] = t - Sin[t]
y[t_] = 1 - Cos[t]

Here is a plot of the orbit of the valve

ParametricPlot[{x[t], y[t]}, {t, 0, 15}]

This picture is of course vertically distorted. A more realistic picture can be obtained using the option AspectRatio->2/15. Do this!

Now we ask the following question. We see that at time t=2π the valve has returned to the ground once more. What is the total length of the orbit of the valve? Considering the fact that the circumference of the wheel is 2π, can you make a guess? In the following we shall explain and do the computation in Mathematica.

The speed vector is simply the derivative of the vector (x(t), y(t)). We denote the speed in the x-direction by vx, the speed in the y-direction by vy.

vx[t_] = D[x[t], t]
vy[t_] = D[y[t], t]

Of course the value of the speed is the length of the vector with components (vx,vy).

speed[t_] = Sqrt[vx[t]^2 + vy[t]^2]

To find the length of the valve's orbit we have to integrate the speed from t=0 to t=2π.

Integrate[speed[t], {t, 0, 2*Pi}]

You probably would not have guessed this answer. Here is a further exercise

Exercise 12.
Determine the length of that part of the orbit which is above the wheel's axis, i.e. the part of the orbit where y>1.

Here is another question. When the bicycle starts moving at time t=0, the valve will leave the origin at a certain angle.

Exercise 13.
You may guess that the valve takes off at a vertical direction. Either give a good argument for this, or use  the speed vector we just found to verify.

And now a completely different curve, the logarithmic spiral, which is already present in the work of Archimedes. The curve is given by the parametrisation

x[t_] = Exp[-t/10]*Cos[t]
y[t_] = Exp[-t/10]*Sin[t]

where t runs from 0 to ∞.

Exercise 14.
a) Make a plot of the logarithmic spiral and compute its length.
b) Instead of  (E^(-t/10)) we take (E^(-a*t))where a is an arbitrary positive number. What is the length of the new spiral in terms of a?

More on Mathematica.

Thoughout this course, we will start every notebook with some comments on how to use Mathematica.

One can also find a nice introduction to Mathematica by clicking on the following link to a Tour of Mathematica. It will take you to a file in the Help Browser.

A bifurcation diagram.

States.

We continue our analysis of the behaviour of the recurrence relation

	\(TraditionalForm\`x\_\(n + 1\) = α\ 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 (note the spaces on the right side in the definition of  f  below):

Clear[α, f, x]

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

Exercise 1.
a) Choose a value for the parameter α, say α = 2. Make the assignment in Mathematica.
b) Verify that we can execute one step of the recurrence relation (1.1)  with the intial value (Subscript[x, 0]) = 1/4 by typing 

	Nest[f, 1/4, 1]
	
Try this in Mathematica.
c) What happens when you increase the value of the last argument to 2?
Hint: Remember that you can type ?Nest for information on this command.
d) Investigate the command NestList[f, 1/4, 3]. Interpret what you get in terms of the recurrence relation (1.1).
e) You will have noticed by now that Mathematica will not round its output (if it is given exact input). However, exact output is sometimes difficult to read! One could use the N command to get a numerical approximation of the result. Another way to get an approximation is to use approximate input. Replace the intitial value 1/4 with .25 in the command NestList[f, 1/4, 3]. Mathematica regards decimal numbers as approximate. Using a single approximate number as input will force Mathematica to use floating point arithmetic (which is also much faster) in all intermediate calculations, and to report an approximate (decimal) answer.
f) Use the NestList command to approximate the values (Subscript[x, 0]),  (Subscript[x, 1]),  (Subscript[x, 2]),   ...  ,  (Subscript[x, 30]) for the recurrence relation (1.1) when α = 2 and (Subscript[x, 0]) = 1/4. Does the sequence appear to converge?

Cobwebbing.

In this section 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?!), and choose some value for α, and do it all in a single input cell. We use a semicolon to finish each line of input.

Clear[α,  f,  x];
α = 2;
f[x_]  :=  α x (E^(-x))

Some start value, like (Subscript[x, 0]) = 0.25, has to be chosen as well. We now use the Cobweb command to get a picture of what happens when we apply the recurrence relation (1.1) repeadedly. The syntax for Cobweb is precisely like the syntax for Nest and NestList. The red dots on the x-axis are the values produced by the NestList command.

Note: Cobweb is not a built-in command. Rather, it was programmed into Mathematica for this course, and will work on the machines in this lab. If you own your own copy of Mathematica, you will need to ask your instructor for a copy of the "package" that defines this command before it will run on your computer.

Cobweb[f,   .25,  1]

Cobweb[f,   .25,  2]

Cobweb[f,   .25,  3]

Cobweb[f,   .25,  30, AspectRatio -> Automatic]

For fun, one may animate the progression with the command CobwebMovie. It will produce one "frame" of a movie for each iteration. Double click on the grouping bracket of the graphics frames when it finishes, then double click on the remaining picture to start the movie. Controls for the animation will appear in the lower left of the notebook window once the animation begins.

CobwebMovie[f,   .25,  30]

Exercise 2.
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 α 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 (Subscript[x, 0]). The sequence of points {(Subscript[x, n])} in the state space is the orbit through (Subscript[x, 0]). Throughout this section we will keep the initial value (Subscript[x, 0]) fixed. It's value is not particularly important, as we will see. We wish to investigate the stability of the orbits as the parameter α varies. For this purpose, we are interested in the values (Subscript[x, n]) when n is large. We can use the command NestList for investigating orbits. Do you agree?

Exercise 3.
a) Enter the cell below, which assigns the value 7.4 to the parameter α.

Clear[α,  f,  x];
α = 7.4;
f[x_]  :=  α x (E^(-x))

b) Use the commands Nest and NestList in a single line of Mathematica code to produce the values (Subscript[x, 100]), (Subscript[x, 101]), (Subscript[x, 102]), ... , (Subscript[x, 120]), for the initial value (Subscript[x, 0]) =  .25.
c) Produce the values (Subscript[x, 10000]), (Subscript[x, 10001]), (Subscript[x, 10002]), ... , (Subscript[x, 10020]), for the initial value (Subscript[x, 0]) =  .25. Comment on the stability of the orbit.
d) Repeat part (c) using four different initial values (Subscript[x, 0]). What do you conclude?
e) Now repeat part (c) again, but this time keep (Subscript[x, 0]) fixed (say at .25, although this is not crucial) and gently vary the parameter α above. Start with the value α = 2.0, and work up to α = 14.0. Report those values of α where the nature of the stability of the orbit through (Subscript[x, 0]) changes, and comment on the nature of the orbit.

Chaotic behaviour.

Exercise 4.
a) Increase α 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 α. For small positive α we find steady states. When we increase α, we find 2-cycles and even 4-cycles. For bigger values of α, 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 α.   The command we use is called BifurcationPlot.

?BifurcationPlot

b) Let α increase from 2 to 30,  fix (Subscript[x, 0]) =  .25, and plot (Subscript[x, n]) for 100 ≤ n ≤ 200 using the BifurcationPlot command. You will need to clear the variable α before using this command, for α will now need to vary.

Clear[α]