Programming with Mathematica

Evaluation

Entering Commands

It is assumed that those taking this course are familiar with the following conventions used to type Mathematica input.

Built-in Mathematica symbol names are typically one or more full-length English words (or a standard abbreviation), where the first letter of each word in the name is capitalized. Mathematica is case sensitive. Typical built-in symbol names are Integrate, FactorInteger, Pi, and Cos.

Function arguments are typed inside square brackets, as in Factor[1 - x^9].

In these notes, formatting functions will be entered on occasion using postfix notation, as in the following.

{{a, b}, {c, d}} // MatrixForm

Parentheses are used only for grouping, as in (a + b)^2 and 1/(2  π).

Collections of objects (such as vectors or arrays of data) are called lists, and are entered as a sequence of elements inside curly brackets {} and separated by commas. The elements of a list can themselves be lists; such nested lists are used to represent matrices, tensors, or arbitrary arrays. Lists need not be “rectangular”. A list of integers ranging from 1 to 3 is entered {1, 2, 3}.

Strings are entered as a sequence of characters inside double quotes, as in "This is a string.".

One way to refer to previous results is to use percent sign (%) notation. A single percent sign % refers to the last expression returned, a double percent sign %% refers to the expression returned two evaluations ago, and so forth; a percent sign followed immediately by a number n (%n) refers to the output line labeled Out[n].

We use the on-line help by entering ?FunctionName, which displays the usage message for FunctionName. The usage message for a function typically describes the function and the types of arguments that the function expects. Here is the usage message for Integrate.

?Integrate

If we do not know the exact name of a function, we can print a listing of functions containing a given string by entering ?*String*. Here is a listing of the functions that contain the word “Plot”.

? *Plot*

Many functions accept options, or named arguments, that control the behavior of the function. Options have the form OptionName->Value, and are entered after the function’s required arguments. Options[FunctionName] returns a list of the options a function accepts and their default values. By default, FactorInteger factors an integer over the real integers.

FactorInteger[12]

Here are the options that FactorInteger accepts.

Options[FactorInteger]

By setting the option GaussianIntegers to True, Mathematica will factor an integer over the Gaussian (complex) integers.

FactorInteger[12, GaussianIntegers -> True]

Many of Mathematica’s built-in mathematical functions and plotting commands are used in these course notes. These are not explained if the meaning can be gathered from the usage message.

How Mathematica Evaluates Input

Every valid input to Mathematica is an expression. While Mathematica allows us to type in many different kinds of questions and commands, all expressions have the same basic internal representation.

For example, we enter a sum by typing the following.

a + b

Internally, however, Mathematica represents the sum in the following way.

Plus[a, b]

Here is a list of two symbols.

{a, b}

Here is the internal form of the list.

List[a, b]

A typical Mathematica expression has the following form (where each of the (Subscript[e, i]) can itself be a general expression).

(Subscript[e, 0])[(Subscript[e, 1]),  (Subscript[e, 2]),  …,  (Subscript[e, n])]

The first expression (Subscript[e, 0]) is called the head of an expression, and can be extracted by evaluating Head[expression].

Head[a + b]

Head[{a, b}]

Head[Sin[1 + a]]

FullForm allows us to view the entire internal form of an expression. Here are the internal forms of some expressions.

{FullForm[a + b + c], FullForm[a == b], FullForm[Sin[2  a]]}

In some cases, use of the function Hold is necessary to prevent evaluation of an expression. In the following example, the argument to FullForm is evaluated before the internal form can be printed.

FullForm[1 + 1]

Hold prevents evaluation of an expression.

FullForm[Hold[1 + 1]]

Some expressions are atomic, and do not have parts that can be extracted. Examples are numbers, strings, and symbols. The head of an atomic expression is a tag that describes its type.

{Head[1234], Head[1.234], Head[(23/45)], Head[1 + I], Head["str"], 
      Head[undefined]}

Note that the head of an atomic expression does not appear in its FullForm representation.

FullForm[7]

Some Mathematica expressions are compound expressions, sequences of complete expressions separated by semicolons. While it is occasionally safe to separate the expressions in a compound expression with newlines instead of semicolons; it is recommended that semicolons be used to avoid ambiguity.

fib[0] = 1; fib[1] = 1

Here is the internal form of the previous compound expression.

FullForm[Hold[fib[0] = 1; fib[1] = 1]]

Mathematica generally evaluates a general expression of the form (Subscript[e, 0])[(Subscript[e, 1]),  (Subscript[e, 2]),  …,  (Subscript[e, n])] from left to right, first evaluating the head (Subscript[e, 0]), then recursively evaluating each of the (Subscript[e, i]) in turn. (That is, Mathematica performs depth-first evaluation.) We can alter the evaluation of an expression by using Hold and its variants, and there are several classes of built-in functions that have nonstandard evaluation sequences.

Mathematica typically returns the value of an expression after it is evaluated, meaning we can refer to the result using % notation. In the case of a compound expression, the value returned is the result of the last expression in the compound expression.

The following input returns a numerical value.

2^50

The value returned can be used in further calculations.

%^(1/25)

In the following compound expression, only the value of the last expression is returned.

(2^25);  (3^100);  (4^50)

%^(1/10)

If the last expression in a line of Mathematica input is followed by a semicolon, the result will not be printed, but can still be referred to using %. In this example, Mathematica computes the value of (2^50), but does not print the result.

2^50;

We still refer to the result using %.

%^(1/25)

Some expressions do not evaluate to a value, but display a result as a side effect; in such cases, the value returned is Null. For example, the command Print prints its argument after evaluating it, but does not return a value.

Print["Print commands return nothing."]

The value returned by the Print statement is Null.

FullForm[%]

Simply entering the string, on the other hand, returns the value of the string.

"This string returns a value."

FullForm[%]

Compound expressions can be used inside other functions, where normally a simple expression would be used. For example, if we plot the compound expression (Print[x]; x) Mathematica creates a graph of the values returned by the compound expression (the values returned are the values of x), but prints each value of x that is used to create the graph.

Plot[(Print[x]; x), {x, 0, 1}, PlotPoints -> 5];

Of course, Mathematica has many rules for simplifying expressions that involve built-in functions. For example, Mathematica recognizes that the expressions Sin[π/4], 2  x - x, Integrate[Exp[-x^2], x], and Plot[Sin[t]/t, {t, -10, 10}] should be replaced with their respective results.

Styles of Programming

Programming Overview

Mathematica allows us to program in several different styles, including procedural, rule-based, mathematical, and functional programming, mirroring the styles of programming used in BASIC, Fortran, C, Lisp, APL, and many other languages. An important point to remember is that, while Mathematica contains programming structures that are similar to those in other programming languages, we must still adhere to Mathematica syntax and punctuation conventions.

A straightforward example for illustrating the different styles of programming is the implementation of the factorial function. The factorial function is built into Mathematica; the factorial of a positive integer n is defined as the product of the integers 1, 2, …, n, and it can be entered into Mathematica as n!.

35!

In the following discussion, we show four different ways to program our own factorial function.

Procedural Programming

One way we can compute a factorial is procedurally, as we would in Fortran or BASIC, using a Do loop. First, we initialize the value of a partial-product variable temp to 1; then we multiply temp by each of the numbers from 1 to the desired value (say, 35) and print the final value of temp. Here is a procedural way to compute 35 factorial.

\(\(temp = 1; \  (*\ initialize\ the\ product\ *) 


    Do[temp = temp*n, {n, 1, 35}]; 
     (*\ loop\ through\ n = 1, 2, …, 35\ *) 

temp\ 
     (*\ return\ the\ result\ *) \)\)

We verify the answer is correct by testing the equality of temp and the built-in factorial function.

temp == 35!

Mathematica contains many functions that control the flow of a procedural program. The functions Do, For, and While set up different types of loops; If, Which, and Switch direct the flow of an evaluation according to conditions that are met; and the pairs Goto-Label and Throw-Catch permit arbitrary branching within a group of expressions.

As a general rule, procedural programming is useful for translating existing programs written in other languages. However, the functional and rule-based styles of programming, discussed below, use operations that directly manipulate Mathematica’s internal representation of expressions, and are therefore more efficient and usually more elegant.

Rule-Based Programming

Instead of defining a procedure, we can define functions to perform computations. For example, we can define a function fac that computes factorials.

It is a property of the factorial function that m! = m(m - 1)! for positive integer values of m. That is, we can express a factorial m! in terms of an easier case, (m - 1)!. This is called a recursive definition, and in Mathematica we can define a function recursively, as in C or Pascal. First we tell Mathematica that fac of a number m is equal to m times fac of m - 1. (Later we will see how to change the definition of fac so that it accepts only positive integer values of m.)

fac[m_] := m fac[m - 1]

The definition is not complete, however, because there is currently no way for the function to stop calling itself. It is necessary to define a stopping condition, or elementary case, for the definition. In this case we define fac of 0 to be 1, corresponding to the definition 0! = 1.

fac[0] = 1;

We are now ready to compute factorial values using fac.

fac[35]

We can compare fac to the built-in factorial function.

fac[35] == 35!

Mathematical Programming

In many cases there is a built-in mathematical function that performs the desired operation. Another way to compute a factorial is to multiply the desired numbers using the built-in function Product. Our function fac2 computes the factorial of m by multiplying values of j as j ranges from 1 to m.

fac2[m_]  :=  Product[j,  {j, 1, m}]

fac2[35]

Functional Programming

The most abstract and most efficient style of programming in Mathematica involves using high-level Mathematica programming constructs, such as Apply and Map. While procedures and recursive functions have to compute many values in the course of a computation, using high-level functions allows us to work with large structures as a whole.

Here is a third factorial function, called fac3, programmed functionally; that is, by applying an operator (in this case, Times) to an entire list of values at once. We generate a list of integers from 1 to m, using the built-in function Range, then effectively insert multiplication signs between the numbers.

fac3[m_]  :=  Apply[Times,  Range[m]]

fac3[35]

We can easily spot-check our answers by testing for equality among the built-in factorial, fac, fac2, and fac3.

35!  ==  fac[35]  ==  fac2[35]  ==  fac3[35]

Exercises: Evaluation

Use FullForm to inspect the internal form of each of the following expressions. (In some cases it will be necessary to use Hold.)

Sin[1/2 + 1/x]

Note that the internal form of an exact rational number a/b is Rational[a, b].

FullForm[Sin[1/2 + 1/x]]

f[x_] := x^x

A function definition made with  :=  returns the value Null, so it is necessary to wrap the definition in Hold.

FullForm[Hold[f[x_]  :=  x^x]]

g'[x]

The function Derivative is a functional, meaning it operates on other functions.

