Some Basic Mathematica(Introduction)

Simple Mathematica calculations

In your notebook, type

   (1 + 2) * 3

and press Enter. You’ll get the result:

   9

You can omit the multiplication sign, and Mathematica is smart enough to know that it is implied:

   (1 + 2) 3
   9

 

But Mathematica can do a lot more than simple integer math. Try

   102 / 9

You get the result

   34 / 3

You see that Mathematica reduced the fraction to its simplest form by cancelling out common factors in the numerator and denominator. Also, notice that the result is still expressed as a fraction, not as a decimal number. This is because the fraction is more exact than any decimal approximation to it, and Mathematica tries to maintain as much accuracy as possible unless you specify otherwise.

 

To get a numerical approximation to this result, type

   N[102 / 9]

You’ll get the result

   11.3333

The N[…] command tells Mathematica to evaluate the quantity in brackets numerically. If we want more decimal places, we can ask Mathematica to calculate the same thing to 20 digits:

   N[102 / 9, 20]
   11.333333333333333333

 

At this point, it’s useful to introduce a Mathematica shortcut that lets us take the previous output cell and include it in the current input cell. Enter

   102 / 9

and you’ll get the standard Mathematica result:

   34 / 3

Then enter in a new input cell

   N[%, 20]

and you’ll get the result

   11.333333333333333333

just as before. The % symbol tells Mathematica to insert the output of the previous calculation.

 

Mathematica knows the basic mathematical functions like sine and cosine. Typically these functions have the name you would expect, but with their first letter capitalized. To take the sine of 34/3 radians, enter

   Sin[%]

and you’ll get the result

   -0.9434996270154848971

Notice that Mathematica has kept 20 significant figures; it remembers that that was the accuracy of the previous output cell and passes that information along to the current cell.

 

If we try to compute the sine of 34/3 radians from scratch, we don’t get the same result:

   Sin[102 / 9]
   Sin[34 / 3]

Here Mathematica has left both the fraction and the sine function unevaluated. To evaluate this numerically, we can use the N[…] command:

   N[%, 20]
   -0.9434996270154848971

Notice that Mathematica uses square brackets to delimit the argument of a function, instead of the more conventional parentheses. This is because parentheses are used exclusively for grouping in Mathematica, as in the calculation we entered at the beginning of this section:

   (1 + 2) * 3

 

Plotting lists of (x, y) points

We’ve seen that Mathematica uses parentheses for grouping and square brackets in functions. Mathematica uses curly braces to delimit lists of numbers. For example, a list of the first ten prime numbers would be

   {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}

If you enter this as Mathematica input, you’ll get the exact same list as output:

   {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}

This is because we haven’t asked Mathematica to do anything to the list. Later on we’ll see how we can use lists of numbers in calculations.

 

The reason we want to discuss lists now is that Mathematica represents (x, y) points as lists of two numbers. For example the point (x = 1, y = 3.7) is represented in Mathematica as

   {1, 3.7}

 

Suppose we have a list of three (x, y) points: (1, 2), (2, 3), and (3, 5). These points are represented as lists of two numbers. The list of the three points is then a nested list, or a list of lists:

   { {1, 2}, {2, 3}, {3, 5} }

To plot a graph of these three points, we use a new Mathematica command:

   ListPlot[ { {1, 2}, {2, 3}, {3, 5} } ]

This tells Mathematica to plot the list of points inside the square brackets. The spaces are optional and are included here mainly for the sake of clarity.

 

To plot ten points representing the first ten prime numbers, we would therefore type

   ListPlot[ { {1, 2}, {2, 3}, {3, 5}, {4, 7},
    {5, 11}, {6, 13}, {7, 17}, {8, 19},
    {9, 23}, {10, 27} } ]

Again, the extra space is optional. If you type the input as it is shown here, however, you’ll notice that your input cell grows large enough to accomodate all three lines of the command. If you forget a curly brace (or add an extra one), Mathematica will inform you of this. One of the points may be hard to see; look closely at the origin of the x and y axes to convince yourself that there is a point there.

 

Now let’s join the points with straight line segments. Type

   ListPlot[ { {1, 2}, {2, 3}, {3, 5}, {4, 7},
    {5, 11}, {6, 13}, {7, 17}, {8, 19},
    {9, 23}, {10, 27} }, PlotJoined -> True ]

The modifier

   PlotJoined -> True

tells Mathematica to connect the points with lines.

 

Here’s another shortcut that can save you a lot of typing. Suppose we wanted to store the list of the first ten primes in memory, so that we don’t have to type it in every time we want to use it. Enter

   primes = { {1, 2}, {2, 3}, {3, 5}, {4, 7},
    {5, 11}, {6, 13}, {7, 17}, {8, 19},
    {9, 23}, {10, 27} }

Mathematica will respond with

   { {1, 2}, {2, 3}, {3, 5}, {4, 7},
     {5, 11}, {6, 13}, {7, 17}, {8, 19},
     {9, 23}, {10, 27} }

This list is now stored in the variable “primes”. We can use this variable in a ListPlot command as follows:

   ListPlot[primes, PlotJoined -> True]

It’s good stylistic practice to give your own variables names that begin with lowercase letters. This prevents us from trying to use the name of a Mathematica command or function as a variable name. Mathematica commands are all case specific, which means that uppercase and lowercase letters are distinct.

 

If you just want to see what is stored in the variable primes, type

   primes

and Mathematica will respond with

   { {1, 2}, {2, 3}, {3, 5}, {4, 7},
     {5, 11}, {6, 13}, {7, 17}, {8, 19},
     {9, 23}, {10, 27} }

 

Plotting functions such as y = f(x)

Suppose we want Mathematica to plot the function y = x**2 over the range of x values 0 < x < 10. Enter

   Plot[x^2, {x, 0, 10}]

and you’ll get such a plot. (Note that Mathematica uses the ^ character to indicate raising a number to a power.)

 

The Plot command can be used to plot virtually any one-dimensional function:

   Plot[Sin[ x + 1/x ], {x, 0.5, 5}]

Let’s take this command apart to see how its components work together.

The command takes two arguments, which are (in order) the function we want to plot and the range of x values we want to use:

   Plot[ <function>, <range> ]

The range is expressed as a list with three elements: the first element is the variable that will be plotted on the horizontal axis; the second element is the lower limit on this variable; and the third element is the upper limit on this variable.

 

The fact that the independent variable is specified as part of the plotting range means that we can call the independent variable anything we like:

   Plot[Sin[time], {time, 0, 4 Pi}]

(Note that Mathematica has memorized the value of pi. You just have to remember that the letter P is capitalized.)

 

Combining two or more plots

Mathematica lets you store plots in variables so that you can combine several individual plots in a composite figure. Enter

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

and you will get two individual plots of the sine and cosine function. To plot both functions on the same set of axes, enter

   Show[splot, cplot]

You can combine different types of plots in this fashion. For example, you might want to combine a plot of experimental data points with a plot of a curve that fits through these points. The data points could be plotted with aListPlot command, and the curve with a Plot command. Then the Show command would combine these two plots.

 

Sometimes you will generate a plot, and after looking at it, decide that you want to save it in a variable so that you can combine it with another plot. To do this, you can use the % shortcut we already mentioned. For example,

   Plot[ Sin[x], {x, 0, 2 Pi} ]
   splot = %

will store the plotted sine curve in the variable splot.

 

The Show command will combine several plots at once:

   splot = Plot[ Sin[x], {x, 0, 2 Pi} ]
   cplot = Plot[ Cos[x], {x, 0, 2 Pi} ]
   tplot = Plot[ Tan[x], {x, 0, 2 Pi} ]
   Show[splot, cplot, tplot]

