A Mathematica Tutorial

Prepared by Dan Schroeder for use by members of the WSU Physics Department, 31 January 2007

This tutorial is intended to quickly introduce some of the handiest features of Mathematica, using examples taken from physics.  To use this tutorial you should actually type all the input into Mathematica and see what output it produces.

1.  Plots and Integrals

Let's start by plotting a function:

Plot[x^3/(Exp[x]-1),{x,0,10}]

This example illustrates several of Mathematica's important syntax conventions:
* Hit shift-return or the keypad "enter" key to send the input to Mathematica
* Regular parentheses ( ) are for grouping
* Square brackets [ ] enclose function arguments
* Curly brackets { } enclose elements of a "list" or array
* Names of built-in functions such as Plot and Exp start with capital letters (Mathematica is case sensitive)
* Arithmetic operators are +, -, *, /, and ^.  The * for multiplication can be omitted when there's no ambiguity, as we'll see in later examples.

Now let's try an integral:

Integrate[x^3/(Exp[x]-1),{x,0,10}]

Yuck!  Mathematica gave an answer in terms of PolyLog functions, which I know nothing about.  So let's just ask for a numerical answer:

N[Integrate[x^3/(Exp[x]-1),{x,0,10}]]

The "N" function says "Give me a numerical answer, with the default working precision."  (If you need higher precision, you can get it.)  But with integrals, it's usually much faster to use the NIntegrate function, which does the calculation numerically from the start:

NIntegrate[x^3/(Exp[x]-1),{x,0,10}]

Instead of typing the same expression over and over, we can give it a name.  Use a lower-case first letter to be sure that your name won't conflict with that of a built-in function:

planck = x^3/(Exp[x]-1)

The = symbol is an assignment operator, like in Fortran or C, but in Mathematica, the values of variables are (in general) algebraic expressions, not numbers.  Now we can use "planckSpectrum" as much as we like in subsequent calculations:

NIntegrate[planck,{x,0,Infinity}]

Notice that "Infinity" is a built-in function so it starts with a capital letter.  Very few built-in function names are abbreviated.  With this upper limit, however, the exact answer is much simpler:

Integrate[planck,{x,0,Infinity}]

The real fun starts when you embed an NIntegrate inside a Plot.  So, for instance, with Mathematica you can calculate and plot the Debye heat capacity in a single line:

Plot[3t^3 NIntegrate[x^4 Exp[x]/(Exp[x]-1)^2,{x,0,1/t}],{t,0,2}]

Notice that Mathematica whines a bit about the singularity at t=0, but displays the graph correctly nevertheless.  (20 years ago it would have taken me about a day to produce this plot.  Now it takes only a few seconds!)

Mathematica documents are called notebooks (.nb).  Here are a few general tips for working with them:  
	* To abort a calculation, use the Abort command under the Kernel menu.  
	* To collapse and expand groups of cells, double-click on the enclosing brackets at the right margin.  
	* To add a comment cell like this, or title or section heading,, click where you want it and before you type anything, choose the desired style from the Format menu.

2.  Root Finding, Maximization, Sums, Derivatives

Besides integration, another common numerical task is root finding.  Here's a plot of a pair of functions that intersect in three places:

Plot[{Tanh[2x],x},{x,-2,2}]

To find the nontrivial intersection near x=1, use the FindRoot function.  Note that the double == sign indicates a test for equality, as in C:

solution = FindRoot[Tanh[2x]==x,{x,1}]

The answer is returned as a list containing a "rule" to replace x with the desired value.  To check the result, we can evaluate tanh(2x) "at the point" specified by the rule:

Tanh[2x] /. solution

The next two instructions will make more sense later.  Here's Mathematica's "full" internal representation of the solution:

FullForm[solution]

And now that we know this representation, we can figure out how to extract the root itself as the 2nd part of the 1st part of the expression, just as if it were an element of an array:

solution[[1,2]]

If you can't remember how to use a built-in function, you can type it with a question mark in front.  The "More" link takes you to the online help, which you can also access via the Help menu:

?FindRoot

(The online help is great for looking up little things, but for extended reading I prefer the hard copy Mathematica book.)

To find the maximum value of a function, use FindMaximum:

FindMaximum[planck,{x,3}]

