An Introduction to Mathematica
			©Sally Otto (2005)

Beginning to use Mathematica

To start, go to Help and click "Tutorial."

Entering commands

Let's first practise using the following commands:

    Factor[f]  - Factors the function f
    D[f, x]   - The derivative of a function f with respect to x
    Integrate[f, x]  - The indefinite integral of f with respect to x
    Integrate[f, {x,a,b}]  - The definite integral of f with respect to x between a and b

For example, let's start by factoring ((1 + x - 6  (x^2))/(1 - 2  x)).  Place the cursor in the following command and hit "enter" (or "shift-return"):

Factor[(1 + x - 6 x^2)/(1 - 2  x)]

In this equation, spaces are used to indicate multiplication and ^ is used to indicate a power (see "Short Guide to Mathematica Functions").  Mathematica should answer 1 + 3 x.  

You will get the same answer if you write the above fraction in a prettier form using the input palette (get this toolbox by clicking File -> Palettes -> BasicInput).  Again, place the cursor in the following command and hit enter (do this each time you see a command in bold):

Factor[((1 + x - 6  (x^2))/(1 - 2  x))]

Notice that these sentences are written in a different font than the input and output commands.  When entering a new line, you can change the format to text by clicking Format -> Style -> Text.  You can change it back to input by Format -> Style -> Input.

Also notice that each equation that Mathematica enters is given a number.  For example, the first command that you entered was labeled In[1], and the answer was labeled Out[1].  This is convenient because it allows us to refer to previous commands and manipulate them without retyping them.  For example, if Out[1] was 1 + 3 x, then we can factor (1+3 x) + 2 (x^2) by typing:

Factor[Out[1]  +  2 (x^2)]

