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"]

Lecture 3: 

Patterns, 
Substitutions, 
Evaluation

Patterns and Substitutions

Wildcards

Patterns represent classes of expressions - they are defined in terms of elementary patterns, called wildcards
_                            any expression
x_                          any expression, given the name x
x_^n_                   any expression (x) raised to any power (n)
a_ + b_              sum of any two expressions
f[x_]                   the function f with any argument
f_[x_]                any function with any argument
f_[x_, y_]      any function with any two arguments
f_[x_, x_]      any function with two identical arguments

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

Some examples

f[x] /. f[_] -> 5

The underscore matched x in the above example. Notice that since the pattern is used nameless, there is no way to refer to what it matches in the rhs of the substitution rule. Compare with

a /.  x_  ->  Sin[x]

Here x_ matched a. By calling the pattern x_, one can use it in the argument of Sin.

x^2 /. x_ -> a

Notice that x_ stands for any (single) expression. Mathematica checks for matches starting with the expression as a whole and then checking its subparts, that's why x_ matched the whole expression above, and not just x or 2.

FullForm[x -> a]

So, -> is the infix form of Rule.  What is the internal representation of the pattern x_?

{FullForm[_], FullForm[x_]}

So, _ is Blank[]. The empty brackets are there because _ can take an argument (see below).

You can use the patterns x_ and n_ on the rhs of a substitution rule (without the underscores)

b^2 /.  x_^n_ ->  {x, n^x}

Pattern matching is purely typographical (well, almost) - Mathematica will not do math to match an expression. In the example below, Mathematica does not expand Log[x y]=Log[x]+Log[y].

Log[x y] /.  a_ + b_ ->  a b

Substitution rules can be given names,

sub = a_ x_  ->  {a,  x}

3 x /. sub

Notice that sub does not match x (thought of as 1 x), because x does not have the structure of the pattern in sub

x /. sub

To see this clearly, look at the FullForm of a_ x_ and x

FullForm[a_ x_]

FullForm[x]

x does not explicitly involve Times so there is no match (but there is a way around this, see below).

Find the complex conjugate of 3+I2

Conjugate[3 + I 2]

The same, with a+Ib

Conjugate[a + I b]

Mathematica does not know if a and b are real, so it prefers to leave the expression as is (instead of expanding it in terms of Re[a], Im[a] etc.). Try I -> -I

a + I b /.  I -> -I

a - I ba + I b /.  I -> -I

Notice

FullForm[3 + I 2]

You can let Mathematica know that a, b are reals

Refine[Conjugate[a + I b], 
      a  ∈  Reals  &&  b  ∈  Reals]

Exercise

What do x_ and y_ match in the following?

a + I b /. Complex[x_, y_] -> Complex[x, -y]

Solution

FullForm[a + I b]

(since a, b are not known to be real), so, x matches 0 and y matches 1.

Often, you need to know the FullForm to make /. work.

1/x /. a_/b_ -> a + b

FullForm[1/x]

FullForm[a/b]

However

1/x /. x^n_ -> {x, n}

since

FullForm[x^n]

which matches that of 1/x.

Mathematica does a minimum of math in pattern matching, e.g., it uses commutativity and associativity of addition and multiplication

f[b + a] /. f[a + x_] -> a f[x]

Constrained Patterns

Head constraints

_head only matches expressions with Head  head

{2, 3.4, 5 + 3  I} /. x_Integer  ->  x^2

Check the heads with Map

Head /@ {2, 3.4, 5 + 3  I}

In _head, head is the argument of _

FullForm[_Integer]

Evaluation

Most built-in functions in Mathematica Map automatically, e.g., Sin[{a,b,c}]={Sin[a],Sin[b],Sin[c]}
and {x,y,z}^2={x^2,y^2,z^2}.
Apply[head,expr], with infix form head@@expr, replaces the head of expr with head.

Define the euclidean norm of a vector

Clear[norm]

norm[list_List] := Plus @@ (list^2)

The definition only applies when the argument is a list.

{norm[{x, y, z}], norm[3.5], norm[200!]}

?norm

norm with any other type of argument is returned unevaluated, unless a definition for that type of argument also exists. Define norm for reals

norm[x_Real] = "Hello"