(There's a similar FindMinimum function.)  Of course, you could also differentiate the function and then use FindRoot.  The derivative function is simply D:

deriv = D[planck,x]

FindRoot[deriv == 0,{x,3}]

Another commonly used function is Sum:

Sum[1/n^2,{n,1,10}]

Sum[1/n^2,{n,1,Infinity}]

Sum[1/n^s,{n,1,Infinity}]

Of course you can embed a Sum inside a Plot.  Here's a truncated Fourier series approximation to a square wave:

Plot[Sum[Sin[k*x]/k,{k,1,25,2}],{x,0,2Pi}]

Notice the step size of 2 in the list of k values to include.  Try increasing the upper limit!

3.  Lists

Mathematica includes all sorts of functions for creating, manipulating, and plotting lists.  The most common way of creating a list is the Table function:

factorials = Table[n!,{n,1,10}]

To extract an element of a list, follow the list's name with the element number in double-square-brackets:

factorials[[5]]

Here's another list, containing the Stirling approximation values to the first 10 factorials:

stirling = N[Table[n^n*Exp[-n]Sqrt[2Pi* n],{n,1,10}]]

(Try out the previous expression without the enclosing N[ ]!)

Basic arithmetic operations act on each element of the list separately.  So, for instance, here's how you can calculate the fractional error in Stirling's approximation for all 10 values at once:

fracError = (factorials-stirling)/factorials

(There's a separate dot product operation, a.b or Dot[a,b], where a and b are lists that are treated as vectors or matrices.)

Now let's construct a two-dimensional list, or a list of lists.  Here each sub-list will be the ordered pair (n,fracError[[n]]):

pairs = Table[{n,fracError[[n]]},{n,1,10}]

Sometimes it's nice to display a two-dimensional list as a matrix:

MatrixForm[pairs]

To extract an element of a matrix, put the row and column numbers in double-square-brackets, separated by a comma:

pairs[[5,2]]

(Equivalently, we could say pairs[[5]][[2]].)

Very often you'll want to plot a set of ordered pairs on a graph.  The function you need is ListPlot:

ListPlot[pairs,Frame->True,PlotRange->{{0,11},{0,0.09`}},PlotStyle->{Hue[0.75`],AbsolutePointSize[5]}]

To make the plot more attractive, I've specified a number of "options":  Frame, PlotRange, and PlotStyle.  These and many other options can be specified in any order, separated by commas, using the syntax "option -> value" (note the hyphen-greaterthan combination to produce the arrow).  Here's how you can learn the names of all the available options, and learn how to use any of them:

Options[ListPlot]

?Frame

4.  Reading and Writing Data

The other thing you often want to do with lists is read and write them from/to external files.  First, it's a good idea to find out what's the current working directory:

Directory[]

On my machine the result is "/Users/dschroeder", which is the Mac-Unix path for my home folder.  Let's change the working directory to my desktop folder, which is immediately within the home folder:

SetDirectory["~/Desktop"]

Now let's write the "pairs" list to an external file called "pairs.txt":

Export["pairs.txt",pairs,"Table"]

Here the final quoted "Table" tells Mathematica what format to use when writing the file.  Open the file in a text editor to take a look!

To read data from a file, use the Import function, again with the "Table" option (assuming the data is in rows and columns, separated by white space):

importedData = Import["data.txt","Table"];

The final semicolon in the previous line suppresses Mathematica's output, which would be long and rather ugly.  A better way to display the data is to use MatrixForm or TableForm:

TableForm[importedData]

More often you'll simply want to plot the data:

dataPlot=ListPlot[importedData,PlotStyle->Hue[0.7`]]

Finally, let's overlay a theoretical curve to try to model some features of the data.  I'll use a sum of a linear "background" function and a Gaussian, with all the constants chosen "by eye" to match what I want to match.  The option DisplayFunction->Identity suppresses the output of the Plot itself.  I then use the Show function to combine dataPlot and curvePlot:

curvePlot=Plot[7.55`-0.0048` x+2.14` Exp[-((((x/18)))^2)],{x,-100,100},PlotStyle->{AbsoluteThickness[2],Hue[0]},DisplayFunction->Identity];
Show[dataPlot,curvePlot]

5.  More Graphics

Another very common task is to plot a two-dimensional function of some parameter.  For example, plotting a pair of sine functions produces what's called a Lissajous figure:

ParametricPlot[{Sin[2 t],Sin[3 t]},{t,0,2 π},AspectRatio->Automatic,GridLines->Automatic,Axes->False,Frame->True,PlotStyle->{Hue[0],AbsoluteThickness[2]}]

(Try removing the various options above, to see what effect each of them has.)  Here's a parametric plot of an elliptical orbit, where e is the eccentricity (which you should try changing):

a=1;epsilon=0.9`;
r=((a (1-(epsilon^2)))/(1+epsilon Cos[theta]));
ParametricPlot[{r Cos[theta],r Sin[theta]},{theta,0,2 π},AspectRatio->Automatic,PlotStyle->{Hue[0.75`],AbsoluteThickness[2]}]