As a short-cut, you can refer to n\(TraditionalForm\`\^th\) output as either Out[n] or %n

Factor[%1 +  2 (x^2)]

If you leave out the number, % refers to the last output.  For example, we can expand the previous result using:

Expand[%]

Mathematica really consists of two pieces of software:  the window (which we see) and the kernel (the brains behind it all).  The numbering of the input/output is done by the kernel.  When you close Mathematica, everything in the kernel is forgotten.  The next time that you open Mathematica, it will have forgotten what Out[1] was.  Mathematica starts counting your input/output all over again.

Now let's turn to derivatives.  First let's find ((d(4 (x^3)))/dx)by entering the following command:

D[4 (x^3), x]

Mathematica should give 12 (x^2), which is the derivative of 4 (x^3)with respect to x.  Similarly, we can find the derivative of Exp[4 (x^3)] with respect to x by entering:

D[Exp[4 (x^3)], x]

If you want to evaluate your answer at a particular point (e.g., x = 1), you can use the replacement command  /. x -> 1 
Try making this replacement in the above derivative.

To write the answer as a real number, use the N[ ] command

N[%]

Anti-derivatives are found using the command Integrate.  For example, we can undo the above differentiation with respect to x using:

Integrate[12 (x^2) Exp[4 (x^3)], x]

The above represents an indefinite integral.  We can use the input palette to write this in a more familiar form:

∫12 (x^2) Exp[4 (x^3)] \[DifferentialD]x

We can also use Mathematica to find definite integrals.  For example, the area under the curve 1+(x^2) from x = 0 to 1 is given by:

Integrate[1 +  (x^2), {x, 0, 1}]

The term in curly braces {} represents the variable, the lower bound, and the upper bound, respectively.  Mathematica typically using curly braces {} to denote lists of information, square brackets [ ] to denote functions, and normal parentheses ( ) to denote terms that should be kept together.

Practise factoring

Use Mathematica to factor the following equations:

	1 + 3  x + 2  (x^2)

	2 - 5 x - 12 x^2
	 
Check all of the answers that Mathematica gives by hand.  Caution:  You have to place the cursor below the area of this text in order to get a new Input line (notice the blue lines down the right-hand side of the window, which tell you what Mathematica treats as a single group of lines).

Practise taking derivatives

Use Mathematica to take the derivatives of the following equations with respect to x:

	1 + 3  x + 2  (x^2)

	 x^3 y^2

	 x^3 y[x]^2
	 
Check all of the answers that Mathematica gives by hand.  In the last two cases, notice that Mathematica gives you different answers.  In the first case, Mathematica treats y as a constant parameter (we didn't tell it otherwise).  If, however, y depends on x, then we have to tell Mathematica that y is a function of x by writing y as y[x].

Practise finding indefinite integrals

Use Mathematica to take the anti-derivatives of the following equations with respect to x (i.e., find the indefinite integrals):

	1 + 3  x + 2  (x^2)

	 x^3 y^2

	 x^3 y[x]^2
	 
Check all of the answers that Mathematica gives by hand (use the fact that (1 - x^2)/(1 - x) factors).  Again, Mathematica gives you different answers in the last two cases.  In the last case, Mathematica knows that y is a function of x, but it does not know what function of x it is, so it cannot evaluate the integral.

Practise finding definite integrals

Use Mathematica to find the integral of the following functions with respect to x between x=0 and x=1:

	1 + 3  x + 2  (x^2)

	(1 - x^2)/(1 - x)

	 x^3 y^2

	 4 (x^3)
	 
Check all of the answers that Mathematica gives by hand.

Identifying the right command

The key to using Mathematica is being comfortable finding out how to do things that you don't already know.

For example, if you can't remember the notation required to integrate a function, you can ask Mathematica for help, either by browsing the help (click Help -> Help Browser) or entering a question mark followed by the command of interest:

?Integrate

If you can't remember exactly what the name of a function was, you can use an * as a wildcard.  For example, the following returns all of the commands that start with Int:

?Int*

Mathematica always capitalizes the first letter of any function that it knows.  Thus, the following request for help doesn't work:

?int*

In the following exercises, try to guess the right name for the command and determine the right notation for use in Mathematica:

Plotting functions

Plot a graph of 1 + 3  x + 2  (x^2)with respect to x.

Plot a 3D graph of 1 + 3  x + 2  x ywith respect to x and y.

Drawing random numbers

Use Mathematica to draw a random real number between 0 and 1.

Use Mathematica to draw a random integer between 0 and 99.

Handling lists

Mathematica is also able to handle lists of information (e.g., vectors and matrices), which are indicated by {}

mylist = {1, 2, 3, 5, 7, 11};

The ;  at the end just tells Mathematica not to output the information to the screen.  Try deleting the ; and reentering.

We can make a list using the Table function (or a For loop):

mynewlist  =  Table[1/i, {i, 1, 10}]

To plot a list, we use ListPlot

ListPlot[mynewlist]

Altering your plots

Figure out how to join the dots in the previous graph.

Figure out how to change the range of your plot.

Make this line dashed by adding the following to the end of your command, right before the last ]
	, PlotStyle -> Dashing[{0.04}]

Making lists

Make a list of 50 random numbers using the Table command.

rantab=

You can find a least-squares fit to the data using Fit:

Fit[rantab, {1, x}, x]

Mathematica actually has the ability to load packages with more specialized commands.  For example, let's load the package on Linear Regression commands:

You can fit a regression line through these points using Regress:

<<Statistics`LinearRegression`

Regress[rantab, {1, x}, x]

There is also a non-linear regression package that you can load (look up information on NonlinearFit).

Writing your own functions

You can define your own functions in Mathematica using the definition command  :=
For example, the following function defines exponential growth:

Clear[expgrow,n0,r,t]

expgrow[n0_,r_,t_] := Exp[r*t]*n0

A few points to emphasize:

(1) It is good to get in the habit of erasing any previous definitions of the parameters and variables using Clear, to avoid accidently using previous values.
(2) It is a good idea to name your own functions using lower case letters to avoid any conflicts with definitions that Mathematica already knows.
(3) The parameters or variables critical to evaluating the function are written inside the [ ], where each one is followed by a _ to indicate that any possible value of the parameter or variable is acceptable.

Enter the following commands to see the difference that the underscore makes:

test[x]  :=  x^2
test2[x_]  :=  x^2

test[x]

test2[x]

test[3]

test2[3]

Defining a function

Define a function, y, that equals sine(x) + x and plot this function from x = 0 to 3.

Solving recursion equations

Mathematica can solve simple recursion equations, like the model of geometric growth:

RSolve[{n[t] == a n[t - 1], n[0] == n0}, n[t], 
      t]

Try altering this code for logistic growth, n[t] = n[t-1] + r n[t-1] (1 - n[t-1]/K), which does NOT have a general solution.

Solving differential equations

Mathematica knows the general analytical solution to some (but by no means all) differential equations.

The procedure is very similar, except we use DSolve rather than RSolve.

The DSolve command we would use is therefore:

DSolve[{D[n[t], t] == a*n[t],n[0]==n0}, {n[t]}, t]

Try to solve the differential equation for the logistic equation dn/dt = r n (1-n/K), which can be solved in continuous time.

Writing simulations in Mathematica

Mathematica can be used as a programming language in which you can run simulations.  However, there are a few tricks you need to learn about writing
simulations in Mathematica.

Say that you have the recursion equation:

	x[t] = a x[t-1]			(exponential growth)
	
where x[t] is a function of one parameter, a, and x[t-1], so you could simply write the function:

	x[t_,a_] = a*x[t-1,a] 
	
THIS WON'T WORK THOUGH!!  It sends Mathematica into an infinite loop because Mathematica tries to find x at time t as a function of x at time t-1, but it doesn't know that either, so it tries to find x at time t-2, etc.

It would be OK to use

	x[t_,a_] := a*x[t-1,a] 
	
The := tells Mathematica not to calculate anything until you tell it to.

The second step would be to set the initial conditions.  Say:

	x[0,a_] = 10

Remember, the underline (_) tells Mathematica to make the definition for all
values of a.

You could then ask Mathematica what x is at any time t for particular values of a:

	x[10, 1.1]
	x[100,1.1]

THIS TURNS OUT TO BE REALLY SLOW!!

The problem is that Mathematica does not remember previous values that it has calculated. For every time t, Mathematica goes back to the recursion to
figure out x all the way up from 0 to t.

To make Mathematica remember past values that it has calculated, we use a trick:

	x[t_,a_] := x[t,a] = a*x[t-1,a] 
	
The first time Mathematica sees x[10, 1.1] it goes to the recursion to calculate it, but the second part (x[t,a] =) then sets x[10, 1.1] forever more to the value on the far right.  Now, every time x[10, 1.1] is called in the future, Mathematica remembers it.

Geometric growth example

Clear[x,a,t]
x[t_,a_] := x[t,a] = a*x[t-1,a]
x[0,a_] := 1

(I always Clear all the parameters and variables first.)

x[10,1.1]

x[100,1.1]

ListPlot[Table[x[t,1.1],{t,0,100}]]

You can't use Plot, which only works for continuous functions (we only know x at the discrete time points of each generation).

Many more graphics commands are available if we input the extra Graphics packages:

<< Graphics`Graphics`

After entering the above, use DisplayTogether to plot a couple of growth trajectories using different parameters.

Write a simulation

Write a recursion equation for the logistic model n[t] = n[t-1] +r n[t-1] (1-n[t-1]/K)
Use your code to graph logistic growth when n[0] = 100, K = 1000, for various values of the intrinsic growth rate, r.

Solving differential equations numerically

The way that Mathematica numerically solves differential equations (e.g., dn/dt = r n (1- n/K)) is a trade secret, but it is said to involve 500 pages of code.  Roughly, Mathematica determines what n must be as a function of time to have the right slope, dn/dt, at each point in time.

To numerically solve a differential equation, we use NDSolve:

NDSolve[{equations you want to hold true},
		{variables you want to solve for}, {t,tmin,tmax}]

where

	tmin = the first time point that you want to know about
	tmax = the last time point that you want to know about
	
For example, consider a model of exponential growth, where the habitat declines in quality over time, so that dn/dt = (r0 - delta t) n.  Given an initial population size, we have two equations we want to hold true:

	D[n[t],t] == (r0-delta*t)*n[t]
and	
	n[0] == n0
	
and we want to solve for one variable (n[t]).

NOTE:  The == is used to denote an equation that we want to be true.  We're not setting the right hand side to the left hand side.

The NDSolve command we would use is:

NDSolve[{D[n[t], t] == (r0 - delta*t)*n[t], n[0] == n0}, {n[t]}, {t, tmin, tmax}]

Mathematic won't solve this yet, since it doesn't know what the parameters are.  What I do is define a new function of all of the parameters and set this function to the above NDSolve:

Clear[habitat, n]
habitat[r0_, delta_, n0_, tmin_, tmax_] := habitat[r0, delta, n0, tmin, tmax] = 
   NDSolve[{D[n[t], t] == (r0 - delta*t)*n[t], n[0] == n0}, {n[t]}, {t, tmin, tmax}]

[NOTE:  habitat[r0_,delta_,n0_,tmin_,tmax_] := 
	habitat[r0,delta,n0,tmin,tmax] =
notation ensures that Mathematica remembers the numerical solution that it finds.]

When we call habitat with all the parameters, Mathematica will return its best numerical guess for the function, n[t], between tmin and tmax (Mathematica calls this an interpolating function).

habitat[0.1, 0.01, 1000, 0, 100]

We can then plot n[t] by making a substitution of the approximate function habitat for n[t].

For example, with r0=0.1, delta =0.01, n0=1000, tmin=0, and tmax=100:

Plot[n[t] /. habitat[0.1, 0.01, 1000, 0, 100], {t, 0, 100}, 
      PlotRange -> {{0, 100}, All}]

Writing more involved code

Simulations and more complicated code can be written together within a Block, which keeps track of parameters and variables locally.  For example, the following code calculates chi-squared statistics for 2-by-2 table of integers in the form (({{a, b}, {c, d}})).

<< "Statistics`ContinuousDistributions`"

contingencyanalysis[{a_, b_, c_, d_}] := 
      Block[{r1, c1, tot = a + b + c + d, expected, 
          chisquared}, 
r1 = (a + b)/
            tot; 
c1 = (a + c)/
            tot; 
expected = {r1 c1 tot, 
            r1 (1 - c1) tot, (1 - r1) 
              c1 tot, (1 - r1) (1 - 
                  c1) tot}; 
chisquared = 
          Plus @@ (({a, b, c, d} - expected)^2/
                expected); 
{chisquared, 
          1 - CDF[ChiSquareDistribution[1], chisquared]}]

First, figure out what Plus @@ does to a list of numbers.

contingencyanalysis[{5, 6, 1, 10}]

To get the actual number we would enter:

N[%]