CHAPTER			6

Functional Programming

Pascal is for building pyramids—imposing, breathtaking, static structures built by armies pushing heavy blocks into place. Lisp is for building organisms—imposing, breathtaking, dynamic structures built by squads fitting fluctuating myriads of simpler organisms into place.
[Abelson]

Some Functional Aspects of Mathematica

What is a functional programming language? Basically, it's a language in which functions can be defined and applied to arguments. It is important that the arguments of a function can themselves be other functions applied to other arguments, etc. There is an abstract, theoretical functional programming language called the lambda calculus (to be discussed in detail in Chapter 11, Section 7) which has exactly three operations: function definition, function application, and substitution—which is essentially a rewrite rule for function application. There are only a few pure functional languages (e.g., Haskel and Miranda). Most so-called functional languages, such as Scheme or ML, are impure in the sense that they have many other features to make programming more convenient. Mathematica belongs to the category of languages that have a pure functional component, but also include many other features. It has the advantage over Scheme and ML of having many built-in mathematical functions to start with that can be combined with each other to create new functions. Later in this chapter we'll look at exactly what constitutes functional programming in Mathematica.

Defining Functions

.  Pure Functions

The function Map is a feature of virtually every functional language. Normally, it takes two arguments (a function and a list), and it applies the function to every entry in the list, producing a list as the output. For instance:

Map[Sin,  {a,  b,  c}]

When using Map and similar operations (see Section 3 for more details), it is convenient not to have to name the function being mapped over a list or other expression. For instance, suppose we have two functions f and g and we want to form the sum f[x] + g[x] for all the entries in a list. One way would be to define a new function h to be the sum of f and g, and then map it down the list.

h[x_] := f[x] + g[x]

Map[h, {a, b, c}]

Instead, without defining a separate function h, the form Function[{x}, f[x] + g[x]] can be used.

Map[Function[{x}, f[x] + g[x]], {a, b, c}]

The expression Function[{x}, body] is a "pure function" with one bound variable x. This notation is essentially the lambda notation of the lambda calculus, discussed in Chapter 11, Section 7, where the form λx . body is written for the same thing. The operation Function[{x}, body] or λx . body is a canonical name for the function that does to any argument whatever body describes as being done to x. Such a function is applied to a value in the same way that built-in functions are. For instance:

Function[{x},  (x^2)][5]

One can of course give a name to such an expression and then use it as a function:

t  =  Function[{x},  Expand[(((1 + x))^3)]];

t[a]

The form Function[{x}, body] is the same as Function[body] or, in postfix notation, body&, where x in body is replaced by # (whose FullForm as we have seen is Slot[1].)  This is a much more convenient form to use in Map and, as will be seen, in many other places. For instance, our original example now becomes

