Mathematica tutorial II, Dynsys 2018

Patterns

Patterns in Mathematica are useful in order to:
define functions
apply functions to classes of expressions => Efficiently transform input data to output data
 make custom-made mathematical operations on expressions
identify and extract/remove subexpressions

Blank patterns (Patterns on the form \!\(\*"\""\)\!\(\*"x_"\)\!\(\*"\""\))

(* Here "_" denotes a pattern that matches anything and the content is assigned to a variable "x" in a transformation rule. If we are not interested in the value of the pattern we can simply write "_" *)

(* Patterns can be used to substitute parts of an expression (here the expression is f[a]) *)
f[a]/.x_->x^2
f[a]/.f[x_]->f[0]
f[a]/.f[_]->f[0]
f[a]/.f[x_]->x^2

(* NB: the pattern variable is local and takes priority over the same symbol on the left-hand side of the substitution rule.  *)
f[x,y]/.f[x,y]->x
f[x,y]/.f[x,x_]->x

(* To avoid confusion, do not use the same symbols for patterns as for your regular variables! *)

(* Patterns can be used in any structure of expressions, and any number of patterns can be used *)
{a+2b,x^c}/.{p1_+p2_,x^p3_}->{p1,p2,p3}

User-defined functions

(* Functions are declared using a delayed assignment (function:=rhs, the right-hande side is evaluated each time the function is used) *)
(* Arguments are entered as patterns *)
f1[x_]:=x^2;
f1[1+a]

(* Functions can take several arguments.  *)
f1[x_,y_]:=(x+y)^2;
f1[a,b]

(* Use parentheses to perform several operations in a function. The return value is the last evaluation in the function (can also use Return[]-function). *)
f1[x_]:=(
y=x+1;
y^2
);
f1[a]

(* OBS: If you terminate the last line by a semicolon, the function will return what comes after the semicolon, i.e. a NULL expression *)
f1[x_]:=(
y=x+1;
y^2;
);
f1[a]
Print[f1[a]]

(* Retrieve information about what we have assigned to f1: *)
?f1

(* By default arguments are local and all variables within the function are global. *)
(* Use Module (or Block) to create environment with local variables *)
a=1;
f[x_]:=Module[{a=2},x+a];
f[3]
a
a=.;

Everything is expressions

(* To better understand how to match patterns it is useful to have some understanding of how Mathematica expressions are built up *)
(* All mathematical formulas, lists, graphics, etc. are expressions on the form f[x,y,...], i.e. as a list {x,y,...} of a type f *)
(* Examples: *)
x+y+z ===Plus[x,y,z]
{a,b,c} ===List[a,b,c]
(a->b)=== Rule[a,b]
a=b ===Set[a,b]

(* Use FullForm to see the form f[x,y,...] of an expression *)
sum=x+y+z;
FullForm[sum]

(* The elements of an expression can be obtained in the same way as elements of a list *)
sum[[1]]
sum[[2;;3]]
(* The zeroth element (Head) gives the type of an expression *)
sum[[0]]
Head[sum]

Head[x]
Head[3.14]

(* We can change elements of the expression *)
sum[[3]]=w;
sum
(* Note that expressions are sorted alphabetically, i.e. changing z to w changes the positions also of x and y *)
FullForm[sum]

(* Changing the head of an expression literally transforms it into another type *)
sum[[0]]=Times;
sum
sum[[0]]=f;
sum

(* Use @@ (Apply) to evaluate expression with a new head, without modifying the expression *)
List@@sum
Apply[List,sum]
sum

(* To get an overview of complicated expressions, use TreeForm *)
A=TreeForm[{{a b,c d^2},{x^3 y^4}}]
TreeForm[A]

(* A more complex tree *)
A=StreamPlot[{y,-Sin[theta]},{theta,-π,π},{y,-4,4},StreamScale->{Automatic,Automatic,0.04}]
TreeForm[A]

Other forms of patterns

(* Pattern with restricted head: x_Head *)
{1,7,c,x}/.p1_->match
{1,7,c,x}/.p1_Integer->match
{1,7,c,x}/.p1_List->match
{1,7,c,x}/.p1_Symbol->match