The two definitions coexist, and Mathematica chooses which one to use according to the type of the argument

{norm[{x, y, z}], norm[3.5], norm[200!]}

Now the first two calls to norm are evaluated, while in the third, only the 200! can be evaluated.

Define the norm of an integer to be a plot of its prime decomposition. For an integer n, FactorInteger[n] returns its prime decomposition in the form 
{{(Subscript[p, 1]), (Subscript[r, 1])}, ...,{(Subscript[p, k]), (Subscript[r, k])}}, where n =(((Subscript[p, 1]))^( (Subscript[r, 1])))...(Subscript[((Subscript[p, k])^r), k]). ListPlot then plots the points {(Subscript[p, i]), (Subscript[r, i])}.

FactorInteger[6!]

norm[n_Integer] := ListPlot[FactorInteger[n]]

Now all three calls are evaluated

{norm[{x, y, z}], norm[3.5], norm[200!]}

Aside (warming up for graphics):

What is the internal representation of a point plot?

FullForm[norm[5!]]

This is just a list of two lists: the first of the points in the plot and the second of the options to be used in the plot (in the form of substitution rules).

norm[2000!]

How about the norm of a point plot? Define it to be the number of points in the plot.

norm[g_Graphics] := Length[g[[1]]]

norm[norm[200!]]

The plot comes from the evaluation of norm[200!], the 46 from the evaluation of norm[norm[200!]]. Check the 46th prime

Prime[46]

Nest[f,x,n] gives f[f[...f[x]...]] (n times). NestList[f, x, n] gives a list of all intermediate results: {x, f[x], f[f[x]],...}

NestList[norm, 200!, 5]

The first element of the list is the evaluation of 200! (norm applied zero times to 200!). The second entry is the norm of 200!, i.e., a plot of its prime decomposition. The third entry is the norm of the plot, i.e., the number of points in the plot (46). The fourth is the norm of 46, i.e., a plot of its prime decomposition (46=2x23). The fifth entry is the number of points in the plot (2) and the sixth is a plot of the prime decomposition of 2.

Check the definitions associated with norm

?? norm

Give a rule for the sum of norms of complex arguments

norm[x_Complex] + norm[y_Complex] := Sqrt[norm[x y]]

Every rule given with Set (=) or SetDelayed (:=) is associated with the head of the lhs. Mathematica prefers not to add any rules to Plus, so it gives it the attribute Protected

?? Plus

One can associate a rule with the heads of the arguments of the lhs using ^= or ^:=

norm[x_Complex] + norm[y_Complex] ^:= norm[Floor[Abs[x y]]]

The rule now is associated with norm, not with Plus

?? norm

norm[5 + I 3] + norm[4 - I 6]

Another example

f2[g[x_], h[y_]] := g[y]^h[x]

The above rule is associated with f2 (the head of the lhs)

?f2

Clear[f2]

The same rule, given in the form

f2[g[x_], h[y_]] ^:= g[y]^h[x]

is no longer associated with f2

?f2

but with the heads of the arguments of the lhs

?g

?h

Clear[f2, g, h]

Values assigned to a symbol via = or := are called downvalues of that symbol. Those assigned via ^= or ^:= are called upvalues.

Define norm[1.1]

norm[1.1] = Pi

?? norm

Which definition does Mathematica use?

norm[1.1]

Mathematica classifies rules according to generality and tries particular cases before general ones.

Any new rule given overwrites any old one with the same lhs

norm[1.1] = E

?? norm

norm[x_Real] := Abs[x]

?? norm

General pattern constraints

pattern /; condition makes pattern match only if condition is met

fact[n_  /;  n > 0] := n!

{fact[5], fact[3.2], fact[-1]}

(Mathematica defines n!=Gamma[n+1]).

Cases[{(Subscript[e, 1]), (Subscript[e, 2]), ...}, pattern] returns a list of the (Subscript[e, i]) that match pattern.

Cases[
{x, 3  x, 3  x + 2, x^2 + 11, 
        Sqrt[x]}, 
a_ x^n_ + 
        c_
]

a_. x matches x too (when a_. is absent, it is taken to have the value 1). 
The default values of a_. depend on the context