Now let's plot a function of two variables.  Here's the formula for one of the hydrogen 3d wavefunctions (actually the probability density), sliced through the xz plane:

h3d = (2z^2-x^2)^2 *Exp[-2Sqrt[x^2+z^2]/3]

My favorite way to display such a function is with a DensityPlot:

DensityPlot[h3d,{x,-20,20},{z,-20,20},Mesh->False,PlotPoints->500,PlotRange->All,ColorFunction->Hue ]

(Try it without the options, then add them back one at a time.  The color is pretty but doesn't really help visualize the function, in my opinion.)  You can also display the same information with a ContourPlot:

ContourPlot[h3d,{x,-20,20},{z,-20,20},PlotPoints->500,PlotRange->All,ContourShading->False,Contours->12]

(Again, try it without the various options.)  Finally, you can display the function as a surface plot using Plot3D:

Plot3D[h3d,{x,-20,20},{z,-20,20},PlotRange->All,Mesh->False,PlotPoints->200,Boxed->False,Axes->False]

If you create a list of related plots, Mathematica lets you animate them as frames in a movie.  As an example, here's the time-dependent solution to the one-dimensional Schrodinger equation for a Gaussian wavepacket moving to the right.  This example also demonstrates that Mathematica can handle complex numbers, with "I" representing the square root of -1, and the Re and Im functions extracting the real and imaginary parts.

k=10; 
wavepacket = Exp[-k^2/4]*Exp[(I*x+k/2)^2/(1+2I*t)]/Sqrt[1+2I*t]

Table[Plot[{Re[wavepacket],Im[wavepacket]},{x,-2,14},PlotPoints->500,PlotRange->{-1,1},PlotStyle->{{Hue[0]},{Hue[.6]}},AspectRatio->.3,Ticks->None,ImageSize->600],{t,0,1,.01}];

Normally you'll want to "collapse" the sequence of plots by double-clicking on the enclosing bracket at the right margin.  To animate the display you can simply double-click on the first graphic, or select them all and choose "Animate Selected Graphics" from the Cell menu.  Notice that when the animation is running, you can control it with the buttons that appear at the lower-left corner of the window frame.

6.  Symbolic Operations

Mathematica is best known for its symbolic manipulation capabilities.  Most of these take a lot of getting used to, so they don't lend themselves to a quick tutorial.  But here's an example of a symbolic calculation in which we start with the partition function for a two-state system and calculate first the average energy and then the heat capacity:

z = 1 + Exp[-beta*epsilon]

Oops!  If you've worked through the examples above, then the symbol "epsilon" already has a numeric value, left over from when we used it for the eccentricity of the elliptical orbit.  Here we want to treat epsilon as simply an undefined symbol.  So we need to erase its previous value with the Clear function.  While we're at it, let's also erase the value of k (which we'll use in a moment):

Clear[epsilon,k]

Now go back and re-enter the expression for z, the partition function.  To get the average energy we differentiate with respect to beta, then multiply by -1/z:

energy = -(1/z)D[z,beta]

This expression can be simplified.  While we're at it, let's replace beta with 1/(k t):

energy = Simplify[-(1/z)D[z,beta]]/. beta->1/(k*t)

Finally, differentiate with respect to t (temperature) to get the heat capacity:

heatcap = D[energy,t]

Now what?  Well, we could plot the result, but first we'd better set epsilon and k equal to 1 (in appropriate units):

Plot[heatcap/.{epsilon->1,k->1},{t,0,3}]

7.  Differential Equations

Many of the laws of physics are expresed as differential equations.  Mathematica can solve differential equations symbolically, but I rarely use this feature.  On the other hand, I frequently use the numerical differential equation solver, NDSolve.  As an example, let's solve the differential equation for a simple pendulum.  The function we want to find is theta[t], and Newton's law says that the second derivative of this function is proportional to minus its sine.  I'll use units in which the constant of proportionality (g/L) is equal to 1.  We also need to specify two boundary conditions for a second-order equation, so I'll take the initial value of theta to be some constant (initially small), and the initial value of the angular velocity to be zero.  The syntax of NDSolve requires that we provide three arguments:  a list of the equation(s) and boundary condition(s), the name of the function to solve for, and the usual list containing the independent variable and its initial and final values:

solution=NDSolve[{theta''[t]==-Sin[theta[t]],theta[0]==0.1` π,theta'[0]==0},theta,{t,0,50}];
Plot[theta[t]/.solution,{t,0,50},PlotStyle->{AbsoluteThickness[2],Hue[0]}]

Notice that I used a semicolon to suppres the output of the NDSolve function itself, and then immediately told Mathematica to plot the solution that it found.  (The solution itself is something called an interpolating function, which stores a bunch of numerical values and rules for interpolating between them.  You can create such objects from lists of data if you ever need to.)

Try increasing the value of theta[0] to see how the solution changes!

If we want to know a particular point on the solution curve, we can use FindRoot:

FindRoot[(theta[t]/.solution)==0,{t,5}]

Here's another differential equation example, this time solving the time-independent Schrodinger equation for a harmonic oscillator potential.  Although the exact solution is known in this case, the method can be used just as easily with other potential energy functions.  In that case, you'd need to guess values of the energy e until the function goes to zero at +infinity.  This is called the "shooting method" of finding eigenvalues.

e=1.5`;xmax=5;
solution=NDSolve[{((psi'))'[x]==(-2 e+(x^2)) psi[x],psi[-xmax]==0,psi'[-xmax]==0.001`},psi,{x,-xmax,xmax}];
Plot[psi[x]/.solution,{x,-xmax,xmax}]

8.  Concluding Remarks

Mathematica can do lots of other things:  arbitrary-precision arithmetic, sophisticated algebraic manipulations, linear algebra, procedural programming with conditions and loops, other types of graphics, image processing, GIS analysis, typesetting of mathematical equations, and interacting with other software.

To learn more about Mathematica, I highly recommend the Mathematica book.  There's a copy of the book in the Physics Majors room.  (It's an old edition, which is actually good in most respects--less cluttered with new features that you don't need.)  I keep a slightly newer edition in my office, which I'll gladly lend out.  The book begins with a quick tour to show you some of the things Mathematica can do.  Then it gets down to business, with four major sections.  The first section, A Practical Introduction, is worth reading sequentially, with a computer at hand to try things out.  Serious users will eventually want to read much of the next section (Principles) as well.  The third and fourth sections (Advanced Mathematics and Reference Guide) are mostly for reference.

Let me end by mentioning some of the things that Mathematica is not.  It is not a replacement for pencil and paper, nor is it a replacement for thinking carefully and using your common sense.  (Just as you shouldn't believe every answer that appears on your pocket calculator, you shouldn't believe every answer that you get out of Mathematica.  In both cases, it's not that the machine makes mistakes--it's that your instructions to the machine may not accomplish what you intended.)

Nor is Mathematica a replacement for all other computational environments.  For some tasks it's easier to use a spreadsheet, or a traditional programming language, or another high-level package like Matlab.  Mathematica's built-in numerical functions are quite speedy, but Mathematica's performance is rarely acceptable for computationally intensive tasks where you need to program your own low-level algorithms.  Although Mathematica is capable of producing beautiful illustrations and typeset equations, there are far better programs for producing publication quality graphics and page layouts.  Like other sophisticated tools, Mathematica is really good at some things and merely adequate at others.  Learn what it can do, and put it to work for you, but don't expect too much of it!

It's also worth mentioning that Mathematica is a commercial product with a rather high price tag.  WSU pays a substantial annual fee for a site license that allows enough copies to be running at once to meet the needs of the Math and Physics Departments.  Faculty and Staff can also obtain home copies under this license.  The WSU Bookstore will sell Mathematica to students at a fairly affordable price, but not to others.  Aside from the high monetary price, there are other costs associated with using a commercial, proprietary system.  You can't see the internal algorithms that Mathematica uses, or change the way it works at a low level to suit your purposes.  The software's copy protection is highly annoying when you need to reinstall it or transfer it to a different computer.  If you create complicated Mathematica documents and want to preserve the capability of modifying them in the future, you're pretty much committed to buying upgrades to retain compatibility with new hardware and operating systems.  If you want to share your work with others, they can read it for free but can't modify it without their own copies of the software.  Sometimes it might be better to use a free, open-source computing environment just to avoid all these hassles.