. Basics of Mathematica Programming

Paul Abbott

Pre-requisites

To improve the appearance of graphics, TraditionalForm will be used as the default format type for text that appears in graphics,

$FormatType = TraditionalForm;

and the default style for text in graphics is set as follows.

$TextStyle = {FontSize -> 10};

In this assignment simple functions are used to illustrate the elements of Mathematica programming style. Aspects of the language including function definitions, transformation and replacement rules, pattern-matching, functional operations, and association are also introduced.

Built-in functions are fast and implemented in the most general way possible. They mesh well with other objects and their properties —  derivatives, integrals, series, special cases, numerical evaluation (for  arbitrary complex parameters), etc.— are incorporated into the kernel.

There are over 1000 built-in functions and often you do not have to write a  program (in the standard sense) but can solve problems by a sequence of  simple built-in operations. Not only are such programs elegant and compact but they are also efficient and modular and can serve as building blocks for  robust large programs.

Mathematica includes a wide array of programming constructs and programming in Mathematica is different to most other languages. Before considering the style of programming most appropriate for a given problem, you should check whether the functions or operations you require already exist in the system. In general, commands are interpreted rather than compiled which is usually convenient.  Moreover, one can do a lot of  computation without writing a collection of functions longer than a couple of lines each. By changing the parameters and selecting a range of cells, one can get Mathematica to update all the cells — similar to working with a spreadsheet — and this is, for many purposes, all the programming needed to get the answer to a wide range of problems.

In these exercises I have included a number of examples which are  applicable to Mathematica programs of all sizes.  Time spent writing few-line programs is rarely wasted because it is more likely to lead to code that is elegant, efficient, robust, maintainable, understandable, verifiable, and extensible.

The range and flavour of Mathematica programming styles are illustrated using simple examples: First, a factorial exposition outlines the different styles available. Two sections on the fundamental tools — functional programming, pattern-matching  and replacement rules — are demonstrated using a number of exercises. The section on association highlights the advantages and pitfalls of association. The last section presents a brief discussion on writing programs by collecting together code fragments.

A number of basic principles are identified which, it is hoped, will help the reader take advantage of the possibilities offered by this high-level, powerful, and flexible, programming language.

. Factorial!

In the following examples, the factorial function is used to present an overview of the range of Mathematica programming styles available. The factorial function is an excellent vehicle for expounding the various Mathematica programming styles because it admits such a large number of points of view:

To some extent each of these viewpoints is naturally implemented using a particular style of programming. An important question to help select a programming style is whether the code will work for both numerical and symbolic parameters (if appropriate) and the ease of generalisation to more parameters, e.g., to two and more dimensions.  In many situations, a specific choice of programming style may present itself as optimal because it most naturally handles the widest range of possibilities; functional programs often satisfy this criteria.

.. Built-in

The factorial function is, of course, built-in:

10!

Alternatively, one can use the Gamma function (Γ(n + 1) ≡ n!) which  generalises factorial to non-integer, negative, and complex values.

Γ(11)

.. Procedural

Fortran, Pascal, or C programmers will probably understand the following procedural code without any explanation:

\(TraditionalForm\`\(facProcedural = 
Function[{n}, 
Module[{i = 1, s = 1}, 
	
            While[i < n, \(i++\); 
              s = s\ i]; 
s]
];\)\)

facProcedural[10] == 10!

The syntactical subleties are the use of Function, denoting a pure or anonymous function, and Module, which delimits a block of code.

Advantages of procedural style include:

Disadvantages include:

.. Recursive

The recursive view of the factorial function leads to the following recursive program which may appear more natural:

\(TraditionalForm
    \`facRecursive[n_] := If[n == 1, 1, n\ facRecursive[n - 1]]\)

facRecursive[10] == 10!

The syntactical subleties here are the use of n_  (read as n blank) to denote a generic variable n whose value is bound and passed to the right-hand side of the expression  and := (delayed assignment or SetDelayed) which delays execution of the right-hand side of the procedure until n is bound to a value.

.. Rule-based

For many situations rule-based programming is often simplest and most elegant. Specifying the initial condition:

facRule[0] = 1;

and the general recursive rule:

facRule[n_] := n facRule[n - 1]

Mathematica automatically applies the recursive definition (up to $RecursionLimit):

facRule[10] == 10!

One advantage of this programming style is that it is often similar to the actual problem specification for a wide range of problems in science and engineering. This fact reduces the likelihood of transcription and algorithmic coding errors.

.. Functional

For the purpose of this overview functional programming is defined as programming using a level of abstraction which emphasises  the essential algorithmic structure.

A functional algorithm for factorial can be stated as multiply together all numbers from 1 to n, and coded using Range, say

Range[10]

and then Applying (using the shorthand infix notation @@) the Times operation to this List to multiply these numbers together:

Times@@%

This is easily converted into a procedure:

facFunctional[n_] := Times@@Range[n]

facFunctional[10] == 10!

FoldList

Another important functional operation is FoldList:

?FoldList

Using FoldList one can efficiently generate, in one pass, a whole table of factorials values (starting from 0!):

FoldList[Times, 1, Range[10]]

.. Series

Truncating (using Normal) the asymptotic series for Γ(n),

Normal[
      Series[Γ(n), {n, ∞, 1}]]

yields, using the relation n!  ≡  
        n Γ(n),

facSeries[n_] = n %

This expression (the first term in Stirling's formula) is suitable for approximate computation of the factorial function. For example, here is the fractional error in this approximation for n from 1 to 10:

Table[
      N[((n! - facSeries[n])/(n!))], {n, 10}]

This idea is quite general: if one has an analytic representation of a function (or a recurrence relation), an asymptotic expansion can be used to generate approximations.

.. Dynamic Programming

As it stands, the code for facRule above is somewhat inefficient. After entering

facRule[10]

one can use ? (i.e., Information) to see the values assigned to facRule:

?facRule

One finds that, after computing facRule[10], intermediate values have not been saved.  This can be addressed by using  dynamic programming, i.e., saving the intermediate results. The syntax for this is always of the  form f[n_]  :=  f[n]  =   ... , and so is easily implemented:

facDynamic[0] = 1;

facDynamic[n_] := facDynamic[n] = 
        n facDynamic[n - 1]

Now, after computing the factorial of 3:

facDynamic[3]

a check on facDynamic reveals that the intermediate values have been recorded:

?facDynamic

The down-side of dynamic programming is that, because all intermediate  values or expressions are being stored, it is possible to run out of memory  using this approach. Another disadvantage, particularly when developing code, is that if you want to change any of the defining functions you need to Remove all rules, make the necessary modification, and then re-enter them.

.. Argument Checking

One defect of facDynamic is that its code only makes sense for positive integer arguments.

Conditions to check the validity of function arguments are easily  incorporated. The argument n to factorial must be Positive and an Integer. After removing the previous definitions for facDynamic:

Remove[facDynamic]

we can implement this in Mathematica syntax as

facDynamic[0] = 1;

\(TraditionalForm
    \`facDynamic[n_Integer?Positive] := 
      \(facDynamic[n] = n\ facDynamic[n - 1]\)\)

Here are four examples. The first evaluates as before,

facDynamic[3]

The second,

facDynamic[y]

fails because y is not an integer. The third

facDynamic[-5]

fails because -5 is not positive and, finally,

facDynamic[1.2]

fails because 1.2 is not an integer.

Mathematica does not require the user to specify the number or type of arguments of a 
function. Here are four reasons why :

An Aside

You can always generate an error message if a function is called with inappropriate arguments.  E.g., since factorial is a function of one argument, the following catch-all rule (for all other cases) issues an error message (using the built-in Message mechanism) if facDynamic is called with an incorrect number of arguments:

facDynamic[args___] := 
      Null /; Length[{args}] ≠ 1 && 
          Message[facDynamic::"argx", facDynamic, 
            Length[{args}]]

As a side-effect of evaluating the Condition — entered using the infix shorthand /;, which can be read as such that:

?/;

the rule on the right-hand side computes the number of arguments using Length and (denoted here by &&) prints a message when it is not equal to one (i.e., Length[{args}]  !=  1):

facDynamic[1, 2]

Note that our incorrect input is returned rather than the (dummy) Null because the conditional test failed. Also, we see that only output  corresponding to valid input has been stored:

?facDynamic

.. Tracing Computations

Mathematica includes a Trace function for tracing steps in a computation.  This can be used to verify  that intermediate factorial values are being saved and not recomputed. For  example, when computing facDynamic[5] only facDynamic[4] has to be computed as facDynamic[3] has already been computed above.

Tracing all computations of facDynamic (selected using the pattern facDynamic[_]) when evaluating facDynamic[5] verifies that this is indeed the case:

Trace[facDynamic[5], facDynamic[_]]

The Trace command can take an arbitrary patterns (see pattern-matching) to focus on particular aspects of a programming problem.  For example, here we only trace those computations of facDynamic[n] for which n > 24:

Trace[facDynamic[30], 
      facDynamic[n_ /; n > 24]]

. Functional Programming

Functional programming is very powerful and general.  Users who have not been exposed to procedural programming often accept and utilise functional programming as though it is second nature. Functional programming operators are powerful programming constructs; this is illustrated here through a number of examples.

.. Horner Form

One can use the Horner representation to convert a list of digits in an arbitrary base into its corresponding  number:

Horner[u_List, x_]:=Fold[Function[{a, b}, a x + b],0,u]

For example, 13247 in base 8 is

Horner[{1, 3, 2, 4, 7}, 8]

This can be checked using the built-in BaseForm function,

BaseForm[5799,  8]

using the notation base^^digits (in any base from 2 to 36),

8^^13247 == 5799

or by applying IntegerDigits in base 8.

IntegerDigits[5799, 8]

Importantly, Horner works for both numbers and symbols — it provides an efficient  represent polynomials (by minimising the number of addition and  multiplication operations):

Horner[{a, b, c, d}, x]

Furthermore, one can verify this program by running it for symbolic arguments. For example,

Expand[Horner[{a, b, c, d}, x]]

is obviously correct. In this way, coding errors can often be located.

As an aside, other applications of Horner representations include speeding up the computation of matrix powers and efficient computation of (x^n).

.. Jacobian Determinant

You've heard of inner products.  Here is an example of the use of the outer product (through Outer) to provide a natural implementation of the Jacobian matrix, required when changing integration variables:

\(TraditionalForm
    \`JacobianMatrix[f_List, vars_List] := Outer[D, f, vars]\)

For example, when changing from cartesian coordinates, {x, y}, to polar coordinates, {r, θ}, the Jacobian matrix is:

\(TraditionalForm
    \`JacobianMatrix[{r\ Cos[θ], r\ Sin[θ]}, {r, θ}]\)

This has determinant

Det[%] // Simplify

In functional programming, generalisation from a scalars to vectors (one dimension) and then to matrices (and even  to tensors) is often straightforward. For example, JacobianMatrix automatically works for an arbitrary number of variables.

.. Digital Roots

In numerology, one adds together all the digits of a number and repeats this process until the result does not change:

\(TraditionalForm
    \`DigitalRoot[n_] := FixedPoint[Plus@@IntegerDigits[#]&, n]\)

DigitalRoot[1234567]

.. Logistic Map

The logistic map is an example of a dynamical system:

(Subscript[Logistic, μ])_[x_] := μ x (1 - x)

Here we iterate (Subscript[Logistic, 3.5]) forty times starting at x = 0.3:

NestList[(Subscript[Logistic, 3.5]), 0.3, 40]

ListPlot[%, PlotJoined -> True];

.. Random Walk

Here is a routine for the two-dimensional, n-step random walk with equal (unit) step lengths. Starting at

z = {0, 0};

construct a routine for generating a uniformly distributed random number 
between {-1, 1}:

rand := Random[Real, {-1, 1}]

Each time you enter rand you get a different answer:

Table[rand, {5}]

Then construct a routine called RandomStep for generating a random 2 dimensional vector of unit length:

\(TraditionalForm
    \`\[LeftDoubleBracketingBar]x_\[RightDoubleBracketingBar] := 
      \@\(x . x\)\)

\(TraditionalForm
    \`Norm[x_] := 
      x\/\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\)

RandomStep := Norm[{rand, rand}]

For n steps form the cumulative sums (see FoldList above):

\(TraditionalForm
    \`RandomWalk[n_Integer] := FoldList[Plus, z, Table[RandomStep, {n}]]\)

E.g.,

RandomWalk[5]

Now plot the steps by putting a line through all the points:

\(TraditionalForm
    \`\(Show[Graphics[Line[RandomWalk[500]]], Axes -> True, 
      AspectRatio -> Automatic, AspectRatio -> Automatic]; \)\)

The heart of RandomWalk is the realisation that

FoldList[Plus, a, {b, c, d}]

forms the cumulative sum (see FoldList above).  This syntax trivially generalises to more dimensions:

FoldList[Plus, {a, b}, {{c, d}, {e, f}}]

. Pattern-matching and Replacement Rules

When simplifying expressions by hand one naturally tends to apply  transformation rules as required.  In addition, humans are very good at  pattern-recognition.  Mathematica, includes replacement rules and pattern-matching operations in a natural way.

Pattern-matching is a simple and elegant way to write programs, often  without requiring functional or procedural methods — although these can be also used. It is probably best explained by giving a number of   examples. Consider the replacement rule:

Clear[x, y, z]

{x, y, z} /. x -> 2 x

The  /.  shorthand denotes the infix form of the ReplaceAll operator, whilst  ->  is the infix form of Rule. In this case there was only one (exact) match.

If we introduce a generic variable x_:

{x, y, z} /. x_ -> 2 x

we see that it matches x, y, and z.

In the following sub-sections a number of examples — ranging from 
elementary to quite sophisticated — are used to illustrate various 
aspects of pattern-matching and replacement rules.

.. Predicates

Part of the power of pattern-matching lies in its ability to compactly   represent more complicated patterns. Also, one can combine pattern-matching with predicate or testing operators using  /;  and ?.

For example, here we convert all numbers less than 3 to 0:

{1, 2, 3, 4, 5} /. x_ /; x < 3 -> 0

By convention, names of built-in testing predicates in Mathematica end in Q (there are a few other built-in functions whose names end in Q that are not predicates).

In general, built-in predicates return the Boolean values True or False. For example,

VectorQ[{a, b}]

Here we use ?PrimeQ (this is not finding out information on PrimeQ — instead it is a syntax for applying a test) to identify all   prime numbers in a list and then use a replacement rule to square them:

{1, 2, 3, 4, 5} /. x_?PrimeQ -> (x^2)

.. ReplaceList

ReplaceList is a very useful pattern-matching tool. Consider:

\(TraditionalForm
    \`ReplaceList[{a, b, c}, {x_, y_} :> {{x}, {y}}]\)

Here, the pattern {x_, y_} cannot match the expression {a, b, c} so we get an empty list as output. However, if the pattern _

?_

is replaced by the pattern __

?__

in one of the pattern variables, then there is one possible match:

\(TraditionalForm
    \`ReplaceList[{a, b, c}, {x__, y_} :> {{x}, {y}}]\)

The general pattern involving ___

?___

takes into account all possible matches:

\(TraditionalForm
    \`ReplaceList[{a, b, c}, {x___, y___} :> {{x}, {y}}]\)

.. Number Puzzle

For the sequence

seq = {6, 44, 30, 15, 24, 12, 33, 23, 18, 8};

it is easy to construct a (user-defined) pattern which uses a test (in  this case Equal) to find (the first) four numbers totalling to 100:

\(TraditionalForm
    \`seq /. {___, a_, ___, b_, ___, c_, ___, d_, ___} :> 
        {a, b, c, d} /; a + b + c + d == 100\)

. Association

The Object-Oriented Programming (OOP) idea of association is very powerful 
 and incorporated into the Mathematica language in a natural way. Often, a  definition or rule should be 
associated with a particular object rather than with a general operation.

If one enters:

g[a_] + g[b_] := g[a + b]

an error message is generated indicating that Plus—the top-level  operator on the left-hand side—is a built-in protected operator and  that this input is effectively trying to attach a rule to Plus.  In this situation, what the user really wants is to be able to associate a specific (user-defined) rule whenever an expression involving a sum of g functions arises.  The mechanism for handling this is the association operation which has the shorthand syntax

?/:

After entering

g /: g[a_] + g[b_] := g[a + b]

the rule is associated with g (instead of being associated with +) using the syntax g /:

?g

This rule is applied repeatedly to expressions until they no longer change.  For example,

g[a] + g[b] + g[c]

This uses pattern-matching — and the associativity (Orderless)  property of +:

Attributes[Plus]

since the rule above was only defined for the sum of two terms.

. Simple Programming

In Notebooks one can do computation without writing programs by forming a  collection of basic functions,  selecting the range of cells, and updating  the cells — similar to working in a spreadsheet.  Then, if required, cutting, pasting, and copying can be used to construct a more robust and  general procedure. Here is an example of this approach.

.. Sequence of Steps

The sequence of steps below displays 10 randomly chosen points in the interval [0, 1] × [0, 1]:

Table[Random[], {10}, {2}];

Point /@ %;

Show[Graphics[{PointSize[0.03], %}], 
        AspectRatio -> Automatic, PlotRange -> All];

This approach is analogous to using a spreadsheet.

.. Converting a Sequence of Steps into a Module

To convert these steps into a program, we collect the above steps together into a Module and introduce a (dummy) local variable (pts):

\(TraditionalForm\`RandomPointPlot[n_, s_:  0.03] := 
      Module[{pts}, 
		pts = Table[Random[], {n}, {2}]; 
		pts = 
          Point /@ pts; 
		Show[Graphics[{PointSize[s], pts}], 
          AspectRatio -> Automatic, PlotRange -> All];]\)

A good approach to programming in Mathematica is to first find the appropriate "sequence of steps" and then, when everything is working well and has been tested in all the cases of interest, convert these steps into a procedure as shown in this example.

. Blocks and Modules — Variable Scoping

In Mathematica there are two special programming constructs — Block and Module — which, like their names suggest, are used to delimit blocks of  code. The variables contained within these delimiters are local — in  the sense illustrated below. The following simple example clearly shows the  difference between Block and Module.

This assignment focuses on the elements of Mathematica and not on how to write complete Mathematica programs. An excellent introduction to programming in Mathematica is found in [Maeder].

.. Block

Define

m = (i^2)

In a Block, the local value of i is used in the body:

Block[{i = a}, i + m]

Trace clearly shows the sequence of steps:

Trace[Block[{i = a}, i + m]]

Block localises the values of a variable but not the symbol itself. This is  called dynamic scoping.

.. Module

In a Module, the i that appears explicitly in the body is treated as a local  variable:

Module[{i = a}, i + m]

In other words, the i inside the Module is different to the external i:

Trace[Module[{i = a}, i + m]]

Module localises the symbols. This is called lexical scoping.

.. Block or Module?

The question arises as to when one should use Block and when one should use Module. The answer is that most of the time one should use Module. A  particular application that requires Block instead of Module is when temporarily changing the value of a system variable. Consider the following recursive procedure for computing factorial:

\(TraditionalForm
    \`fac[n_] := 

      Module[{h, m}, 
		h[0] = 1; 
		h[m_] := \(h[m] = m\ h[m - 1]\); 

        		h[n]
]\)

Everything seems to work fine:

fac[20] == 20!

Note also that h is "hidden":

?h

However, when we try a larger number we run into a problem:

fac[256]

It is apparent that we have hit a recursion limit.  This is controlled by  the global variable, $RecursionLimit:

$RecursionLimit

Here is a modified routine (overwriting the previous one) which uses Block to temporarily change the recursion limit and Module to do the computation:

\(TraditionalForm
    \`\(fac[n_] := 

      Block[{$RecursionLimit = ∞}, 
		
        Module[{h, m}, 
			h[0] = 1; 
			
          h[m_] := \(h[m] = m\ h[m - 1]\); 
			h[n]
		]
]\ \)\)

Now we have no problems with a much larger calculation:

fac[500] == 500!

and the global $RecursionLimit keeps its usual value:

$RecursionLimit

. Summary

To summarise, here are some basic points relevant to Mathematica programming: