Mathematica Programming:
Introduction for Numerical Analysis

Jason Cantarella

Mathematica as a calculator.

A good place to start with Mathematica is with entering expressions and evaluating them. Here are some examples:

22!

Sin[π/2]

E^0.2

N[π,22]

N[1.0/3.0,15]

1.0/3.0 // InputForm

Notice that:

You hit Shift-Return to evaluate the expression. Return just moves you down in the text window.

To type π, you enter Esc-pi-Esc (this works for other greek letters, too). Mathematica knows what π is automatically.

e is written with a capital ‘E’ (lowercase ‘e’) is just another variable name to Mathematica.

Numbers with a decimal point are treated as machine floating point numbers (actually, they are doubles). Numbers without are kept in fractional form, which is very often NOT what you want in this class.

You can convert anything to decimal form with “N”. The results aren’t what you might expect if you “N” a floating point number, but if you put it in “InputForm”, you can see all the digits.

Literate Mathematica

Your Mathematica work is a document called a “Notebook”. Your work will make a lot more sense (to both of us) if you use the different types of styling available in notebooks. The brackets on the right hand side of the notebook denote “cells”, each of which has a type. For instance, this is a “Text” cell. Here, return and shift-return do nothing but move to the next line. The cell above is a “Section” cell.

This is an (item) cell. You can select the cell type from the Format->Style menu, or with the keys Command-1 throught Command-9. You’ll mostly use the following:

Command-1 Title. This goes at the top of your notebook. Example title: “Homework #2, Math 4500”.

Command-2 Subtitle. This goes next. Example subtitle: “by Rachel Whitaker”.

Command-7 Text. This is where you explain what you’re doing. Put in LOTS of text cells. Mathematica doesn’t really have a good commenting mechanism for your code, so the preferred method is to make your functions really short, and intersperse a lot of text cells in the middle.

Command-4 Section. When you switch problems, or parts of a computation, a good idea is to use a “Section” heading to keep the reader oriented.

Command-9 Input. This is where regular Mathematica input goes.

You can change fonts (or the look of the entire document) with Format->Stylesheet. You should explore; you might like some of them.

The Palettes->Writing Assistant has an interface for “Typesetting” mathematical formulas, such as (∑_(n=1)^k)(n^2). I use it, but personally I hate this kind of thing.

Getting Help.

You will do a LOT of reading documentation and getting help with Mathematica. Luckily, the built-in help functions are excellent.

Help->Documentation Center  .

Help->Virtual Book.

Getting Help.

A help function that you might not guess is the question mark:

?Sin

Or the double question mark:

??Sin

Functions. The short form.

You can write down functions right away without worrying too much about what the syntax actually means. (We’ll get to that shortly.)

f[x_] := Tan[x] Exp[x];

f[7.0]

f[π]

There are several things to notice about the function syntax.

On the left hand side, the variable appears with an underscore, like x_. No underscore, no live function. Also, no error message, so you’ll make this mistake at least once.

The assignment is “:=” and not “=”. Again, no colon, no live function (but again, no error message).

Note the color coding on the right hand side:

Black. This is already defined.

Blue. This is not defined yet.

Green. This is a local variable for this function.

Functions take anything.

Absolutely anything.

f["My name is Ishmael"]

f[Graphics[Circle[]]]

f[Sound[SoundNote["C"],2]]

Plots. The short form.

Mathematica graphics and plotting is probably something that we could have a semester course on. However, most of what you’ll actually want to do is easy.

Plot[f[x],{x,0,1}]

Plots. Some options.

The most important generic options for Plot are probably AspectRatio->Automatic, Filling->Axis.

Plot[f[x],{x,0,1},AspectRatio->Automatic,Filling->Axis]

Plots. PlotRange.

You can set the range of the plot with PlotRange.

Plot[f[x],{x,0,1},AspectRatio->Automatic,Filling->Axis,PlotRange->{{0,1},{0,1}}]

Plots. Plotting data.

If you have a list of data to plot, you can do that with “ListPlot” or (usually better) “ListLinePlot”, which connects the data points.