Or you can combine the plots in stages:

   splot = Plot[ Sin[x], {x, 0, 2 Pi} ]
   cplot = Plot[ Cos[x], {x, 0, 2 Pi} ]
   sandc = Show[splot, cplot]
   tplot = Plot[ Tan[x], {x, 0, 2 Pi} ]
   Show[sandc, tplot]

 

When you combine two or more plots, you may want to adjust the limits of the x and y axes to focus attention on a particular region of the plot. The PlotRange modifier to the Show command lets you do this. Enter

   Show[splot, cplot, tplot,
    PlotRange -> { {0, 10}, {-10, 10} }]

This tells Mathematica to extend the horizontal axis so that it includes the range 0 < x < 10 and the vertical axis so that it includes the range -10 < y < 10. The general format for the PlotRange modifier is

   PlotRange -> { { <Xmin>, <Xmax> }, { <Ymin>, <Ymax> } }

You can see that the range is specified as a nested list, but not as a list of two (x, y) points.

 

Note that when we changed the axis limits, Mathematica did not extend the plots across the entire x axis. If you want the curves to extend all the way to x = 10, you need to specify this in the original Plot commands:

   splot = Plot[ Sin[x], {x, 0, 10} ]
   cplot = Plot[ Cos[x], {x, 0, 10} ]
   tplot = Plot[ Tan[x], {x, 0, 10} ]
   Show[splot, cplot, tplot,
    PlotRange -> { {0, 10}, {-10, 10} }]

 

Fitting data points to polynomial expressions

Suppose we have a set of experimental data that we want to fit with a straight line. For specificity, imagine that we have measured absorbance as a function of concentration in a test of the Beer-Lambert law:

   concentration (M)     absorbance

        0.01                0.095
        0.02                0.202
        0.04                0.413
        0.06                0.629
        0.10                1.101
        0.15                1.735

First let’s put this data into a list of (x, y) points:

   data = { {0.01, 0.095}, {0.02, 0.202}, {0.04, 0.413},
    {0.06, 0.629}, {0.10, 1.101}, {0.15, 1.735} }

We can plot the data with

   ListPlot[data]

from which we see that the data looks fairly linear at low concentrations, but starts to curve upward as the concentration increases. (Well, it’s actually a bit hard to see this, but it will become more evident once we try to fit the data with a straight line.)

 

To fit the data to a straight line, type

   Fit[data, {1, conc}, {conc}]

This tells Mathematica to fit the data in terms of two functions of the variable conc; these functions are 1 and conc. You should get the result

   -0.044375 + 11.6875 conc

from which we see that the slope and y-intercept of the line are 11.6875 and -0.044375, respectively. (Note that Mathematica uses its default of six significant figures even though our data has fewer significant figures.)

 

The general form of the Fit command is:

   Fit[ <data-list>, <function-list>, <variable-list> ]

The first argument is a list of (x, y) data points, in the standard Mathematica form of a nested list. The second argument is a list of functions that Mathematica will use to try to fit the points. The third argument is a list of the independent variables in the data set. In this example, we have only one independent variable, so the variable list is a list with one element: the variable conc. (Of course, we could name the independent variable anything we want, as long as the same name appears in the list of functions and the list of variables.) Later we’ll see how we can fit data as a function of two or more independent variables.

 

You might be asking yourself, “Why is there a 1 in the list of functions?” Good question! Let’s try leaving the 1 out:

   Fit[data, {conc}, {conc}]
   11.2461 conc

Here we have instructed Mathematica to fit the data to a straight line that goes through the origin. The 1 in the variable list is used to fit the y-intercept in our original example. Think of the fit as finding the best constants A and B such that the data is described by the line

   A * 1 + B * conc

When we leave out the 1, we are telling Mathematica not to look for a constant term in the fitting equation. This is equivalent to finding the best constant B such that the data is described by the line

   B * conc

with A implicitly set to zero.

 

