Tutorial 1: Introduction To Mathematica

Introduction To Calculus, Carleton College Math 111 01-02, Fall 2007, Prof. Joshua R. Davis

revised Thursday 2007 September 13, 2:58 PM

Introduction

Mathematica is a computer algebra system, which means that it's a computer program capable of sophisticated mathematical computation and symbolic manipulation. There are several similar products on the market with varying capabilities. We will be using Mathematica frequently this semester, instead of hand calculators, so you are expected to develop a basic familiarity with it. On the other hand, do not rely on it too much; you are not allowed to use it on exams.

This short tutorial interweaves explanatory text with actual Mathematica commands for you to execute. I strongly recommend that you execute every command as directed; also try some commands on your own, to experiment. Nothing in this tutorial needs to be handed in, but you will have to use similar techniques on your homework.

Mathematica As A Calculator

Mathematica is essentially a sophisticated calculator. You type in a mathematical expression and hit the Enter key; then the program evaluates the expression for you and returns the answer. For an example, I have typed in

3 + 5

Click on that expression with your mouse and then hit Enter (not Return). Mathematica pauses a little bit to start up its computational system ("kernel") and then promptly spits out an answer of 8. It also labels the "3 + 5" as an input and the "8" as an output, to help you keep things clear. In fact, you can enter several inputs at once, and it will give you the output for each one. For example, click on this next chunk of input and press Enter:

-Pi * 37
3^4

You can see that * indicates multiplication and ^ indicates an exponent. The symbol Pi is the mathematical constant 3.14159... that you've seen before. Notice that Mathematica does not give you  - Pi * 37 as a decimal number, because the decimal expansion is infinitely long and messy, and Mathematica keeps everything exact by default. To get a decimal number, use the N function:

N[-Pi * 37]

When you use the N function, you are saying to Mathematica, "Compute the stuff inside the square brackets, and then give me a numerical approximation." Whenever you use functions in Mathematica, you always put the input in square brackets. Here is sin(0):

Sin[0]

The authors of Mathematica chose to use square brackets for function input instead of parentheses because parentheses are used for so many different things in mathematics; requiring square brackets for functions makes the notation a little easier to process in software. Also, Sin is capitalized because all built-in Mathematica functions (such as N and Sin) and constants (such as Pi and E) are capitalized. Now try

Sin[1]

Mathematica is not very helpful here, because the exact value of sin(1) is not a nice number that it can print out explicitly. Perhaps you want a numerical approximation:

N[Sin[1]]

Notice how I've nested the Sin function inside the N function. Mathematica can nest functions as deeply as you like. In mathematics, this nesting is usually called composition.

By the way, Mathematica operates in radians --- and so do we in Math 111 --- because radians are so much better-behaved, mathematically, than degrees, gradians, etc. We'll talk about that in class.

Mathematica can also do algebra. Try

(1 + x)^3

It formats the expression nicely, but it doesn't know that I really want the thing multiplied out. For that, I need to use the Expand function:

Expand[(1 + x)^3]

Mathematica can also simplify complicated algebraic expressions.

Simplify[x^2 + (1 - Cos[x]^2)/Cos[x]^2 - x^5 * x^(-3)]

Naming Things

On most calculators there is a memory function that lets you store the result of a computation for later use. Mathematica has a very friendly storage facility that lets you give each result a descriptive name. For example, here is the number of seconds in a year (how did I compute it?):

secondsperyear = 60 * 60 * 24 * 365

That command does two things: it computes the number 60 * 60 * 24 * 365, and then stores the answer under the name "secondsperyear". This is handy if I'm doing a problem that uses this number more than one time. For example, the stars Proxima Centauri and Sirius are 4.2 and 8.6 light-years from Earth, respectively. How many light-seconds away are they? (A light-year is the distance that light travels in a year. A light-second is the distance that light travels in a second.) To figure it out, we multiply the light-years by the number of seconds in a year:

4.2 * secondsperyear
8.6 * secondsperyear

Mathematica reports the two answers in scientific notation, since they are so big. The name "secondsperyear" that I chose was entirely up to me; I could have used the name "x" or "jimmy", but I didn't because those are not as descriptive. Furthermore, simple names like "x" are a bad idea because we often want to use "x" as a variable. Suppose that we carried out the preceding example using "x" instead of "secondsperyear":

x = 60 * 60 * 24 * 365
4.2 * x
8.6 * x

