Mathematica workshop - worksheet 1

The basics

Please don't send people copies of this worksheet. Instead, direct them to
 http://j-star.org/mathematica_course.html 
 where the latest version will always be available.
© Jony Hudson, Imperial College London, 2006-2011
 [I don't like legal crap normally, but some dude copied the last Mathematica notebook I put on the web nearly verbatim into a textbook. Without asking me!]

This worksheet aims to introduce you to the Mathematica way of doing things through some simple exercises. By doing things in a way that fits naturally with Mathematica  you will hopefully find that subsequent, more realistic, excersises are more easily solved. But first, what do I mean by "the Mathematica way" ?

The Mathematica way

Mathematica isn't like most of the programming languages that many people know. And, at first, this can make it difficult to write structured, efficient Mathematica code. Hopefully this section will make some of the key differences clear.

There are several old families of programming languages that most modern languages are descended from. Almost all of the programming languages that are taught at school or university (to scientists at least), and most that are in everyday use are from the "Algol/Fortran" family. C, Pascal, Basic, Fortran, Java, C# are all from this family. Mathematica is not descended from Algol, but draws its core ideas from other language families (LISP, ML, rewrite systems). These ideas are very different to those behind the Algol family. To make matters more confusing, Wolfram (the makers of Mathematica) have tried very hard to make the system feel familiar to average users by presenting an Algol-like face. You might call Mathematica a rewrite system in LISP clothing, in Algol clothing. Or you might not. But either way keep in mind that when you're writing Mathematica code that the ideas behind the system might be different from what you're used to. Enough waffle - what are these ideas? Here are some of the key differences between Mathematica and the Algol family:

1) Uniform representation for everything. This is one of the most significant differences, but possibly the most abstract, so don't worry if it doesn't make too much sense now. Eventually it will. In the Algol family there are usually strong distinctions between source code, compiled object code, executable code, program data - these different types of data are usually stored in different formats and processed by different tools (the compiler, the linker, runtime system etc). Not so in Mathematica. Everything from the notebook files (the source code), to the functions in memory, to program data - graphics, sounds, tables, lists, vectors - are all represented by the same universal format and all processed by one tool, the Mathematica kernel. This has profound implications for the way code can be written. In particular, writing code that generates and modifies code, not data is very easy (data, code, who cares if it's all the same !?) and often very useful. 

2) Functions are first class objects. (In fact, because of the uniform representation, everything can be and is a first class object). In most languages in common use there is a distinction between data (variables) and functions. And usually functions cannot be manipulated as easily or as thoroughly as variables. Hence the first-class (variables, data, objects in OO languages) and second-class (functions) monikers. In Mathematica you can assign functions to variables, write functions that take functions as parameters (functionals), operators that act on functions - anything you like really. [You can sort of do most of this in most Algol languages - using pointers in C, interfaces/anonymous inner classes in Java, delegates in C# etc - but it's often inelegant and quite limited. Mathematica really treats functions as just another kind of data (of course it does, because there's a uniform representation for everything!)]

3) List based. Everything in Mathematica is represented as an "expression". And these expressions are fundamentally list based (this is what Mathematica gets from LISP - the LIS in LISP stands for list). As a result it almost always pays to structure your code around lists and use higher order functions on these lists to achieve your aims. I'll explain what this means below as it's much easier to show by example.

4) Pattern matching and rules. The core of Mathematica is a powerful pattern-matcher. This is used behind the scenes to do an awful lot of things and it's also available for you to use. Rules allow you to use the pattern-matcher to effect transformation of expressions. They make it very easy to program some applications that are traditionally quite complex, but they do take a bit of getting used to.

5) Interactive, interpreted. In Algol-like languages one typically writes the code in a text editor, compiles the code to an executable and then runs the code, possibly with a debugger attached. This cycle is repeated over and over until the program does what you want. Not so in Mathematica: expressions are evaluated as soon as you enter them. This changes the way that you work - the program is built while it's running, debugging is done by writing code to inspect the running system, code just kind of grows by itself. It takes some discipline to write in this style - the main trick is that if something seems even slightly unsatisfactory then re-write it straight away - but it does allow you to be very productive once you get used to it.