Data = Table[{x,f[x]},{x,0,1,0.1}]

ListPlot[Data,PlotMarkers->{Automatic,Medium}]

Plots. ListLinePlot.

You can get these connected by line segments with “ListLinePlot”.

ListLinePlot[Data,PlotMarkers->{Automatic,Medium}]

Plotting Multiple Things Together.

To plot several functions on the same axes, use “Show”.

Show[Plot[f[x],{x,0,1},Filling->Axis],ListPlot[Data,PlotMarkers->{Automatic,Medium}]]

Multiple plots in Plot.

If you’re interested in plotting multiple functions, you can do this directly in Plot:

Plot[{Sin[x],f[x]},{x,0,1},Filling->Axis]

3d Graphics.

Somehow, you can’t resist making cool 3d plots of stuff, even when it doesn’t really help that much. But it’s cool.

Plot3D[Sin[x y],{x,0,π},{y,0,π}]

A somewhat more useful 3d plot is the ContourPlot.

ContourPlot[Sin[x y],{x,0,π},{y,0,π}]

Algebra and Calculus.

We will do a lot of calculus and algebra with Mathematica. So let’s see some examples.
First, you can take derivatives with D.

D[f[x],x]

D[f[x],x,x]

D[f[x],{x,5}]

This is a mess, and Mathematica may be able to help. You can always try to improve an expression with “Simplify”.

Simplify[D[f[x],{x,5}]]

Calculus continued.

You can do surprisingly complicated integrals with Integrate.

Integrate[f[x],x]

value = Integrate[f[x],{x,0,1}]

Since this answer doesn’t really mean anything to us, let’s evaluate it numerically.

N[value,15]

Taylor’s theorem and Series.

We’ll do a lot of Taylor series in this class, so let’s learn how to do that in Mathematica.

Series[f[x],{x,0,5}]

Notice the syntax: Series[function,{variable,point to expand around,highest order term}]. You can convert this to a polynomial with Normal

poly = Normal[Series[f[x],{x,0,5}]]

Now poly isn’t a function yet. Why not? Because x isn’t a variable. To convert poly into a function automatically, we’ll need some trickery. In the meantime, we can always put it into InputForm and copy and paste to use it to define a function.

poly // InputForm

Taylor’s theorem and series (continued).

p[x_] := x+x^2+(5*x^3)/6+x^4/2+(41*x^5)/120;

Plot[{f[x],p[x]},{x,0,1},Filling->Axis]

The other thing we probably want to do is think about the coefficients as list of numbers.

CoefficientList[poly,x]

Two kinds of equals. “=”

It’s not obvious from the surface, but Mathematica is really just the world’s most complicated search-and-replace function. That’s really what’s happening behind the scenes in essentially every Mathematica operation. So let’s start with the very basics. There are two kinds of “equals” signs in Mathematica.

bob = π + 7

bob

The “=” means “the left-hand side is forever after whatever the right hand side evaluates to right now”. This is a “dead assignment”, because the result is always the same.

fred = bob + 2

bob = 3

fred

Two kinds of equals. “:=”

On the other hand, the “:=” means “whenever you’re asked to evaluate the left hand side, replace it with the right hand side and evaluate that instead”. This is a “live assignment” because it changes as the stuff on the right-hand side gets updated. This is why it’s so important to use the “:=” when defining a function.

steve := bob + 2

steve

bob = 9

steve

Patterns.

Search and replace gets a lot more interesting if you’re allowed to use variables to match and replace stuff with a certain structure. The most basic pattern (you’ve already seen this) is the function definition:

newf[x_] := x^2;

newf[steve]

newf[π]

Here, the x_ means “a pattern (blank) matching any single object, which I will refer to as x”. Since a list is a “single object”, you can pass a list of numbers

newf[{1,2,3,4,5}]

But you can’t pass a pair of objects, since that doesn’t match the pattern.

newf[1,2]

Patterns, continued.

You can define a new function with the SAME NAME, but a different pattern.

newf[x_,y_] := x*y;

newf[7,8]

newf[{7,8}]