We get the same answers as before, which is fine. But now we go off to work on some other problems, during the course of which we want to expand out the algebraic expression (1 + x)^3. We've already seen how to do that:

Expand[(1 + x)^3]

What went wrong? Why didn't we get 1 + 3 x + 3 x^2 + x^3, as in the preceding section? Well, Mathematica no longer treats x as a variable, because we've declared that x has the specific value of 60 * 60 * 24 * 365. So Mathematica is computing (1 + 60 * 60 * 24 * 365)^3, which is apparently 31363163502213982608001. We want Mathematica to forget that x was ever assigned a value --- to go back to having it as a variable. To this end, use the Clear function:

Clear[x]
Expand[(1 + x)^3]

If Mathematica is ever behaving strangely for you, it's probably because some of the symbols you're using have been assigned values that you've forgotten about. Use the Clear function liberally to avoid such problems (or save your work, quit Mathematica, open your notebook again, and activate only the inputs you need).

The ability to store things is especially useful when dealing with data. Suppose I have a bunch of student test scores and I want to compute the mean (that is, the average) and the median (the middle score). I can enter the scores in a list by enclosing them in curly brackets; then I give that list the name "grades":

grades = {76, 80, 64, 42, 91, 56, 54, 59, 49, 76, 75, 60, 76, 61, 54, 90, 50, 51, 64}

There are built-in functions to compute the mean and median:

Mean[grades]
Median[grades]

Notice that Mean pops out the exact answer. If you want a numerical approximation, what function do you use?

All names in Mathematica are case-sensitive. The command

expand[(1 + x)^3]

does nothing useful, because Mathematica has no idea what the "expand" function is. Remember that all built-in Mathematica functions and constants begin with upper-case letters. To avoid confusion, when making up your own names you should always begin them with lower-case letters. That's why we chose the name "secondsperyear" instead of "Secondsperyear" or "SecondsPerYear", and the name "x" instead of "X".

Similarly, if you mess up and use parentheses for the function input, then Mathematica doesn't know what you're talking about:

Expand((1 + x)^3)

In fact, in the new version (version 6) of the software, it colors the parentheses red to indicate a probable error.

At this point you may be thinking that Mathematica is a rather finicky program. That's true; you have to tell it exactly what you want. The upside of this is that you can tell it surprisingly subtle things. Mathematica is in fact a fully capable computer programming language, with a vast library of built-in functions. To learn more about its capabilities, check out the Help menu or the web. We're going to talk a little about graphing, and then get on to stuff like the homework.

Two Kinds Of Graphing

Before running any other commands in this section, please first enter these two. They just make our graphs easier to see; you aren't expected to know how they work.

SetOptions[Plot, PlotStyle->Thickness[Large]];
SetOptions[ListPlot, PlotStyle->PointSize[Large]];

There are two basic graphing functions that we're going to use a lot. The first one, Plot, is used to graph functions defined by a formula. For example, the following command plots the function y = (x / 14 - 130)^2, for values of x between 1800 and 2000.

Plot[(x/14 - 130)^2, {x, 1800, 2000}]