In the example presented here, we might want to force the y-intercept to be zero; after all, that is the sensible Beer-Lambert law prediction. But if our apparatus is dirty or malfunctioning, it may register a “baseline” absorbance even when there is nothing in our sample cell. This would show up as a constant error in the absorbance, independent of concentration. So the fit to the line

   A * 1 + B * conc

lets us estimate this error. Notice that the slopes of the two straight lines are slightly different.

 

As we mentioned in class, the Beer-Lambert law can break down at high concentrations. Let’s add a quadratic term to our fitting equation to see if the data is described better that way:

   linfit = Fit[data, {1, conc}, {conc}]
   -0.044375 + 11.6875 conc
   quadfit = Fit[data, {1, conc, conc^2}, {conc}]
   -0.00325599 + 9.91865 conc + 11.1374 conc^2

Now let’s plot both fits, along with the experimental data:

   points = ListPlot[data]
   curve1 = Plot[linfit, {conc, 0, 0.15}]
   curve2 = Plot[quadfit, {conc, 0, 0.15}]
   Show[points, curve1]
   Show[points, curve2]

We see that the quadratic curve does seem to fit the data better. In addition, the quadratic fit leads to a y-intercept which is much closer to zero, in accord with our intuition.

 

Generating and manipulating lists of numbers

The Table command can be used to generate a list of numbers using a predefined mathematical expression. Suppose, for example, that we want to make a list of the squares of the numbers from 1 to 10. This command will generate such a list:

   Table[ n^2, {n, 1, 10} ]
   { 1, 4, 9, 16, 25, 36, 49, 64, 81, 100 }

 

The general format of the Table command is

   Table[ <expression>, <iterator> ]

where <iterator> tells Mathematica what sequence of numbers to use in generating the list from <expression>. In this example, the iterator

   {n, 1, 10}

indicates that the variable n will go from 1 to 10.

 

The loop variable does not have to be an integer. A list of evenly spaced numbers in the interval between 0 and 1 can be generated by

   Table[ x, {x, 0, 1, 0.01} ]
   { 0, 0.01, 0.02, ..., 0.99, 1 }

Here the iterator

   {x, 0, 1, 0.01}

indicates that the lower limit on the variable x is 0, the upper limit is 1, and the interval between successive values of x is 0.01.

 

By using more complicated expressions in the Table command, we can generate lists of (x, y) points to be used in ListPlot. Try

   sine = Table[ {x, Sin[x]}, {x, 0, 2 Pi, Pi/100} ]
   ListPlot[sine]

In this example, the expression used to construct the list is itself a list. So the end result is a nested list of (x, y) points suitable for plotting with ListPlot.

 

So now we know how to build lists. What about taking lists apart? Define a simple list of five numbers:

   list = Table[2^n, {n, 0, 4}]
   { 1, 2, 4, 8, 16 }

To pick out the fourth item in this list, we type

   list[[4]]
   8

The double brackets are Mathematica’s way of extracting a single element from a list of items.

 

Now suppose we have a list of three (x, y) points:

   data = Table[ {n, 3^n}, {n, 0, 2} ]
   { {0, 1}, {1, 3}, {2, 9} }

This is a nested list of three two-element lists, each representing a single (x, y) data point. To get the third point, we type

   data[[3]]
   {2, 9}

Mathematica gives us the third element of the list, which is itself a list of two numbers. To get just the y coordinate of this point, type

   data[[3, 2]]
   9

Here we have two numbers in brackets, which we can think of as indices into the nested list. Basically we are asking Mathematica to take the third element of the original list, and then take the second element of this sub-list.

 

Complex numbers in Mathematica

Mathematica uses the capital letter I to represent the square root of -1. Type

   Sqrt[-1]

and you’ll get the answer

   I

You can use I in expressions: the complex number 2 + 3i is represented as

   2 + 3 I

in Mathematica. The generic complex number (x + y i) is written as

   x + y I