(* Condition on pattern: pattern/;condition->rule; Takes current pattern hit, checks whether condition is true and if so, substitutes the pattern using the rule *)
{1,4,7,9,13}/.p1_/;PrimeQ[p1]->match

(* Alternative pat1 | pat2 applies rule if either pat1 or pat2 matches *)
{1,4,7,9,13}/.(p1_/;PrimeQ[p1])|(p1_/;EvenQ[p1])->2p1

(* We can use the ":" operator to assign a pattern (parenthesis to the right of ":") to a new pattern variable (p3, to the left of ":")  *)
{1,4,7,9,13}/.p3:((p1_/;PrimeQ[p1])|(p2_/;EvenQ[p2]))->2p3

(* One should be careful not to use unmatched patterns. In the example below the unmatched pattern seems to evaluate to unity... *)
{2,4,7,9,13}/.(p1_/;(Print["Cond 1: p1=",p1,"; p2=",p2,"; Evaluation:",PrimeQ[p1]];PrimeQ[p1]))|(p2_/;(Print["Cond 2: p1=",p1,"; p2=",p2,"; Evaluation:",EvenQ[p1]];EvenQ[p2])):>(Print["p1=",p1,"; p2=",p2];a p1+b p2)

(* The patterns mentioned above are the ones I use most, but in some cases it is good to know that there exists additional pattern types: *)
(* Blank sequence p1__ matches any sequence of one or more expressions (useful to match arguments in a function  ) *)
f[x,y]/.f[p1__]->{p1}
(* Repeated p1.. gives a sequence of expressions that each matches the pattern p1 *)
(* Longest[p1] and Shortest[p1] matches the longest/shortest sequence of patterns consistent with p1 *)

(* Use of patterns to create function headings. Functions are matched in order given by ?f  *)
f[x_]:=x^2;
f[x_Integer]:=x^3;
f[x_/;x>2]:=x;

?f

f[a]
f[3]
f[3.2]
f[1.5]

Patterns in built-in functions

(* In the examples above patterns were used in substitution rules *)
(* There are many additional functions that make use of patterns, some examples being: *)