a_. x       default 1 (e.g., matches x)
x+a_.       default 0 (e.g., matches x)
x^a_.       default 1 (e.g., matches x)

Cases[
{x, 3  x, 3  x + 2, x^2 + 11, 
        Sqrt[x]}, 
a_.  x^n_ + 
        c_
]

Cases[
{x, 3  x, 3  x + 2, x^2 + 11, 
        Sqrt[x]}, 
a_.  x^n_.  + 
        c_
]

Cases[
{x, 3  x, 3  x + 2, x^2 + 11, 
        Sqrt[x]}, 
a_.  x^n_.  + 
        c_. 
]

FullForm[Sqrt[x]]

Cases[
{x, 3  x, 3  x + 2, x^2 + 11, 
        Sqrt[x]}, 
a_.  x^n_.  + 
          c_.  /; (! PrimeQ[n]) &&  
          n ≥ 1
]

Logical operators

&&          logical And
||          logical Or
!             logical Not

PrimeQ returns True if its argument is a prime and False otherwise.
The first argument of Cases can have any head, the result retains that head

Cases[
x + 3  x + 3  x + 2 + x^2 + 11 + 
        Sqrt[x], 
 
      a_.  x^n_.  /; (! PrimeQ[n]) &&  
          n ≥ 1
]

What happens if the dot is taken out from a_.?

Cases[
x + 3  x + 3  x + 2 + x^2 + 11 + 
        Sqrt[x], 
 
      a_ x^n_.  /; ! PrimeQ[n] &&  
          n ≥ 1
]

FullForm[x + 3  x + 3  x + 2 + x^2 + 11 + Sqrt[x]]

Select

Select[list,test] returns a list of those elements of list for which the boolean pure function test returns True. PrimeQ returns True if its argument is a prime and False otherwise.

Select[{-1.2, 3, 199, 37, -5, 2}, PrimeQ]