Map[(f[#] + g[#])&, {a, b, c}]

The three expressions

all produce the same value when applied to an argument. The third is clearly the most concise form.

\(\(#\^2 &\)[5]\)

To map it down a list, use Map, or /@.

\(Map[#\^2 &, \ {a, \ b, \ c}]\)

\(\(#\^2 &\) /@ {a, \ b, \ c}\)

A combination of symbols like \(\(#\^2 &\)\(/@\)\)  can be hard to read. It is somewhat improved by extra spaces in the form #^2& /@ {a, b, c}, but in general we will avoid such combinations (although I am personally very fond of them) unless they force themselves on us.

m[#1,  #2]/n[#1,  #2] &[a,  b]

One can also operate on pure functions and get pure functions as the output. For instance, the derivative of a function f can be written as f'.

Sin'

Notice that Sin with no argument is a pure function and the output of Sin' is written explicitly as a pure function with a # and an &, although presumably, a simple Cos would have been sufficient. The following also work:

\({\(#\^2 &\)', \ \(\(\(\(\(#\^7 &\)'\)'\)'\)'\)'}\)

Pure functions written in this form with # and & are a distinctive and very attractive feature of the Mathematica programming language.

.  Four Kinds of Function Definition

Functions are such an important feature of Mathematica that they are represented in (at least) four different ways.

..  Expressions

An expression like x^2 + b x + c contains a symbol x that is intended to be interpreted (by us) as a variable; that means, we are intended to interpret the whole expression as describing a function of the variable x. However, it is only we who know this; there is no way for Mathematica to know it unless we somehow tell Mathematica what the variable is. Thus, in commands like

Solve[ x^2 + b x + c == 0, x], 
          Plot[Sin[Cos[x] + Tan[x]], {x, 0, Pi}],
          Sum[ (i^3), {i, 0, 10}] 
          Product[(x + i), {i, 1, 4}], etc.,

the first argument is such an expression and the second argument includes a description of the appropriate variable, either by just naming it or by including it as the first argument of an iterator.

2  x + 5  /.  x -> 2

..  Named pure functions

In a function definition such as f[x_] := (x^2), the thing being defined is f, which should be thought of as the "function in itself". The x in this definition is a dummy variable which is not really there at all. (In Chapter 7, Section 3.2, we'll see that x_ is a pattern.) If we insist on thinking of x as the variable in the definition of f, then it is a bound variable in the sense of logic. The f defined here is in fact a name for a pure function and can be used wherever pure functions are appropriate, just like the names of built-in functions. For example:

f[x_]  :=  (x^2);

Map[f, {a, b, c}]

Thus, defining a function using SetDelayed (i.e., :=) gives a name to a pure function. The definitions

square[x_] := (x^2),
	square = Function[{x}, (x^2)] and
	square = \(TraditionalForm\`#\^2\)&

are essentially equivalent. However, there are subtle differences described in Chapter 11, Section 7.

..  Nameless pure functions with bound variables

Expressions like Function[{x}, (x^2)] define  functions using a syntax that is essentially the same as the lambda calculus. Such a definition involves a bound variable, in this case x, whose name clearly doesn't matter; i.e., Function[{y}, (y^2)] describes the same function. The function itself has no other name attached to it, so it is a nameless pure function. Functions of more than one variable are allowed, for example

Function[{x, y}, x + y][2, 3]

but they need to be given the proper number of values as arguments (two here). This is to be distinguished from the following "curried" version, which is a function of one variable returning as value another function of one variable.

Function[{x}, Function[{y}, x + y]][2]

%[3]

..  Anonymous functions

The point of the syntax using # and & is that it is possible to construct a nameless pure function with no variables, bound or otherwise, i.e., an anonymous function. Anonymous pure functions are functions named by a canonical variable-free name, i.e., a # - & expression. When there is more than one slot, enough arguments have  to be given to fill all slots. "Currying" as above is not possible with anonymous functions, except by rather grotesque contortions.

Evaluate[(#1 + #2)&[2]]&/. #2 -> #1

%[3]

..  Conversion between forms of functions

Each of these notions is appropriate in its own place. It is possible to convert from one to the other.

Function[{x}, Evaluate[(\(TraditionalForm\`#\^2\) + 2#)&[x]]]

For instance, the anonymous function (\(TraditionalForm\`#\^2\) + 2#)& can be plotted using iv) a) above:

\(\(Plot[\(\((#\^2 + 2  #)\) &\)[x], \ {x, \ 0, \ 10}];\)\)

.  The fundamental dictum of Mathematica programming

The purpose of all of these operations based on Map is to make it possible to treat lists as a whole. For instance, one really poor way to square the entries in a list is as follows:

list = {a, b, c};

squares  =  Table[list[[i]]^2,  {i,  Length[list]}]

What this does is to tear apart the original list by extracting its parts one at a time, carrying out the squaring operation on each part, and then reassembling the new parts into a new list. Another way, which is only marginally better, is to use the typical Lisp style of recursive programming.

sqEls[{}] :=  {};
sqEls[list_]  :=  
      Join[{First[list]^2},  sqEls[Rest[list]]]

sqEls[{a,  b,  c}]

Here, we have ripped out the first entry and squared it. Note that the procedure also has to be initialized by saying what it does to the empty list. This kind of disassembly of lists is forbidden in Mathematica-style functional programming. In a generalized form, here is the fundamental dictum of Mathematica programming.

Treat mathematical structures as wholes.
Never tear them apart and rebuild them again.

Mathematica allows one to bypass both of the preceding strategies and think of a mathematical structure as a whole in the Platonic sense of mathematics. Mathematical structures are complete entities, even though they may be made up of parts or elements, and in many cases one can carry out constructions on these entities directly without having to be concerned with their parts.

Applying Functions to Values

We have already discussed and used the Mathematica facilities for defining functions and applying them to arguments. For instance, f[x_] := (x^2) defines the squaring function and f[2] or f@2 or 2//f applies it to the value 2. But there is more to practical functional programming than this. There are a number of built-in operations that take arbitrary functions as arguments and do something with them. For instance, there are several built-in commands that take a function as first argument and an arbitrary expression as second and apply the function to various parts of the expression. The first of these, Map, as we saw in Section 2.1, normally has two arguments, a function and a list, and it applies the function to every entry in the list, producing a list as the output. As usual in Mathematica, the list argument can be replaced by an arbitrary expression and levelspecs can be used to specify the level of the expression at which the function is applied.

.  Map and Its Relations

..  Map

The operation Map[function, expr] or in infix form function /@ expr applies function to each subexpression at level 1 in expr. The operation Map applied to a list, therefore, just applies the function to each entry in the list; e.g.,

Map[Sin, {a, b, c}]

Furthermore, Mathematica can map a function down the arguments of an arbitrary expression, not just a list. For instance, continuing with exp1 as in the preceding chapter,

exp1  =  (x^3) + (((1 + z))^2);

Map[Sin, exp1]

Be sure you understand the output here. Even more, Map takes a third argument, with a new effect. Map[function, expr, levelspec] applies function to the parts of expr described by levelspec. For instance, we can apply Sin to all of the leaves (described by the levelspec {-1}), or to all of the subexpressions whose depth is at least 2 (described by the levelspec -2).

Map[Sin, exp1, {-1}]

Map[Sin, exp1, -2]

This is a very powerful facility and is one of our main tools in manipulating expressions. The next  sections describe some variations on Map.

..  MapAll

MapAll[function, expr] or function //@ expr applies function to every proper subexpression of expr.  It is the recursive form of function /@ expr; i.e., it applies function to all of the arguments, then to all of the arguments of the arguments, etc.

MapAll[Sin, exp1]

This operation is almost the same as using Map with the levelspec Infinity or -1.

Map[Sin, exp1, Infinity]

The difference is that the levelspec Infinity does not include the whole expression itself.

..  MapAt

MapAt[function, expr, positionlist] applies function to the parts of expr described by the list of partspecs in positionlist. Here a partspec, as usual, is described by the list of edges in the tree form of the expression from the root to the given subtree. We'll use this to apply Sin just to the variables x and z in exp1 by giving a list of two partspecs.

MapAt[Sin, exp1, {{1,1}, {2,1,2}}]

See also MapThread and MapIndexed in The Mathematica Book.

..  Apply

The various versions of Map act on the arguments of a function. Apply acts only on the head of an expression. What Apply[head, expr] or head @@ expr does is to replace the head of expr by head (and possibly carries out a subsequent simplification or computation.) For instance:

Apply[Plus, {2, 3, 4}]

Here, the head of {2, 3, 4} is List, which does nothing to its arguments. When List is replaced by Plus, then something happens, since Plus calculates the sum of its arguments. Here is a slightly trickier example.

Apply[Plus, 2 3 x]

This time the head of 2 3 x is Times, which multiplies out as much of the expression as it can, yielding Times[6, x]. When this head is replaced by Plus, we get the indicated result.

Apply[Sin,  {{f[a],  f[b]},  {f[c],  f[d]}},  {2}]

Apply is frequently used if one wants to first prepare a number of ingredients and then apply some operation to them. The ingredients can be held in a list until they are ready and then the head of the list is changed to the appropriate operation by Apply.

..  Through

Once it is brought to your attention, Map seems like an obvious operation that ought to be available in any programming language. But it has an asymmetrical aspect: one function is applied to a list of values. What about applying a list of functions to a single value? That can be done too, using the operation Through. For instance:

Through[{Sin,  Cos,  Tan}[a]]

If the functions in the list are themselves listable, then they can be applied to a list of values.

Through[{Sin,  Cos,  Tan}[{a,  b,  c}]]

.  Thread, Inner and Outer

..  Thread

The built-in operation Thread does the same thing for an arbitrary head that the attribute of listability does for heads with this property.

Thread[ff[{2, 3, 4}, {a, b, c}, {x, y, z}]]

Actually, Thread works with either lists of the same length or individual arguments.

Thread[ff[{2, 3, 4}, {a, b, c}, 2, x]]

Furthermore, for Thread the arguments don't even have to be lists; they can also have an arbitrary head, which is then included as a second argument to Thread. (Note that Listable works only with lists.)

Thread[ff[hh[2, 3, 4], hh[a, b, c], hh[x, y, z]], hh]

Thread can be used, for instance, to construct lists used in substitutions. Consider the following fragment of Mathematica code.

var  =  {x,  y,  z};
point  =  {1,  2, 3};
expr  =  (x^2) y + 2  y z - 4  x (z^3);

Thread can be used to produce a list of substitutions.

Thread[Rule[var, point]]

This can then be used to substitute the point into the expression.

expr /. Thread[Rule[var, point]]

..  Outer

Outer was discussed briefly earlier. It basically does for functions of several variables what Map does for functions of one variable. Outer[function, lists ...] takes any number of lists as second through nth arguments and outputs the function applied to all choices of arguments from the lists arranged in a nested form to make this suitable for tensor computations. Thus, Outer of a function with two simple lists produces a matrix. Interestingly, this works with heads other than List. Note that the numbers of arguments in the lists do not have to be the same.

{Outer[Times, {a, b, c}, {1, 2}],
 Outer[Times, set[a, b, c], set[1, 2, 3]]}//MatrixForm

Outer with three simple lists produces what could be regarded as a list of three matrices.

Outer[Times, {a, b, c}, {1, 2, 3}, {u, v, w}]

Outer with two matrices as arguments produces something of depth 4.

tensor = Outer[Times, {{a, b}, {c, d}}, 
             {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}]

Transpose can be used to rearrange this in many ways.

Transpose[tensor, {1, 4, 2, 3}]

See also Distribute and Through.

..  Inner

There are two more complicated ways of applying operations to arguments consisting of lists. Inner is a generalization of Dot. The first example shows how to write Dot in terms of Inner, while the fourth example here shows that, in fact, Inner[f, list1, list2, g] is the same as Apply[g, Thread[f[list1, list2]]]. (See the next chapter for Apply.)  The fifth example shows that second and third arguments don't have to have head List. All that matters is that the heads be the same and that they have the same number of arguments, which are then extracted to be used by f and g.

{Inner[Times, {a, b, c}, {1, 2, 3}, Plus],
 Inner[Plus, {a, b, c}, {1, 2, 3}, Times],
 Inner[f, {a, b, c}, {1, 2, 3}, g],
 Apply[g, Thread[f[{a, b, c}, {1, 2, 3}]]],
 Inner[f, hello[a, b, c], hello[1, 2, 3], g]}//MatrixForm

..  Distribute

In particular, Distribute can sometimes be used to do the same things as Outer, without going inside deeper list structures. The syntax is somewhat different.

Distribute[{{a, b}, {c, d}}, List]

Distribute[{{{a, b}, {c, d}}, {{e, f}, {g, h}}}, List]

However, this fails to do the expected thing if the head of the first argument evaluates its arguments.

Distribute[Plus[{a, b}, {c, d}], List]

.  Nest and Fold

..  Nest, NestList, and FixedPoint

There are two more pairs of useful operations that fit the discussion here. They are common ingredients of functional programming languages. The first is Nest and its related operations. Nest[function, x, n] applies function to x and repeats it n times;  e.g., Nest[f, x, 3] returns f[f[f[x]]], while NestList[function, x, n] makes a list of these repeated operations a total of n + 1 times (since it starts from 0). For instance:

Nest[(# 2)&, a, 3]

NestList[f, a, 3]

NestList[(# 2)&, a, 3]

The following works because, as we have seen, D is Listable in its first argument.

NestList[D[#, y]&, r[y] == Sin[y] Cos[y], 4]//TableForm

Here is an example producing a simple continued fraction.

Nest[(1/(1 + #))&, x, 3]

An operation that is closely related to Nest is FixedPoint, which nests its operation until there is no change. For instance, everyone is familiar with what happens if the Cos key on a pocket calculator is pushed repeatedly. In principle, FixedPoint is what happens if it is pushed forever.

{Nest[Cos, 0.5, 6], 
 Nest[Cos, 0.5, 12], 
 FixedPoint[Cos, 0.5]}

Actually, FixedPoint stops after machine accuracy is achieved. Look up the options to FixedPoint to see how to change this. There is also an operation FixedPointList. The Practice section gives some examples. See Chapter 11, Section 6, for more serious uses of FixedPoint.

..  Fold and FoldList

The second pair of functions is Fold and FoldList, which do something similar to Nest, but for functions of two variables.

takes a function f of two variables, a starting seed value, and a list of subsequent values and returns

For instance:

Fold[f, a, {b, c, d}]

Similarly, FoldList produces a list of the successive values of this procedure.

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

Here are a couple of examples.

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

FoldList[Power, 2, {2, 3, 4, 5}]

It is usually easy to see when it is appropriate to use the function Nest; namely, there is some top level operation that is to be repeated a number of times. (This is sometimes called "tail" recursion, although in Mathematica it might be better to call it "head" recursion.) However, it is not so easy to see when it is appropriate to use Fold. What happens is that, at each repetition of the operation, new information is fed in from the list of values in the third argument. That is, instead of the third argument being a number saying how many times the operation is to be performed, it is a list of values to be used in repeating the operation. In the rest of the book we shall see a number of non-trivial uses of Fold, each of which is a triumph of human ingenuity.

.  Substitution

Function application in a functional programming language usually means substitution of a value for a variable. Thus, we expect that defining the squaring function by f[x_] := x^2 and then applying f to 2 should be the same as evaluating the substitution x^2 /. x -> 2. Of course it is, but as will be seen in Chapter 11, Section 4, this form of substitution sometimes doesn't work correctly. Furthermore, as will be discussed in the next chapter, this form of substitution is really an application of a local rewrite rule and should not be thought of as a substitution at all. However, there is another operation in Mathematica that exactly implements the idea of substitution in functional programming languages; namely, With. For instance,

With[{x  =  2},  (x^2)]

In many functional languages, this would be written in the form let x = 2 in (x^2). Instead of giving the value where the function is to be applied in the first argument of With, one can also specify the function in this  position.

\(With[{square\  = \ #\^2 &}, \ square[2]]\)

In the theory of functional programming languages, based on the lambda calculus, an expression of the form let x = 2 in (x^2), as above, is synonymous with applying the pure function λx.(x^2) to the value 2. Since applying pure functions to values is the only thing that is done in functional programming, such programs consist mainly of let expressions. This style of programming can be adopted in Mathematica and often has attractive results. For instance, consider the following method to calculate improper integrals with possible singularities at the end points.

improperIntegrate[expr_, {x_, a_, b_}] := 
   With[{integral = Integrate[expr, x]},
     Limit[integral, x -> b, Direction -> 1] -
     Limit[integral, x -> a, Direction -> -1]]

This works nicely on typical examples.

improperIntegrate[
      (1/(((2  x - 1))^((2/3)))),  {x,  (1/2),  2}]

Functional Programs

Functional programs are nested sequences of "button pushes"; i.e., they are single expressions made up solely from built-in commands and built-in constants. Sometimes these are called "one-liners". It is possible to do many intricate operations just using such one-liners. There is a column devoted just to them in the Mathematica Journal. The basic rule for a "strict" one-liner, as it will be called here, is that the only ingredients allowed on the right-hand side are built-in operations and constants or argument names that occur in the left-hand side. This rules out nearly all Table constructions, since they require a (bound) variable in the iterator argument which does not occur in the left-hand side. It also rules out use of Integrate and Solve unless they are parts of function definitions which include the variable specification on the left-hand side. It rules out expressions of the form Function[{x}, body], since that also includes the bound variable x, but anonymous pure functions do exactly the same things without introducing any bound variables, and this is what makes it possible to construct strict one-liners. Non-strict, or ordinary, one-liners have no such restrictions. The only condition for them is that, in theory at least, they should be written on one, possibly very long, line.

.  Simple Examples of Functional Programs

Some of the things we did interactively or in more than one step in the first three chapters can be put together to make simple functional programs. For instance, the interactive sequence of operations

Integrate[(x/(1 - (x^3))),  x]
D[%,  x]
Simplify[%]

can be put together into a single nested operation.

Simplify[D[Integrate[(x/(1 - (x^3))),  x],  x]]

This violates the rules for a strict one-liner since it involves the bound variable x, but otherwise it is just a nested sequence of built-in commands. Here are the results of the same process applied to a number of other interactive constructions from Chapters 1 and 3. Note that we have replaced Table constructions by mapping a pure function down an index list constructed by Range or by using listability.

ListPlot[ N[ Log[ Map[ #!&, Range[20] ] ] ] ];

Fit[N[Log[Map[#!&, Range[20]]]], {1, x, x^2}, x]

Consider the construction N[Log[Map[#!&, Range[20]]]] here. There are, in fact, at least three ways to construct this list of numbers. First, just build a table as we did in Chapter 1.

Table[N[Log[n]], {n, 20}]

Second, map a pure function (the factorial function) down the list of desired numbers, constructed by the Range operation.

N[Log[Map[#!&, Range[20]]]];

Third, use the fact that !, Log, and N are Listable to get the result from a very brief command.

N[Log[Range[20]!]];

In the next three examples, we have to use Map because ToString and Rationalize are not Listable.

Map[{ToExpression[ToString[N[Pi, #]]], 
     N[Sin[ToExpression[ToString[N[Pi, #]]]], 11] }&, 
     Range[5]] // TableForm

\(Map[Rationalize[N[π], \ 0.1\^#] &, \ Range[10]]\)

\(Union[Map[Rationalize[N[π], \ 0.1\^#] &, \ Range[20]]]\)

This way of using Mathematica constitutes functional programming in Mathematica. It views the basic entities of Mathematica as functions (possibly of many variables) and the basic operation as composition of functions, or rather, the iterated application of functions to arguments. One-liners either can just carry out some specific calculation or can be used as definitions of functions whose arguments can then be given values to do something interesting.

.  Developing a Functional Program

The large stock of built-in functions make it possible to solve rather intricate problems in a straightforward way. For instance, in Section 3.1.3 above, starting with the expression

exp1  =  (x^3) + (((1 + z))^2);

we applied Sin to the variables x and z using MapAt.

MapAt[Sin, exp1, {{1,1}, {2,1,2}}]

You may have wondered how one would know what partspecs to give without carefully analysing the expression. But Mathematica will do this for you by itself, using Position.

Position[exp1, x]

Thus, the following one-liner does it all.

MapAt[Sin, 
      exp1, 
      Join[Position[exp1, x], Position[exp1, z]]]

However, this violates the strict rule by referring to the variables x and z. The idea here can be developed farther, by having Mathematica do more of the work. Perhaps you feel that Mathematica should also find the variables without our having to tell it what they are and that would turn this into a strict one-liner. In this kind of expression, the variables are particular leaves in the tree form of the expression. We would get all the leaves by the following:

Level[exp1, {-1}]

We want to select x and z from this list. The predicate Not[NumberQ[#]]& is True just for them.

Select[Level[exp1, {-1}], Not[NumberQ[#]]&]

Next we need the position in exp1 of each entry of this list, which we find by mapping the pure function Position[exp1, #]& down this list.

Map[Position[exp1, #]&, %]

This has too many brackets, but Flatten with a levelspec will get rid of them.

Flatten[%, 1]

Now we can put this all together, using the mouse to replace each % by its construction in the previous line, which yields a one-liner function definition that will take any expression expr (instead of exp1) and apply some given function fun (instead of Sin) just to the variables in it.

mapVarsOnly[fun_, expr_] :=
	MapAt[fun, expr,
		  Flatten[
            Map[Position[expr, #]&, 
	 			Select[Level[expr, {-1}], 
	 	               Not[NumberQ[#]]&]], 1]]

This is a true one-liner since the only ingredients on the right-hand side are built-in functions and constants, together with expr and fun, which occur on the left-hand side. Try this out with Sin and exp1 to check that it works as expected.

mapVarsOnly[Sin, exp1]

Now, try some other examples.

\(mapVarsOnly[\@# &, \ exp1]\)

mapVarsOnly[
      ArcTan,  ((((x - (y^2) + 3))^w)/(Sqrt[(u^3) + 3  v]))]

As an example of a one-liner, this is OK, but in fact it would fail on an expression that has what we would regard as a symbolic constant, e.g., a, because it would treat that as a variable, too. In the exercises we ask you to fix the definition so that mapVarsOnly only treats letters between p and z as variables.

.  Frequencies

List manipulations are important in functional programming. Suppose we want to write a function that takes a list as its argument and returns a list of the number of times each entry occurs in the original list. Start with a list to use as an example.

list = 
   {a, d, s, f, d, a, s, a, d, f, d, f, g, d, a, f, g};

Union will give us the "set" of distinct entries in this list, written in canonical Mathematica order.

Union[list]

Our problem is to determine how many times each entry here occurs in list. There is a built-in function that will do that for a single entry.

Count[list, f]

We don't just want the number 4, but we want it associated with the symbol f so we know what it means. Thus, we want the pair {f, 4} as output. This is easily constructed by a pure function.

{#, Count[list, #]}&[f]

So all we have to do is put these together correctly in order to design a function that gives each distinct element and the number of times it occurs.

frequencies[list_] := 
    Map[{#, Count[list, #]}&, Union[list] ]

Try this out on our example.

frequencies[list]

This is read as saying that a occurs 4 times, d 5 times, etc. We will use frequencies later in constructing our own BarChart graphics function.

.  Newton's Method

..  One variable

Newton's method is a procedure for finding a zero of a function. There is of course the built-in function FindRoot, but we want to construct our own version to see how it works. Given an expression representing a function of x, e.g., expr = (x^2) - 3, and some starting value x0, then the procedure calculates a sequence of values given by the iterative formula

(Subscript[x, 0]) =  x0;
(Subscript[x, n + 1]) = x - (expr/D)[expr,  x] /. x -> (Subscript[x, n])

This is repeated until the results change by less than some specified error. We will use FixedPoint and let Mathematica decide when to stop. What we have to do is to take the right-hand side of the iterative formula and turn it into a pure function, which is then repeatedly applied to x0. It is clearest to do this as a separate function called oneNewtonStep as follows:

\(oneNewtonStep[
        expr_, \ {x_, \ x0_}]\  := 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x - 
          expr\/D[expr, \ x] /. x -> N[x0]\)

To force the computation to be done numerically rather than exactly, we use N[x0] instead of x0 as the last argument here. The value of oneNewtonStep at stage n is what is to be used as the starting point for stage n + 1. That means we want to consider it as a pure function of the initial point x0. We can either use Nest some given number of times, or let Mathematica decide how often to iterate this procedure by using FixedPoint.

newton[expr_, {x_, x0_}] :=
  FixedPoint[oneNewtonStep[expr, {x, #}]&, N[x0]]

Here are a couple of examples.

newton[x^2 - 3, {x, 1.0}]

newton[x - Cos[x], {x, 0.5}]

..  Several variables

Essentially the same formula works for n functions of n variables. Newton's method then finds values of all the variables so that all of the functions are zero. We just imagine that x means an n-dimensional vector and expr means n functions of n variables. The derivative D becomes the Jacobian matrix, and division means multiplication in the sense of the Dot product by the inverse. The formula  becomes

(Subscript[x, n + 1]) = (x - Inverse[jacobian[expr,  x]] . expr[x] /. 
          x -> (Subscript[x, n])

Recall that the Jacobian is given by the operation

jacobian[exprs_, vars_] := N[Outer[D, exprs, vars]]

Our intention now is that exprs is to be a list of expressions and vars a list of variables. The initial value will be a list vars0 of values. We have to change the notation slightly so that newton won't become confused about being given one function or a list of functions. (For a better way to handle this, see the next chapter.)

\(oneNewtonStep[exprs_, \ vars_, \ 
        vars0_] := 
	\((vars - Inverse[jacobian[exprs, \ vars]] . exprs)\) /. 
        Thread[vars -> N[vars0]]\)

In the single-variable case, we just said x -> (Subscript[x, 0]), but now we need a list of rules. Thread does exactly the right thing; e.g.,

Thread[{x,  y,  z} -> {1,  2,  3}]

The final function is almost the same as the one variable case.

newton[exprs_, vars_, vars0_] :=
  FixedPoint[oneNewtonStep[exprs, vars, #]&, N[vars0]]

Here is a simple example.

exprs1 = {x^2 + y^2 - 13, x^3 - y^3 - 19};

newton[exprs1, {x, y}, {2, 1}]

In the exercises, you are asked to restructure this program so that the answer is a list of substitutions.

Practice

1.  MapThread[Rule, {x, y, z}, {1, 2, 3}]

2.  MapIndexed[Nest[Sin, #1, Sequence@@#2]&, {a, b, c}]

3.    \(#\^2\)&[anything]

4.    #&[anything]

5.  1&[anything]

6.  something&[anything]

7.    polys = Table[1 - (x^n),  {n,  10}];

8.  Factor /@ polys // ColumnForm

9.    \(Expand[#\^2]\)& /@ polys // ColumnForm

10. FixedPointList[Cos, .5]

11.  \(FixedPointList[
        Cos, \  .5, \ 
		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ SameTest\  -> \((Abs[#1 - #2]\  < \ 10\^\(-6\) &)\)]\)

The following are taken from the One-Liners column of the Mathematica Journal, Vol. 1, 1991. Try to understand what they do and how they work. In some cases, minor or major changes have been made to comply with the canon that one-liners should not introduce any extraneous variables on their right-hand sides.

12.	rootPlot[poly_, z_] :=
      ListPlot[{Re[z], Im[z]}/. Solve[N[poly == 0], z],
               Prolog -> PointSize[0.04]]

13.  poly[n_,  z_] :=
	          (z^n) - 
        Plus @@ (Power[z, #] &  /@  Range[0,  n - 1])

14.	rootPlot[poly[20, z], z]

15.  newtonRoot[f_,  x0_] := FixedPoint[(# - f[#]/(f')[#]) &,  x0]

16.	newtonRoot[(# - Cos[#])&, 0.5]

17.  \(\(ListDensityPlot\)\([\)\(
\)\(		\ \ \ \ \ \ \ \)\(Outer[
          If[IntegerQ[\@\(#1\^2 + #2\^2\), \ 0, \ 
                1] &, 
				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Range[50], \ Range[50]]]\)\)

18.	phasePlot[f_, {x_, xmin_, xmax_}] :=
       ParametricPlot[
           Evaluate[{f, D[f, x]}], {x, xmin, xmax}]

19.  phasePlot[Sin[(x^2)],  {x,  0,  2  π}]

20.	reverseInteger[n_] :=
      Dot[Power[10, #]& /@ Range[0, Floor[N[Log[10, n]]]], 
          IntegerDigits[n] ]

21.	reverseInteger[123456789]

Exercises

Observe the fundamental dictum of Mathematica programming in working these exercises.

displayed in a nice form. Hint: try Fold.

power1[x_, n_, base_] :=
	Fold[(#1^2 #2)&, 1, x^IntegerDigits[n, base]]

8	   3	   0	     0
			                 3	   2	   1	     2
			                 0	   1	   2	     2
		                        	   0	   2            2	   14

\(\(\(				\)\(legendre \((f, \ g)\)\  = \ ∫\_\(-1\)\%1 
          f \((x)\) g \((x)\) \[DifferentialD]x\)\)\)