You can see from this that this is a) very convenient, and b) can get you into a lot of trouble. Mathematica will give you no warning when you redefine a variable name you’ve already used, reassign something you needed elsewhere, or otherwise mess yourself up. So I recommend:

If you’re just confused, restart Mathematica and re-evaluate your notebook from the top. If you’ve tried a lot of different things (and erased them), some previous definition is probably messing you up.

Don’t do too much in one notebook. It’s not unreasonable to answer different homework questions in different notebooks, for example.

Don’t do a lot of editing and erasing, or evaluating your notebook out of sequence. The results can be very weird.

Transformation Rules and Advanced Patterns.

Sometimes, Simplify won’t do what you want it to and you want to transform an expression in a more controlled way. You do this with a transformation rule, written x->y and “applied” with /.

poly

poly /. x->Cos[y]

poly /. x-> 5.0

More advanced transformation rules.

You can get much more subtle than this! Suppose I have

ratF =(Sqrt[2 x + 1] - 3)/(2x - 8)

Simplify[ratF]

I took calculus, so I know that I really want to multiply by the conjugate of the stuff on top, or make a more general transformation:

ratF /. (a_ + Sqrt[b_]) :> (a^2 - b)/(a-Sqrt[b])

Simplify[%]

Transformation Rules, Continued.

Sometimes, one application of a transformation rule doesn’t get the job done...

Sin[Sin[x]] /. Sin[x_] :> Cos[x]

In this case, you can use the //. operation to apply the transformation rules until the expression stops changing.

Sin[Sin[x]] //. Sin[x_] :> Cos[x]

If you give Mathematica a list of transformation rules, it will try them in order. You can do some pretty sophisticated things this way.

DerivativeRules = { MyD[x_ + y_] :> MyD[x] + MyD[y], MyD[x_ y_] :> MyD[x] y + x MyD[y], MyD[Sin[x_]] :> Cos[x] MyD[x], MyD[x_^n_] :> n x^(n-1), MyD[x] -> 1};

MyD[x^3 Sin[x^ 2 + 4]] /. DerivativeRules

MyD[x^3 Sin[x^ 2 + 4]] //. DerivativeRules

When Mathematica gives you transformation rules.

You can get transformation rules as output from Mathematica as well. For instance, when you use Solve to solve an equation, you get list of transformation rules which substitute in the various solutions.

sols = Solve[x^2 - 4 x + 7 == 0,x]

x /. sols[[1]]

x /. sols[[2]]

Heads.

One of the really kind of weird things about Mathematica is that it doesn’t really have “structure” data types. What we have is elementary data types and lists. Lots of lists. The data type of an object is called its “Head” in Mathematica.

Head[1]

Head[1.0]

Head[π]

When we pass stuff to a function in Mathematica, we are just wrapping it in a new Head to something new (and checking to see if there’s a transformation which applies to stuff with that Head).

Head[f[2,3,5]]

f[2,3,5][[1]]

As you’ve probably guessed, you can address parts of an expression with [[1]], [[2]], etc. Something to watch out for is [[0]], which stands for the head.

f[2,3,5][[0]]

You can change the head directly with the Apply function, which is abbreviated @@. As in

Times @@ f[2,3,5]

Lists.

The default expression type has the head “List”. This is what you’ll usually operate with (and lists of lists, and lists of lists of lists, and so forth).

Head[{1,2,3}]

Lists have “levels”, which correspond to the “level” of nesting in the list. So

nestedlist = {{1,2},{3,4},{5,6,7}}

has two levels. We can address these as

nestedlist[[1,2]]

Or get heads

nestedlist[[1,0]]

Most functions will “thread” over a list, meaning that they will replace a list with a new list contained the function evaluated on each member of the list.

Sin[{1,2,3}]

Table.

A standard way to generate lists is to use the Table function.

sinTable = Table[{x,Sin[x]},{x,0,2π,π/20.0}]

Tables, continued.

No sooner do you generate a table than you want to start messing with it.

sinTable // TableView

Operations on Lists.

