Principia Matematicae

Dr. C. Chryssomalakos, ICN - UNAM

This is a set of notes on some general aspects of Mathematica. Some of the material is scattered in the references cited in the bibliography, mostly Refs. [2], [3], [5] ---the rest are things unearthed while struggling with the subject.

SetDirectory["E:\chryss\math\workshop"];

Functions

There are four basic types of function definition in Mathematica.

Expressions

Define a function f via

f = x^2 + 3  x + 2

Evaluate f at x=2

f /.  x -> 2

The /. means to substitute according to the rules that follow, in this case x-> 2.

FullForm[f /.  x -> 2]

FullForm[Hold[f /.  x -> 2]]

/. is infix of ReplaceAll, -> is infix of Rule

Named pure functions

This defines h[anything]=anything^2

h[x_] := x^2

Notice the use of :=. The rhs is evaluated when the function is called, not when defined.

Compare

h2[x_] := Expand[(x + a)^2]

h3[x_] = Expand[(x + a)^2]

h2[b + c]

h3[b + c]

Consider another example

h4[f_, x_] := D[f, x]

h4[Sin[x], x]

Clear[f]

h5[f_, x_] = D[f, x]

?? f

?? e3

h5[Sin[x], x]

Nameless pure functions with bound variables

You can map a function down a list with Map

Map[Tan, {x, 1, a + b}]

You can define the function without giving it a name

Map[
Function[{t}, Sin[t^2]], 
{x, 
        1, a + b}
]

Define f1 as a pure function with one bound variable

f1 = Function[{s}, Cos[(s + 1)^2]]

f1[2]

D[f1[r], r]

Plot[
f1[z], 
{z, 0, 
        10}
]

Integrate[f1[x], x]

General form is Function[{(Subscript[x, 1]), (Subscript[x, 2]),  ... , (Subscript[x, n])}, body].

Anonymous pure functions

A pure function can be defined without the use of a bound variable

f2 = #^2 &

f2[2]

The # denotes the argument, the & the end of the definition.

Map[
Sin[#^2] &, 
{a, b, 
        c + d}
]

Comment

Functions can be iterated with Nest

g = Nest[Sin, x, 100]

Plot[
Evaluate[
        Nest[Sin, x, 5]], 
{x, 0, 
        10  Pi}
]

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

Plot[Nest[1/(1 + #) &, x, 17], {x, - .9, 10}, 
      PlotRange -> All]

Anonymous pure functions can have any number of arguments

m[#1, #2]/f[#1 #2] &[x, y]

Conversion between forms of functions

Starting with an expression

expr5 = x + Log[x]

get a named pure function,

f4[x_] = expr5

f4[y]

?? f4

Compare with SetDelayed

f5[x_] := expr5

f5[y]

?? f5

You can also get a pure function with a bound variable

Function[{x}, Evaluate[expr5]][a]

Evaluate forces evaluation of the body of Function. Compare with

Function[{x}, expr5][a]

(Notice: Wrong in "Mastering Mathematica"!)

Finally, you can get an anonymous function

Evaluate[expr5 /. x ->  #] &

Evaluate[expr5 /. x ->  #] &[a]

The order of evaluation is:
1) Evaluate the body
    a) Evaluate expr5: x+Log[x] /. x--> #&[a]
    b) Effect the substitution: # +Log[#]&[a]
2) Substitute the argument for #: a+Log[a]

Compare with

expr5 /. x ->  # &

expr5 /. x ->  # &[a]

(Notice: Wrong in "Mastering Mathematica"!)

The order of evaluation is:
1) Substitute the argument for #: expr5 /. x--> a
2) Evaluate expr5: x+Log[x] /. x--> a
3) Effect the substitution: a+Log[a]

Starting with a named pure function

f[x_] := x + Log[x]

get an expression

expr6 = f[x]

a pure function with a bound variable

Function[{x}, f[x]]

Function[{x}, f[x]][a]