The second graphing function is used to plot data. You can give it a list of points, where each point is a pair of numbers (x, y), and it will show all of them in a scatter plot. To demonstrate this, I first set up a list giving the U.S. population (in millions) for each decade from 1900 to 1990 (adapted from http://www.npg.org/facts/us_historical_pops.htm):

mydata = {{1900,76},{1910,92},{1920,106},{1930,123},{1940,132},{1950,152},{1960,181},{1970,205},{1980,227},{1990,249}}

I then use the ListPlot function to plot the data:

ListPlot[mydata]

Notice that Mathematica draws the x-axis at y = 100, not y = 0. Mathematica automatically makes choices like this in an attempt to produce attractive figures. You can turn it off, but let's not worry about that right now; just keep an eye on it.

Frequently in this tutorial we will want to plot some data, then fit a curve to the data, then plot the data and the curve together. To show two plots together, store each plot in a symbol of your choosing, and then use the Show function:

mycurveplot = Plot[(x/14 - 130)^2, {x, 1800, 2000}]
mydataplot = ListPlot[mydata]
Show[mycurveplot, mydataplot]

Together, these plots look pretty silly, because they have nothing to do with each other. In the following sections, however, we choose our curves to fit the data nicely.

Crickets

Problem 1.2 #22 in your textbook concerns crickets chirping at various temperatures (Fahrenheit). As in the preceding section, we enter the data as a list of points. Then we plot them with the ListPlot function, and store the result in dataPlot.

tempchirps = {{50, 20}, {55, 46}, {60, 79}, {65, 91}, {70, 113}, {75, 140}, {80, 173}, {85, 198}, {90, 211}}
dataplot = ListPlot[tempchirps]

The data look like they lie on a line (roughly), so we're going to fit a line y = m x + b to them and plot the line on top of the data. Here's one way to do it, for m = 4 and b = -100. First we clear the names m, b, and x, just to be safe. Then we use some weird "/." syntax. Don't worry about how it works; just know that its effect is to plug the given values of m and b into m x + b. Then we plot that line and show it with the data.

Clear[m, b, x]
myline = m * x + b /. {m -> 4, b -> -100}
lineplot = Plot[myline, {x, 0, 100}]
Show[dataplot, lineplot]

The line y = 4 x - 100 doesn't seem to pass through the points very well. Before going on any further, edit that last input to try different values for m and b, and see whether the line improves. Keep trying different values until you get a line that fits the data pretty well.

Have you played around until you have good values for m and b? No? Then try some more.

Now have you found yourself good values for m and b? Okay, then let's learn how to make Mathematica itself find m and b. It takes several steps.

First we clear the symbols m, b, and x, just to be safe.

Then we invoke the FindFit function, telling it (1) what data to fit, (2) what sort of function to fit to them, (3) what parameters it's supposed to play with, and (4) what the independent variable is. This function computes the best possible choices for m and b. We give meaning those choices the name "bestmb".

(By the way, the best choices is the ones that minimize the least-squares error, which is the most common measure of error in statistics. What math lies at the heart of this error-minimization process? Calculus! More on that later in the term.)

In the next command, we plug that m and b into m x + b using the same /. syntax as before. In the final two commands, we plot the fit line and display it along with the data, as before.

Clear[m, b, x]
bestmb = FindFit[tempchirps, m * x + b, {m, b}, x]
myline = m * x + b /. bestmb
lineplot = Plot[myline, {x, 0, 100}]
Show[dataplot, lineplot]

By inspecting the output, we can see that FindFit chose m to be 4.85667 and b to be -220.967. The plot shows that this m and b describe a really good fit line. We conclude that if y is the number of chirps and x is the temperature, then the data are well-described by y = 4.86 x - 221. Then we can ask follow-up questions such as, "At what temperature does chirping stop?" or "How much chirping is there when the temperature is 1000000 F?"

Fluoxetine

The drug Fluoxetine is marketed under the brand name Prozac. On Day 0 a patient is injected with this drug. The concentration of the drug in her blood is measured over the next week, as it is slowly metabolized in the brain and liver. Here are the recorded concentrations for days 0 to 7, in ng/ml:

drugdata = {{0,21.6},{1,13.2},{2,8.64},{3,5.28},{4,3.24},{5,1.98},{6,1.29},{7,0.91}}
dataplot = ListPlot[drugdata]

This is certainly not a line, but it might be an exponential curve. To assess that possibility, we're going to take the (natural) logarithm of the concentration values. We do it using that weird "/." syntax again. Don't worry too much about what it means. Just read the command and believe me that it takes the drug data, replaces each concentration with its logarithm, and stores the new data list under the name "semilogdrugdata".

semilogdrugdata = drugdata /. {day_, concentration_} -> {day, Log[concentration]}

The term "semilog" is used when we take the log of the y-values but not of the x-values. A plot of the semilog data is called, in standard statistics terminology, a semilog plot.

semilogdataplot = ListPlot[semilogdrugdata]

On the semilog plot, the data seem to lie on a straight line. So let's have Mathematica fit a line to the these data, just as we did earlier for crickets. Notice that I've adjusted x to run from 0 to 7. The result is quite good.

Clear[m, b, x]
bestmb =  FindFit[semilogdrugdata, m * x + b, {m, b}, x]
myline = m * x + b /. bestmb
lineplot = Plot[myline, {x, 0, 7}]
Show[semilogdataplot, lineplot]

So far, we've fitted a line to the semilog data. What does this tell us about the original data? Let's think about this carefully.

Remember that the semilog plot doesn't show how x relates to y; it shows how x relates to log y. We found a line with slope m = -0.460776 and intercept b = 3.04925. But this does not mean that y = -0.460776 x + 3.04925; it means that log y = -0.460776 x + 3.04925, because the semilog plot relates x to log y, not to y!

So what we know is that log y = -0.460776 x + 3.04925. Then we just do a little algebra to solve for y. (This is done in class; make sure you know how to do it.) We get

  y = e^3.04925 * e^(-0.460776 x).

This is an exponential function. Let's plot it, and show the plot with the original data plot (not the semilog plot).

curveplot = Plot[y = E^3.04925 * E^(-0.460776 * x), {x, 0, 7}]
Show[dataplot, curveplot]

The point of this exercise is that a line on a semilog plot corresponds to an exponential curve on the original plot. To fit an exponential curve to a set of data, all we must do is fit a line to their semilog data.

The exponential fit curve lets us answer many follow-up questions. For example, perhaps you are familiar with the concept of half-life --- the amount of time needed for half of a substance to decay away. (This concept exists only for exponential decay.) In the fluoxetine data, the concentration on day 0 is about 22. Half of that is about 11. Inspecting the exponential fit curve, we see that the concentration is about 11 after about 1.5 days. So the half-life of fluoxetine is about 1.5 days.

AIDS Cases

Here is data about the number of AIDS cases in the U.S. for the time period 1981-1992 (from E. K. Yeargers, R. W. Shonkwiler, and J. V. Herod, An Introduction to the Mathematics of Biology: with Computer Algebra Models, Birkhäuser, Boston, 1996). Each datum consists of a year (the year of the epidemic, where year 1 is 1981) and the number of cases up to and including that year.

data = {{1,97}, {2,709}, {3, 2698},{4,6928},{5,15242},{6,29944},{7,52902},{8,83903},{9,120612},{10,161711},{11,206247},{12,257085}}
dataplot = ListPlot[data]

The data don't plot on a straight line. Let's try a semilog plot. Make sure you understand these commands; they are very similar to those used in our Fluoxetine example.

semilogdata = data /. {year_, cases_} -> {year, Log[cases]}
semilogdataplot = ListPlot[semilogdata]

That doesn't look like a straight line either. So let's try something different. Let's take the log of the year and the cases, to produce what's called a log-log plot:

loglogdata = data /. {year_, cases_} -> {Log[year], Log[cases]}
loglogdataplot = ListPlot[loglogdata]

That's pretty close to a straight line. So let's fit a line to it and graph the fit line. It is important to remember that a log-log plot does not show a relationship between the original x and y, or between x and log y; it shows a relationship between log x and log y. To make sure I remember this, I call my independent variable logx, not x. Notice that x ranges from 1 to 12, so log x ranges from Log[1] to Log[12]. I have adjusted my Plot command to use this range of values.

Clear[m, b, logx]
bestmb = FindFit[loglogdata, m * logx + b, {m, b}, logx]
myline = m * logx + b /. bestmb
lineplot = Plot[myline, {logx, Log[1], Log[12]}]
Show[loglogdataplot, lineplot]

We got m = 3.27409 and b = 4.41522. Remember that this is a log-log plot, so what we've found is that log y = 3.27409 (log x) + 4.41522, not y = 3.27409 x + 4.41522. Again doing some algebra to isolate y (this is done in class), we get

  y = e^4.41522 * x^3.27409.

This is a power function. Let's plot it against the original data.

curveplot = Plot[E^4.41522 * x^3.27409, {x, 1, 12}]
Show[dataplot, curveplot]

The point of this exercises is that a line on a log-log plot corresponds to a power function graph on the original plot. To fit a power function to a set of data, all we must do is fit a line to its log-log data.

Other Fit Curves

We have seen three kinds of fits: linear, exponential (which looks linear in a semilog plot), and power law (which looks linear in a log-log plot). These are just examples; there are many data sets that don't fit into any of these three categories. Consider Problem 1.2 #24. If anything, those data look piecewise-linear. There are a best linear fit to them, a best exponential fit to them, and a best power law fit to them, but none of these models is excellent.

When a scientist is presented with data, the most appropriate kind of model is often not obvious. The scientist chooses a model based on her other knowledge of the system that produced the data, or based on a desire for a simple model, or based on trial and error. After she publishes her data with a particular model, another scientist may attempt to model the same data in a more sophisticated way. There are always better models to be had. When you model data in this class, try to choose a model that's better than others, but don't worry about finding a model that captures every detail; then you may never finish your homework.