Here the best I can do for you is to give you some examples and recommend to you the Virtual Book tutorial “Operations on Lists”. However, this is really important stuff, and you’ll probably want to do a lot of it. Given sinTable, a common question would be “Can I get separate lists of the x and y coordinates?” Sure! Use Transpose:

Transpose[sinTable] // TableView

yCoords = Transpose[sinTable][[2]]

Now I want to see the positive entries in this list. To do this, you use Select. This requires a function (to be called on each list member). You get the list entries for which this function returns “True”.

IsPositive[x_] := TrueQ[x > 0];

IsPositive[3]

IsPositive[-5]

Select[yCoords,IsPositive]

Operations on Lists II.

There are weirder things to do with lists, which you often want to do. We start with the “range” operator x ;; y which is kind of like Matlab. You can use this as a general submatrix operation on a nested list if you need to.

lessyCoords =yCoords[[2;;10]]

You can join lists with Join:

Join[{1,2,3},{"a","b"}]

It’s convenient to split things up into sublists sometimes, which you do with Partition.

Partition[lessyCoords,2]

And you can eliminate all that sometimes annoying nesting with Flatten:

sinTable[[1;;5]]

Flatten[sinTable[[1;;5]]]

A last move that comes in handy is the “riffle”, which combines two lists:

Riffle[{1,2,3,4,5},{"a","b"}]

The classic “riffle” is used for commas or spaces in text output, as in

StringJoin[Riffle[{"My","name","is","Ishmael"},{" "}]]

Iteration.

This is a really important topic, especially because the pseudocode used in the book and the lectures doesn’t always translate into the Very High Level syntax of a programming language like Mathematica. If you’ve taken programming, you’ve seen “the loop”:

for(i=0;i<10;i++) {

   x += Sin[i];
   
 }

Don’t ever do this in Mathematica. I’m so serious about this that even though it’s perfectly possible to code this in Mathematica, I am not going to tell you how.   

The point is that we need to make a distinction between two kinds of iteration. In the first, the same operation is being done on large number of different inputs and the results are listed. Like

for(i=0;i<10;i++) { 

    y[i] = sin[i];
    
 };
 
 This is called a “Map” in Mathematica, and the point is that this computation can be parallelized. It is very easy to crash your computer trying to demo this, because any simple computation is quick. But let’s try some harder ones.

Map.

We use Map when the results of the iteration don’t depend on each other, and we are applying the same function to a list of inputs, producing a list of outputs. I have to invent something nontrivial to do or we wouldn’t see any benefit from it, so let’s take an image and apply a lot of meaningless filters to it.

Picture = ;

StupidEnhance[v_] := ColorQuantize[ImageAdjust[Sharpen[Blur[Picture,v]]],25];

StupidEnhance[3.0]

Map[StupidEnhance,{1,2,3}]

Map! Map!

If we’re doing a lot of mapping, it can take a while.