FullForm[g'[x]]

Notice that the head of a derivative is an expression, and not just a symbol.

FullForm[Head[g'[x]]]

{(x^y)^z,  x^(y^z)}

This exercise shows the effect of grouping expressions in different ways.

FullForm[{(x^y)^z, x^(y^z)}]

{x == y, x === y}

In this case, it is necessary to wrap the expression in Hold because x === y evaluates to the value False.

FullForm[{x == y, x === y}]

Here we use Hold.

FullForm[Hold[{x == y, x === y}]]

GaussianIntegers -> True

FullForm[GaussianIntegers -> True]

Functions

Defining Variables

We define variables in Mathematica by entering name = value, using a single equal sign. (The internal form of name = value is Set[name, value].) Variable names can be as long as desired, and can be any combination of letters (both upper- and lower-case) and numbers, with the restriction that a variable name cannot begin with a number. Note that a variable name cannot contain an underscore character, so var_name is not a valid variable name.

To set the variable called newvar equal to 15, we enter the following.

newvar = 15

Wherever newvar is used, it is replaced with its assigned value.

Sqrt[newvar]

Variables can take symbolic values as well, so the following is a valid assignment. (We keep Mathematica from printing the result of a variable assignment by ending the input line with a semicolon.)

zzz = (xxx^2) + (yyy^2);

The right-hand side of a variable assignment can be a function or program, in which case Mathematica sets the value of the variable equal to the value returned after evaluating the right-hand side. Here we set the value of solutionset to be the result of a Solve command.

solutionset = Solve[(x^2) == x + 1, x]

We can even set a variable equal to a function name. For instance, we can set the variable int equal to the built-in Integrate function.

int = Integrate;

We can now use int in places we would use Integrate.

int[Exp[-(x^2)], x]

It is important to realize that variable assignments are permanent. The values of newvar, zzz, solutionset, and int will remain in memory until we quit Mathematica, or until we use Clear to tell Mathematica to clear the value of the variable. Here is a calculation that uses the value of newvar defined above.

(newvar^2) + newvar + 200

Clear[newvar] clears the value of newvar.

Clear[newvar]

The symbol newvar no longer has a value.

(newvar^2) + newvar + 200

Here we clear the values of the other variables defined earlier.

Clear[zzz, solutionset, int]

(To clear the value of every variable and function defined in the current session without having to enter a complete list, type Clear["Global`*"]. The form Global` is explained in the chapter “Packages”.)

Mathematica applies the standard rules of algebra to symbols that have not been assigned values. One common use of symbols without definitions is to treat them as units.

30  year (365  day/year) (24  hour/day)

In some cases we want to reevaluate the right side of a variable definition each time the variable is used. For example, we may wish to define a variable rnd that returns a random integer from 1 to 10, so that it returns a different number each time. To tell Mathematica to reevaluate the right side of an assignment each time the variable is used, we use the colon-equal sign :=. The internal form of name := value is SetDelayed[name, value].

The following definition of rnd works as desired. Notice that no value is returned when we use :=.

rnd := Random[Integer, {1, 10}]

Every time we use the name rnd, Mathematica reevaluates the right-hand side of the assignment, and returns a (most likely) different random number.

{rnd, rnd, rnd}

If we were to use the equal sign, Mathematica would evaluate the right-hand side of the assignment immediately, and return the same value every time the variable is used.

samernd = Random[Integer, {1, 10}]

{samernd, samernd, samernd}

Assign the symbol const the value of an approximation to (1/(2  π)) with 30 significant digits.

We use N to obtain the numeric approximation.

const = N[1/(2 Pi), 30]

The constant can now be used inside other expressions.

const (E^(-(x^2)))

Define a variable named AorB that returns either the symbol A or the symbol B, with a 50 percent chance of each.

If we want to choose randomly between A and B, and have the same value returned each time AorB is used, we use the equal sign = in the definition.

AorB = If[Random[] < 0.5, A, B];

{AorB, AorB, AorB, AorB}

If we want a (possibly) different value to be returned each time we use AorB, we use the colon-equal sign := in the definition, which tells Mathematica to reevaluate the right-hand side of the definition each time we use AorB.

Clear[AorB]

AorB := If[Random[] < 0.5, A, B];

{AorB, AorB, AorB, AorB}

Introduction to Functions and Patterns

We define a function by giving it a name, a template that describes its arguments (called a pattern), and an expression that describes the operations to be performed on the arguments. On the left side of the following definition, sinc is the name of the function, x is the name of the pattern, and the underscore _ is a pattern, or blank, that will match any single argument; on the right side of the definition is an expression that describes the operations to be performed on an expression that matches the pattern.

sinc[x_] := Sin[x]/x

Notice that we do not include underscore characters on the right-hand side of the definition, and that nothing is returned after making the function definition.

When we enter the following expression, Mathematica looks in its rule table for sinc and sees that it contains a rule for evaluating sinc of a single expression; the single expression matches the pattern x_, and Mathematica replaces every occurrence of x on the right-hand side of the definition with the argument.

sinc[1.23456]

It is important not to leave a space between a pattern’s name and the underscore, though we can separate the pattern name and the pattern structure with a colon, if desired.

sinc[x : _] := Sin[x]/x

sinc[(π/3)]

The pattern x_ matches any single expression, so the expression a + b is a valid argument to sinc.

sinc[a + b]

The pattern x_ will not match zero arguments or a sequence of more than one argument. Mathematica returns an expression unevaluated when it has no definitions to apply, as in the case of sinc called with zero arguments.

sinc[]

Similarly, sinc does not evaluate if given a sequence of two arguments.

sinc[1, 2]

A list is considered a single argument, so sinc will evaluate if given a list of arguments.

sinc[{1, 2, 3}]

The name we give to a pattern is unimportant, and it is treated as a dummy variable. The pattern x_ on the left-hand side of the definition of sinc could have expressed as num_, z_, or anything else.

Notice that we used a colon-equal sign := (SetDelayed) to separate the left and right sides of the definition of sinc. It is usually safe, as a rule of thumb, to use the colon-equal sign whenever we define functions.

The expression on the right side of a function definition can be as complicated as desired. Functions we define can call other built-in and user-defined functions.

\(complexPlot[z_] := 
	
      ListPlot[
		Transpose[{Re[z], Im[z]}], 
        PlotStyle -> PointSize[0.03], 
		PlotRange -> All]\)

complexPlot[N[Fourier[{1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1}]]];

The right side of a function definition can be a compound expression, in which case the value returned by the function is the value of the last expression in the definition. We must enclose the entire compound expression in parentheses to make it clear that the entire expression is part of the definition.

sinc[x_] := (Print["The argument to sinc is ", x]; Sin[x]/x)

sinc[-1.1]

Evaluating sinc[a, b] does not work because we have not entered a definition for sinc when it is given any number of arguments other than one.

sinc[a, b]

To give a definition to sinc when it is applied to two arguments, we add the following definition.

sinc[x_, y_] := Sin[x]/y

Now, when we give sinc two arguments, the first matches the pattern x_ and the second matches the pattern y_. The arguments that match these patterns are substituted wherever x and y appear in the right-hand side of the definition.

sinc[a, 7]

Notice that we can have several definitions attached to the same symbol. We can find all the definitions Mathematica has for a function by entering ? FunctionName.

? sinc

The value a function returns does not have to be a number or a simple expression. The following function returns a list of graphics primitives, with the arguments substituted for the appropriate dummy variables.

CircledDigit[d_, {x_, y_}] := {GrayLevel[0], Disk[{x, y}, 0.05], 
        GrayLevel[1], Text[ToString[d], {x, y}]}

Here is the result of one use of CircledDigit.

CircledDigit[8, {0, 1}]

The function can be used in conjunction with other functions, in this case inside a Plot statement.

\(\(Plot[Sin[x\^2], {x, \(-1\), 1}, 
	
      Epilog -> 

        		{CircledDigit[1, {0.50, 0.25}], 
          CircledDigit[2, {\(-0.75\), 0.53}]}, 
      AspectRatio -> Automatic]; \)\)

Sometimes we may wish to define a function to act only on arguments with a particular form. For example, we can define a function called recip that works only on symbolic rational expressions, and returns the reciprocal of its argument.

recip[a_/b_] := b/a

If recip is given a symbolic rational expression, the numerator matches the pattern a_, the denominator matches the pattern b_, and Mathematica returns b/a.

recip[c/d]

recip[((u + v)/w)]

The argument u + v does not match the pattern a_/b_, so recip[u + v] returns unevaluated.

recip[u + v]

Note that the argument 1/2 does not match a_/b_ because the internal form of 1/2 is Rational[1, 2].

recip[(1/2)]

It is acceptable to leave a pattern unnamed if we do not refer to the corresponding argument on the right side of the definition. The following definition accepts any expression as an argument, but does not use the expression when computing a value.

ignore[_] := "Ignored."

ignore[10]

A useful function for investigating which expressions match a given pattern is MatchQ, which takes an expression and a pattern, and returns True or False to indicate if the given expression matches the given pattern.

Here we test if several expressions match the pattern a_ + b_.

MatchQ[c + d, a_ + b_]

MatchQ[(u + v) + w, a_ + b_]

The expression 3  z + z does not match the pattern a_ + b_ because Mathematica evaluates the expression to 4  z, which is not in the form of a sum.

MatchQ[3 z + z, a_ + b_]

Define a function dist that takes two arguments, x and y, and returns Abs[x - y].

The definition is straightforward: we give two patterns, x_ and y_, to match the two arguments.

dist[x_, y_] := Abs[x - y]

Here we verify that the function works as desired.

dist[-1, 6]

Define a function avg that returns the average of its two arguments.

Again, the definition is straightforward.

avg[a_, b_] := (a + b)/2

Here is a test of the function.

avg[10, 10000]

Create a function dx that accepts a single argument, and returns the derivative of its argument with respect to x.

Because x is not the name of a pattern on the left side of the definition, Mathematica computes the derivative with respect to the literal symbol x.

dx[expr_] := D[expr, x]

Here is a test of the function.

dx[Sin[Exp[x]]]

Limiting Arguments to Certain Types

We can limit the type of arguments that a function will accept by putting the head of an acceptable argument after the underscore in the pattern. The following function mersenne accepts only one argument n with the head Integer, and returns the Mersenne number corresponding to the (n^th) prime number (M((Subscript[p, n]))).

mersenne[n_Integer] := 2^Prime[n] - 1

The argument 10 has the head Integer, so it matches the pattern n_Integer.

mersenne[10]

Because mersenne has not been defined for noninteger arguments, an attempt to use mersenne with a noninteger argument returns unevaluated.

{mersenne[1.5], mersenne[π], mersenne[var]}

The head of an equation is Equal.

Head[(x^2) == 1]

We can create a function that recognizes equations by using the pattern _Equal on the left side of the definition. The following function solveX solves a given equation with respect to x. (Because x is not the name of a pattern on the left side of the definition, solveX solves with respect to the literal symbol x.)

solveX[eqn_Equal] := Solve[eqn, x]

solveX[(x^2) == 1]

Built-in mathematical constants such as Pi and E have the head Symbol.

If a function can accept any of several types of arguments, we use alternative patterns of the form (Subscript[patname_head, 1]) | (Subscript[patname_head, 2]). Any number of alternatives are allowed.

IsASymbolOrInteger[x_Integer | x_Symbol] := True

IsASymbolOrInteger[10]

IsASymbolOrInteger[Pi]

IsASymbolOrInteger[22/7]

In many cases we define a function that accepts a list as an argument. If the length of the list is unimportant, we can use the pattern _List on the left side of the definition.

reverseAnyList[L_List] := Reverse[L]

reverseAnyList[{1, 2, 3}]

reverseAnyList[{1, 2, 3, 4}]

If the length of the list is significant, we explicitly enter the allowable form. The following function accepts only a list with three elements.

reverseTriple[{x_, y_, z_}] := {z, y, x}

reverseTriple[{1, 2, 3}]

The argument here has four elements, so it does not match the given pattern, and the function returns unevaluated.

reverseTriple[{1, 2, 3, 4}]

It is important to keep in mind the difference between a function that accepts a list of several arguments and a function that accepts a sequence of several arguments. The pattern inside reverseTriple does not match the following input because the arguments are in the form of a sequence of three arguments, not a list of three arguments.

reverseTriple[1, 2, 3]

Define a function add3list that accepts a list of three integers and returns the sum of the arguments. Define a function add3seq that accepts a sequence of three integers and returns their sum.

The functions add3list and add3seq differ only on their left sides: add3list contains list brackets surrounding the patterns, while add3seq does not.

add3list[{x_Integer, y_Integer, z_Integer}] := x + y + z

add3seq[x_Integer, y_Integer, z_Integer] := x + y + z

Here are tests of the functions.

add3list[{1, 100, 10000}]

add3seq[1, 100, 10000]

Neither function accepts an argument intended for the other.

add3list[1, 2, 3]

add3seq[{1, 2, 3}]

Create a function doubleString that accepts a string and returns two copies of the string joined together. (Use the function StringJoin on the right side of the definition.)

To limit arguments to those with the head String, we use the pattern s_String.

doubleString[s_String] := StringJoin[s, s]

doubleString["Twice"]

Conditions

Beyond limiting arguments to certain types, we can limit them according to given conditions by using the conditional operator slash-semicolon /;. Expressions of the form expr /; cond have the internal form Condition[expr, cond], and we read them as “expr provided that cond”. Functions commonly used in conditions are the tests for equality == and !=, which respectively test if two expressions are equal or unequal; the numerical tests <, <=, >, >=; and the built-in predicate functions TestQ, which are explained in the next section.

An example of a function requiring a condition is the function mersenne, which should accept only positive integers; nonpositive arguments cause a warning message.

mersenne[ -2 ]

To limit the pattern n_Integer to positive integers, we tell Mathematica to apply the definition only if n is greater than zero by stating the condition after the /; operator. (It is first necessary to clear the previous definition for mersenne.)

Clear[mersenne]

mersenne[n_Integer] := (2^Prime)[n] - 1 /; n > 0

Mathematica no longer tries to evaluate mersenne for nonpositive arguments.

{mersenne[15], mersenne[-15]}

Alternatively, we can include the condition on the left side of the definition, directly following the pattern.

Clear[mersenne]

mersenne[n_Integer /; n > 0] := (2^Prime)[n] - 1

Again, the function works as intended.

{mersenne[15], mersenne[-15]}

Using the same technique, we define a square-wave function. We define squarewave[x] to be -1 when the greatest integer less than x is even, and 1 when the greatest integer less than x is odd.

squarewave[x_] := -1 /; EvenQ[Floor[x]]

squarewave[x_] := 1 /; OddQ[Floor[x]]

Here is a plot of squarewave.

Plot[squarewave[x], {x, -2.5, 2.5}, 
      PlotStyle -> Thickness[0.02]];

Note that we can use the conditional function If to achieve the same results.

squarewave2[x_] := If[EvenQ[Floor[x]], -1, 1]

The graph of squarewave2 is identical to the graph of squarewave.

Plot[squarewave2[x], {x, -2.5, 2.5}, 
      PlotStyle -> Thickness[0.02]];

We define compound predicate functions by using the logical functions And (&&) and Or (||). The following is one rule in a user-defined integration function (the power rule). Following the body of the definition is the compound condition that the exponent n must be free of the variable of integration v, and the exponent n must not equal -1.

\(myIntegrate[v_\^n_, v_] := 
      v\^\(n + 1\)\/\(n + 1\) /; FreeQ[n, v] && n != \(-1\)\)

Here is a list of two results from myIntegrate, the first with an argument that satisfies the compound condition, and the second with an argument that does not.

{myIntegrate[(x^13), x], myIntegrate[(1/x), x]}

Define a function uniform that accepts a single real argument and returns the value 1 if its argument is between 0 and 1, and returns 0 otherwise.

Here is the first case, when 0 < x < 1.

uniform[x_Real] := 1 /; 0 < x < 1

Here is the second case, applicable everywhere else. It is not necessary to put a condition testing if x is outside the interval (0, 1) because Mathematica will first attempt to apply the more specific rule defined above.

uniform[x_Real] := 0

Here is a plot of uniform.

Plot[uniform[t], {t, -0.5, 1.5}, 
      PlotStyle -> Thickness[0.03], AspectRatio -> Automatic];

Note that the pattern x_Real prevents uniform from returning a value when given an integer or rational argument.

uniform[(1/2)]

Create a piecewise function of two variables that acts differently inside the unit circle in the x-y plane than outside.

Here is one such function, first defining the case where x and y fall inside the unit circle with the condition Sqrt[x^2  +  y^2]  <  1.0.

pw[x_, y_] := 1 - ((x^2) + (y^2)) /; Sqrt[(x^2) + (y^2)] < 1.0

Here is the second case.

pw[x_, y_] := 0

Here is a graph of our function.

Plot3D[pw[t, u], {t, -2, 2}, {u, -2, 2}, 
      PlotPoints -> 35, Mesh -> False, PlotRange -> All];

Define the power rule for differentiation by defining a function that accepts two arguments, an expression of the form (v^n) and a variable v, and returns n (v^(n - 1)), provided that the exponent n is free of the variable v.

To test that the exponent is free of the variable of differentiation, we use FreeQ.

\(deriv[v_\^n_, v_] := n\ v\^\(n - 1\) /; FreeQ[n, v]\)

The following case satisfies the condition.

deriv[(x^7), x]

The following case does not.

deriv[(w^w), w]

Conditional Patterns

In the previous section we tested conditions after the patterns in a function were matched; using conditional patterns, we make a condition part of a pattern. A conditional pattern has the form patname_?condition (or patname_head?condition), and an argument matches a conditional pattern only if condition returns the value True.

Mathematica has a function called Positive that returns True if its argument is positive.

{Positive[-15], Positive[0], Positive[15]}

To include the condition Positive in the pattern of the function mersenne, we make the following definition (after first clearing the previous definition).

Clear[mersenne]

mersenne[n_Integer?Positive] := (2^Prime)[n] - 1

The argument 20 matches the pattern n_Integer?Positive, and the function evaluates to a value.

mersenne[20]

Nonpositive arguments now return unevaluated.

{mersenne[0], mersenne[-20]}

Mathematica has many built-in predicate functions, functions that return True if their condition is met, and False otherwise. By convention, predicate functions end with the letter Q. Here are examples of some built-in predicate functions. (A complete listing of functions that end with Q can be found by entering ?*Q or Names["*Q"].)

EvenQ[12]

PrimeQ[(2^((2^4))) + 1]

VectorQ[{1, 2, 3, 4}]

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

IntervalMemberQ[Interval[{3, 4}], π]

StringQ["Is this a string?"]

Any predicate function can be included in a pattern. Here is a custom predicate function that returns True if its argument is inside the unit circle in the complex plane.

pred[z_] := True /; Abs[z] <= 1.0

We can now use pred as part of a conditional pattern.

invert[x_?pred] := 1/x

The following argument satisfies pred, so the function returns a value.

invert[0.02 - 0.03 I]

The following argument does not satisfy pred, so the function returns unevaluated.

invert[3 + 4 I]

Using a conditional pattern, create a function that returns the symbol Yes if its argument is a prime number and No otherwise.

The following conditional pattern tests if its argument is prime.

isItPrime[n_?PrimeQ] := Yes

This pattern matches any argument that is not a prime number.

isItPrime[n_] := No

Here is a test of the function.

{isItPrime[17], isItPrime[28]}

Define a function called vectorLength that returns the length of a vector.

Here is a definition of vectorLength that accepts only vectors by using the conditional pattern v_?VectorQ. It computes the length using the formula (Sqrt[v.v]).

vectorLength[v_?VectorQ] := (Sqrt[v.v])

Here is the length of a vector.

vectorLength[{1, 2, n, 4}]

Here the argument to vectorLength is not a vector, so the function returns unevaluated.

vectorLength[1 + n]

Define a function called collatz that has the following recursive definition.

collatz(n)={({{collatz(n/2),n even}, {collatz(3n + 1), n odd,  n ≠ 1}, {1,  n = 1}})

Here we use the predicate functions EvenQ and OddQ as part of a conditional pattern.

collatz[n_?EvenQ] := collatz[n/2]

collatz[n_?OddQ] := collatz[3 n + 1]

collatz[1] = 1;

The specific case collatz[1] gets matched before the general case collatz[n_?OddQ], so it is not necessary to add a test n ≠ 1 to the latter rule.

It is conjectured that for any positive integer argument, collatz as defined here will return 1 (as opposed to going into an endless loop). Here is some sample input.

collatz[10]

collatz[100]

Here is a modification that prints all of the numbers that are encountered during the evaluation of collatz.

Clear[collatz]

collatz[n_?EvenQ] := (Print[n]; collatz[n/2])

collatz[n_?OddQ] := (Print[n]; collatz[3 n + 1])

collatz[1] = 1;

Here are the numbers encountered during the evaluation of collatz[10].

collatz[10]

The function Count returns the number of elements in an expression (typically a list) that match a given pattern. Use Count, CoefficientList, and a conditional pattern to count the number of odd coefficients in the expansion of (((1 + x))^11).

Here is the list of coefficients.

CoefficientList[(((1 + x))^11), x]

Here is the number of odd coefficients.

Count[%, _?OddQ]

Default and Optional Patterns

Here is a function pow that accepts two arguments, and returns the first argument raised to the power of the second.

pow[a_, b_] := (a^b)

The first argument matches a_ and the second argument matches b_; and Mathematica makes the appropriate substitutions on the right-hand side of the definition.

pow[-6, 3]

One way to assign a default value for the exponent (say, 2) is to give pow a second definition for the cases in which it is called with one argument.

pow[a_] := (a^2)

The newer rule is applied when pow is called with one argument, effectively defining a default value of 2 for the exponent.

pow[256]

The on-line help displays the two definitions for pow.

?pow

A more concise way to do the same thing involves using default patterns. A default pattern looks like patname_:defaultvalue, and defaultvalue is returned as the value of patname if the corresponding argument is omitted. Here is a function powtwo that combines the two previous definitions into one statement.

powtwo[a_, b_: 2] := (a^b)

If we give powtwo two arguments, both are used when computing the function.

powtwo[-6, 3]

If only one argument is given, Mathematica uses the default value 2 for the exponent’s pattern b_.

powtwo[256]

Some operators such as addition and exponentiation have built-in default values. Here is a definition of the power rule for a differentiation function. (For this definition, we ignore the possibility that the exponent depends on the variable of differentiation.)

\(deriv[var_\^n_, var_] := n\ var\^\(n - 1\)\)

The pattern inside the deriv function matches expressions of the form (x^n).

deriv[(z^5), z]

Mathematica matches patterns according to form, and not necessarily mathematical equality. For this reason, Mathematica does not recognize x as (x^1), and the following case returns unevaluated.

deriv[x, x]

One way to define deriv for this case is to give the exponent pattern a default value of 1, using the pattern n_:1.

Clear[deriv]

\(deriv[var_\^\(n_: 1\), var_] := n\ var\^\(n - 1\)\)

The function now handles this case correctly.

deriv[x, x]

To use exponentiation’s default (or identity) value in a pattern, we include a period (.) directly after the underscore. The default value of an exponent is 1, and we tell Mathematica to allow the case (x^1) as follows.

\(deriv2[var_\^n_., var_] := n\ var\^\(n - 1\)\)

The case (x^1) now works as intended.

deriv2[x, x]

Other defaults are x_ + y_., in which y defaults to 0; and x_*y_., in which y defaults to 1.

Define a function ann that computes the net present value of an annuity-immediate using the formula (∑_(k = 1)^n)(((1 + i))^(-k)), where n (the term of the annuity) and i (the interest rate) are arguments to ann, and i has a default value of 5%.

We assign a default value to the argument i by including the pattern i_:0.05.

ann[n_, i_:0.05] := Sum[(1 + i)^-k, {k, 1, n}]

We can include the interest rate as a second argument to ann.

ann[10, 0.04]

If we omit the second argument, the default value is used.

ann[10]

Variable Numbers of Arguments

In some cases we will define functions that should accept any number of arguments. For example, the built-in operator Plus accepts a variable number of arguments.

{Plus[a, b], Plus[a, b, c, d]} // ColumnForm

To handle the possibility of any number of arguments, the pattern __ (two underscores) matches one or more arguments, and ___ (three underscores) matches zero or more arguments.

We can define a function called HasConsecutivePairQ that turns True if it contains the same argument twice in a row. One way to detect such pairs is to use the same pattern name twice in a row on the left side of the function definition, as in the following simplified example.

PairOnlyQ[z_, z_] := True

PairOnlyQ[3, 3]

To allow a matching consecutive pair anywhere among the function’s arguments, we insert triple-underscore patterns before and after the matching-pair patterns, which tells Mathematica to allow any number (including zero) of arbitrary arguments before and after the matching pair. Any sequence of arguments before or after the matching pair is ignored, so we do not name the corresponding patterns.

HasConsecutivePairQ[___, z_, z_, ___] := True

We insert a second rule that matches every other case and returns False.

HasConsecutivePairQ[___] := False

The following use of HasConsecutivePairQ has a pair of consecutive 3s, and returns the value True.

HasConsecutivePairQ[0, 1, 2, 3, 3, 4, 5]

The following has no consecutive matching pairs, so it returns False.

HasConsecutivePairQ[0, 1, 2, 0, 1, 2]

To define a HasPairAnywhereQ function that returns True if any two of its arguments match, we create the following definition (followed by a definition that matches anything else and returns False).

HasPairAnywhereQ[___, z_, ___, z_, ___] := True; HasPairAnywhereQ[___] := False

There can now be any number of elements between the matching pair of elements.

HasPairAnywhereQ[0, 1, 2, 1, 3]

Here there are no elements between the matching items.

HasPairAnywhereQ[0, 1, 1, 2]

To define a function HasPairAtLeastOneSpaceApartQ, we use a double-underscore between the patterns with the same name. The double-underscore matches any number of elements except zero.

HasPairAtLeastOneSpaceApartQ[___, z_, __, z_, ___] := True

HasPairAtLeastOneSpaceApartQ[___] := False

Now the function matches pairs with at least one element between the matching items.

HasPairAtLeastOneSpaceApartQ[0, 1, 2, 3, 1, 4]

However, the function no longer recognizes consecutive pairs.

HasPairAtLeastOneSpaceApartQ[0, 1, 1, 2]

Double- and triple-blank patterns are often used to create error-catching definitions, or “fall-through” cases, the idea being that they catch input that does not match any of the “correct” forms of a function.

Define a function HasTwoConsecutivePairsQ that returns True if given two consecutive matching pairs of arguments (with any number of arguments between the pairs).

Each consecutive pair will match a pattern of the form (x_,  x_); that is, two consecutive patterns with the same name. Because there can be any number of arguments before and after each pair, we surround the matching patterns with triple-blank patterns that match zero or more arguments.

HasTwoConsecutivePairsQ[___, x_, x_, ___, y_, y_, ___] := True

In the following test, there are two consecutive 3s and two consecutive 7s.

HasTwoConsecutivePairsQ[1, 2, 3, 3, 4, 5, 6, 7, 7, 8]

In the following test, there is only one matching consecutive pair, so the arguments do not match the pattern and the function returns unevaluated.

HasTwoConsecutivePairsQ[1, 2, 3, 3, 4, 5, 6, 7, 8]

Define a function called SortList that sorts a list of numbers into descending order, using double- or triple-blank patterns.

A simpler function SortTwo that sorts a list of two elements into descending order might have the following definition: SortTwo of a list {x, y} returns {y, x} provided that y > x.

SortTwo[{x_, y_}] := {y, x} /; y > x

A recursive definition of SortList switches any two consecutive elements of the list for which the first element is greater than the second element. To allow any number of elements before and after the pair under consideration, we use triple-blank patterns (bef___ and aft___) before and after the pair.

Clear[SortList]

SortList[{bef___, x_, y_, aft___}] := 
      SortList[{bef, y, x, aft}] /; y > x

When we give SortList a list of numbers, the definition calls itself until no consecutive pairs of arguments have the property that the second is greater than the first.

SortList[Table[Random[Integer, {1, 10}], {15}]]

To define SortList to return a list of sorted values, we define a second “fall-through” case that applies to any case that does not match the first pattern, and simply returns the argument.

SortList[L_] := L

Now SortList returns a sorted list.

SortList[Table[Random[Integer, {1, 10}], {15}]]

Define a function RemoveDuplicates that removes duplicates from a list of elements.

To recognize duplicates we use the same pattern name twice in a function definition.

AreTheseTheSame[x_, x_] := True

Here is a fall-through case that applies to nonmatching pairs.

AreTheseTheSame[x_, y_] := False

Here is a test of our function.

{AreTheseTheSame[1, 1], AreTheseTheSame[1, 2]}

Here is a recursive definition that recognizes duplicate elements in a list, with any number of elements (including zero) before, between, and after the duplicate items, using triple-blank patterns. (The function is defined so that for any duplicate pair, it discards the element that occurs later in the list.)

RemoveDuplicates[{bef___, x_, betw___, x_, aft___}] := 
      RemoveDuplicates[{bef, x, betw, aft}]

Here is a list that contains some duplicate elements.

listWithDuplicates = Table[Random[Integer, {1, 10}], {15}]

Here is a test of the function.

RemoveDuplicates[listWithDuplicates]

For RemoveDuplicates to return a list, we define a fall-through case that matches any list that does not match the first case we defined, and returns just the list.

RemoveDuplicates[L_] := L

The function now works as intended.

RemoveDuplicates[listWithDuplicates]

Repeated Patterns

If we give the same name to more than one pattern inside a function, the arguments corresponding to the patterns must be the same. For example, the function deriv, defined previously, requires the same symbol for both the base of the expression we wish to differentiate and the variable of differentiation; this is because the pattern for the base and the pattern for the variable of differentiation have the same name var.

\(deriv[var_\^n_, var_] := n\ var\^\(n - 1\)\)

Here is a function called merge that accepts a pair of matching arguments, and returns one copy of the argument.

merge[x_, x_] := x

Here is a case that matches the given patterns.

merge[1, 1]

This case does not match the patterns.

merge[0, 1]

Unlike merge, which was defined only for a sequence of two identical arguments, there are cases in which we do not know beforehand how many arguments will be given, and we require that the arguments be the same. For example, the following function requires a list of pairs as its argument.

NPV[args : {{_, _}..}] := 
      (∑_(i = 1)^(Length[args]))args[[i, 1]] 
          (1.05^(-args[[i, 2]]))

The double period .. in the pattern args : {{_, _}..} tells Mathematica that the pattern {_, _} will be repeated one or more times; the entire structure is called args. NPV returns the net present value of payments whose amounts are defined by the first element in each pair in args, to be paid at times defined by the second argument in each pair (assuming 5% interest compounded each period).

NPV[{{p1, 1}, {p2, 4}, {p3, 10}}]

If we give NPV an argument that is not a list of pairs, it returns unevaluated.

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

Create a pattern that matches a list containing one or more lists of length 3.

The pattern {_, _, _} matches a single list of length 3.

MatchQ[{1, 2, 3}, {_, _, _}]

To match a list containing one or more sublists of length 3, we use the repeated pattern {{_, _, _}..}. Here are two tests of the pattern.

MatchQ[{{1, 2, 3}, {2, 3, 4}, {3, 4, 5}}, {{_, _, _}..}]

MatchQ[{{1, 2, 3, 4}, {2, 3}, {3, 4, 5}}, {{_, _, _}..}]

Here is a simple function definition using the repeated pattern.

ScatterPlot3D[L : {{_, _, _}..}] := 
      Show[Graphics3D[{PointSize[0.03], Map[Point, L]}]]

ScatterPlot3D[{{1, 1, 0}, {1, 0, 1}, {0, 1, 1}, {0, 0, 0}}];

Upvalues

The typical function definitions we have seen have the following form. The left side of the definition has the head mylog, and consequently Mathematica attaches the definition to the symbol mylog.

mylog[x_ y_] := mylog[x] + mylog[y]

In some cases, however, we want to associate a function definition with a symbol that is not the head of the left side of the definition. For example, if we want a product of myexp functions to simplify to myexp of the sum of the arguments, we might enter the following definition.

myexp[x_] myexp[y_] := myexp[x + y]

Mathematica is unable to create this rule because a product of myexp functions has the head Times, and Times is protected against accidental modification, meaning we cannot by default attach new definitions to it. Here we see that the head of the left side of the previous definition is Times.

Head[myexp[x_] myexp[y_]]

One way to teach Mathematica this definition is to define an upvalue for myexp; that is, a definition whose left side has the head Times, but which should be attached to the symbol myexp.

myexp /: myexp[x_] myexp[y_] := myexp[x + y]

Now Mathematica attaches the simplification rule to myexp and not Times, and carries out the desired simplification.

myexp[a b] myexp[c] myexp[d]

The on-line help displays the upvalue for myexp using the abbreviation ^:=.

?myexp

Create a symbol Undefined that has the property that if it is used in a sum or product, the result is the symbol Undefined.

Because we cannot attach new rules to the operators Plus and Times, we must define upvalues for the symbol Undefined. The following two rules state that if Undefined is used in a sum or product, the result is the symbol Undefined.

Undefined /: Undefined + _ := Undefined

Undefined /: Undefined*_ := Undefined

2 Undefined + 1/2

The Ordering of Definitions

The principle Mathematica uses to order multiple definitions for a function is to apply specific rules before applying more general rules, in order to keep specific rules from being shadowed by general rules. For example, here are two definitions for the factorial function fac (after clearing any previous definitions).

Clear[fac]; fac[n_] := n fac[n - 1]; fac[0] = 1;

Asking for information about fac reveals the rules Mathematica uses when computing values of fac, and the order in which the rules are applied. We see in the following listing that the specific definition fac[0] = 1 is tried before the more general definition fac[n_] := n fac[n - 1]; if this were not the case, the definition attached to fac[0] would never be applied.

?fac

The rules attached to a function are called its downvalues, and we can get a list of a function’s downvalues by entering DownValues[FunctionName]. Items found earlier in the list of downvalues are applied before items found later in the list.

DownValues[fac]

In some cases Mathematica cannot determine which rule among a set of rules is more general. In this case, Mathematica evaluates the rules according to the order in which they were entered. Here are two equally general definitions.

ambig[x_Integer, y_] := "rule one"

ambig[x_, y_Integer] := "rule two"

In this case, Mathematica applies the first definition to arguments that match both definitions.

ambig[1, 2]

To change the order in which definitions are applied, we manipulate the list of downvalues. Here is the current ordering.

DownValues[ambig]

Here we reverse the order in which the definitions are applied.

DownValues[ambig] = Reverse[DownValues[ambig]]

Mathematica now applies the second definition.

ambig[1, 2]

Transformation Rules

Functions work by replacing the dummy variables in the right-hand side of a definition with specific values, which are entered as arguments to the function. Another general method for substituting particular values into an expression is to use transformation rules.

One way to substitute particular values for the variables x and y in the formula Sqrt[x^2  +  y^2] is to assign particular values to x and y, then evaluate the formula. Here we set y equal to 7.

y = 7;

Evaluating Sqrt[x^2  +  y^2] now reflects the new value of y.

Sqrt[(x^2) + (y^2)]

One problem with this approach is that the value 7 is substituted wherever y appears in an expression evaluated after the assignment, at least until y is cleared of its value. For instance, if we try to solve an equation with respect to y, Mathematica substitutes the value 7 into y, in effect assuming we want to solve the equation with respect to 7.

Solve[(y^10) + (x^10) == 1, y]

To clear the value of y, we use Clear.

Clear[y]

A better method for replacing variables (or any subpart of an expression) with particular values is to use transformation rules, or simply rules. A rule is an object of the form before -> after (with internal form Rule[before, after]), and we carry out a substitution by using the slash-period /. operator in the form expr /. rules, where expr is the expression to which the rules apply, and rules is a single rule or list of rules. (The internal form of expr /. rules is ReplaceAll[expr, rules].)

For example, to replace y in the expression Sqrt[x^2  +  y^2] with 7, we type the following.

Sqrt[(x^2) + (y^2)] /. y -> 7

Although we replaced y with 7 in the expression, the global value for y remains undefined.

y

We are not limited to numerical replacements. Here we replace y with 1 + a.

Sqrt[(x^2) + (y^2)] /. y -> 1 + a

We can replace any number of variables in an expression with other values by giving a list of replacement rules.

Sqrt[(x^2) + (y^2)] /. {x -> 12, y -> 5}

We can also create a list of expressions corresponding to different sets of replacement rules by telling Mathematica to substitute more than one set of values for an expression. In the following we compute three values for Sqrt[x^2 + y^2] by substituting into it a list containing three lists of replacement rules.

Sqrt[(x^2) + (y^2)]
         /. {{x -> 3, y -> 4}, {x -> 12, y -> 5}, {
          x -> 1, y -> 1}}

A subtle case arises when we perform more than one substitution on an expression with only one variable. For instance, to substitute two values for x into an expression we might try the following.

ex = Log[x (1 - x)];

ex /. {x -> 0.35, x -> 0.55}

While we may have expected two values in return, we see that Mathematica returned only one value. The reason is that Mathematica performs the substitutions from left to right: First the rule x -> 0.35 was applied to ex, replacing all occurrences of x with 0.35; then the rule x -> 0.55 was applied to that result, which was free of any occurrences of x.

To avoid this difficulty, we must put each set of replacement rules inside its own list, even if each set contains only one rule: this indicates to Mathematica that each replacement rule should be considered a separate solution set.

ex /. {{x -> 0.35}, {x -> 0.55}}

In general, Mathematica will return as many values from a substitution as there are sets of replacement rules.

Replacement rules are important because functions such as Solve return results in the form of a list of lists of replacement rules.

Solve[(x^2) - 1 == 0, x]

Here we verify the solutions by substituting each one into the original equation.

(x^2) - 1 == 0 /. x -> -1

(x^2) - 1 == 0 /. x -> 1

Instead of testing the solutions separately, we can test both solutions at the same time by substituting a list of solution sets into the equation.

(x^2) - 1 == 0 /. {{x -> -1}, {x -> 1}}

A more efficient way to do the same thing is to name the list of replacement rules. Here we solve the same equation, this time naming the solution set of replacement rules solset.

solset = Solve[(x^2) - 1 == 0, x]

Instead of retyping the solutions, we can directly substitute solset into the equation (note that the entire rule is contained in solset, so it is important not to type x -> solset).

(x^2) - 1 == 0 /. solset

Similarly, we can generate a list containing only the solutions by substituting the rules containing the solutions into the variable(s) of the equation.

x /. solset

Any part of an expression can be replaced with another symbol or expression.

Sin[x] /. Sin -> Cos

Sin[x] /. {Sin -> Cos, x -> 1 + (y^2)}

Here are rules that perform unit conversion.

24 inch /. inch -> 2.54 cm

10000.0 feet /. feet -> (1/5280)miles

Even the head of an expression can be replaced. Below we replace the head List with the head Plus, which adds the elements of the list.

{a, b, c} /. List -> Plus

In addition to using rules to replace literal symbols, we can use patterns with rules. Here we replace every occurrence of x raised to an even power n with zero.

1 + x + (x^2) + (x^3) + (x^4) /. x^n_?EvenQ -> 0

Here we replace every occurrence of (Sqrt[(z^2)]) in an expression with z.

Sqrt[(((1 + x))^2)] - (y^2) /. Sqrt[z_^2] -> z

By using a pattern on the left side of a rule, we can use functions without making explicit definitions.

funct[10] /. funct[x_] -> 1 + (x^2)

We can combine function definitions and rules. The following function rootPlot plots the roots of an equation in the complex plane, using the ideas developed earlier.

Solve and NSolve return solutions to systems of equations as lists of replacement rules.

NSolve[(x^5) == 1 + I, x]

To isolate the solutions, we substitute the values in the replacement rules for the variable x.

x /. %

Here is one way to convert a list of complex and real numbers into a list of ordered pairs to be plotted with ListPlot, using rules and a conditional pattern that matches any numeric quantity: the rule converts a numeric quantity x into an ordered pair of the form {Re[x], Im[x]}.

% /. x_?NumericQ -> {Re[x], Im[x]}

The pairs can now be plotted using ListPlot.

ListPlot[%, PlotStyle -> PointSize[0.03]];

We combine the previous steps into one function called rootPlot.

\(\(rootPlot[eqn_Equal, var_] := 
	
      ListPlot[\(var /. NSolve[eqn, var]\) /. 
          x_?NumericQ -> {Re[x], Im[x]}, 
        PlotStyle -> PointSize[0.03], PlotRange -> All, 
        AspectRatio -> Automatic]; \)\)

Here is a plot of the seven roots of \(TraditionalForm
      \`\(w\^7\)i = 
        w\^5 - \(w\^3\)i + w\^2 - w + 1\).

rootPlot[I (w^7) == (w^5) - I (w^3) + (w^2) - w + 1, w];

Variable Scoping

Many function definitions can be simplified by the use of auxiliary functions and variables. The following function returns a list containing the minimum and maximum x  values from a list of {x, y} pairs.

findminmax[data : {{_, _}..}] := {Min[First[Transpose[data]]], 
        Max[First[Transpose[data]]]}

findminmax[{{2, 5}, {7, 1}, {1, 3}}]

Changing the right side of the definition into a compound expression in which we assign the auxiliary variable tdata the value First[Transpose[data]] will simplify the function somewhat.

Clear[findminmax]

findminmax[data : {{_, _}..}] := 
      (tdata = First[Transpose[data]]; {Min[tdata], Max[tdata]})

findminmax[{{2, 5}, {7, 1}, {1, 3}}]

A problem with this approach is that the value of tdata remains in memory until we explicitly clear it.

tdata

Clear[tdata]

A better approach is to use temporary, or local, variables on the right-hand side of a definition, which is done by enclosing the body of the definition in a Module, and listing the variables that are local to the routine. In the following definition of findminmax, we put the right side of the definition inside Module, declaring tdata to be a local variable.

Clear[findminmax]

\(findminmax[data : {{_, _}..}] := 
	
      Module[{tdata}, 
		
        tdata = First[Transpose[data]]; {Min[tdata], Max[tdata]}]\)

The function works as before.

findminmax[{{2, 5}, {7, 1}, {1, 3}}]

However, tdata has not been assigned a permanent value.

tdata

The following function accepts a list of ordered pairs that lie in the unit square, then creates a graph of the ordered pairs, along with a red rectangle that highlights the middle third of the vertical range of the data. In order to compute the corners of the rectangle, it is necessary to use several local variables.

Clear[bandplot]

\(bandplot[data : {{_, _}..}] := 
      Module[{tdata, ymin, ymax, ymed, yrange}, tdata = Transpose[data]; 

        		{ymin, ymax} = {
            Min[tdata[[2]]], 
            Max[tdata[[2]]]}; 
        ymed = \(1\/2\) \((ymin + ymax)\); yrange = ymax - ymin; 
		
        ListPlot[data, 
          Prolog -> {Hue[0], 
              Rectangle[{0, ymed - yrange/6}, {1, ymed + yrange/6}], 
              GrayLevel[0]}, PlotRange -> All, 
          PlotStyle -> PointSize[0.02], AxesOrigin -> {0, 0}]]\)

bandplot[Table[{Random[], Random[]}, {30}]];

Module performs lexical scoping; that is, it renames local variables so that they do not conflict with global variable names. Other scoping functions are Block, which performs dynamic scoping, and With, which creates local constants.

Side Effects and Dynamic Programming

In addition to returning a value, many functions have side effects such as printing intermediate results or showing a graph. The result of a Plot command is a collection of the graphics primitives and directives necessary to reproduce the graph (the collection has the abbreviated output form “⁃Graphics⁃”); the graph itself is merely a side effect of the computation.

Plot[1/Abs[(t^2) - 3 t - π], {t, -2, 5}]

Other possible side effects are Print statements, and functions that change values of global symbols. The following function fib computes Fibonacci numbers recursively, using the formula (Subscript[F, n]) = (Subscript[F, n - 1]) + (Subscript[F, n - 2]), with the starting values (Subscript[F, 1]) = (Subscript[F, 2]) = 1. First, we define the starting values.

Clear[fib]; fib[1] = fib[2] = 1;

The right-hand side of the following definition stores, or caches, intermediate values of fib by immediately setting each intermediate value fib[n] equal to its value.

fib[n_] := (fib[n] = fib[n - 1] + fib[n - 2])

Here is the seventh Fibonacci number.

fib[7]

We see that the definition of fib has the side effect of storing intermediate results.

? fib

The style of programming that involves a function modifying its own definition is called dynamic programming.

Lists

Creating and Interpreting Lists

The basic data structure in Mathematica is the list. The elements of a list can be any combination of numbers, variables, strings, functions, and other lists; enclosed within curly brackets { } and separated by commas.

Mathematica lists often have natural interpretations as vectors, matrices, or tensors, depending on how deeply they are nested. A one-dimensional list acts as a vector, a two-dimensional list as a matrix, and higher-dimensional lists as tensors.

One way to enter an arbitrary list is to type a sequence of elements inside curly brackets and separated by commas. The following one-dimensional list vec can be considered a vector with the three components a, b, and c.

vec = {a, b, c};

An easier way to create a vector of evenly spaced numbers is to use Range. Range[max] creates a list of integers from 1 through max; Range[min, max] creates a list of numbers from min through max in increments of 1; and Range[min, max, inc] returns a list of numbers ranging from min through max with increment inc.

vec2 = Range[2, 5, 1/2]

A more general function for creating lists is Table, which generates lists (nested to any depth) of elements that may depend on the given iterators. For example, to generate a list of values of (k^2) as k takes on values from 1 to 10, we evaluate the following.

Table[(k^2), {k, 1, 10}]

We create tables of more complicated functions in a similar way.

vec3 = 
      Table[N[Exp[π I t]], {t, 0, 2 π, π/4}]

Many mathematical and plotting functions are designed to operate on lists. Here is a plot of the real parts of the elements of vec3.

ListPlot[Re[vec3], PlotStyle -> PointSize[0.03]];

Fourier expects a list of values.

Fourier[vec3]

We generate two-dimensional tables by including two iterators.

rect = Table[(i/j), {i, 1, 5}, {j, 1, 2}]

It may be useful to use TableForm to display lists of two or more dimensions.

TableForm[rect]

The options TableAlignments, TableSpacing, and TableHeadings allow us to customize the appearance of a table.

\(TableForm[rect, 
	TableAlignments -> {Center, Center}, 
	
      TableSpacing -> {1, 10}, 
      TableHeadings -> {None, {"Before", "After"}}]\)

How we interpret an array such as rect depends on the context in which it arises. In some cases, it may be considered a 5 × 2 matrix, in others a list of five ordered pairs. Here we interpret rect as a matrix, and matrix-multiply it by its transpose using Dot (.).

rect.Transpose[rect] // MatrixForm

Here we interpret rect as a list of ordered pairs to be plotted.

ListPlot[rect, PlotStyle -> PointSize[0.02]];

The first argument to Table can be a conditional statement such as If or Which. An element of the following matrix is 1 if its indices i and j satisfy \(TraditionalForm
      \`|i - j| ≤ 2\), and is otherwise 0.

TableForm[Table[If[Abs[i - j] ≤ 2, 1, 0], {i, 7}, {j, 7}]]

There is no need for Mathematica arrays to be rectangular. Here the upper limit of the inner loop j depends on the value of the outer loop i.

trimat = Table[i + j, {i, 1, 4}, {j, 1, i}]

The result is a triangular array.

TableForm[trimat]

There is no limit to the number of dimensions that can be contained in a Mathematica list. Here we generate a four-dimensional tensor (alternatively, a list of lists of lists of lists, or a matrix of matrices) by including four iterators in the Table command.

tensor = Table[i + j + k - l, {i, 2}, {j, 3}, {k, 4}, {l, 5}];

TableForm[tensor]

It is important to make a distinction between list operations and matrix operations, particularly with respect to matrix multiplication. The function Dot (which we enter as a period .) is used for matrix and vector multiplication. Here is the product of a matrix and a vector.

{{a, b}, {c, d}}.{x, y}

There is no difference in Mathematica between a row vector and a column vector, as matrix and vector functions automatically interpret vectors in the appropriate way. Here we multiply a matrix by a vector on both the left- and right-hand sides. Notice that the vector is typed the same way in each case.

{x, y}.{{a, b}, {c, d}}.{x, y} // Simplify

Here is the matrix product of two matrices.

{{1, 2}, {3, 4}}.{{a, b}, {c, d}} // MatrixForm

Here we use ordinary multiplication *, and Mathematica performs termwise multiplication.

{{1, 2}, {3, 4}}*{{a, b}, {c, d}} // MatrixForm

As a general rule, if two lists have the same shape, Mathematica will perform termwise arithmetic, and return a list with the same shape. For instance, the difference of two vectors of equal length returns a vector of the differences of the elements.

{1, 2, 3} - {a, b, c}

Notice also that Mathematica will combine all elements of a list with a single expression. If we add a constant to a vector, Mathematica returns a vector where the constant has been added to each element.

{a, b, c} + 2

Similarly, multiplying or adding a single expression (a scalar) to a matrix results in the expression being multiplied or added to each element of the matrix.

{{1, 2}, {3, 4}} 𝒞 //  MatrixForm

{{1, 2}, {3, 4}}  + 𝓀 //  MatrixForm

Mathematica produces a warning if two lists cannot be combined because they have different lengths.

{1, 2, 3} + {a, b, c, d}

Exercises: Creating and Interpreting Lists

Create a list of the first ten prime numbers. (Prime[n] returns the (n^th) prime number.)

We use Table to generate a list of the (k^th) prime numbers as k goes from 1 to 10. When we omit a starting value for the iterator k, it is assumed to be 1.

Table[Prime[k], {k, 10}]

Alternatively, we can avoid using the Table command by recognizing that the mathematical function Prime automatically maps itself over lists.

Prime[Range[10]]

Define a function called zeroVector that creates a vector containing n zeros, where n is an argument to zeroVector.

Here is one definition using Table. We can omit the name of the iterator if it does not occur in the expression given to Table.

zeroVector[n_Integer?Positive] := Table[0, {n}]

Here is a vector of ten zeros.

zeroVector[10]

We can define a zeroMatrix function in a similar way.

zeroMatrix[m_, n_] := Table[0, {m}, {n}]

Here is a 3 × 5 matrix of zeros.

zeroMatrix[3, 5] // MatrixForm

Generate a 5 × 5 matrix, in which each element is the row number plus the column number, or the number 5, whichever is greater.

We use Max to determine if 5 or i + j (the row number plus the column number) is greater.

Table[Max[{5, i + j}], {i, 5}, {j, 5}] // MatrixForm

Beginning with the matrix {{a, b}, {c, d}} and the vector {x, y}, create the vector {a x + b y, c x + d y}. Beginning with the same matrix and vector, create the matrix {{a x, b x}, {c y, d y}}.

This exercise illustrates the difference between list and matrix operations. Here are definitions of the matrix and vector.

exermat = {{a, b}, {c, d}}; exervec = {x, y};

For matrix-vector multiplication we use Dot.

exermat.exervec

To perform termwise multiplication we use Times.

exermat*exervec // MatrixForm

Compute the eigenvectors of a numeric approximation to the matrix

and see if they are normalized (that is, see if each eigenvector has unit length).

Here are the eigenvectors of the approximate matrix.

evecs=Eigenvectors[N[(({{1, 2, 3}, {3, 2, 2}, {7, 1, 2}}))]];TableForm[evecs]

The following function determines the length of a vector.

veclength[z_] := (Sqrt[z.z])

Here we use Table to test if each eigenvector has unit length.

Table[veclength[evecs[[i]]], {i, 3}]

Alternatively, we use the high-level programming function Map to map veclength over each vector in evecs.

Map[veclength, evecs]

Extracting Parts of Lists

Given a one-dimensional list, we extract elements according to their position by using Part. In Mathematica the first element of a list is at position 1, not 0. (Part 0 of a list is the head List.)

Part[{a, b, c}, 1]

A compact, equivalent way to express the same command is to enclosed the desired position between double brackets [[ and ]].

{a, b, c}[[1]]

In a multidimensional array, an element of a list can be another list. For example, here is a table of hurricane names, wind speeds (in knots), and pressures (in millibars),where each triple is collected in a sublist.

hurr =
	{{"Allison",   65,  987},
	 {"Barry",     50,  995},
	 {"Chantal",   60,  992},
	 {"Dean",      40,  999},
	 {"Erin",      75,  975},
	 {"Felix",    115,  950},
	 {"Gabrielle", 60,  998},
	 {"Humberto",  90,  970},
	 {"Iris",      75,  970},
	 {"Jerry",     35, 1003}};

TableForm[hurr, TableSpacing -> {0, 3}]

To extract an entire sublist, we specify its position using Part (or the part-specification brackets [[ ]]).

Part[hurr, 2]

To extract a single element of a two-dimensional list, we refer to its coordinate, its row and column. For instance, to extract the wind speed of hurricane Erin, we first select Erin’s whole record, which is at position 5, and then select the second part of the record.

Part[hurr, 5, 2]

Here are two more ways to enter the same expression.

{hurr[[5, 2]], hurr[[5]][[2]]}

Similarly, we extract the pressure of the storm Gabrielle by extracting the sublist at position 7, then extracting the third part of Gabrielle’s record.

hurr[[7, 3]]

To extract elements from deeply nested lists, it may be necessary to specify many coordinates.

In general, specifying one position in an n-dimensional (“rectangular”) list results in a list with n - 1 dimensions. Extracting the first element of the 4-dimensional list called tensor results in a 3-dimensional list.

tensor[[1]] // TableForm

Similarly, extracting an expression specified by two coordinates from an n-dimensional list returns a list with n - 2 dimensions. Part (1, 1) of tensor is a two-dimensional matrix.

tensor[[1,1]] // MatrixForm

Specifying a part of tensor with three coordinates returns a vector.

tensor[[1,1,1]]

To extract an atomic element from tensor, we specify four coordinates.

tensor[[1,1,1,1]]

It is important to be aware that there is a difference between the part specification hurr[[1,2]] and hurr[[{1,2}]] (the latter has curly brackets inside the double square brackets). The former returns the element at position (1, 2) of the list hurr.

hurr[[1,2]]

The latter is equivalent to {hurr[[1]], hurr[[2]]}.

hurr[[{1,2}]]

Note that [[{ }]] notation provides a quick way to manipulate the elements of a list. Here is an example of using this notation to rearrange the elements of a list.

plist = {a, b, c, d, e, f};

plist[[{1, 4, 2, 5, 3, 6}]]

Using Table and Part, we can define a function called sampleWithReplacement that takes a sample of a given size from a list of data.

sampleWithReplacement[L_List, n_Integer?Positive] := 
      Module[{ℓ = Length[L]}, 
        Table[L[[Random[Integer, {1, ℓ}]]], {n}]]

sampleWithReplacement[{a, b, c, d, e}, 10]

We can extract a range of elements using Take. Here we extract the first five elements of a list.

Take[{a, b, c, d, e, f, g}, 5]

Here we extract the last five elements.

Take[{a, b, c, d, e, f, g}, -5]

We can extract a sublist from the middle of the list, as well. Here we extract the third through sixth elements of a list.

Take[{a, b, c, d, e, f, g}, {3, 6}]

Drop works similarly, removing elements from a list. Here we drop the first three elements from a list.

Drop[{a, b, c, d, e, f, g}, 3]

Note that Take and Drop do not alter the value of a list. Here is a vector of values.

alpha
       = {𝒶, 𝒷, 𝒸, 𝒹, ℯ, 
        𝒻};

If we drop some elements from alpha, Mathematica returns the remaining elements.

Drop[alpha, {3, 5}]

However, the original list remains unchanged.

alpha

To alter the value of alpha, we must set the value of alpha to the desired subset.

alpha = Drop[alpha, {3, 5}];

alpha

Here is a function called sampleWithoutReplacement that takes a given number of samples from a list, without choosing the same element more than once.

sampleWithoutReplacement[L_List, n_Integer?Positive] := 
      Module[{ℒ = L, ℓ = Length[L], pos}, 
        Table[pos = Random[Integer, {1, ℓ}]; 
          ℒ = Drop[ℒ, {pos}]; 
          ℓ = ℓ - 1; L[[pos]], {n}]]

Here is a sample of fifteen data points from a list of integers from 1 to 50.

sampleWithoutReplacement[Range[50], 15]

Exercises: Extracting Parts of Lists

Using Part, extract the b from {{a, b, c}, {d, e, f}}.

The b is in the first row and second column of the matrix, so we use Part as follows.

Part[{{a, b, c}, {d, e, f}}, 1, 2]

We can use the double-bracket shorthand for Part as follows.

{{a, b, c}, {d, e, f}}[[1, 2]]

Switch the third and fourth elements of the list {a, b, c, d, e}.

We can use the part-selection notation [[{ }]] to switch the elements.

{a, b, c, d, e}[[{1, 2, 4, 3, 5}]]

Use Table and Part to define a function called EveryOtherElement that extracts every other element of a list.

The key is to set up an iterator inside Table that extracts every second element from a list.

EveryOtherElement[L_List] := Table[Part[L, n], {n, 1, Length[L], 2}]

EveryOtherElement[{a, b, c, d, e, f, g, h, i, j, k, l}]

Define a function MatrixTrace that computes the sum of each of the elements along the diagonal of a matrix.

Here is a matrix.

testmat=(({{a, b, c}, {d, e, f}, {g, h, i}}));

The elements along the diagonal have equal coordinates in the horizontal and vertical directions.

Table[testmat[[n, n]], {n, 1, 3}]

To add the elements we use Sum.

Sum[testmat[[n, n]], {n, 1, 3}]

Here is a general definition of MatrixTrace.

MatrixTrace[m_] := Sum[m[[n, n]], {n, Length[m]}]

Here is a square matrix.

\(primemat = Partition[Prime[Range[16]], 4]; 
MatrixForm[primemat]\)

The function works as expected.

MatrixTrace[primemat]

A hand calculation shows that the answer is correct.

2 + 13 + 31 + 53

Here is a modification to allow only square matrices.

Clear[MatrixTrace]

MatrixTrace[m_] := 
      Sum[m[[n, n]], {n, Length[m]}] /; 
        MatrixQ[m] && First[Dimensions[m]] == Last[Dimensions[m]]

Characteristics of Elements

Often it is necessary to select elements from a list according to certain criteria, rather than according to position. The Select function returns the elements of a given list that meet a specified criterion.

For example, one way to extract all of the prime numbers less than 100 is first to generate a list of integers from 1 to 100, then select from that list all the elements that satisfy the predicate function PrimeQ.

Select[Range[100], PrimeQ]

Similarly, we can choose all pairs that satisfy some specified criteria from a list of ordered pairs. Here is a list of points along a spiral.

pts = Table[{r Cos[r], r Sin[r]}, {r, 0, 25, 0.25}];

Here is a predicate function that acts on pairs of numbers and returns True if its first element times its second element is greater than zero.

quad1or3[{x_, y_}] := (x y > 0)

Here is a second predicate that returns True if its second element is greater than or equal to its first element.

ygex[{x_, y_}] := (y >= x)

We want each of the following graphs to use the same options, and SetOptions allows us to change a function’s default option settings. (We will change the option settings back to their default values after this example.)

\(\(SetOptions[ListPlot, AspectRatio -> 1, 
	
      PlotRange -> {{\(-25\), 25}, {\(-25\), 25}}, 
      PlotStyle -> PointSize[0.02], DisplayFunction -> Identity, 

      	Ticks -> None]; \)\)

The first plot is of all the points in the list pts.

ptsplot = ListPlot[pts, PlotLabel -> "All points"];

The second plot is of the points in the first or third quadrant (those that satisfy the predicate function quad1or3).

q1or3pts = Select[pts, quad1or3];

q1or3plot = ListPlot[q1or3pts, PlotLabel -> "Quad. 1 or 3"];

The last plot is of the points for which the y value is greater than or equal to the x value.

ygexpts = Select[pts, ygex];

ygexplot = 
      ListPlot[ygexpts, PlotLabel -> "y ≥ x"];

Here are the three plots.

Show[GraphicsArray[{ptsplot, q1or3plot, ygexplot}]];

Here we set the options to ListPlot back to the default settings.

\(\(SetOptions[ListPlot, 
	AspectRatio -> 1/GoldenRatio, 
	
      PlotRange -> Automatic, 
	PlotStyle -> Automatic, 
      DisplayFunction :> $DisplayFunction, 
      Ticks -> Automatic]; \)\)

We can test whether an element is present in a list by using MemberQ and FreeQ.

MemberQ[{a, b, c, d, e}, c]

FreeQ[{a, (a^b), (a^((b^c)))}, c]

Exercises: Characteristics of Elements

Use Select to extract all the numbers from {a, 1, b, 2, c, 3}.

We use the built-in predicate function NumberQ.

Select[{a, 1, b, 2, c, 3}, NumberQ]

Create the following list of x, y pairs.

pts = 
      Flatten[Table[{x, y}, {x, -3, 3, 0.3}, {y, -3, 3, 0.3}], 1];

Select from it all points whose coordinates satisfy sin(x) < sin(y).

Here is a plot of all the points.

ListPlot[pts, AspectRatio -> 1, PlotStyle -> PointSize[0.02]];

Here is a predicate function that accepts an ordered pair and returns True if the sine of the horizontal coordinate is less than the sine of the vertical coordinate.

inRegionsQ[{x_, y_}] := Sin[x] < Sin[y]

We use Select to extract the points that satisfy inRegionsQ. The following is a shortened version of the output.

Short[Select[pts, inRegionsQ], 3]

Here is a plot of the points that satisfy inRegionsQ.

ListPlot[Select[pts, inRegionsQ], AspectRatio -> 1, 
      PlotStyle -> PointSize[0.02]];

Clear[pts, inRegionsQ]

Using Names and StringLength, find all Mathematica function names that contain more than twenty-five letters. Find all function names that contain fewer than three letters.

Here is a shortened list of all built-in Mathematica function names.

Short[Names["System`*"], 3]

Here is a predicate function that returns True if its argument is a function whose name contains more than twenty-five letters.

longNameQ[s_String] := StringLength[s] > 25

We use Select to extract the desired function names.

Select[Names["System`*"], longNameQ]

The following function returns True if its argument contains fewer than three letters.

shortNameQ[s_String] := StringLength[s] < 3

Here are all the names with one or two letters.

Select[Names["System`*"], shortNameQ]

List Structure

Characteristics of Structure

There are many functions to determine characteristics of a list, such as its length, depth, or dimensions; further, we can test if a list is a vector or matrix.

The length of a vector is the number of elements it contains.

vec = {a, b, c};

Length[vec]

Because the vector vec is one-dimensional, Dimensions[vec] returns a list containing the size of vec in its one dimension.

Dimensions[vec]

The length of a matrix is the number of rows it contains; that is, the number of sublists it contains.

Length[rect]

The dimensions of a matrix are its number of rows and columns.

Dimensions[rect]

In general, the length of an n-dimensional list is the number of sublists it contains that have n - 1 dimensions. The list tensor contains two three-dimensional sublists, so its length is 2.

Length[tensor]

Dimensions returns the length in each “direction” of a nested list.

Dimensions[tensor]

If we have an irregularly shaped array of data, we can find the length of each sublist by mapping the function Length over the array.

irreg = {{1, 2, 3, 4}, {2, 4, 6}, {3, 6, 9, 12, 15, 18}, {4, 8}}; 
    TableForm[irreg]

Map[Length, irreg]

Here is a function that pads (with zeros) an irregularly shaped array, so that each sublist has the same length.

padArray[𝒜_List] := 
      Module[{𝓂𝒶𝓍 = 
            Max[Map[Length, 𝒜]]}, 
        Map[Take[Join[#, Table[0, {𝓂𝒶𝓍}]], 
              𝓂𝒶𝓍]&, 𝒜]]

padArray[irreg] // TableForm

Changing List Structure

Flatten is used to remove all but the outermost list brackets in a nested list of lists.

Flatten[hurr]

Here is a nested list called prods.

prods = Table[m n, {m, 2, 5}, {n, 1, 5}]

Flatten removes all of the inner list brackets from prods, returning a one-dimensional list of values.

Flatten[prods]

The flattened list can now be sorted, have its duplicates removed, and so forth.

Sort[Flatten[prods]]

Union[Flatten[prods]]

Lists can also be partitioned into sublists with a given length. Here a list of integers from 1 through 12 is partitioned into sublists of length 2, 3, and 4.

Partition[Range[12], 2]

Partition[Range[12], 3]

Partition[Range[12], 4]

Notice that all sublists created by Partition must be of the same length: if we partition a length-12 list into sublists of length 5, the two final elements are dropped.

Partition[Range[12], 5]

It is also possible to partition a list into “overlapping” sublists by giving the offset between subsequent lists as a third argument to Partition.

Partition[{a, b, c, d, e}, 3, 1]

Partition[{a, b, c, d, e}, 3, 2]

Often we have lists of data to be manipulated. Manipulating vectors (one-dimensional lists) is straightforward, but to manipulate multidimensional lists may require manipulation of the list’s structure.

For example, say we have a list of ordered pairs called mydata.

mydata = Table[{n, Sin[n + 0.1 Random[]]}, {n, 0, 9, 0.3}]

Transposing mydata results in a list containing two lists: the first being a list of x-coordinates and the second a list of y-coordinates.

Transpose[mydata]

With mydata in its transposed form, we can select the first sublist (the list of all x-coordinates) and manipulate it separately, and so forth. For example, here are the original points.

ListPlot[mydata, PlotJoined -> True];

Here we give the name tdata to the transpose of mydata.

tdata = Transpose[mydata]

Extracting the first part of the transpose of mydata returns a list containing all the x-values in mydata.

tdata[[1]]

Once the x-values have been extracted, we can manipulate them. Here we take the natural logarithm of one plus each of the elements.

Log[1 + %]

To get the transformed data back into the form of ordered pairs, we create a list containing a list of the transformed x-values, then a list of the original y-values.

transdata = {%, tdata[[2]]}

Transposing the previous list returns a list of ordered pairs.

final = Transpose[transdata]

Here is a graph of the transformed data.

ListPlot[final, PlotJoined -> True];

Note that one effect of Transpose is to reverse the dimensions of a matrix.

Dimensions[mydata]

Dimensions[Transpose[mydata]]

Here are the dimensions of our array of hurricane data.

Dimensions[hurr]

Transposing the array returns a list containing a list of all the hurricane names, a list of all the wind speeds, and a list of all the pressures.

transhurr = Transpose[hurr]

Given the transformed list, we can obtain ordered pairs of wind speeds and pressures by first extracting the respective lists.

transhurr[[{2, 3}]]

Next we transpose the previous list to obtain the desired pairs.

hurrpairs = Transpose[ transhurr[[{2, 3}]] ]

The data are now in a form to be plotted using ListPlot.

ListPlot[ hurrpairs,

	FrameLabel->{"Wind Speed", "Pressure"},
	PlotLabel->
		StyleForm["1995 Tropical Storm Data",
			FontFamily->"Helvetica", FontSize-> 10],
	PlotStyle->PointSize[0.02],
	Frame->True, RotateLabel->True,
	TextStyle->{FontSize->9} ];

Exercises: List Structure

Define a function percentageOfPrimes that returns the percentage of prime numbers in a list.

To determine the percentage of primes in a list, we divide the number of primes by the length of the list. One way to determine the number of primes in a list is to use Length and Select.

Length[Select[Range[100], PrimeQ]]

A more straightforward method is to use Count and the conditional pattern _?PrimeQ.

Count[Range[100], _?PrimeQ]

Here is a definition of percentageOfPrimes.

percentageOfPrimes[𝒫_List] := 
      Length[Select[𝒫, PrimeQ]]/Length[𝒫]

Here is the percentage of primes in a list of integers from 1 to 100.

percentageOfPrimes[Range[100]]

Here is the percentage of primes in a list of odd numbers from 1 to 100.

percentageOfPrimes[Range[1, 100, 2]]

Generate the following list.

pts = Table[{r Cos[r], r Sin[r]}, {r, 0.0, 25.0, 0.25}];

Find the maximum x-value. Find the minimum y-value.

To extract the list of x-values, we take the first part of the transpose of pts. Once the list of x-values has been extracted, we use Max to determine the maximum.

Max[ Transpose[pts][[1]] ]

We follow a similar process to determine the minimum of the y-values.

Min[ Transpose[pts][[2]] ]

Sorting and Grouping Data

Sort sorts a list of numbers into increasing order.

Sort[{1, 3, 5, 7, 9, 2, 4, 6, 8}]

We can specify a function that determines if two adjacent elements are sorted. Here we tell Mathematica to sort a list of numbers in decreasing order, by specifying that two elements are sorted if they satisfy the function Greater.

Sort[{1, 3, 5, 7, 9, 2, 4, 6, 8}, Greater]

When working with multidimensional data, it is often necessary to explicitly define a sorting function. Here is our list of hurricane data.

hurr

To sort the list in ascending order of wind speed, we first define a sorting function as follows. The function returns True if the second element (the wind speed) of one record is greater than the second element of another.

hurrSort[rec1_, rec2_] := rec1[[2]] < rec2[[2]]

Sort[hurr, hurrSort] // TableForm

Here is a pure sorting function equivalent to hurrSort.

Sort[hurr, #1[[2]] < #2[[2]]&] // TableForm

Split partitions a list into sublists that have identical elements. Here is a list of results from tossing a coin twenty times.

tosses = Table[If[Random[] < 0.5, "H", "T"], {20}]

Split divides the list into lists of runs.

runs = Split[tosses]

By mapping the function Length over the results of Split, we can compute the length of each run of heads or tails.

Map[Length, runs]

As an extension of the previous example, we generate a list containing pairs of the form (result, length of run).

Map[{First[#], Length[#]}&, runs]

This is called run length encoding.

Exercises: Sorting and Grouping Data

Using RotateRight, compute the first differences of a list of data.

Here is a list of data generated from a quadratic curve.

quaddata = Table[(n^2) - 5 n + 1, {n, 0, 10}]

We compute the first difference by rotating each element of the list of data one position to the right, and subtracting the rotated list from the original list. Here is the list of first differences.

quaddata - RotateRight[quaddata]

The first element of the subtracted list is meaningless (the first element of the original list minus the last element of the original list), so we use Drop to remove it.

Drop[%, 1]

Extract the five largest elements from a list of random numbers.

Here is a list of random numbers.

randomdata = Table[Random[], {15}]

Sort arranges a list into ascending order.

Sort[randomdata]

To extract the five largest elements in the original list, we extract the last five elements from the sorted list.

Take[%, -5]

Define a function called runEncode that takes a list of discrete random data and converts it into lists of the form (element, length of run). Define a function called runDecode that reconstructs a list resulting from runEncode.

Here is some coin-tossing data taken from the discussion.

tosses = Table[If[Random[] < 0.5, "H", "T"], {20}]

runEncode[ℒ_List] := 
      Map[{First[#], Length[#]}&, Split[ℒ]]

td = runEncode[tosses]

This definition of runDecode reverses the encoding of runEncode.

runDecode[ℒ_List] := 
      Flatten[Map[Table[First[#], {Last[#]}]&, ℒ]]

Here are the original data.

runDecode[td]

% == tosses

Expressions

Lists are just a special case of Mathematica expressions, and it is the uniformity among Mathematica expressions that gives us great freedom in manipulating them.

Recall that the most general Mathematica expression has the following form.

(Subscript[e, 0])[(Subscript[e, 1]), (Subscript[e, 2]), …, (Subscript[e, n])]

Recall that the first element (Subscript[e, 0]) is what we have called the head of the expression, and can be determined by evaluating Head[expression]. In the list {a, b, c}, the value of (Subscript[e, 0]) is List, the value of (Subscript[e, 1]) is a, the value of (Subscript[e, 2]) is b, and the value of (Subscript[e, 3]) is c. FullForm displays the internal form of an expression.

FullForm[{a, b, c}]

Just as we extract sublists from lists, we extract subexpressions from expressions.

Compare the FullForm of {a, b, c} and a + b + c.

FullForm[a + b + c]

We use Part to extract an element from a list.

Part[{a, b, c}, 1]

In the same way, we use Part to extract an element from a sum.

Part[a + b + c, 1]

Recall that part 0 of an expression is its head.

Part[a + b + c, 0]

Similarly, here is part of a more complicated sum.

Part[1 + x + x^2 + x^3, 3]

To extract the exponent in the third term of the sum, we refer to its coordinate. The internal representation of exponentiation is Power.

FullForm[x^2]

The exponent is the second part of the expression (x^2), so we can extract the 2 as follows.

Part[1 + x + x^2 + x^3, 3, 2]

Alternatively, we can use [[ ]] notation.

(1 + x + x^2 + x^3)[[3, 2]]

We can use replacement rules with expressions. In this case we convert a list of rules of the form OptionName -> Value to a list of option names.

Options[TableForm]

% /. Rule[name_, value_] -> name

Complicated expressions have complicated internal forms, and we can always examine the internal form of an expression by using FullForm. In some cases we must tell Mathematica not to evaluate an expression by wrapping it inside Hold.

FullForm[Hold[f[x_] := Sin[1/x] /; x ≠ 0 && -3 < x < 3]]

Exercises: Expressions

How many terms are in the expansion of (((x + y + z))^10)?

Length gives the number of elements in any expression.

Length[Expand[(((x + y + z))^10)]]

Extract the right-hand side of the rule x -> 1 + a.

Here is the internal form of the rule.

FullForm[x -> 1 + a]

This internal form is similar to a list containing the left- and right-hand sides of the rule.

FullForm[{x, 1 + a}]

We extract the second element of a list using Part.

Part[{x, 1 + a}, 2]

In the same way we extract the right-hand side of the rule by using Part.

Part[x -> 1 + a, 2]

Strings

In many cases it is convenient to keep data in the form of a string, and then manipulate the string directly. Mathematica contains many string manipulation functions that work similarly to list manipulation commands.

Reverse[{b, a, c, k, w, a, r, d, s}]

StringReverse["backwards"]

Join[{a, b, c}, {x, y, z}]

StringJoin["abc", "xyz"]

Mathematica does not simplify an expression inside quotation marks.

"1 + 1"

ToExpression changes a string into a Mathematica expression, which is then evaluated.

ToExpression[%]

Here is a data set in the form of a list of strings.

register = {
"123456789012345678901234",
"Alabama 2785103100413816",
"Alaska  0306017700110483",
"Arizona 2039083900782878",
"Arkansas0940051700141457"};

The first eight columns of each line contain the name of a state, the ninth through twelfth columns contain the number of automobiles registered (in thousands), the thirteen through sixteenth columns contain the number of thousands of trucks and buses registered, the seventeenth through twentieth columns contain the number of thousands of motorcycles registered, and the remaining columns contain the total for each state. (The first line of data is used to number the columns.)

Here we drop the first element of the list of data.

register = Rest[register]

The following function dissects each string into its respective list of data.

dissect[𝒮_String] := {StringTake[𝒮, 8], 
        StringTake[𝒮, {9, 12}], 
        StringTake[𝒮, {13, 16}], 
        StringTake[𝒮, {17, 20}], 
        StringTake[𝒮, {21, 24}]}

Map[dissect, register] // TableForm

As defined above, dissect leaves each number in the data set in the form of a string. The following addition to the definition of dissect changes the numerical strings into numbers.

dissect[𝒮_String] := 
      Flatten[{StringTake[𝒮, 8], 
          ToExpression[{StringTake[𝒮, {9, 12}], 
              StringTake[𝒮, {13, 16}], 
              StringTake[𝒮, {17, 20}], 
              StringTake[𝒮, {21, 24}]}]}]

Map[dissect, register] // TableForm

Just as parts of lists and expressions can be replaced using ReplaceAll (/.), parts of strings can be replaced using StringReplace. In the following string we replace every occurrence of the word “Cos” with “cos” and replace every square bracket with a matching parenthesis.

StringReplace[ "Cos[Cos[x]]",
	{"Cos" -> "cos", "[" -> "(", "]" -> ")"} ]

Exercises: Strings

How many characters are in the word Mathematica?

StringLength returns the number of characters in a string.

StringLength["Mathematica"]

Use Characters and Permutations to find all anagrams (rearrangements of letters) of the word math.

Here are the characters in the string “math”.

Characters["math"]

Permutations returns all possible permutations of the elements in a list.

anagrams = Permutations[Characters["math"]]

There are 24 (4!) permutations.

Length[%]

To combine each list of characters into a string, we use StringJoin.

Table[StringJoin[anagrams[[i]]], {i, 24}]

Alternatively, we can use the high-level programming function Map to map StringJoin over each list of characters.

Map[StringJoin, anagrams]

High-Level Programming

Introduction to High-Level Programming

Traditional (procedural) programming involves setting up loops and branching structures, such as Do, If, and While, typically repeating a calculation (such as incrementing a variable) for a series of values.  High-level programming allows us to work on lists and other expressions as a whole, rather than looping over their elements.

In traditional programming, we might add the numbers from 1 to 100 by initializing a partial-sums variable, then looping over each value of the sum and adding it to the partial-sums variable. Following is one such procedural way to add the numbers from 1 to 100 using Do.

psum = 0; Do[psum = psum + i, {i, 1, 100}]; psum

High-level programming allows us to add the numbers by creating a list of integers from 1 to 100, then directly changing the list into a sum, using Apply, defined below.

Apply[Plus, Range[100]]

High-level programs allow us to take advantage of knowing Mathematica’s internal representation of an expression, and consequently create programs that execute much more efficiently.

Apply

We often collect lists of elements to be given as arguments to some operator, and the Mathematica function that allows us to do this is Apply. If we collect arguments to be given to a function f in a list called args, we might use Apply as follows. Apply effectively replaces the head of an expression, in this case replacing the head of args (List) with the head f.

args = {a, b, c, d};

Apply[f, args]

Notice that we give Apply (and the other functions in this chapter) only the function name (f), not the name and argument (f[x]).

Recall that the internal form of a + b is Plus[a, b], the internal form of a+b+c is Plus[a, b, c], and so forth. Considering this fact, we can give the list of numbers from 1 to 100 directly as an argument to the Plus function, by using Apply to change the head of a list of numbers to Plus.

Apply[Plus, Range[100]]

In this way we can create an elegant function for computing the mean of a list of numbers by adding all the elements of the list, then dividing by the length of the list.

mean[L_List] := Apply[Plus, L]/Length[L]

mean[{0.2, 2.0, 20.0, 200.0}]

Similarly, the following method of computing 10 factorial (10!) is much faster than a recursive definition. Rather than having to call Times ten times, this definition uses it only once.

Apply[Times, Range[10]]

Apply is useful when giving a list of arguments to a function that expects a sequence of arguments. For example the function ArcTan expects a sequence of two arguments, while a point in the plane (such as {3.0,4.0}) is represented by a list of two arguments. To compute the arc tangent of a point in the plane, we apply ArcTan to the list of coordinates.

Apply[ArcTan, {3.0, 4.0}]

Mathematica’s high-level programming constructs operate not only on lists, but on general expressions. For instance, Apply replaces the head of any expression with another expression.

Apply[List, a + b + c]

Map

Map accepts a function and a list of arguments, and returns a list of function values. Here is a function that adds one to the reciprocal of its integer argument.

fun[t_Integer] := 1 + (1/t)

By mapping fun over {1, 2, 3} we obtain the list of values {fun[1], fun[2], fun[3]}.

Map[fun, Range[3]]

Given a list of ordered triples, for instance, we can use Map to map the graphics primitive Cuboid over the list, obtaining a list of cubes at the specified coordinates. Here is a list of ordered triples.

mydata = 
      5.0 Table[{Sin[j], Cos[j], Cos[j]}, {j, 0, 2 π, (π/12)}];

Here we map Cuboid over the list of triples.

cubes = Map[Cuboid, mydata];

Here are the first five cuboids in the list cubes.

Take[cubes, 5]

The resulting list of graphics primitives can then be displayed.

Show[Graphics3D[cubes]];

Here is a list called elts that contains several records, where each record contains the name of a chemical element, its abbreviation, its heat of fusion (in kilojoules per mole), and its atomic number.

elts =

{{Hydrogen, "H", 0.12, 1},
 {Lithium, "Li", 4.6, 3}, 
 {Sodium, "Na", 2.64, 11},
 {Potassium, "K", 2.4, 19}, 
 {Rubidium, "Rb", 2.2, 37},
 {Cesium, "Cs", 2.09, 55}, 
 {Lanthanum, "La", 10.04, 57}, 
 {Terbium, "Tb", 16.3, 65}};

Here we define a function placeabbr that accepts a list of four items, then returns a graphics primitive of the abbreviation placed at the {atomic number, heat of fusion} coordinate. (The element name is discarded.)

placeabbr[{_, ab_String, hof_, atno_}] := Text[ab, {atno, hof}]

Mapping placeabbr over the list elts generates a list of the text graphics primitives corresponding to each record in elts.

graphelts = Map[placeabbr, elts]

Here is a plot of the resulting graphics primitives.

Show[Graphics[graphelts], Frame -> True, 
      PlotRange -> {{-1, 67}, {-1, 18}}, 
      TextStyle -> {FontFamily -> "Times", FontSize -> 9}, 
      FrameLabel -> {"Atomic number", "Heat of fusion"}];

If we map a function onto an array, Mathematica maps the function onto each row of the array.

Map[f, {{1, 2}, {3, 4}}]

To map a function to the elements of an array, we specify a level at which the function should be mapped. The elements of a matrix are at the second level of the matrix, so we enter the following.

Map[f, {{1, 2}, {3, 4}}, {2}]

Like Apply, Map operates over general expressions.

Map[Sqrt, a + b + c]

Iteration

Whereas recursion involves explicitly defining a function in terms of itself (as with the factorial definition fac[n_] := n fac[n - 1], in which the function name appears on both the left- and right-hand sides of the definition), iteration uses high-level functions defined specially for having a function call itself.

NestList accepts a function to iterate, a starting value, and a number of iterations, and returns a list of nested values of the function. (Nest returns the last value of NestList.)

NestList[f, x, 4] // ColumnForm

Here are the first ten iterates of x → 1 + (1/x), starting with x = 1.0.

flip[x_] := 1 + (1/x)

NestList[flip, 1.0, 10] // ColumnForm

Here is one way to generate a two-dimensional random walk. We start with the ordered pair {0, 0}, and iterate a function that adds a random integer from -1 to 1 to each coordinate. The result is a list of ordered pairs.

\(onestep[{x_, y_}] := 

      	{x + Random[Integer, {\(-1\), 1}], y + Random[Integer, {\(-1\), 1}]}
        \)

onestep[{0,0}]

The following command generates a list of iterates of onestep. The output is suppressed by ending the line with a semicolon.

walkdata = NestList[onestep, {0, 0}, 999];

Here is an abbreviated version of the list of points.

Short[walkdata]

To convert walkdata into a graph, we connect and plot the ordered pairs using ListPlot with the option PlotJoined set to True.

ListPlot[walkdata, PlotJoined -> True, Axes -> False, 
      AspectRatio -> Automatic];

FoldList and Fold are used to iterate functions of two variables. Here we use FoldList with symbolic arguments.

FoldList[f, x1, {x2, x3, x4}] // TableForm

Here we generate a list of partial products of elements in a list using FoldList with the operator Times.

FoldList[Times, 1, {2, a, 3, b, 5, c}]

Here are two ways to get the final result.

Apply[Times, {2, a, 3, b, 5, c}]

Fold[Times, 1, {2, a, 3, b, 5, c}]

Mathematica contains a function called IntegerDigits that returns a list of digits of a given number in a given base. Here are the base-2 digits of the decimal number 714.

IntegerDigits[714, 2]

To reverse the process, here is one way to convert a list of binary digits into a decimal number using Fold.

bitlist = {1, 0, 1, 1, 0, 0, 1, 0, 1, 0};

Using Fold, we iterate a function that multiplies its first argument by 2, then adds the second argument.

shiftbits[x_, y_] := 2 x + y

Folding shiftbits over the list of bits returns the decimal value.

Fold[shiftbits, 0, bitlist]

Our function reproduces the built-in function FromDigits.

FromDigits[bitlist,  2]

FixedPoint iterates a function until two successive iterates are equal, or the maximum number of iterations are met. Here we iterate x → 1 + (1/x), starting with 1.0, until a fixed point is reached.

FixedPoint[flip, 1.0]

The final value is the golden ratio.

N[GoldenRatio]

Here is part of a fractal Julia set calculation, one that follows the iterations of z → (z^2) + i, starting with the value 0.1 + 0.1  i. In this case we specify that a maximum of 35 iterations should be carried out, and that the iteration should stop when an iterate reaches a distance greater than two units from the origin.

iter[z_] := (z^2) + I

stop[z1_, z2_] := Abs[z2] > 2.0

FixedPointList[iter, 0.1 + 0.1 I, 35, SameTest -> stop]

To plot the path of the iteration, we again convert the complex numbers into ordered pairs, then use ListPlot to plot the ordered pairs.

ListPlot[% /. Complex[a_, b_] -> {a, b}, 
      PlotJoined -> True, AspectRatio -> Automatic];

Using this technique, we create an elegant program to generate the fractal Mandelbrot set.

iter2[z_] := (z^2) + a + b I

DensityPlot[
      -Length[FixedPointList[iter2, a + b I, 35, SameTest -> stop]]
          , {a, -2, 1/2}, {b, -1, 1}, Mesh -> False, 
      AspectRatio -> Automatic, PlotPoints -> 75];

Other High-Level Programming Functions

Thread

Thread constructs lists by collecting “parallel” elements from a series of lists, then giving each group of elements a specified head.

Thread[f[{a, b}, {1, 2}, {A, B}]]

If one of the elements in the series of lists is an atom, the atom appears in each expression created by Thread.

Thread[f[{a, b}, 1, {A, B}]]

Mathematica effectively uses Thread when performing list arithmetic, such as adding two lists or adding a scalar to each element of a list.

{a, b, c} + E + {1, 2, 3}

Thread[Plus[{a, b, c}, E, {1, 2, 3}]]

Thread is useful in constructing tables of rules or equations, given a list of names and a list of values.

Thread[{a, b, c} -> {1, 2, 3}]

Thread[{a, b, c} == {1, 2, 3}]

Thread duplicates Transpose if given a list of lists with equal lengths.

Thread[{{a, b, c}, {1, 2, 3}}]

Transpose[{{a, b, c}, {1, 2, 3}}]

Through

Where Map applies a function to a list of arguments, Through applies a list of functions to an argument.

Through[{Sin, Cos, Log, Exp}[x]]

Here is a list of complex numbers.

fdata = Fourier[Table[Random[], {10}]]

In order to plot the real parts of the data against the imaginary parts, we need to convert the complex numbers into a list of real-imaginary pairs. Using Through, we can create a list containing a list of all the real parts followed by a list of all the imaginary parts.

fpairs = Through[{Re, Im}[fdata]]

To convert the previous list into a list of real-imaginary pairs, we use Transpose.

fpairs = Transpose[fpairs]

Here is a plot of the real parts versus the imaginary parts of the transformed data.

ListPlot[fpairs, PlotStyle -> PointSize[0.02]];

Outer

Outer[head, (Subscript[list, 1]), (Subscript[list, 2]), …] constructs lists of all possible combinations of elements from more than one list, even if they have unequal lengths, then applies the specified head to each of the lists.

The following diagram illustrates what happens when we enter Outer[List, {a, b}, {c, d, e}, {f, g}]. Mathematica constructs triples created by following each of the downward paths to the bottom of the diagram, then applies the head List to each of the triples.

Here is the output.

Outer[List, {a, b}, {c, d, e}, {f, g}] // TableForm

Here we compute all possible pairs of elements, one each from {1, 2, 3} and {a, b, c, d}, and then give each pair the head Plus.

Outer[Plus, {1, 2, 3}, {a, b, c, d}] // TableForm

Here is a concise way to construct the Jacobian of two functions f and g using Outer. The following command generates a table of all pairings of one function (f or g) and one variable (x or y), and then for each pair takes the derivative D of the function with respect to the variable.

Outer[D, {f[x, y], g[x, y]}, {x, y}] // MatrixForm

Here is a similar method of constructing the Wronskian matrix for a list of functions (used to verify fundamental solution sets of differential equations). The function WD sets up the differentiation function D for use with pairs consisting of a function of x and the number of times it is to be differentiated.

WD[fn_, n_] := D[fn, {x, n}]

Here is some sample output of WD.

{WD[f[x], 1], WD[g[x], 3]}

Here are the lists of the functions f, g, and h each differentiated 0, 1, and 2 times.

Outer[WD, {f[x], g[x], h[x]}, {0, 1, 2}]

In order to generate the Wronskian matrix of f, g, and h, we need to group the same-order derivatives, which we do by transposing the previous matrix.

Transpose[%] // MatrixForm

Exercises: High-level Programming

Define a function scatterPlot3D that accepts a list of ordered triples, and generates a three-dimensional plot of the triples.

To generate a scatter plot of three-dimensional points, we map the three-dimensional graphics primitive Point over each ordered triple.

Map[Point, {{1, 2, 3}, {2, 3, 4}, {3, 4, 5}}]

To turn this procedure into a function, we must define a function that accepts a list of triples.

scatterPlot3D[triples : {{_, _, _}..}] := 
      Show[Graphics3D[Map[Point, triples]], Axes -> True]

Here is a list of 100 random ordered triples.

data3D = Table[Random[], {100}, {3}];

We see that the data are of the desired shape.

Short[data3D, 3]

Here is a test of the function.

scatterPlot3D[data3D];

Clear[data3D]

Define a function grayScatterPlot3D that accepts a list of 4-tuples (a 4-tuple is an ordered list of four numbers), where the fourth argument is a number between 0 and 1 that describes the gray level of an ordered triple.

Here is a function that accepts a 4-tuple and returns a list of graphics primitives, the gray level followed by the point at the proper location.

makePoint[{x_, y_, z_, g_}] := {GrayLevel[g], Point[{x, y, z}]}

Here is a test of the function.

makePoint[{1, 2, 3, 0.5}]

Here is the function definition. Notice that we include the PointSize[0.02] directive inside the list given to Graphics3D to change the size of all points in the graph.

\(\(grayScatterPlot3D[tuples : {{_, _, _, _}..}] := 
	
      Show[Graphics3D[{PointSize[0.02], Map[makePoint, tuples]}], 
        Axes -> True]; \)\)

Here is a list of 4-tuples.

data4D = Table[Random[], {100}, {4}];

Here is the first element in the list of data.

First[data4D]

Here is a test of the function.

grayScatterPlot3D[data4D];

Clear[data4D]

Define a function that accepts a list of lists, and returns a list of the averages of each of the sublists.

To compute the average of one list, we use Apply.

thelist = {a, b, c};

Apply[Plus, thelist]/Length[thelist]

To compute the average of each sublist in a list, we map the previous procedure over a list of lists.

avgone[onelist_] := Apply[Plus, onelist]/Length[onelist]

The repeated pattern L:{_List..} matches a list containing one or more lists.

averages[L : {_List..}] := Map[avgone, L]

averages[{{1, 1, 2, 3, 5, 8, 13, 21}, {a, b, c}, {1, a, 2, b, 3}}]

Use Outer to generate every possible pair of elements (for simplicity, allowing repeated pairs) from the list {a,b,c,d,e}.

The following command takes one element from each copy of the list, giving each pair the head List.

allpairs = Outer[List, {a, b, c, d, e}, {a, b, c, d, e}]

One way to remove duplicates is to use Flatten and DeleteCases. If we flatten the list (remove all the inner brackets) and then partition the result into lists of length two, we get a list of all the pairs.

Partition[Flatten[allpairs], 2]

Flatten[allpairs, 1] also removes the brackets from the first level of allpairs, leaving the deeper list brackets intact.

Flatten[allpairs, 1]

Once we have a list containing just the collection of pairs (with no extra list brackets), we can remove duplicates using DeleteCases.

nonmatchingpairs = DeleteCases[%, {x_, x_}]

If the order of elements in a pair is unimportant, we can remove pairs that contain the same elements regardless of order by first sorting each pair, then removing duplicates with Union. Here we map Sort over the list of pairs, which sorts each pair.

Map[Sort, nonmatchingpairs]

Taking the union of the resulting list with itself removes matching pairs.

Union[%]

As expected, there are (({{5}, {2}})) distinct pairs.

Length[%] == Binomial[5, 2]

Use Outer to create a list of strings that look like the binary numbers 000, 001, 010, 011, ..., 111.

We first consider the easier problem of creating the four pairs 00, 01, 10, and 11. We want to create a list of all pairs formed by choosing one element from a list {0, 1} followed by one element from a second list {0, 1}. To begin, we give each pair the generic head f.

Outer[f, {0, 1}, {0, 1}]

Because we want to create strings, we replace the digits 0 and 1 with the strings “0” and “1”.

Outer[f, {"0", "1"}, {"0", "1"}]

To concatenate each pair of strings, we replace the head f with the function StringJoin, which returns the desired list of strings.

Outer[StringJoin, {"0", "1"}, {"0", "1"}]

To generate the list of three-digit binary numbers, we include in Outer a third copy of the list {"0", "1"}.

Outer[StringJoin, {"0", "1"}, {"0", "1"}, {"0", "1"}]

Removing the inner braces from the previous list will give us a flat list of the desired strings.

Flatten[%]

Beginning with the string "F", use Nest and StringReplace to iterate the rule that replaces every occurrence of "F" in a string with "F+F--F+F", three times.

Here is one call to the function StringReplace.

StringReplace["F", "F" -> "F+F--F+F"]

Here is a function that carries out the replacement for any string.

srep[s_String] := StringReplace[s, "F" -> "F+F--F+F"]

To iterate the replacement we use Nest.

Nest[srep, "F", 3]

Create a one-dimensional random walk using NestList.

We begin with a height value of zero.

height = 0.0;

Here is a function that adds a random number between -1 and 1 to its argument.

step1d[x_] := x + Random[Real, {-1, 1}]

If we use Nest to iterate this function, we find the final position of the random walk.

Nest[step1d, 0.0, 500]

By using NestList, we collect all of the intermediate iterations.

randomwalk = NestList[step1d, 0.0, 500];

Here is a shortened version of the list of points on the random walk.

Short[randomwalk]

To plot the random walk, we use ListPlot.

ListPlot[randomwalk, PlotJoined -> True];

Clear[randomwalk]

The expected value of a discrete random variable Y is defined to be the sum (over all possible values of Y) of the product of the value of an event  y and the probability of y occurring ((∑_y) y p(y)). Compute the expected value corresponding to the list {{1, 0.5}, {2, 0.25}, {3, 0.2}, {4, 0.05}}, where each sublist is a {y, p(y)} pair.

To multiply two elements of a list, we use Apply and Times.

Apply[Times, {2, 0.25}]

Here is an auxiliary function that multiplies the two elements of its argument.

multPairs[L_] := Apply[Times, L]

An identical definition is the following.

multPairs[{x_, y_}] := x y

To multiply each pair in the list, we map either definition of multPairs over the list.

Map[multPairs, {{1, 0.5}, {2, 0.5}, {3, 0.2}, {4, 0.05}}]

To add the elements in the resulting list, we use Apply and Plus.

Apply[Plus, Map[multPairs, {{1, 0.5}, {2, 0.5}, {3, 0.2}, {4, 0.05}}]]

The expected value of Y is 2.3.

Pure Functions

Motivation

The functions we have so far examined have had the following form.

sqrtPlusOne[x_] := 1 + Sqrt[x]

Once such a function is defined, the symbol name sqrtPlusOne remains in memory until we remove it with Remove, which deletes both the value and the name of a symbol. (Clear deletes only the value of a symbol.)

Pure, or anonymous, functions are functions that do not have names; by using pure functions in appropriate places, we decrease the number of symbol names in memory, and consequently decrease the possibility of a name conflict.

Pure functions are useful inside operators that take a function as an argument, as with the high-level programming constructs like Map and Apply, particularly when the argument function will not be used again.

Defining Pure Functions

Here is a simple definition of a function named func with a single dummy variable called x.

func[x_] := 1 + (x^2)

We rewrite the definition of func[x_] as the following anonymous function. The symbol Function is a container for an anonymous function, and in this case it contains a dummy variable name x (and not a pattern such as x_) and a function body that describes what to do with the variable.

Function[x, 1 + (x^2)]

We give a pure function an argument inside square brackets, just as with a named function.

Function[x, 1 + (x^2)][10]

func[10]

One advantage to using a pure function is that there is no symbol name (func in the case of the named function) consuming memory.

We define a pure function of two arguments similarly.

Function[{x, y}, 1 + (x^y)]

Here we give the pure function two arguments, again inside square brackets.

Function[{x, y}, 1 + (x^y)][10, 3]

A shorter form of a pure function equivalent to func is the following. Each occurrence of the hash mark # represents a slot into which the function’s argument will be substituted. This format does not require that we name the dummy variables used in the definition.

\(\(Function[1 + #\^2]\)[10]\)

We define pure functions that accept more than one argument by numbering the slots #1, #2, and so forth.

\(\(Function[1 + #1\^#2]\)[10, 3]\)

The shortest form for a pure function that is equivalent to func is the following. Notice that the pure function definition looks like the right-hand side of the definition for func, but with # taking the place of the dummy variable x. The ampersand & marks the end of the pure function definition, and it is important not to omit it.

\(\((1 + #\^2)\)&\)

We give this pure function an argument inside square brackets, just as with the other forms.

\(\(\((1 + #\^2)\)&\)[10]\)

Again we number the slots if the function accepts more than one argument. The following function raises the first argument to the power of the second argument.

\(\((#1\^#2&)\)[E, \(-π\)]\)

A pure function may have as many arguments as desired. Here is a pure function that accepts four arguments.

\(\((#1\^#2 < #3\^#4&)\)[2, 3, 3, 2]\)

Here is a pure predicate function. The following function returns True if its argument is between 3 and 10. (Note that # is an abbreviation for #1, representing the first argument of a function.)

(3 < # < 10)&

This is the equivalent long form of the pure function.

Function[x, 3 < x < 10]

Here are two tests of the function.

(3 < # < 10)&[7]

Function[x,  3 < x < 10][-7]

Interpolating functions in Mathematica behave as pure functions. Here is an interpolating function that goes through a list of the first ten prime numbers.

ip = Interpolation[Prime[Range[10]]]

We give the resulting InterpolatingFunction object an argument inside square brackets, just as with a pure function.

ip[5.25]

Here is a plot of the interpolating function. (Note the use of a pure function in the value given to Epilog.)

Plot[ip[x], {x, 1, 10}, 
      Epilog -> {PointSize[0.02], 
          Map[Point[{#, Prime[#]}]&, Range[10]]}];

Similarly, the derivative of a built-in function without any arguments is returned in the form of a pure function. Here is the derivative of the Sin function.

Sin'

Using Pure Functions as Arguments

Pure functions are useful as criteria inside functions like Select. Here we select from the list of integers from 1 to 100 all the elements whose square root is an integer (that is, all the perfect squares) using the pure predicate function IntegerQ[Sqrt[#]]&, which returns True if the square root of its argument is an integer.

Select[Range[100], IntegerQ[Sqrt[#]]&]

Here we select integers between 3 and 10 from the integers from 1 to 15.

Select[Range[15], 3 < # < 10 &]

Pure predicate functions are often used in conditional patterns. The following function returns True if its argument is between 3 and 10, and otherwise returns False.

between3and10[x_?(3 < #1 < 10&)] := True

between3and10[x_] := False

{between3and10[7], between3and10[17]}

Sort accepts an optional second argument, a function that determines if one element should be listed before another. Here we use a pure function as the ordering function inside Sort, sorting our array of hurricane data in order of descending wind speed.

hurr =
	{{"Allison",   65,  987},
	 {"Barry",     50,  995},
	 {"Chantal",   60,  992},
	 {"Dean",      40,  999},
	 {"Erin",      75,  975},
	 {"Felix",    115,  950},
	 {"Gabrielle", 60,  998},
	 {"Humberto",  90,  970},
	 {"Iris",      75,  970},
	 {"Jerry",     35, 1003}};

The pure sorting function expects two arguments, which we label #1 and #2; and we arrange the sublists of hurr so that the second part of the first argument is greater than the second part of the second argument.

TableForm[Sort[hurr, #1[[2]] > #2[[2]]&]]

Pure functions are useful inside high-level functions such as Map and Apply. Here we map a pure function that computes a numeric approximation to the base-2 logarithm of its argument over a list of numbers.

Map[Function[N[Log[2, #]]], Range[2, 10]]

Here is the same calculation with the pure function written more compactly.

Map[N[Log[2, #1]]&, Range[2, 10]]

Exercises: Pure Functions

Map a pure function that squares its argument over a list of numbers ranging from 1 to 10.

Here is a pure function that squares its argument.

Function[x, (x^2)]

Here is the function applied to a single argument.

Function[x, (x^2)][16]

We give the pure function as the first argument to Map.

Map[Function[x, (x^2)], Range[10]]

Here is a more compact form.

\(Map[#\^2&, Range[10]]\)

The function Power automatically maps itself over a list, so using Map is unnecessary.

\(\(#\^2&\)[Range[10]]\)

Using a pure function as an argument to FixedPoint, find the fixed point of the iteration x → 1 + (1/x), where x has a starting value of 1.0. Find the fixed point of x → (Sqrt[1 + x]), again with x starting at 1.0.

Here is a pure function that adds 1 to its reciprocal.

Function[x, 1 + 1/x]

Here is the same function given as an argument to FixedPoint.

FixedPoint[Function[x, 1 + 1/x], 1.0]

Here is a more compact version of the same expression.

FixedPoint[1 + 1/#&, 1.0]

Here is the fixed point of x → (Sqrt[1 + x]).

FixedPoint[(Sqrt[1 + #])&, 1.0]

Both iterations converge to the golden ratio.

N[GoldenRatio]

Create a one-dimensional random walk using NestList and a pure function.

Inside NestList we use a pure function that adds a random number (in this case an integer between -1 and 1) to its argument.

NestList[# + Random[Integer, {-1, 1}]&, 0, 10]

Here is a longer list of points on the random walk.

NestList[# + Random[Integer, {-1, 1}]&, 0, 200];

ListPlot[%, PlotJoined -> True];

Use Table and FindRoot to generate a list of roots of the equation sin((x^3)) = cos(3 (x^2)), with starting values for FindRoot ranging from -3.1 to 3.1 in increments of 0.25. Remove duplicates in the list of roots using Union by creating a pure function as the value of the option SameTest, considering two roots duplicates if they differ by less than (1/10).

Here is a list of the roots.

Table[
      FindRoot[Sin[(x^3)] == Cos[3 (x^2)], {x, st}], {st, -3.1, 3.1, 0.25}]

To extract only the roots, we substitute the list of rules into the variable x.

theroots = x /. %

Union removes duplicates of a single list.

Union[theroots]

By default, Union considers two elements duplicates if they are equal. However, because FindRoot uses numerical methods, uncertainty can appear in the roots it finds, and the uncertainty can cause Mathematica not to recognize that different approximations to the same root are the same number.

To change the default criterion that determines which elements are duplicates, we change the value of the option SameTest. SameTest expects a function of two variables that returns True if its two arguments are “the same”. Here is a named function that returns True if the distance between its two arguments is less than (1/10).

theSameQ[x_, y_] := Abs[x - y] < (1/10)

theSameQ[9.75, 9.82]

Here we use Union with a pure function equivalent to theSameQ as the value of SameTest. Here is the list of the roots with the duplicates (as we have defined them) removed.

Union[theroots, SameTest -> (Abs[#1 - #2] < (1/10)&)]

Packages

Contexts

In addition to the notion of variables and constants that are local to a routine (as created by Module, Block, and With), Mathematica supports the organization of symbols with the notion of contexts.

The idea of contexts is similar to the idea of different directories in a file system. Just as there is no conflict between two files with the name ReadMe.txt if the two files are contained in different directories, there is no conflict between two symbols called newvar if they are defined within different contexts.

Every symbol defined in Mathematica is defined within a context, and the full name of a symbol x is context`x, where context is the name of the context in which the symbol is defined and x is called the short name of the variable. We call the backquote character ` a context mark, similar to the directory separator character for full file pathnames (“:” on a Macintosh file system,  “/”on a Unix file system, and “\” on a Microsoft Windows file system). The default context for new symbols is Global`, and we create symbols in other contexts at any time by simply specifying the full name of the variable, in particular including the context name.

Here we define a new variable called newvar. Mathematica places the symbol in the context Global`, as revealed by entering ?newvar.

newvar = 13; ?newvar

In every Mathematica session the system variable $Context stores the current context. When referring to a symbol within this context we may use its short name.

$Context

To create another symbol with the short name newvar, we specify a new context in which to place it, in this case mycontext`.

mycontext`newvar = (A^2) + (B^2);

When we refer to the short name newvar, Mathematica looks first in the current context, Global`.

newvar

We can always refer to a symbol in another context by using its full name.

mycontext`newvar

We can tell Mathematica to look first in the context mycontext` by changing the value of $Context.

$Context = "mycontext`"

When we type the short name newvar, Mathematica now looks first in the context mycontext`.

newvar

However, we can refer to the variable newvar in the context Global` by using its full name.

Global`newvar

We can make Global` the current context again by resetting the value of $Context.

$Context = "Global`"

The global variable $ContextPath contains a context search path, which is similar to the directory search path used when working with the local file system. When we use the short name of a symbol that is not defined in the current context, Mathematica will search this list of contexts for a symbol with the given short name.

$ContextPath

Contexts may be hierarchical, so the full name of a symbol may look like (Subscript[c, 1])`(Subscript[c, 2])`(Subscript[c, 3])`name, where the (Subscript[c, i]) are context names.

hidden`deeply`in`a`new`context`newvar = (1/(Sqrt[2  π]));

A Mathematica package typically defines its own context. When we begin a Mathematica session, the default context is Global` and symbols we enter will be in this context unless we specify otherwise. $ContextPath will, however, contain several other contexts, including System`, where built-in symbols such as Pi reside, and the contexts of Mathematica packages that have been loaded. Context[symbol] gives the first context in which a symbol appears.

Context[Pi]

Context[newvar]

When a variable does not appear in the current context, Mathematica looks in other contexts in the order defined in the list $ContextPath. Here is the current value of $ContextPath.

$ContextPath

The previous list indicates that Mathematica first looks for symbol names in the Global` context, where user-defined variables and functions are created, then in the System` context, where built-in functions and objects are defined. We can assign a new value to the symbol Pi by creating a symbol with the short name Pi in the Global` context, with the value 3.14159.

Global`Pi = 3.14159;

Because Mathematica looks for symbols in the context Global` before looking in System`, the value of Global`Pi is used when the short name Pi is used. When the same short name appears in more than one context, we say the first occurrence of the name shadows any other occurrence.

Pi (r^2)

We can refer to the built-in value of Pi by using its full name.

System`Pi (r^2)

In order to remove the symbol Global`Pi, we cannot simply use Clear because Clear erases the value of Global`Pi but leaves the symbol name in memory.

Clear[Pi]

Mathematica still understands the short name Pi is in the context Global`, even though it has no value.

Context[Pi]

To remove the symbol name from memory, we must use Remove.

Remove[Pi]

Mathematica now recognizes that the short name Pi refers to System`Pi.

Context[Pi]

Shadowing problems arise when we try to use a function defined in a package before loading the package. Here we attempt to use the function ImplicitPlot before loading the package Graphics`ImplicitPlot`. Because Mathematica has no definition of ImplicitPlot, the input is returned unevaluated.

ImplicitPlot[(x^2) - (y^2) == 1, {x, -2, 2}]

If we now load the package Graphics`ImplicitPlot`, Mathematica warns us that the symbol ImplicitPlot is defined in more than one context.

Needs["Graphics`ImplicitPlot`"]

If we now attempt to use ImplicitPlot, Mathematica defers to the empty definition of Global`ImplicitPlot.

ImplicitPlot[(x^2) - (y^2) == 1, {x, -2, 2}]

Context[ImplicitPlot]

To correct this conflict, we must remove the symbol Global`ImplicitPlot.

Remove[ImplicitPlot]

The function now works as expected.

ImplicitPlot[(x^2) - (y^2) == 1, {x, -2, 2}];

The short name ImplicitPlot now refers to the symbol with the full name Graphics`ImplicitPlot`ImplicitPlot.

Context[ImplicitPlot]

The Structure of Packages

A Mathematica package typically defines a new context and several new symbols that are intended to be used outside of the package. The convention is to create a new context whose name is related to the name of the package, define all new symbols within that context, and add the new context to the beginning of the the context search path $ContextPath. In the case of Mathematica’s standard add-on packages, the context created is DirectoryName`PackageName`. On platforms that allow long filenames, a package’s context corresponds exactly to the file pathname (where the file pathname separator is replaced with the context separator `). On platforms that allow only short file names, the file names are shortened versions of the context names. The command ContextToFileName[context] returns the string specifying the file name appropriate to the local file system.

Here is the file name corresponding to the context Graphics`ImplicitPlot`. (On the computer used to create this notebook, “:” is the file pathname separator.)

ContextToFileName["Graphics`ImplicitPlot`"]

There are four basic functions for creating and moving among contexts that are used when creating a Mathematica package. BeginPackage["context`"] makes context ` and System` the only active contexts. BeginPackage["context`", {"(Subscript[need, 1])", "(Subscript[need, 2])", …}] also loads the packages (Subscript[need, i]) with Needs["(Subscript[need, i])"]. EndPackage[] restores $Context and $ContextPath to their values before BeginPackage was invoked, and prepends the new context to $ContextPath. Begin["context`"] makes context` the current context; End[ ] returns the present context, and reverts to the previous one.

With these functions we construct the basic outline of a well-designed Mathematica package.

Here are the steps we follow to create a package. First we use BeginPackage to make the package context the current working context, if necessary loading any auxiliary packages.

BeginPackage["PackageName`","(Subscript[Package, 1])`","(Subscript[Package, 2])`\(\*"","\)…]

Next we define usage messages for the functions that will be visible to the user. (Mathematica’s message-printing mechanism is explained in the chapter “Special Topics”.)

pkgfn::usage = "Usage message for pkgfn."

After defining messages for the symbols that are to be exported, we create a new working context for symbols that are not to be exported by using Begin. By convention, the context for hidden symbols is PackageName`Private`. Notice that the argument to Begin starts with a context mark `, which indicates that we are creating a subcontext of the current context.

Begin["`Private`"]

In the `Private` context we enter any function and variable definitions. Symbols created here that do not have usage messages will be hidden from the user. The function aux defined here will not be visible to a user of the package because it does not have a usage message.

aux[n_] := defn

pkgfn[n_] := n + aux[n] + …

When we have finished making definitions, we end the `Private` context using End.

End[];

Finally, we end the package context.

EndPackage[];

Here is an example of a small package that defines some combinatorial functions.

Because each function has a usage message defined before Begin["`Private`"], it will be visible to the user.

\(BeginPackage["DiscreteMath`CombinatorialFunctions`"]; 


    Subfactorial::usage\  = 
	
      "Subfactorial[n] gives the number of permutations of n objects which leave no object fixed."; 


    CatalanNumber::usage\  = 
	
      "CatalanNumber[n] gives the nth Catalan number."; 


    Hofstadter::usage\  = 
	
      "Hofstadter[n] gives Hofstadter's 		function for positive integer n."; 

Begin["`Private`"]; 


    Subfactorial[n_Integer?Positive]\  := \ 
	
      \(n!\)\ Sum[\((\(-1\))\)^k/\(k!\), \ {k, \ 0, \ n}]; 


    CatalanNumber[n_Integer]\  := \ Binomial[2  n, n]/\((n + 1)\); 


    Hofstadter[n_Integer?Positive]\  := 
	
      \(Hofstadter[n]\  = 
	
        Hofstadter[n\  - \ Hofstadter[n - 1]]\  + \ 
          Hofstadter[n\  - \ Hofstadter[n - 2]]\); 


    Hofstadter[1]\  = \ \(Hofstadter[2]\  = \ 1\); 

End[\ ]; 

    EndPackage[\ ]; \)

Before loading the package, the function Subfactorial is undefined.

?Subfactorial

Using Needs is the best way to load packages into Mathematica.

Needs["DiscreteMath`CombinatorialFunctions`"]

The package context DiscreteMath`CombinatorialFunctions` is now at the beginning of $ContextPath.

$ContextPath

The symbol Subfactorial is defined in the context DiscreteMath`CombinatorialFunctions`.

Context[ Subfactorial ]

We now use Subfactorial as if it were a built-in function.

?Subfactorial

Subfactorial[10]

ListPlot[Log[Map[Subfactorial, Range[2, 10]]], 
      PlotStyle -> PointSize[0.02]];

Typically packages are created in the form of Mathematica notebooks, and are stored in the AddOns/StandardPackages subdirectory inside the main Mathematica directory. When creating notebooks that contain packages, it is necessary to give the cells in the notebook the initialization cell property, by selecting the cells in the notebook, pulling down the Cell menu, and choosing Initialization Cell from the Cell Properties submenu. (Alternatively, with the package’s notebook open we can pull down the File menu and choose Package Format from the Save As Special submenu.) It may also be necessary to append the directory that contains the package to the global variable $Path if Mathematica cannot find the appropriate file.

When learning to create packages it is helpful to examine and modify existing, working packages.

Special Topics

Usage and Warning Messages

The mechanism Mathematica uses to return usage and warning messages can be used to pass usage and warning messages for user-defined functions. A message has the general form symbol::tag, where symbol is the function or object that the message is attached to, and tag is a unique identifier attached to the message. (The internal form of the form symbol::tag is MessageName[symbol, "tag"].)

To define a message, we set the appropriate symbol::tag object equal to the string containing the text of the message. If we want to create a message for the function func, we assign it as follows, arbitrarily giving the message the tag generalmsg. (We first clear any previous definitions of func.)

Clear[func]

func::generalmsg = "This is a message for func.";

We use Message to print a message using the same mechanism as built-in functions.

Message[func::generalmsg]

We use Message on the right-hand side of a function definition to generate warning or error messages for users of func. Here the message is printed whenever func is called with a single argument.

func[_] := Message[func::generalmsg]

func[0]

A special tag for messages is “usage”; symbol::usage is defined to be the usage message for symbol, which is printed when we enter ? symbol. Here is a usage message for func.

func::usage = "func is a function that prints a message.";

The usage message is now returned when we enter ? func.

? func

Here is a more complete definition of the factorial function fac that includes message definitions. First we clear any previous definitions of fac.

Clear[fac]

Here is a usage message for fac, which is printed when we enter ? fac.

fac::usage = "fac[n] returns n!";

Here are the two rules for computing fac of a positive integer argument.

fac[n_Integer?Positive] := n fac[n - 1];

fac[0] = 1;

Here is an error message that is printed when fac is given the wrong number of arguments, or an argument that is not a positive integer.

fac::wrongargs :=
  "fac expects a single positive integer argument."

fac[___] := Message[fac::wrongargs]

Here is an example of fac given two arguments.

fac[2, 3]

We can define more sophisticated messages that print helpful information about errors. The following commands define a message that is returned when fac is called with a single argument that is not a positive integer. Inside the string portion of the message we include two backquotes ``, a form that marks the position into which an expression will be spliced.

fac::error :=
 "Positive integer argument expected in fac[``]."

We then add a definition which states that when fac is called with a single numerical argument that is not a positive integer, the message fac::error will be printed with the argument spliced into the string of the message.

fac[wrong_?NumericQ] := Message[fac::error, wrong]

Mathematica now splices the incorrect input into the warning message.

fac[-6]

Options

Suppose we wish to define a function printvec that prints a vector in either row or column form, where the direction is controlled by the setting of an option called Direction.

Here is the left side of a function printvec that takes no options.

printvec[vec_?VectorQ] := body

Options are entered after the required arguments of a function. Because a user can enter any number of options, including zero, we put a triple-blank pattern after the required arguments to printvec, which matches the entire sequence of options. Options are entered as rules of the form OptionName -> value, so each option must have the head Rule. (There may be any number of options, including zero, so we use the triple-underscore pattern opts___Rule.)The left side of the definition now looks like the following.

printvec[vec_?VectorQ, opts___Rule] := body

Somewhere in the body of the definition we must extract the value, if any, of the option Direction. To replace the symbol Direction with the value typed by the user, we carry out the substitution defined by the replacement rule.

Direction /. {opts}

To complete the definition, we must have a way to substitute the default value for Direction if the user does not enter a value. First we create a list of the options that printvec accepts, including the default value on the right side of the replacement rule.

Options[printvec] = {Direction -> Horizontal};

Finally we can substitute the default value for Direction by substituting the list Options[printvec] after any user-chosen value has been substituted. The following line will first attempt to substitute a user-chosen value for Direction, and if no value has been specified, the default value will be substituted.

Direction /. {opts} /. Options[printvec]

Here is a complete definition of printvec, which uses Switch to test which of the various possible settings for Direction is in effect.

\(printvec[vec_?VectorQ, opts___Rule] := 
	
      Module[{dir = \(Direction /. {opts}\) /. Options[printvec]}, 
		
        Switch[dir, 
			Horizontal, vec, 
			Vertical, 
          ColumnForm[vec], 
			_, Message[Direction::error]]]\)

For completeness we define an error message that Direction must be set to either Horizontal or Vertical.

Direction::error =
  "Direction must be Horizontal or Vertical.";

The default setting for Direction is Horizontal, which is used if the option is omitted.

printvec[{1, 2, 3, 4}]

We can explicitly set Direction to Horizontal, if desired.

printvec[{1, 2, 3, 4}, Direction -> Horizontal]

If we set Direction to Vertical, the vector is printed using ColumnForm.

printvec[{1, 2, 3, 4}, Direction -> Vertical]

If the setting for Direction is neither Horizontal nor Vertical, the error message is printed.

printvec[{1, 2, 3, 4}, Direction -> Diagonal]

Mathematica will ignore options not defined for printvec.

printvec[{1, 2, 3, 4}, Color -> Red]

The add-on package Utilities`FilterOptions` assists in checking if the options entered are valid.

Formatting Results

The Format command allows us to define special printing forms for our functions. The following statement tells Mathematica to use the form R in output for the symbol myConstant.

Format[myConstant] := "R"

Solve[(x^2) - myConstant == 0, x]

Internally, the symbol is still represented as myConstant.

InputForm[%]

NumberForm causes numbers to be printed in a special format.

NumberForm[N[(π^10), 20], DigitBlock -> 3, 
      NumberSeparator -> {",", " "}]

Wrapping SequenceForm around a sequence of expressions prints the concatenated form of strings and other expressions. Here is a form appropriate for printing amounts of money in dollars and cents.

\(DollarForm[n_] := 
	
      SequenceForm["$", 
        PaddedForm[n, {10, 2}, NumberPadding -> {"", "0"}]]\)

DollarForm[1234.5678]

Using Mathematica’s box representation of typeset expressions, we can define a two-dimensional print form for an expression.

Format[Water] := 
      DisplayForm[
        RowBox[{SubscriptBox["H", "2"], "", "O"}]]

Format[Sugar] := 
      DisplayForm[
        RowBox[{SubscriptBox["C", "12"], "", 
            SubscriptBox["H", "22"], "", 
            SubscriptBox["O", "11"]}]]

Sugar + 3 Water

Here we change the format of the function ArcSin so that it displays as “arc sin”. (The function Unprotect allows us to modify a built-in function.)

Unprotect[ArcSin];

Format[ArcSin[z_]] := 
      DisplayForm[RowBox[{"arc sin(", ToBoxes[z], ")"}]]

ArcSin[(2^z) + (z^2)]

ArcSin[(2^z) + (z^2)] // TraditionalForm

Here we remove the new formatting rule and restore protection for ArcSin.

Clear[ArcSin];  Protect[ArcSin];

Debugging and Performance

Trace generates a list of subexpressions encountered in the evaluation of an expression. The following is a list of the expressions encountered during evaluation of the input 1 + 2(3 + 4).

Trace[1 + 2 (3 + 4)]

Here is our recursive factorial function fac.

Clear[fac]

fac[n_] := n fac[n - 1]; fac[1] = 1;

Trace[fac[3]] generates a list of the expressions encountered while evaluating fac[3].

Trace[fac[3]]

Timing tells us the amount of time the Mathematica kernel spends performing a calculation. (The value returned by Timing does not include the time taken to print the result on the screen.)

Timing[fac[20]]

By ending each statement inside the Timing command with a semicolon, we see the time spent during each calculation without displaying the result. Here we compare the timings for the recursive, functional, and built-in definitions of the factorial function.

Timing[fac[250]; ]

Timing[Apply[Times, Range[250]]; ]

Timing[250!; ]

We can speed up numeric functions by using Compile. Compile effectively assumes that the arguments we give to a function and the results of the function are machine-precision real, integer, or complex numbers, which eliminates the extra time it typically takes to determine the type of an argument.

For example, here is the timing for mapping the built-in function Sqrt over a list of integers.

(Subscript[t, 1]) = First[Timing[Map[Sqrt, Range[1, 1000]]; ]]

Here we create a compiled version of Sqrt optimized for machine-precision numerical arguments. Note that the definition of a compiled function is similar to that of a pure function.

cSqrt = Compile[{x}, Sqrt[x]]

Here we use the compiled square root function with a single argument.

cSqrt[100]

Here is the timing for mapping the compiled version of Sqrt over a list of integers. Because Mathematica assumes that the arguments to cSqrt are machine-precision numbers, we achieve a great decrease in the amount of time the function takes to execute.

(Subscript[t, 2]) = First[Timing[Map[cSqrt, Range[1, 1000]]; ]]

Here is the ratio of the two timings. On the machine used to execute these examples, the uncompiled version of Sqrt took over twenty times as long to evaluate.

(Subscript[t, 1])/(Subscript[t, 2])

A great many mathematical functions and programming commands can be compiled with Compile, which allows great increases in evaluation speed.

Attributes

Built-in and user-defined functions may have various attributes set. Attributes are switches that affect some characteristics of a function, such as how it interacts with other expressions and how its arguments are evaluated. A list of the attributes currently set for a function can be found by typing Attributes[FunctionName]. For example, most built-in functions and symbols have the attribute Protected, which prevents them from being accidentally modified.

Attributes[Pi]

If we try to alter the definition of Pi, Mathematica warns us that the function has the attribute Protected.

Pi = (22/7);

Unprotect allows us to alter a symbol that has the attribute Protected.

Unprotect[Pi]

Pi = (22/7);

Mathematica now uses the new value for Pi.

(1/(2 Pi))

To clear the value we defined for Pi, we use Unset.

Unset[Pi]; (1/(2 Pi))

We give a function or object the Protected attribute by using Protect. Here we restore the protection for Pi.

Protect[Pi]

Another common attribute of many Mathematica functions is Listable. A Listable function is automatically mapped over all the elements in a list given as an argument to the function. For instance, Sin[list] returns the sine of each element of list.

Sin[{0, Pi/2, Pi, 3 Pi/2}]

Here is a function recip that returns the reciprocal of an integer argument.

recip[n_Integer] := (1/n)

If we give recip a list as an argument, the input is returned unevaluated, because the argument has the head List, which does not match the pattern n_Integer.

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

However, if we give recip the attribute Listable, it automatically applies itself to each element in a list given as an argument. First we see that recip has no attributes set.

Attributes[recip]

We can set and clear the attributes of a function using SetAttributes and ClearAttributes. Here we set the Listable attribute for recip.

SetAttributes[recip, Listable]

We see that recip now has the attribute Listable.

Attributes[recip]

Now recip automatically maps itself over all of the elements in a list given as an argument.

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