(compare with

Function[{x}, Evaluate[f[x]]][a]

try to trace it...)

and finally, an anonymous function

f[x] /. x -> # &

f[x] /. x -> # &[a]

(compare with

Evaluate[f[x] /. x -> #] &

Evaluate[f[x] /. x -> #] &[a]

trace it...)

Starting with a pure function with a bound variable

f6 = Function[{x}, x + Log[x]]

you can get an expression,

f6[x]

a named pure function,

f6

or an anonymous pure function

f6[x] /. x ->  # &

f6[x] /. x ->  # &[a]

(compare with

Evaluate[f6[x] /. x ->  #] &

Evaluate[f6[x] /. x ->  #] &[a]

)

Exercise: anonymous--> rest

Functional programming

Example 1: Frequencies

Consider a list

list1 = {a, a, c, b, a, c, b, a, c}

Problem: Given list1,construct {{a, (Subscript[n, a])}, {b, (Subscript[n, b])}, {c, (Subscript[n, c])}},
where (Subscript[n, i]) 
        counts the appearances of i.

First solution

list1

Sort[list1]

{#, 1} &[a]

/@ is the infix form of Map

{#, 1} &  /@  Sort[list1]

Try merging neighbors

{#, 1} & /@ 
        Sort[list1] /. {f_, {x_, m_}, {x_, n_}, g_} ->  {f, {x, m + n}, 
          g}

Does not work because rule in rhs only applies to lists with four elements (f_ matches any single expression).

 f__  (two underscores) matches one or more expressions

{#, 1} & /@ 
        Sort[list1] /. {f__, {x_, m_}, {x_, n_}, 
          g__} ->  {f, {x, m + n}, g}

/. only applies the substitution once (but tricky... mostly wrong in books). 
//. keeps applying it until it can no more (dangerous!, control with MaxIterations).

{#, 1} & /@ 
        Sort[list1] //. {f__, {x_, m_}, {x_, n_}, 
          g__} ->  {f, {x, m + n}, g}

First and last entry have only one neighbor. f___ (three underscores) matches zero or more expressions.

{#, 1} & /@ 
        Sort[list1] //. {f___, {x_, m_}, {x_, n_}, 
          g___} ->  {f, {x, m + n}, g}

Putting everything together

frec[
        ls_] := {#, 1} & /@ 
          Sort[ls] //. {f___, {x_, m_}, {x_, n_}, 
            g___} ->  {f, {x, m + n}, g}

frec[list1]

Try frec on another list

list2 = {x, y, z, w, x, w, z, y, y, w, x, z, w, y, x, w}

frec[list2]

Second solution

Union[list1] returns a list of distinct elements of list1

Union[list1]

\(Count[list, \ 
        element]\ returns\ the\ number\ of\ appearences\ of\ element\ in\ \(\(list\)\(.\)\)\)

Count[list1, a]

Construct an anonymous function

{#, Count[list1, #]} &

{#, Count[list1, #]} &[a]

{#, Count[list1, #]} &[b]

Map[f,{x,y,z,...}] returns {f[x],f[y],f[z],...}

Map[
{#, Count[list1, #]} &, 
Union[
        list1]
]

Put everything together

frequencies[
        ls1_] := 
Map[
{#, 
            Count[ls1, #]} &, 
Union[
          ls1]
]

frequencies works with any argument

frequencies[list2]

Comparison of the two solutions

Generate some long sequences

ListPlot[
      frequencies[
        IntegerDigits[10000!]], 
PlotJoined -> 
        True, 
AxesOrigin -> {0, 
          0}, 
ImageSize ->  
        500
]

Length[IntegerDigits[10000!]]

ListPlot[
      frequencies[
        IntegerDigits[
          Floor[10^20000 Pi]]], 
PlotJoined -> 
        True, 
ImageSize ->  
        300
]

frequencies[IntegerDigits[Floor[10^10000 Sqrt[2]]]]

Try Timing frec

Timing[frec[IntegerDigits[5000!]]]

Length[IntegerDigits[5000!]]

Compare with frequencies

Timing[
frequencies[
        IntegerDigits[20000!]]
]

Example 2: Reverse Integer

Write a function that reverses the digits of its integer argument.

rI[n_] := 
Dot[
Map[
Power[10, #] &, 
Range[0, 
            Floor[N[Log[10, 
                  n]]]]
], 
IntegerDigits[
          n]
]

rI[1897]

Exercise

What does power1 do?

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

Example 3: Flowing Cats

Read and Plot Data

\(\(data = 
        ReadList[
"image.dat", 
Number, 
RecordLists -> 
            True
];\)\)

ListDensityPlot[
data, 
Mesh -> False, 
FrameTicks -> 
          None
];

FullForm[data]

ListPlot[Sort[Flatten[data]]];

ListPlot3D[
data, 
Lighting -> False, 
ColorFunction -> 
          Automatic, 
Mesh -> 
          False, 
ViewPoint  ->  {10, 0, 
            5}
];

ListPlot3D[
Abs[
          Fourier[data]], 
Lighting -> 
          False, 
ColorFunction -> 
          Hue, 
Mesh -> 
          False, 
ViewPoint  ->  {0.2, -2, 
            3}
];

Kitten flow on a torus

Define a hamiltonian

ham = -x^2 + p^2;

xpsub = {x -> x[t], p -> p[t]};

Solve Hamilton's equations... a more refined procedure would be needed to handle critical points of the potential --- as is it works fine for -x^2 but not for, e.g., (x^2-a^2)^2.

sol[{x0_, p0_}] := 
{x, 
          p} /. NDSolve[
{
x'[
                  t] == (D[ham, p] /. 
                    xpsub), 
p'[
                  t] == -(D[ham, x] /. 
                      xpsub), 
x[0] == 
                x0, 
p[0] == 
                p0
}, 
{x, 
              p}, 
{t, -1, 
              0}
][[1]]

xpRange stores the dimensions of the image

xpRange = Length /@ {data[[1]], data}

Find the flow {x0, p0, t} on a torus: intra[{x0, p0}, t0] is the position of the initial (t=0) point {x0, p0} at time t0.

For a function f of two variables, 
Thread[f[{a, b, c, ...}, {x, y, z, ...}]] 
gives 
{f[a, x], f[b, y], f[c, z], ...}, 
with similar action on a function of n variables. Below it is used to put the image on a torus via Mod. Through[{f, g, h, ...}[x]] gives
{f[x], g[x], h[x], ...} (but it is a more general command --- check the Help Browser).

intra[{x0_, p0_}, 
        time_] := 
Thread[
(Mod[#1, #2] + 1) &[
          Floor /@ Through[sol[{x0, p0}][
                time]], 
xpRange
]
]

Off[General::spell, General::spell1]

datap[t] is the image at time t (check MapIndexed and Apply (with infix form @@) in the Help Browser).

datap[t_] := 
      MapIndexed[
(Part @@ 
Prepend[
Reverse[
                  intra[Reverse[#2], -t]], 
data
]) &, 
data, 
{2}
]

Table[
ListDensityPlot[
datap[ .025 n], 
Mesh -> 
            False, 
FrameTicks -> 
            None
], 
{n, 0, 
          20}
];

Bibliography

[1] S. Wolfram, "The Mathematica Book", Cambridge University Press (1999), fourth edition.

[2] J. W. Gray, "Mastering Mathematica", Academic Press (1997).

[3] H. Ruskeepaa, "Mathematica Navigator", Academic Press (1998).

[4] R. Maeder, "Programming in Mathematica", Addison Wesley (1997), third edition.

[5] T. W. Gray, J. Glynn, "Exploring Mathematics with Mathematica", Addison Wesley (1991)

[6] A. Gray, "Modern Differential Geometry of Curves and Surfaces with Mathematica", CRC Press (1998).