(Map[StupidEnhance,1.0*Range[1,100]] // AbsoluteTiming)[[1]]

This can be easily parallelized by replacing “Map” with “ParallelMap”. Nothing else changes, but there’s some overhead, particularly if your computer has only a few cores, or is doing something else at the time.xn+1 = xn(2 − xnR)
is 1/R. This iteration was actually used to implement the division operation on some old computers which had ci

(ParallelMap[StupidEnhance,1.0*Range[1,100]] // AbsoluteTiming)[[1]]

Abbreviating Map.

Map can be abbreviated /@ and usually is:

Sin /@ Range[0,3,0.1]

You say this as “Sin mapped onto (stuff)”. It’s equivalent to

Map[Sin,Range[0,3,0.1]]

Nest.

The alternative to a “Map” is the kind of iteration where each iteration applies a transformation to some output. This is really common in numerical algorithms, which tend to iterate on a solution in the hopes of improving it. Here’s a problem from Homework #3.

This means that the iteration

myG[x_] := x(2 - x * 4);

should converge to 1/4 if we evaluate g[g[g[g[g[x]]]]] for any x close enough to 1/4. We do this with “Nest”:

Nest[myG,0.02,15]

We can see the intermediate evaluations with NestList:

NestList[myG,0.02,15]

This shows the iteration converging as hoped. Note that unlike “Map”, “Nest” can’t be parallelized.

Iterations that you shouldn’t code yourself.

Very often, you’ll want to do something that seems to require a loop, but is already a builtin function if you do a little digging. It’s definitely worth it to use the Mathematica function instead of writing your own code. Not only does it make for more readable code, it’s essentially always a lot quicker than your own code. This includes code that searches for or counts stuff in a list.

For example:

Range[1,5,0.1]

Mean[Sin /@ Range[1,5,0.1]]

Total[Sin /@ Range[1,5,0.1]]

Iterations that you shouldn’t code yourself. II.

Max[Sin /@ Range[1,5,0.1]]

Min[Sin /@ Range[1,5,0.1]]

Sin /@ Range[1,5,0.1]

Counting stuff in lists is one of the few places that I know that you want to use a “conditional pattern”. This is written x_ /; (something), where /; means “such that”. So to count the positive values of Sin in this list, we’d use

Count[Sin /@ Range[1,5,0.1],x_ /; x > 0]

Text Processing.

It’s worth noting, though maybe not spending a lot of time on right now, that Mathematica has an entire regular expression engine which is just as powerful as Perl or Python (or grep) if you want to use it. So for instance, if you want to count the number of times the letter “e” appears in the phrase “She sells seashells by the seashore” this can be done with

StringCount["She sells seashells by the seashore","e"]

You can construct more complicated patterns, like “e or h” (with the vertical bar):

StringCount["She sells seashells by the seashore","e"|"h"]

Or any letter character

StringCount["She sells seashells by the seashore",LetterCharacter]

Or any character at all

StringCount["She sells seashells by the seashore",_]

Similarly, you can do string transformations using StringReplace and ordinary transformation rules:

CheesyBadGuyAccent = {"s" -> "z","S" -> "Z"};

StringReplace["She sells seashells by the seashore",CheesyBadGuyAccent]

Speak[%]

You’ll do a lot more with this in homework, and you’ll see that one of the advantages of Mathematica’s string processing functions is that they are stupidly fast (in fact, you can comfortably deal with strings which are millions of characters long).

Interactivity. Manipulate.

Again, this is an immense topic. But the basics are well worth understanding. The idea is that you often don’t just want to display math or calculate something, you want to interact with the process by changing some parameter and seeing how the results change in real time. This is very easy to do in Mathematica.

Manipulate[Plot[Sin[k x],{x,0,2π},Filling->Axis],{k,0,20}]

Manipulate can display anything you like, including numbers, formulae, images, and expressions. So you can check out Pascal’s triangle with

Manipulate[Expand[(1+x)^n],{n,1,100,1}]

Manipulate basics.

And expand Sin[x] in a Taylor series with

Manipulate[Series[Sin[x],{x,0,n}],{n,0,35,1}]

Note the extra “by 1” step in the iterator above “{n,0,35,1}”. If we didn’t have it, the iterator would try to iterate over real values. Hilarity ensues.

Manipulate[Series[Sin[x],{x,0,n}],{n,0,35}]

Manipulate advanced topics.

Actually, you can do a LOT with Manipulate, including having input which isn’t sliders. This is just an example:

Manipulate[Plot[InterpolatingPolynomial[points,x],{x,-2,2},PlotRange->2,Filling->Axis],
{{points,RandomReal[{-2,2},{5,2}]},Locator}]

Manipulate advanced topics.

What makes this work is a special kind of control called a “Locator” which evaluates to its position on the screen.

Manipulate[Graphics[Polygon[pts],PlotRange->1],
{{pts,{{0,0},{1,0},{0,1}}},Locator}]

You can carry this as far as you want to generate little web applets in the “Computable Document Format” which run in “MathematicaPlayer”. In general, these ONLY get Manipulate inputs because otherwise they would BE Mathematica.

Mathematica in teaching.

Well, this entire class is sort of a demo for this. But let me point out for the graduate students, that the “Wolfram Demonstrations Project” is a curated collection of hundreds of Mathematica applets for demonstrating various things. They are awesome.