Lecture 8 : Functions and Loops in Mathematica

Today we' ll cover how to create and use functions in Mathematica. We'll cover when you should use Set vs. SetDelayed, recursion in Mathematica, and how to memoize functions.  We'll also cover functions which only apply to certain arguments, anonymous functions and the associated functional programming commands: Map and Apply, as well as other instances where it is natural and useful to use anonymous functions: Sort and Ordering, for example. If we have the time, we'll cover the looping structures: For, Do, and While loops. That's it for Mathematica programming. From now on, we'll be learning how to put together the stuff we learned in earlier lectures to do interesting things-- lectures will consist mostly of examples.

Please ask questions if anything is unclear.

Notices: 
- On prob 2 of the homework due on Wednesday, θsamples should be 4000 evenly spaced points between (4000^(-1)) and 2π, not random numbers in that range; also, on part a, ii) the idea is to use Sequence and Apply to avoid typing in multiple iterators, the way it is done in part a, i). In prob 3, be careful when defining lhs to use SetDelayed
- The problem set assigned next week (due Wed Nov. 26) will be the last; that brings you to a total of 4 Mathematica problem sets. You'll have an extra week to work on your Mathematica project.

Functions

Functions are useful in Mathematica even when you' re simply using it for calculations, because they allow you to avoid retyping commonly used expressions. Functions are simply expressions whose values are to be determined at a later date. To let Mathematica know what the parameters to the function are, we use patterns, and to tell it to delay evaluation of the expression, we use delayed assignment := (SetDelayed) instead of = (Set)

(* Simple function of one variable *)
sinc[x_] := Sin[x]/x;

sinc[0]
sinc[0] =1; (* note that a more specific set will override the setdelayed *)
(* sinc[x_] := If[x ==0, 1, ((Sin[x])/x)] (* alternative definition *) *)
sinc[0]
?sinc

You may wonder why we use SetDelayed instead of Set to define functions. Indeed, in some cases it makes no difference, but consider these examples:

Clear[x,y,a]

f1[x_] = FullSimplify[x^2];
g1[x_] = Expand[x^2];
h1[n_] = Table[a[i], {i,1,n}];

f2[x_] := FullSimplify[x^2];
g2[x_] := Expand[x^2];
h2[n_] := Table[a[i], {i,1,n}];

f1[Cos[x]+Sin[x]]
g1[x+y]
h1[3] (* still works, despite the earlier warning *)

f2[Cos[x] + Sin[x]]
g2[x+y]
h2[3]

The problem here is that when f1 and g1 are defined, Mathematica applies FullSimplify and Expand immediately to x^2, so these definitions are equivalent to f1[x] = g1[x] = (x^2). Using SetDelayed, Mathematica waits until x has been replaced by an expression before evaluating the calls to FullSimplify and Expand. Likewise, when you create a Table, Mathematica expects to know how many entries it should contain-- it won't accept symbolic expressions-- so to avoid errors, you must use SetDelayed.

If a function is defined using a compound expression (grouped with parentheses, or using Module, Block, or With), its value is that of the last expression. If the last expression has no value, the value is Null.

examplefun[x_] := ( Print["The argument was " <> ToString@x]; x^2; x^6); (* <> is the infix notation for StringJoin *)
examplefun[x] + 2

Because functions are defined with patterns, we can write functions that give you different types of output depending on the type of inputs. Only one pattern will match a particular expression, so only that function definition will apply.

RandomMat[n_] := RandomInteger[{-5,5}, {n,n}] (* random square matrix *)
RandomMat[n_, m_] := RandomInteger[{-5,5}, {n,m}] (* random rectangular matrix *)

RandomMat[4] // MatrixForm
RandomMat[4,5] // MatrixForm

(* function of 4 variables *)
distance[x1_, y1_, x2_, y2_] := (Sqrt[(((x1-x2))^2)+(((y1-y2))^2)]);

(* more elegant *)
distance[p1_, p2_] := (Sqrt[(((p1[[1]] - p2[[1]]))^2)+(((p1[[2]] - p2[[2]]))^2)]);

(* another alternative *)
distance[{x1_, y1_}, {x2_, y2_}] := (Sqrt[(((x1-x2))^2)+(((y1-y2))^2)]);