Select[
{-1.2, 3, 199, 37, -5, 
        2}, 
PrimeQ[#] &&  # > 0 && 
          OddQ[#] &
]

Exercise

Find the powers of 2 that start with 7.
Solution

findD[m_, d_, 
        M_] := 
Select[
{#, m^#} & /@ 
          Range[0, M], 
        First[IntegerDigits[#[[2]]]] == 
            d &
]

findD[2, 7, 500] // TableForm

Delayed Substitution Rules

Notice

x[a, b[r, g]] /. 
      x[z_, s_] -> (z s /. 
            s[[1]] -> s[[2]])

The rhs of the rule was evaluated immediatelly (before s was replaced by b[r,g]) so that s[[1]] and s[[2]] could not be defined. Using :> instead of -> (i.e., 
RuleDelayed instead of Rule) causes the rhs to be evaluated after the substitutions have been made

x[a, b[r, g]] /. 
      x[z_, s_] :> (z s /. 
            s[[1]] -> s[[2]])

Encryption (again)

contFract[list_List] := 
        Fold[(#2 + 1/#1 &), First[Reverse[list]], 
          Rest[Reverse[list]]];

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

str = "Hello.";

ToCharacterCode[str]

N[contFract[ToCharacterCode[str]], 1000]

contFractSin[list_List] := 
        Fold[(Sin[#2] + 1/Sin[#1] &), First[Reverse[list]], 
          Rest[Reverse[list]]];

str2 = N[contFractSin[ToCharacterCode[str]], 1000]

contFractSinList[list_List] := 
      FoldList[(Sin[#2] + 1/Sin[#1]  &), First[Reverse[list]], 
        Rest[Reverse[list]]]

ListPlot[N[contFractSinList[ToCharacterCode[str]], 1000]]

str2 = "Hello. I am an experimental text. I have to be rather long, \
so, please, bear with me."

ListPlot[N[contFractSinList[ToCharacterCode[str2]], 1000]]

str3 = Table[Random[Integer, {1, 100}], {n, 1, 100}]

ListPlot[N[contFractSinList[str3], 1000]]

Newton´s method

lfpl[x_] := 
      Length[FixedPointList[Function[z, (1/z^4 + 4  z)/5], x]]

pp = 401;
data = 
        Table[lfpl[x  +  I y], {y, -1. , 1. , 2/pp}, {x, -1. , 
            1. , 2/pp}];

ListDensityPlot[
data, 
Mesh -> False, 
ColorFunction -> (Hue[
                3  #] &), 
MeshRange -> {{-2, 
              2}, {-2, 2}},  
ImageSize -> 
          600
];

ListPlot3D[
data, 
Lighting -> 
        True, 
ColorFunction -> 
        Automatic, 
Mesh -> 
        False, 
ImageSize -> 
        600, 
PlotRange -> 
        All, 
ViewPoint -> {3, 4, 
          9}
]

Symbolic Manipulation in Quantum Mechanics: an Example

I give here a realistic example of symbolic manipulation using Mathematica. The problem described arises in the study of the motion of a spin-1/2 particle on a cylinder, in the presence of a radial magnetic field --- details can be found in [7] (hopefully available on the web by the end of November). Here the focus will be on the algebraic manipulation problem itself. 

Consider a family of functions f[n](x) satisfying the recursion relation

f[n](x) = (e^(-(b^2)/
              4)) (∫_0^1){f[
                    n - 1] (-b/(Sqrt[2])) - (e^(-b/(Sqrt[2])  (Subscript[xs, n]))) 
                  f[n - 1] (-(Subscript[xs, n]) - 
                      b/(Sqrt[2])) + (∑_(k = 1)^(n - 1))f[
                      n - k - 1] (-b/(Sqrt[2])) 
                    f[k] ((Subscript[xs, n]))} d(Subscript[s, n])/(Subscript[s, n])

(f[0](x)=1). The first part of the problem is to calculate f[n](x). The second part is to evaluate the wavefunction <z| f[n]((a^†)) |0>, where |0> is the ground state of a simple harmonic oscillator (SHO)  and (a^†) is the standard SHO creation operator.

Calculation of f[n](x).

r2 = Sqrt[2];

Given a variable v, generate a new variable vp

p[v_] := 
ToExpression[
        StringInsert[ToString[v], "p", -1]]

{p[s], p[f][3]}

Inert integration

The integrals that appear in the recursion relation for f[n](x) above cannot be evaluated in closed form, in general. Using Integrate would force Mathematica to try to evaluate them. Use instead an operation int and define its properties in the form of substitution rules.

It is convenient to group the rules in various lists (instead of a single one), so that one can use particular groups of rules and/or control the order of their application. 

int[f,{(Subscript[s, 1]),...,(Subscript[s, k])}] represents the multiple integral  (∫_0^1)... (∫_0^1)f \(TraditionalForm\`\(\(ds\_1 ... \) ds\_k\)\/\(\(s\_1 ... \) s\_k\)\)
(notice the inclusion of (Subscript[s, 1])...(Subscript[s, k]) in the measure). The aim is to express f[n](x) as such a multiple integral. 

intRules1 is a list of two substitution rules. The first one is 
c int[f,v]=int[c f,v], where c is a constant (i.e., independent of the integration variables contained in the list v).
FreeQ[expr,var]is True if expr does not contain var. Map FreeQ[expr,#]& down a list of variables to see if any of the variables in the list is contained in expr.

FreeQ[x + y + z, #] & /@ {a, b, z}

Change List to And

And @@ (FreeQ[x + y + z, #] & /@ {a, b, z})

(since x+y+z contains z).

The second rule in intRules1 adds up integrals with the same integration variables

intRules1 = 
{
c_ int[f_, v_] /; 
            And @@ (FreeQ[c, #] & /@ v) :>  
          int[Expand[c f], v], 
int[f_, v_] + 
            int[g_, v_] ->  
          int[f + g, v]
}

The first rule uses delayed substitution so that Expand has something to operate on. Notice that the rule can be applied to a product of integrals, as long as one of them does not involve the variables of the other.

There are two more rules to implement:

∫f(x) dx + ∫g(y) dy = ∫(f(x) + g(x)) dx 
(addition of integrals with renaming of integration variables) and

∫(∫f(x, 
                  y) dx) dy = ∫f(x, 
            y) dx dy.

intRules2 = 
{
(intX___ + 
                int[f_, v_] + int[g_, vp_]) /; 
            Length[v] == Length[vp] :> 
          intX + int[f + g /. Thread[vp -> v], v], 
        int[int[f_, v1_], v2_] :> 
          int[f, Flatten[{v1, v2}]]
}

The first rule in intRules2 checks for equality of the number of integration variables and effects the renaming with Thread. The delayed substitution is necessary so that Thread sees lists as arguments (similarly for Flatten in the second rule).

Computing f[n](x)

Define the recursion for f[n](x). Use primed variables in the product of integrals

Clear[fI]

fI[0, x_, v_] = 1;

fI[n_Integer /; n > 0, x_, 
          s_[n_]] := (int[
E^(-c^2/
                            2) (
fI[n - 1, -c, 
                            s[n - 1]]
 - 
                          E^(-c x s[n])*
                            fI[n - 1, -x s[n] - c, s[n - 1]] + 
                          Sum[
fI[n - k - 1, -c, 
                                s[n - k - 1]] fI[k, 
                                x s[n], p[s][k]], 
{k, 
                              1, n - 
                                1}
]
) // 
                    ExpandAll, 
{s[n]}
] //. 
                intRules1) //. intRules2 // ExpandAll;

Notice the repeated application of the rules via //.. There is no danger of infinite looping because the rules consistently "bring in" everything inside int. A set of rules some of which "bring in" while others "take out" can lead to trouble.

?? fI

fI[1, x, s[1]]

fI[2, x, s[2]]

fI[3, x, s[3]]

fI[4, x, s[4]]

Collect powers of x (in fact, just the linear term) in the exponents (verify that they live in level 4)

fIC[n_Integer /; n > -1, x_] := 
      Map[Collect[#, x] &, fI[n, x, s[n]], 4]

fIC[1, x]

fIC[2, x]

fIC[3, x]

Notice how x is pulled out as a common factor whenever possible.

Calculation of < z| f[n] ((a^†))|0>

Define wF[n,z] = <z| f[n]((a^†)) |0>. 

f[0](x) = 1, so that wF[0, z] = <z| 0> =(Subscript[N, 0]) (e^(-(z^2)/2))

wF[0, z] = (Subscript[N, 0]) E^(-z^2/2);

Similarly, any constant term in f[n](x) gives rise to a gaussian centered at the origin. Terms of the form (e^(cx + b)) give rise to translated gaussians

<z| (e^((ca^†) + b))|0> = ... = (Subscript[N, 0]) (e^(b + (c^2)/2)) (e^(-(((z - (Sqrt[2]) c))^2)/2))

(use the fact that \(TraditionalForm\`e\^\(ib \( p\&^\)\)\) is a translation operator and the Baker-Cambell-Hausdorff formula to rewrite it as a product of exponentials).

If[condition, action1, action2] performs action1 if condition is True and action2 otherwise. This is mapped down the sum of exponentials that constitutes the first part of fIC[n,x] (there is an implied integration as in fIC[n,x] which is suppressed for simplicity).

wF[n_Integer  /;  n > -1, 
        z_] := 
Map[
If[
FreeQ[#, x], 
#  (Subscript[N, 0]) 
              E^(-z^2/2), 
# /. 
              E^(c_. *x +  b_. ) ->  (Subscript[N, 0]) E^(b + c^2/2) 
                  E^(-(z - r2*c)^2/
                        2)
] &, 
fIC[
            n, x][[1]], 
1
]

wF[0, z]

wF[1, z]

(remember to integrate with d(Subscript[s, 1])/(Subscript[s, 1])).

wF[2, z]

Take out as common factor (Subscript[N, 0]) (e^(-(nb^2)/4))

wFN[0, z] = wF[0, z];

wFN[n_Integer  /;  n > -1, 
        z_] := 
(Subscript[N, 0]) E^(-n b^2/4) 
          Expand[E^(n b^2/4)/(Subscript[N, 0])  wF[n, z] /. c -> b/Sqrt[2] //. 
              a_/Sqrt[2] + b_/Sqrt[2] -> (a + b)/
                  Sqrt[2]
] // TraditionalForm

wFN[1, z]

wFN[2, z]

wFN[3, z]

wFN[4, z]

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).
[7] C. Chryssomalakos, A. Franco, A. Reyes, "Spin 1/2 Particle on a Cylinder with Radial Magnetic Field", in preparation.