But enough talk! Time to do some coding. Don't worry if the above stuff doesn't make a lot of sense right now (you might re-read it in a few weeks when you're on the way to becoming an expert). The main things to get used to to be productive in Mathematica are its list-orientation and the importance of defining functions. We'll go through some exercises now that hopefully illustrate what this is about.

Exercises

Code your solutions to the exercises in the notebook just after the solution label (maybe now's a good time to save a copy). You need to evaluate any input cells in the problem definition.

Mapping over lists

One of the most useful functions in Mathematica is Map. It's an example of what computer-scientists call a "higher order list function". With a fancy name like that it must be pretty useful, right? (the name comes from the functional programming world - I don't think Wolfram use this term themselves (they tend to downplay the influence of LISP and ML on Mathematica - probably a Stephen Wolfram ego thing!)). It's actually very simple: all it does is apply a function to every element of a list.

Read the help page on Map (highlight Map and hit F1) and then try the following exercises.

Reversing sub-lists

(Evaluate the cell below before writing your solution)

l={{a,b},{c,d},{e,f},{g,h}};

Swap each pair of elements of the sub-lists of l (i.e. the result should be {{b,a},{d,c},{f,e},{h,g}} ).

Solution

Longer sub-lists

l={{a,b,c},{d,e,f},{g,h,k}};

Reverse the order of the elements in the sub-lists.
Hint: you should be able to use exactly the same line of code to solve this problem as the previous problem. After all you're doing the same thing really.
Hint 2: this section is about Map - does your answer have a Map in it ? You rarely need to write a loop in Mathematica.

Solution

Ragged sub-lists

l={{a,b},{c,d,e},{f,g,h,k}};

Reverse the sub-lists.
Hint: once again, the same line of code should work, as it's the same idea.
Note: compare how much easier this is than to doing it in C, for instance.

Solution

Defining functions, Module

Functions are one of the key organisational structures in Mathematica. You'd better get used to defining them! Pretty much as soon as you spot a pattern emerging in your code you should abstract that pattern, often into a function. Don't be lazy - always define a function if you feel you should.

The Mathematica command Module lets you clump together several lines of code. These lines of code can have their own local variables. It's very useful for defining self contained functions. Here's an example of defining a function using Module:

plotNoisySine[max_]:=Module[{data},
data=Table[{j,5 Sin[(j/10)]+Random[]},{j,max}];
ListPlot[data]
]

plotNoisySine[100]

plotNoisySine[1000]

Notice that the Module has lumped several commands together into one function (you don't need a Module to do this, but for now this is a safe way to do it). The Module has also stopped the data expression from leaking out - it's local to the Module. If I try and see its value:

data

It has none. Modules help keep your functions self-contained and avoid nasty side-effects (like when two functions have supposedly internal parameters that share the same name - the Modules stop these from interfering). It's worth reading about Module in the help pages.

Generating and plotting data

Generate a table of the values of Sin[x] for x from 0 to 50 in steps of 0.1. Plot it.
Hint: you'll need the functions Table, Sin, and ListPlot. Read about them in the help pages.
Hint 2: don't worry about the x-axis for now. ListPlot will assume x values {1,2,3,…} if you do not provide any. That's ok for these exercise. Perhaps you can come back later and fix the x-axis?

Solution

A different function

Repeat the above exercise but this time generate data for a Riemann-Siegel function(the Mathematica function is RiemannSiegelZ), square and plot.
Hint: you're just doing the same thing as above but a bit different, so see if you can re-use that code. And by re-use I don't mean copy and paste. You might have to modify your code above to make it work for this question too. To be clear, I want you to use the function you defined above to answer this part of the question - go change it if you need to. Do not define a new plotting function to solve this problem!
Hint 2: remember, functions can be passed as arguments to functions.
Hint 3: this section is about functions, you have defined a function, haven't you?

Solution

Transforming the data

Now plot the modulus of the Fourier transform of Sin[x]. Let the argument run from 0 to 1000 in steps of 1.
Hint: look up the functions Abs and Fourier in the help pages.
Hint 2: go back and generalize your solution to the first problem of this section if need be. Again, I don’t want you to define a new function to solve this problem - go and modify the first solution so that it’s general enough to solve all three parts.
Hint 3: maybe defining a function to do the Fourier transform and Abs will help.
Hint 4: you might find the function Identity useful when you revise your code above.
Hint 5: the answer might look a bit funny to you. Mathematica returns twice as much data as you might expect for Fourier - the second half is the same as the first, reflected.

Solution

Pure functions

OK, enough of named functions. Time for something new. Often you want to write a "throwaway" function, one that you're only going to use once. Mathematica provides the pure function for this purpose (usually called a λ-expression in other languages that have them, after Church's λ-calculus) [and confusingly, “pure function” means something else in almost every other functional programming language - again a case of Wolfram being insular]. The function that is used for defining pure functions is bogglingly called Function. An example:

Function[2((#)^2)]

The # symbol is used to represent the argument of the function (use #1, #2 etc if you want more that one argument). This defines a function (the one we've just defined is now pretty much useless as it doesn't have a name so we can't refer to it (except with % in a pinch)). The function is called the same way as any other function - by following it with square brackets containing the arguments:

Function[2 ((#)^2)][4]

This might seem very clunky, having to rewrite the function each time, but remember that these pure functions are meant to be throwaway. If you find yourself writing the same pure function over and over then perhaps it ought to have a name. Writing pure functions is another thing that's so common in Mathematica that it has a shorthand (in fact, you can see it above in Mathematica's output): simply write the function body followed by an &.

2((#)^2)&[4]

Here, 2((#)^2)& is the anonymous, pure function, and I'm calling it with argument 4, giving 2 (4^2)=32.

A very common idiom is mapping a pure function over a list. For example, if I want to take a list and calculate 3x+5, for each element x in the list:

l={1,2,99,1000};
(3#+5)&/@l

[Note: the parentheses around the function are not needed, but I think they make things clearer.]
[Advanced note: I don't actually need to do the Map here, as the simple algebraic operations that make up the functon are marked Listable, which means they automatically Map over lists. Ignore this for now, but keep it in the back of your mind as it will probably trip you up one day.]

A very common use of mapping a pure function over a list is manipulating datasets. Imagine I've got a list of the form {{(Subscript[x, 1]),(Subscript[y, 1])},{(Subscript[x, 2]),(Subscript[y, 2])},...} and want to normalise and offset the x-data and take the logarithm of the y-data. One way to do it is to write a pure function that transforms an {x,y} pair and Map it over the list:

dat={{0,1},{1,3},{2,9},{3,16},{4,25}};
{((#[[1]]-2)/2),Log[#[[2]]]}&/@dat

Just to be clear what's going on here: my pure function (the bit before the &) takes a single argument (there's only a #, no #1, #2 etc). This argument must be a list with (at least) two elements, because I'm taking the first and second parts of the argument. (Here #[[1]] represents the first element of the argument list (the x-value) and similarly for #[[2]]. (You get the double-bracket by typing esc [ [ esc - it's a shorthand for the function Part which, as always, you can find in the help pages.) ) The pure function returns a list of two elements (it's defined with curly brackets). I then Map the pure function over the data list.

This might be a bit confusing, so let's go through some exercises to get used to these ideas.

A simple pure function

Write a pure function that calculates 5 to the (x^th)power, where x is the function's argument. Evaluate the function with the argument 99.

Solution

Evaluating to a List

Write a pure function that evaluates to a two-element list. The first element should take the value 5 to the (x^th) power and the second x to the (5^th) power. Evalutate the function with the argument 5. [The answer should be {3125,3125}.]

Solution

Mapping over a list

l={1,8,75,300};

Write a pure function that calculates log(x)+sin(x) for argument x. Evaluate the function for each element of the list l.

Solution

Generating a List of Lists

l={1,2,4,8};

Write a pure function that generates a two element list {x,(x^2)} for argument x. Map it over the list l. [Your answer should be {{1,1},{2,4},{4,16},{8,64}} ].

Solution

Manipulating nested Lists

l={{1,5},{7,89},{43,60}};

Write a pure function that takes a two element list, {x,y}, and evaluates to a two element list that contains the sum and difference, {x+y,x-y}. Map this function over the list l. [You should get {{6,-4},{96,-82},{103,-17}} ].

Solution

Apply

Now that you've mastered pure functions we'll move onto the final concept for this worksheet, the function Apply. Apply can be rather confusing at first as it manipulates Mathematica's expression structures at a low-level. It's pretty useful though, so it's worth figuring out how it works. Let's start by looking a little deeper into the way Mathematica represents things. The command FullForm will let us gain an insight into the expression structure; it shows us how the kernel represents an expression:

a+b//FullForm

a b//FullForm

{a,b}//FullForm

(a^x)+(b^y)+5x//FullForm

You see that Mathematica expressions are of the form Head[ sequence of expressions or atoms ] , where Head is a symbol like List, Plus, Power and an atom is a thing like x, a, b, 5. This is the universal representation that I was talking about above. It's a simple but very powerful idea that everything can be represented as one of these nested-list expressions. It's very instructive to look at the FullForms (and TreeForms) of various Mathematica expressions. Here's an example to get you thinking:

{FullForm[#],TreeForm[#]}&[]

But we shouldn' t get too distracted by all of this elegance right now. Back to the Apply function. You notice that the expressions above for {a,b} and a+b only differ in their heads. The apply function lets us change the head of an expression - effectively mutating it into a different kind of thing. We can change a list into a sum like so :

Apply[Plus,{a,b}]

Again, there's an infix shorthand for this function, @@ :

Plus@@{a,b}

Here's a more complicated example, transforming a list of pairs into a list of exponentials by mapping an Apply across the list (I make a pure function to Apply the Power (!)):

l={{a,3},{b,6},{d,11}};
Power@@#&/@l

Mapping Apply comes up pretty often, so there's a shorthand for it too, @@@:

Power@@@l

We can generate a polynomial by swapping the List head of the generated List for a Plus:

Plus@@Power@@@l

You can begin to get a feeling for how this uniform expression representation makes easy to write code that generates code. But again, I'm getting distracted by the magic of it all. On to some exercises.

Simple Apply

s=a+b+c+d;

Change the above sum into a list. [Your answer should be {a,b,c,d} ].

Solution

Mapped Apply

p=(a+b)(c+d)(e+f);

Extract the variables from p into a nested list. [Your answer should be {{a,b},{c,d},{e,f}} ].
Hint: use FullForm to see what the expression looks like at input, and after your manipulations.

Solution

A worked example

Here's an example to pull together the concepts that we've just covered. It shows how you can do pretty sophisticated stuff very concisely and clearly in Mathematica once you get used to dealing with functions and lists.

Power series expansions of an exponential

I want to plot the truncated power series expansion of Sin[x] for different highest orders and compare it to the function itself. I'll show all the plots on one graph in several fetching colors.

Here's the code. We'll break it down below. Play with the slider to change the maximum order of the expansion.

Note: there are built-in Mathematica functions to do a number of the things I do by hand in this example. In general you should use a function if it exists, but for the sake of demonstration I'll not use them here.

Manipulate[
Show[MapThread[
Plot[Tooltip[#1,#1],{x,0,4π},PlotStyle->#2,PlotRange->{All,{-3,3}}]&,{
Join[{Sin[x]},Plus@@Table[((((-1))^(n+1))/((2 n-1)!))(x^(2n-1)),{n,1,#}]&/@Range[highestOrder]],
Hue[1-((#-1)/( highestOrder+1))]&/@Range[highestOrder+1]
}]],
{highestOrder,0,16,1}]

Step-by-step

This is how I built the function up, step-by-step. This follows the order in which I wrote it. 

First I want to generate a power series expansion. The first ten terms are:

Table[((((-1))^(n+1))/((2 n-1)!))(x^(2n-1)),{n,1,10}]

I can convert this to a sum by swapping the List head for a Plus (there's a function Sum which does all this in one step, but this way's more fun! (or your could use Series and Normal)).

Plus@@Table[((((-1))^(n+1))/((2 n-1)!))(x^(2n-1)),{n,1,10}]

Now let's make this a pure function of the highest exponent.

Plus@@Table[((((-1))^(n+1))/((2 n-1)!))(x^(2n-1)),{n,1,#}]&[10]

Generate a list of expansions from highest order 1 to 10 by mapping this pure function over a Range.

Plus@@Table[((((-1))^(n+1))/((2 n-1)!))(x^(2n-1)),{n,1,#}]&/@Range[10]

Tack the real function onto the start of this List with a Join.

Join[{Sin[x]},Plus@@Table[((((-1))^(n+1))/((2 n-1)!))(x^(2n-1)),{n,1,#}]&/@Range[10]]

Write a pure function that takes two arguments - a function to plot and a color. I'll wrap the argument to Plot in a Tooltip which will annotate the resulting plot if you hover the mouse over it.

Plot[Tooltip[#1,#1],{x,0,10},PlotStyle->#2]&

We now want to Map this plotting function over the list of functions and a list of colors. The function that does this is a generalization of Map, MapThread.

highestOrder=10;
MapThread[
Plot[Tooltip[#1,#1],{x,0,10},PlotStyle->#2,PlotRange->{All,{0,1000}}]&,
{
Join[{Sin[x]},Plus@@Table[((((-1))^(n+1))/((2 n-1)!))(x^(2n-1)),{n,1,#}]&/@Range[10]],Hue[1-((#-1)/( highestOrder+1))]&/@Range[highestOrder+1]
}];

We've generated a list of Graphics expressions, now we just need to display them all together. Show is the function to do this.

Show[
MapThread[
Plot[Tooltip[#1,#1],{x,0,10},PlotStyle->#2,PlotRange->{All,{-3,3}}]&,{
Join[{Sin[x]},Plus@@Table[((((-1))^(n+1))/((2 n-1)!))(x^(2n-1)),{n,1,#}]&/@Range[10]],Hue[1-((#-1)/( highestOrder+1))]&/@Range[10+1]
}
]
]

All that remains is to pull out the magic number 10 and wrap it in a Manipulate.

Playtime

Here are some more free-form exercises to practice the above techniques. When you solve them try and think about doing things in the Mathematica style. Don't worry if you can't finish them all - they aren't required for understanding the subsequent material.

Expression building

Write a function which takes two arguments, a list of numbers {(Subscript[c, 0]),(Subscript[c, 1]),(Subscript[c, 2]),…} and a single number x and evaluates (Subscript[c, 0])+(Subscript[c, 1])x+(Subscript[c, 2])(x^2)+… at the value of x provided.

Solution

Radical!

Write a function that calculates the value of (Sqrt[6+(Sqrt[6+(Sqrt[6+…])])])for an arbitrary number of 6s. Produce a list of the results from one 6 to ten 6s.
Hint: you might want to familiarise yourself with one of the other higher-order functions, Nest.
Hint 2: Mathematica likes to keep track of everything exactly if it can. You might need to use the function N to turn your answers into numerical results.

Solution

π

Calculate π using Leibniz's formula, π=4(∑_(j=0)^∞)((((-1))^j)/(2 j+1)). Write a function that takes an integer n and calculates the sum of the first n terms in this series. Plot the approximate value versus n for n up to 20. Show the exact value of π on the same plot.
Hint: Mathematica has a built in Sum function which does just what it says.
Note: Mathematica can do the infinite sum analytically - try it!

Solution

Monte Carlo π

You can approximate π by generating points uniformly distributed within the unit square, and then counting what fraction of them lie within the unit circle.Write a function to do this for n points.
Advanced: can you plot the points and the circle? Can you color the points which fall in the circle red and the points outside blue?

Solution