or, equivalently,

   x + I y

Mathematica has a tendency to “alphabetize” things, so it will usually print out (x + y i) in the second form.

 

Mathematica uses the function Conjugate to take the complex conjugate of a number. Try it:

   a = 2 + 3 I
   Conjugate[a]
   2 - 3 I

We know that the complex conjugate of (x + y i) is (x – y i). But Mathematica gives us

   Conjugate[x + y I]
   Conjugate[x + I y]

which, needless to say, is not very informative.

 

The problem is that we have not specified whether x and y are real or complex numbers, and Mathematica won’t make any assumptions about x and y without our help. If x and y are themselves complex numbers, then the conjugate of (x + y i) is not simply (x – i y). To tell Mathematica that x and y are real numbers, use the ComplexExpand command:

   ComplexExpand[ Conjugate[x + y I] ]
   x - I y

Now we get the expected result. Remember to use ComplexExpand when working with complex functions like wavefunctions.

 

Derivatives and integrals in Mathematica

Mathematica uses the command D[…] to take derivatives of variables and functions. As an example, type

   D[x^3, x]
   3 x^2

The general format of the D[…] command is

   D[ <function>, <variable> ]

which tells Mathematica to take the derivative of <function> with respect to <variable>. Mathematica treats all derivatives as partial derivatives, so we have

   D[x y^2, x]
   y^2
   D[x y^2, y]
   2 x y

 

To take the second derivative, we can just use the D[…] command twice in a row:

   D[ D[x^3, x], x ]
   6 x

 

The command to compute integrals is Integrate. It works much like D[…] in that we specify an integrand and a variable of integration:

   Integrate[x^2, x]
   x^3 / 3
   Integrate[x y, y]
   x y^2 / 2

Note that Mathematica omits the constant of integration that is so common in calculus textbooks.

 

We can also use Integrate to compute definite integrals by specifying a lower and upper limit of integration:

   Integrate[x^2, {x, 1, 2}]
   7 / 3

computes the integral of x^2 from x = 1 to x = 2. The special keyword Infinity can be used to specify a lower or upper limit of integration:

   Integrate[(1/x)^2, {x, 1/2, Infinity}]
   2

 

Some definite integrals are too hard for Mathematica to compute analytically. In this case, we can estimate the integral numerically using NIntegrate. (Note the double capital letter at the beginning of this command.) As an example, try

   Integrate[Sin[x^4], {x, 0, 1}]
   Integrate[Sin[x^4], {x, 0, 1}]

You see that when an integral is too hard for Mathematica, it simply “spits it out”. But we can compute the integral numerically:

   NIntegrate[Sin[x^4], {x, 0, 1}]
   0.18757

 

Another way to integrate numerically is to nest the Integrate command within the N[…] command:

   N[ Integrate[Sin[x^4], {x, 0, 1}] ]
   0.18757

This has the advantage that we can ask for more significant figures:

   N[ Integrate[Sin[x^4], {x, 0, 1}], 20]
   0.187569544685

Note that we don’t get 20 significant figures because numerical integration is inherently approximate. In this case we only got 18 significant figures. When integrating numerically, you usually need to ask for a few more significant figures than you really want.

 

Special functions in Mathematica

Here is a list of various functions that Mathematica has already defined for you. This page will grow as the semester continues.

Trigonometric functions:

  • Sin[x] gives the sine of x in radians.
  • Cos[x] gives the cosine of x in radians.
  • Tan[x] gives the tangent of x in radians.

 

Exponentiation and logarithms:

  • Exp[x] gives e to the power x.
  • Log[x] gives the natural logarithm of x.

 

Square roots:

  • Sqrt[x] gives the square root of x.

Contour and surface plots in Mathematica

Mathematica can make very nice contour and surface plots of three-dimensional functions such as z = z(x, y). The commands that generate these plots are ContourPlot and Plot3D for contour and surface plots, respectively.

