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. (You can also indicate multiplication by just putting a space between two items.) The symbol Pi is the mathematical constant 3.14159... that you know. 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

Here, 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]

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

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 "jimmy", but I didn't because that's not descriptive.

Now try this:

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

Why didn't the second command produce 1 + 3 x + 3 x^2 + x^3, as it did above? Well, Mathematica no longer treats x as a variable, because I have declared that x has the specific value of 3. So Mathematica is computing 4^3, which is 64. If you want to go back to having x be a variable, then 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, restart Mathematica, 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", 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.

Graphing

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, using a thick line style.

Plot[y = (x/14 - 130)^2, {x, 1800, 2000}, PlotStyle -> Thick]

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):

popData = {{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, using a large point size:

ListPlot[popData, PlotStyle -> PointSize[Large]]

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 do this, store each plot in a symbol of your choosing, and then use the Show function to show them together:

curvePlot = Plot[y = (x/14 - 130)^2, {x, 1800, 2000}, PlotStyle -> Thick]
dataPlot = ListPlot[popData, PlotStyle -> PointSize[Large]]
Show[curvePlot, dataPlot]

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, PlotStyle ->  PointSize[Large]]

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:

fitLine = 4 x - 100
linePlot = Plot[y =fitLine, {x, 0, 100}, PlotStyle -> Thick]
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. First we clear the symbols m, b, and x.

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 returns the best possible choices for m and b, meaning those 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 these values for m and b into the formula m x + b (using the weird /. syntax). We call the result fitLine.

In the final two commands, we plot the fit line and display it along with the data, as before.

Clear[m, b, x]
fitMAndB = FindFit[tempChirps, m x + b, {m, b}, x]
fitLine = m x + b /. fitMAndB
linePlot = Plot[y = fitLine, {x, 0, 100}, PlotStyle -> Thick]
Show[dataPlot, linePlot]

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 you 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 by the brain and liver. Here are the recorded concentrations for days 0 to 7, in ng/ml:

fluoxData = {{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[fluoxData, PlotStyle ->  PointSize[Large]]

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. Read the first command below, and try to figure out what it's doing. What it's doing is replacing each datum (day, concentration) with a new datum consisting of the same day and the log of the concentration. We store this log-data in logFluoxData. Then we plot that log-data, and name the plot logDataPlot, since it's the plot of the log-data.

logFluoxData = fluoxData /. {day_, concen_} -> {day, Log[concen]}
logDataPlot = ListPlot[logFluoxData, PlotStyle -> PointSize[Large]]

By the way, this kind of plot, where the x-data is plotted against the log of the y-data, is usually called a semilog plot. In this case, the data seem to lie on a straight line. So let's have Mathematica fit a line to it and plot the line. The result is quite good.

Clear[m, b, x]
fitMAndB =  FindFit[logFluoxData, m x + b, {m, b}, x]
fitLine = m x + b /. fitMAndB
linePlot = Plot[y = fitLine, {x, 0, 7}, PlotStyle -> Thick]
Show[logDataPlot, linePlot]

So if y is the log of the concentration and x is time in days, then the data are modeled by y = -0.46 x + 3.05. This means that the concentration itself is e raised to -0.46 x + 3.05. Let's graph that against the original data:

curvePlot = Plot[y = E^(-0.46 x+3.05), {x, 0, 7}, PlotStyle -> Thick]
Show[dataPlot, curvePlot]

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 only exists 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.

aIDSData = {{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[aIDSData, PlotStyle -> PointSize[Large]]

The data don't plot on a straight line. Let's try a semilog plot:

logAIDSData = aIDSData /. {year_, cases_} -> {year, Log[cases]}
logDataPlot = ListPlot[logAIDSData, PlotStyle -> PointSize[Large]]

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

logLogAIDSData = aIDSData /. {year_, cases_} -> {Log[year], Log[cases]}
logLogDataPlot = ListPlot[logLogAIDSData, PlotStyle -> PointSize[Large]]

That's pretty close to a straight line. So let's fit a line to it and graph the fit line. Notice that the range of x-values we use to graph the line is the range of x-values in the log-log data.

Clear[m, b, x]
fitMAndB = FindFit[logLogAIDSData, m x + b, {m, b}, x]
fitLine = m x + b /. fitMAndB
linePlot = Plot[y =fitLine, {x, Log[1], Log[12]}, PlotStyle -> Thick]
Show[logLogDataPlot, linePlot]

If y is the log of the number of cases and x is the log of the year of the epidemic, then we have found y = 3.27 x + 4.42. Exponentiating both sides we get that the number of cases is about e^(3.27 x + 4.42) = e^4.42 (e^x)^3.27 = e^4.42 t^3.27, where t = e^x is the year of the epidemic. So let's graph this power function against the original AIDS data:

curvePlot = Plot[y = E^4.42 * x^3.27, {x, 1, 12}, PlotStyle -> Thick]
Show[dataPlot, curvePlot]

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.