(* ReplaceRepeated[] (//.) works like the substition rule, but it is applied until expression no longer changes () *)
(* This is useful when we want to apply several transformation rules to an expression *)
(* For example contraction rules on Kronecker deltas and tensors: *)
rules={Delta[p1_ p2_]Delta[p2_ p3_]->Delta[p1 p3],Delta[p1_ p2_]A_[p2_]-> A[p1]}

(* Standard substitution only applies the rule once *)
Delta[i j]Delta[j k]f[k]/.rules
(* Use repeated substitution to apply the rule until convergence *)
Delta[i j]Delta[j k]f[k]//.rules

(* Alternative (suitable for this course), use FixedPoint[function,expression] which applies a function on an expression until nothing changes *)
FixedPoint[#/.rules&,Delta[i j]Delta[j k]f[k]]
(* Here & denotes the end of an anynomous function with argumemt # (initially the argument is given by the input expression, in subsequent steps the argument is the output from the previous step) *)
(* Alternative notation for the anonomus function: *)
FixedPoint[Function[x,x/.rules],Delta[i j]Delta[j k]f[k]]

(* Use cases to extract all occurences of a pattern in an expression, Cases[expression,pattern] *)
list={a,4,x^2,9,y^2 z}
(* By default Cases traverses the first tree-level in the expression *)
Cases[list,p1_^2]
TreeForm[list]
(* Use level specification {0,∞} to make sure to find all patterns at all tree-levels in the expression *)
Cases[list,p1_^2,{0,∞}]

(* What can we do if we want to match integers on the right-hand side of rules? *)
rules={a->3,b->x,c->17};
(* The code below does not do because Case interpret it as: match pattern p1 and substitute the rule to the right of \!\(\*"\""\)\!\(\*"->"\)\!\(\*"\""\). This is sometimes useful, but not here! *)
Cases[rules,p1_->p2_Integer]
(* Instead use a condition to check we have a rule with right-hand side being an integer *)
Cases[rules,p1_/;(p1[[0]]==Rule)&&(IntegerQ[p1[[2]]])]
(* Finally add a substitution rule to extract the integers *)
Cases[rules,p1_/;(p1[[0]]==Rule)&&(IntegerQ[p1[[2]]]):>p1[[2]]]

(* Similarly MemberQ determines if a pattern exists in the expression and DeleteCases removes all elements in an expression that matches the pattern *)

Programming in Mathematica

How you are supposed to program in Mathematica: Functional
Use rules and pattern matching to transform input data to output data
Written as a sequence of functions applied to an initial expression
Avoid loops at all costs
Avoid changing values of variables , functions do not have side effects, applying the function to a variable always give the same result
Generally slower execution time than procedural (but in Mathematica it is often faster!)
If you know what you do it is often quicker to develop and modify code compared to a procedural approach. However, I find it hard to use for larger projects;

How I secretly program: Procedural
C++ and Java like layout
Written as a list of instructions executed in sequence that changes values of variables:
expr;expr;expr;...
functions may have side effects, applying the function to a variable may give different results because the state has changed from last time the function was applied
Fast because of the straightforward translation from code to machine instructions
Big lists of instructions are easily split into smaller entities (functions)

Some examples of the two approaches

(* Start with two lists *)
list1={1,6,9};
list2={a,b,c};

(* Procedural: Apply a function to the list elements using a loop (Do[] runs the loop variable from 1 to the list length Length[list1])  *)
result={0,0,0};
Do[
result[[i]]=f[list1[[i]],list2[[i]]];
,{i,Length[list1]}];
result

(* Functional: Thread the function f over any lists appearing in the argument *)
result=Thread[f[list1,list2]]

Example functions used in functional programming

(* For Functional programming Mathematica has many built-in functions to apply functions to expressions (lists). Some examples: *)
(* Map[], Applies function to all elements in first level in expression *)
f /@ {a,b,c}
(* Apply[], discussed above, changes head of expression *)
f @@ {a,b,c}
(* Apply function n times to expression *)
Nest[f,{a,b,c},2]
(* For a simple list, thread does the same as f /@ {a,b,c}  *)
Thread[ f[{a,b,c}]]

(* These functions in combination with Fold and FoldList are very powerful *)
FoldList[f,x,{a,b,c}] 
(* Fold gives the last element obtained by FoldList[] *)
Fold[f,x,{a,b,c}] 

(* Example: Create a nested polynomial *)
Fold[x #1+#2&,0,{a,b,c,d,e}]

Example functions used in procedural programming

(* Loops *)
(* Do, While, For, Table *)

(* Conditional expressions *)
(* If,Switch,And(&&),Equal(==),Less(<) *)

(* NB: You would expect that If[expression] do someting if expression is true and something else if expression is false, but the If[] in Mathematica has three conditional results: True, False and otherwise *)
(* For example try to check if a is equal to 3 *)
If[a==3,True,False]
(* Since a is undefined this does not evaluate to false, it evaluates to undefined, the third option in If: *)
If[a==3,true,false,undefined]
(* Changing the equaltiy check to check if expressions are the same evaluates to False *)
If[a===3,true,false,undefined]

(* Modular environments *)
(* Module, With, Block *)

Other Topics

Numerical precision (for high-precision calculations)

(* Use Precision[] to retrieve number of digits of precision in an expression: *)
(* Integers are exact *)
Precision[4]
(* Floats are Machine precision by default (16 digits) *)
Precision[4.0]

(* Use SetPrecision[x,accuracy] to set number of relevant digits in an expression *)
x=SetPrecision[4.000000000000000000000000000000000001000000000000000000000000000000000000000000,100]
y=SetPrecision[4.000000000000000000000000000000000000000000000000000000000000000000000000001000,100]

(* Evaluation of x-y would not work in machine precision  *)
(* We loose accuracy due to cancellation in the subtraction *)
x-y
(* It turns out we lost 36.9 digits of accuracy... *)
Precision[x]
100-Precision[x-y]

(* Note that the high precision is lost if we mix it with low-precision numbers *)
x+1.1
Precision[x+1.1]

Manipulate
Allows for interactive changing of a parameter in an expression
Example for dynamic billiards in an elliptic pool table: http://demonstrations.wolfram.com/DynamicBilliardsInEllipse/