MATHEMATICA PROGRAMMING
5 Feb 2003

Rule Based Programming

Recall a definition for the factorial function:

fac[n_Integer]  :=  Infinity  /;  n  <  0

fac[-1]

fac[0]  =  1

fac[n_Integer]  :=  n * fac[n - 1]

In Mathematica, when multiple rules are defined, the more specific rule is applied first, and general rules are tried after the more specific ones.  ?fac lists in the order of rules.

Definition for the greatest common divisor:  the greatest common divisor d is the largest number that divides both a and b.   The definition below is based on a theorem in number theory which says, if for integer a,b,c such that |c| < b, and 
                      a = k b + c,
then gcd[a,b] = gcd[b,c].  Doing this (division) recursively until remainder c becomes zero. This is called Euclidean algorithm.

gcd[a_,  0]  :=  a
gcd[a_,  b_] :=  gcd[b,  Mod[a, b]]

gcd[3, 9]

What about gcd[...] with one, three or more arguments?

Define our own differentiation algorithm

Let's agree that we do only partial differentiation, so that dif[expr, x] will be zero unless expr contains x explicitly.  To expression this fact, we can say:

dif[f_,  x_]  :=  0  /;  FreeQ[f, x]

dif[a, x]

The differentiation operation is a linear operator.  We program in the linearity of differentiation:

dif[ f_  +  g_,  x_]  :=  dif[f, x]  +  dif[g, x]

dif[ a_ * g_,  x_]  :=  
      a * dif[g, x]  /;  FreeQ[a, x]

We know the differential of polynomials (actually x^n).

dif[x, x]  :=  1

dif[x_^n_,  x_]  :=  
      n*x^(n - 1)  /;  (FreeQ[n, x])

After these definitions, dif[ ] can work with any polynomials.

Rule-Based Programming: operator overloading

Suppose we like to use the notation  S ∨ T to denote union of two sets S and T, and S ∧ T to mean the intersection of the two sets.  Since the two special symbol ∧ and ∨ is used for logical operations on the truth values True and False, can we use them for set operation?   We'll use a list to represent a set.

{a}  ∧  {a, b, c}

Since the operation is not defined, we get back the input.  Here is a definition that "overload" the operator ∧ and ∨.

Unprotect[And,  Or]

List[x___]  ∧  List[y___]  :=  
      Intersection[List[x], List[y]]

{a, b, c}  ∧  {a, e}

List[x___]  ∨  List[y___]  :=  
      Union[List[x],  List[y]]

{a, b, c}  ||  {a, e, r}

{a,  b}  &&  {}  ||  {c, d}

Functional Programming Style

In functional programming style we apply functions to expression rather than through control of execution flow to achieve the computational result. [That is, in functional programing style, we apply functions that do not have side effect, and we do not use explicit loops.]

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

@@ is the Apply[ ] operator; Times @@ expr is short form for Apply[Times, expr].

Range[4]

Apply[Times,  {1, 2, 3, 4}]

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

Observe that:

Times  @@  (a + b + c)

Map[f, expr] is another function of similar style.

Map[f,  h[a, b, c]]

Map[Sin,  {1,  Pi/2,  2  Pi}]

Concept of Pure Function in Mathematica

The first argument in Apply[ ], or Map[ ] requires just the name of the function, e.g., Sin, not Sin[x].  We call such objects pure functions.  You can think of a pure function as a mathematical operator, ready to act on an expression.   Here is a general way of defining a pure function.  Note that Function[x,x^2] is equivalent to #^2 &.

Function[x,  x^2]

Map[ #^2   &,  {1, 2, 3}]

A standard way of writing a pure function without given the function a name is to use # (or #1, #2, etc ) as the argument of the function, construct arbitrary expression, and end the expression with &, such as  #*Sin[#] &.

Let us look at a program that finds zero of a function, f(x)=0, by applying FixedPoint[]; this is an example of functional programming style.  First the built-in function:

?FixedPoint

We use Newton iteration method to find the solution to f(x) = 0 numerically, which is based on a Taylor series expansion to first order:
                   f(x) = f((Subscript[x, 0]))+ f' ((Subscript[x, 0])) (x - (Subscript[x, 0])) + O((((x - (Subscript[x, 0])))^2))
 If f(x)  =  0, then we can solve for x as
 \(TraditionalForm\`\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(x\  = \ 
        x\_0\  - \ \(f(x\_0)\)\/\(f' \((x\_0)\)\)\)\)\),
 thus the solution is a fixed point of the function g(x)  = x  -  f(x)/f' (x),  
      i . e . ,  the zero of f has the property g(x)  =  x.  FixedPoint[] function gets the solution by iterating (Subscript[x, n + 1])= g((Subscript[x, n])),  
                n  ->  ∞.

newtonzero[f_,  x0_]  :=  
      Module[ {x,  fp = f'},  
x  =  
          FixedPoint[(# - f[#]/fp[#]) &,  
            x0]; 
x
]

newtonzero[(Cos[#] - #) &,  0.0]

FindRoot[ Cos[x] == x,  {x,  0}]

We can give newtonzero another definition so that it accept normal function (not pure function) as input.

newtonzero[expr_,  x_Symbol,  x0_Real]  := 
      newtonzero[Function[x, expr],  x0]

newtonzero[Cos[x]  -  x,  x,  0.3]

Here is the way usage message is entered for newtonzero, you get the message when enter ?newtonzero:

newtonzero::usage  = "newtonzero[f,x0] finds f[x]=0, with trial \
starting value x0.  f must be a pure function while x0 should be Real.  \
newtonzero[expr, x, x0] finds x such that expr=0, where expr is an expression \
in the variable x (a symbol), and x0 is the starting value (a Real). N.B. the \
iterations need not always converge."

More single-line Mathematica definitions:
Vector calculus

Grad[f_,  arg_]  :=  D[f, #] &  /@  arg

Div[v_,  arg_]  :=  Inner[D,  v,  arg,  Plus]

JocobianMatrix[f_,  v_]  :=  Outer[D,  f,  v]

TripleProduct[a_,  b_,  c_]  :=  Det[{a, b, c}]