These commands have the same basic form, although there are more “options” for Plot3D. To make a plot of some function of x and y in which x ranges from xmin to xmax and y ranges from ymin to ymax, enter

   ContourPlot[ <function>, {x, <xmin>, <xmax>}, {y, <ymin>, <ymax> ]

or

   Plot3D[ <function>, {x, <xmin>, <xmax>}, {y, <ymin>, <ymax> ]

For example, a contour plot of the function Exp[-(x^2 + 3 y^2)] near the origin can be generated as follows:

   ContourPlot[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2} ]

A three-dimensional surface plot of the same function is generated with the command

   Plot3D[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2} ]

 

If you try these commands, you might find that the surface plot seems “chopped off” near the origin. We can rectify this by including a larger range of z coordinates in the plot, using the PlotRange option that we first saw in two-dimensional plots:

   Plot3D[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2},
    PlotRange -> {0, 2} ]

This command tells Mathematica to include all parts of the function with z values between 0 and 2. We can also use the PlotRange option to specify boundaries on all three coordinates, for example to zoom into a portion of the surface:

   Plot3D[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2} ]
   Show[ %, PlotRange -> {{0, 1}, {0, 1}, {0, 2}} ]

 

You might notice that the resolution in this “zoomed” plot is fairly poor. Mathematica creates surface plots by scanning over a rectangular grid of points and calculating the height of the surface at each point. When we zoom in, Mathematica continues to use the original grid of points. We can make the grid finer by using the PlotPoints option in Plot3D. (Note that this option does not work in the Show command!) Here is how we might increase the resolution of this zoomed plot:

   Plot3D[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2}, PlotPoints-> 30 ]
   Show[ %, PlotRange -> {{0, 1}, {0, 1}, {0, 2}} ]

 

Of course, if you decided to zoom in after viewing a surface plot, you can combine the PlotPoints and PlotRange options in the Plot3D command to do everything in one step:

   Plot3D[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2}, PlotPoints-> 30,
    PlotRange -> {{0, 1}, {0, 1}, {0, 2}} ]

Even this is wasteful in a sense; we are generating a 30-by-30 grid of points in which x and y both range from -2 to 2, and then throwing away all of the points except those where x and y are between 0 and 1. Instead, we might try

   Plot3D[ Exp[-(x^2 + 3 y^2)], {x, 0, 1}, {y, 0, 1}, PlotPoints-> 30,
    PlotRange -> {0, 2} ]

which gives us higher resolution in the region of interest. 

The PlotPoints option applies to contour plots as well. Compare these two plots:

   ContourPlot[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2} ]

   ContourPlot[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2},
    PlotPoints -> 30 ]

The plot with more points is much smoother, and represents the oval shape of the surface more faithfully. But how can we get rid of the shades of gray so that we can see the contours better? Try

   ContourPlot[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2},
    PlotPoints -> 30, ContourShading -> False ]

Nice, isn’t it? The only problem is that we don’t know what z values the contour lines correspond to. With the Contours option we can specify exactly where we want the contours. Try entering

   ContourPlot[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2},
    PlotPoints -> 30, ContourShading -> False,
    Contours -> {0.2, 0.4, 0.6, 0.8},
    PlotRange -> {0, 2} ]

to see how you can control the placement of contour lines. Note that we need to include both the Contours and the PlotRange options to make sure that all of the contour values fit inside the range of z values that Mathematica will plot. Remember that when you only give one pair of numbers for PlotRange, Mathematica applies this range to the z axis. 

Finally, note that the variables x and y could be called anything we want. These two commands produce exactly the same results:

   ContourPlot[ Exp[-(x^2 + 3 y^2)], {x, -2, 2}, {y, -2, 2} ]

   ContourPlot[ Exp[-(u^2 + 3 v^2)], {u, -2, 2}, {v, -2, 2} ]

 


 


 

 


 

 

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s