(* dimension independent, variable metric *)
distance[p1_, p2_, D_] := Norm[(p1-p2).D.(p1-p2)];

distance[{1,2,3},{4,5,6}, DiagonalMatrix[{1,2,3}]]
distance[{1,2,3},{4,5,6}]
distance[{1,2},{3,4}]
distance[1,2,3,4]

As an important aside, note that DSolve and RSolve return their solutions as lists of rules, which you can then turn into functions using what we' ve learned. Let' s revisit an example from the first lecture

Quiet[rules = DSolve[ {r''[x] + r'[x]r[x]==6, r[0]==1, r'[0]==1}, r[x], x]] (* returns a list of rules for r[x]*)

soln[z_] := (r[x] /. rules[[1,1]]) /. x->z; (* one way to turn the result into a function; we'll see an easier method when we cover anonymous functions *)
Plot[soln[z], {z,0,1}]

When determining the value of a function call, Mathematica goes with the most specific pattern match

WeirdConvert[a_ + b_] := WeirdConvert[a] + WeirdConvert[b];
WeirdConvert[a_ b_] := WeirdConvert[a] + WeirdConvert[b];
WeirdConvert[ Power[a_,p_]] := WeirdConvert[a]+p;
WeirdConvert[a_] := a;
WeirdConvert[Sqrt[a_]] := WeirdConvert[(a^p)];

WeirdConvert[(Sqrt[a b +c (d/(e+f))])]

By default, all symbols in Mathematica are global. To use symbols local to a particular function, use the Block or Module commands. Module is safest, because it creates completely new local symbols, whereas Block uses the same symbol, but gives it a local value. When in doubt, use Module. Block and Module are compound expressions, so the value of the function is that of the last statement executed (if the last statement ends in a semicolon, the value is Null)

(* example of why you might need to use local variables *)
f[g_, x_] := Integrate[g[x+s], {s,0,1}];
f[Cos, y] (* this gives the expected value *)
f[Cos,s] (* this doesn't work *)

(* possible remedies *)
f[g_, x_] := Module[{s}, Integrate[g[x+s], {s,0,1}]]; (* internally Mathematica replaces s with a new symbol of the form s$nn so there won't be any conflicts *)
f[Cos, s] (* this works *)

f[g_,x_] := Block[{s}, Integrate[g[x+s], {s,0,1}]]; (* Mathematica temporarily replaces the value of s *)
f[Cos, s] (* this doesn't work-- hence the caution to use Block when you're not clear which to use *)

(* further examples *)
Clear[x,y]
y=x
Block[{x}, x=a; x+y] (* Mathematica replaces the value of x locally *)
x
Module[{x}, x = a; x +y] (* Mathematica replaces x with x$nn *)
x

Let' s write a function which solves systems of equations

(* without using Thread (uglier) *)
solveSystem[sysmat_, rhs_] := Module[{lhs = Table[x[i], {i, 1, Length@Transpose@A}]},
lhs /. First@Solve[ Table[ sysmat[[row]] . rhs ==lhs[[row]], {row, 1, Length[A]}], lhs] (* assume the system is well-posed *)
]

(* using Thread *)
solveSystem[sysmat_, rhs_] := Module[{lhs = Table[x[i], {i, 1, Length@Transpose@A}]},
lhs /. First@Solve[ Thread[ sysmat . rhs ==lhs], lhs] (* assume the system is well-posed *)
]

A = RandomReal[{-1,1}, {10,10}];
v = RandomReal[{-1,1}, 10];

solveSystem[A,v]

Recursion and Memoization

The idea behind recursion in Mathematica is the same as that in any other programming language-- you write a function that calls itself. As usual, you have to be careful that you don't inadvertently cause an endless loop; that is, you have to provide an exit condition. The special variable $RecursionLimit specifies how deep Mathematica will allow recursion to go before giving an error message. Use Block to temporarily change the value of $RecursionLimit (note that you cannot use Module, because it internally changes the name of its local variables)

(* Fibonnaci sequence *)
fib[0] = 1;
fib[k_] := fib[k-1]+k;

fib[n] (* since this expression is symbolic, Mathematica will continue recursion until it has reached $RecursionLimit-- it's usually not a good idea to call a recursive function with symbolic arguments, unless you've specially designed it to handle them *)

(* Optional: a safer way to define fib-- only match expressions that are Integers *)
Clear[fib]
fib[0] = 1;
fib[n_Integer] := fib[n-1] + n;
fib[n]

fib[300]
Block[{$RecursionLimit = 4000}, fib[300]] (* temporarily increase $RecursionLimit *)

The idea behind memoization is to store the values of a function as they are computed, so they don' t have to be recomputed. For example, if I memoized my fib definition, I wouldn't have the problem of nested recursion. This would increase the speed with which any fib value is computed. To memoize in Mathematica, define your function as, e.g.
f[x_] := f[x] = ...
The first time f is called with some argument, the expression f[x] = ... will be evaluated, assigning the value of ... to f[x], and returning it as the value of the function call. The next time you call f with that particular expression, Mathematica will use the value for f[x] you just assigned. For example:

Clear[fib];
fib[0] =1;
fib[n_] := fib[n] = fib[n-1] +n;

Block[{$RecursionLimit = 4000}, fib[300]]
??fib

As a practical example of memoization, consider generating a sample path of Brownian motion using a random Fourier series

X[k_]:=RandomReal[NormalDistribution[0,1]]; (* random coefficients *)
W[t_]:=(t/(Sqrt[π]))X[0]+(Sqrt[(2/π)])(∑_(k=1)^n)((X[k])/k)Sin[k t];(*Fourier-Wiener sine series *)
n =300;

Plot[W[t],{t,0,π}] //Timing (* the problem is that each time Plot evalues (Subscript[W, t]) it uses new coefficients *)

 X[k_]=RandomReal[NormalDistribution[0,1]]; (* attempt at a cheap fix; this doesn't work *) 
Plot[W[t], {t,0,π}]

X[k_] := X[k] = RandomReal[NormalDistribution[0,1]]; (* the correct way to do this sensibly *)
Plot[W[t],{t,0,π}] // Timing (* also much quicker *)
(* if you want a new sample path of W, you need to clear and redefine X *)

Anonymous Functions and Operators to use them with

You use anonymous functions in Mathematica quite often. Unlike Matlab, where you might find yourself assigning them to variables to use them, you usually define and use them directly.
There are two ways to define an anonymous function: using the # & syntax, or the Function syntax.
The Function syntax is: Function[arg, body] or Function[arglist, body]. 
The # & syntax is: (body which uses #1, ... #n as arguments) &
Function[body] is another syntax which is equivalent to the (body) & syntax
Note that # is equivalent to #1 inside of the body
Examples will clarify:

f = Function[#^2] ;
f[1]
f[{1,x}]
??f

g = Function[Norm[#1, #2]];
g[a,b]
??g

You can use anonymous functions directly as the head of expressions

(# &) @@ {1,2} (* remember Apply replaces the head of expressions; in this case, it replaces the List with an anonymous function *)
(#1 &) @@ {1,2}
(#2 &) @@ {1,2}
Function[Log[#]][3,4]
Function[Log[#2]][3,4]

But by far, you most often use anonymous functions in a transient manner, when you need to use a function only temporarily, so it doesn't make sense to store their definitions. In this connection, it is useful to know that the Map function, which has the syntax Map[f, list] or f /@ list, returns a list constructed by applying f to every entry in list. For example, let's revisit the Monte Carlo calculation for pi:

(* incrementally *)
insideCircleQ = Function[pt, If[pt[[1]]^2 + pt[[2]]^2 ≤ 1, 1, 0]];
randomPoints = Table[{RandomReal[], RandomReal[]}, {10^5}];
4 Mean[insideCircleQ /@ randomPoints] //N

(* quick and dirty *)
4Mean[Map[ If[#[[1]]^2 + #[[2]]^2 ≤ 1, 1, 0] &, Table[{RandomReal[], RandomReal[]}, {10^5}]]] //N

As another example of the application of Map and anonymous functions, let' s calculate the (∞, 1) norm of a matrix A. This is a combinatorial quantity, so the list manipulation functions in Mathematica make calculating it a breeze :

InftyBallSamples[A_] := Function[point, Norm[A . point, 1]] /@ Tuples[{1,-1}, Length[Transpose[A]]]; (* get the 1-norm of A at all the points on the ball *)
FindIndexOfLargestSample[A_] := First@Ordering[InftyBallSamples[A],-1]; (* find the index of the (1,-1) tuple where the maximum 1-norm of A occurs *)

InftyOneNorm[A_] := Norm[A . Part[Tuples[{1, -1},Length[Transpose[A]]], FindIndexOfLargestSample[A]], 1]; (* put it all together to get the infty-1 norm of A; how would you do this in Matlab? *)

A = (({{1, 3, 4}, {0, 5, 6}, {0, 7, 8}}));
InftyOneNorm[A]

Other functions that may be useful include Sort and Ordering. You saw Ordering used above (read the Help for details). Sort takes a list and returns it sorted into a canonical order

Sort[{4,3,1,5}] (* for numbers, that is increasing *)
Sort[{z,y,h,a}] (* for symbols, it is lexographically increasing *)

For our purposes, it' s useful to note that you can provide an additional argument to Sort if you want to control how the list is ordered. For example, if I wanted a list of numbers to be sorted in decreasing order, I could call Sort as Sort[list, Greater]. Sort[list, pred] orders the list by evaluating pred[x, y] for every pair of x and y in the list-- this implies that pred must accept two arguments. If pred is true, then x comes before y in the list.

Sort[RandomInteger[{-4,5}, 6], Greater]

As an example, here' s a much simplified definition of InftyOneNorm[A]

InftyOneNorm[A_] :=Norm[ A .(First@Sort[Tuples[{-1,1}, Length@Transpose@A], Norm[A . #1, 1] > Norm[A.#2, 1] &]), 1];
InftyOneNorm[A]

When you use DSolve, RSolve, etc. you can get it to return an anonymous function by asking it to solve for r instead of r[x] (see example below for clarification)

Quiet[s = First[r /. DSolve[ {r''[x] + r'[x]r[x]==6, r[0]==1, r'[0]==1}, r, x]]] (* returns an anonymous function*)

 Plot[s[x], {x,0,1}] (* simpler than not using anonymous functions *)

Looping Structures in Mathematica

The most commonly used control structures in Mathematica are the For, Do, and While loops. Their most useful syntaxes are
For[startexpr, testexpr, increxpr, body]
Do[expr, {num}]
Do[expr,{i, num}]
Do[expr,{i, min, max}]
Do[expr,{i, min,max,stepsize}]
While[testexpr, body]
The Break command exits the inner most loop while the Continue command exits the current iteration of the inner loop.
Each expr may contain multiple expressions, using a Block, Module, or simply by separating multiple statements by semicolons

Note that the value of a control expression is Null, so if for example you want to build up a list using a loop, you need to assign the values to a variable to store them. The usage is straightforward-- pretty much the same way you'd use For, Do, or While loops in other languages.

3 For[i=2, i<5, i++, i] (* this is not 15, as you might expect *)

(* this doesn't work-- you can't access parts of lists unless they exist *)
For[ i=1, i < 40, i++, a[[i]] = If[Mod[i,3]==0, 1, 0]]

(* instead use the idiom of storing the values in a[i] *)
For[ i=1, i <40, i++, a[i] = If[Mod[i,3] ==0, 1, 0]]
a = Table[a[i], {i,1,39}] (* if you need them in list form *)

Another example of a For loop : find the first factorial that is also a Fibonacci number

FibQ[x_] := Or[IntegerQ[Sqrt@(5x^2+4)], IntegerQ[Sqrt@(5z^2-4)]];
TimeConstrained[ For[i = 2, !FibQ[i!], i++, i], 30] (* the principle is sound, but this might run forever, so constrain it to run for at most half a minute *)

An otherwise pointless example that illustrates two important things : the use of compound expression in loops, and how Module implements local variables

For[a = 1; Block[{a}, a = 2], a < 5, a++, Module[{a}, Print[a]]; Print[a]]

If you use For loops in functions, you can return a value from inside the loop by using the Return command

f[n_] := For[i=n, i>0, i--, If[IntegerQ@Log[2, i], Return[i]]] (* find the closest power of two less than n *)
ListPlot[f /@ Range[500]]