More Advanced Mathematica

Expressions

Clear[x]

We have dealt quite a bit with expressions

e.g.

x

x+x

a*x^2+b*x+c

Lists...

{x,0,1}

Plots...

Plot[x,{x,0,1}]

Equations...

a*x^2+b*x+c==0

Functions

addOne[x_]:=x+1

Rules

n->10

What exactly are expressions?

Expressions are built by combining more basic expressions into more complicated expressions

e.g.

a+b

Plus[a,b]

can be thought of as being made up of expressions for a, b, and +

So a complicated expression can be broken down into simpler expressions

At a certain point, we obtain simple, irreducible expressions that cannot be broken down anymore

These irreducible expressions are called atoms

There are three main kinds of atoms:

Numbers

Strings

Symbols

We have seen many examples of numbers: these can be integers or inexact (floating point) numbers

Strings are more-or-less self-explanatory:

"This is a string"

Symbols are what distinguish Mathematica from other kinds of languages

In the example above, a and b are symbols

The + sign is also a notation for the Plus symbol

We can decompose an expression into its atoms using FullForm

FullForm[a*(x^2)+b*x+c]

More complicated expressions

Let’s return to our quadratic equation

quadratic=a*x^2+b*x+c

We can decompose it into atoms:

FullForm[quadratic]

We see that this expression is created using nested symbols and numbers

The nested structure can be visualized intuitively using TreeForm

TreeForm[quadratic]

Extracting parts of expressions

We can programmatically extract parts of expressions:

Recall:

quadratic//TreeForm

The root of the expression tree can be obtained by calling Head

Head[quadratic]

This is always the same thing as getting the 0th index:

quadratic[[0]]

The “levels” of this tree can be extracted using Level

Level[quadratic,{0}]

Level[quadratic,{1}]

Level[quadratic,{2}]

Level[quadratic,{3}]

Level[quadratic,{4}]

More sophisticated expressions

plt=Plot[x,{x,0,1}]

FullForm[plt]

Holding expressions

Consider the expression

expr=x+x

FullForm[expr]

Mathematica evaluates our input expression and converts it to a canonical form

If we want to maintain the original expression in unevaluated form, use Hold

uneval=Hold[x+x]

TreeForm[uneval]

Hold can be released using ReleaseHold

ReleaseHold[uneval]

How are symbols evaluated

Suppose we start with an undefined symbol

x

We can assign a value to this symbol using = (Set)

x=Random[]

How does Mathematica know to replace every instance of x with its value?

It uses a global look-up table. We can find the rules for evaluating symbols by calling OwnValues

OwnValues[x]

The symbol :> (typed :>) means replace the lhs with the rhs, and only then evaluate the rhs

HoldPattern is a way of expressing a rule for a symbol without having the symbol transformed by other existing rules

y:=Random[]

OwnValues[y]

Since we used SetDelayed, we see that in this case y is replaced by a new random number each time

How functions are evaluated

Recall our Fibonacci sequence

fib[0]:=0
fib[1]:=1
fib[n_]:=fib[n-1]+fib[n-2]

OwnValues[fib]

fib has no OwnValues. That’s because if we evaluate

fib

no rule is applied. Rules are only applied if we give an argument to fib

fib[5]

These kinds of rules (that require arguments) are known as DownValues

DownValues[fib]

The order of these rules shows the order in which Mathematica will attempt to apply them

There are also UpValues and SubValues, but they are a bit obscure...

Understanding function evaluation

Trace helps us understand the order in which Mathematica evaluates expressions

Sin[Log[Exp[1],10/3]]

Trace[Sin[Log[Exp[1],10/3]]]

Plot[Sin[x]+Cos[x],{x,0,3}]

Trace doesn’t work with all functions, but it works with many

Trace[FindRoot[Sin[x]+Cos[x],{x,2}]]

However, FindRoot supports EvaluationMonitor

FindRoot[Sin[x]+Cos[x],{x,2},EvaluationMonitor:>Print[x]]

Patterns

We have seen one type of pattern before: x_

The underscore will match any expression, and then the result will be called x

We can test if a pattern matches an expression using MatchQ

MatchQ[1,x_]

MatchQ["Anything will match the pattern x_",x_]

There are other kinds of patterns.

Except will match anything except what’s in the brackets

MatchQ[2,Except[1]]

MatchQ[1,Except[1]]

You can match for specific heads of expressions with _h

MatchQ[1,n_Integer]

MatchQ[2/3,n_Integer]

Head[2/3]

Conditions can be added to a pattern with /;

MatchQ[-1.5,n_/;n<0]

MatchQ[4,n_/;n<0]

Boolean functions can also be added with the ? operator

MatchQ[4,n_?EvenQ]

MatchQ[5,n_?EvenQ]

Patterns can be combined using | for or

MatchQ[5,_?EvenQ|5]

Checking for specific heads allows us to perform type checking

fib[0]:=0
fib[1]:=1
fib[n_Integer]:=fib[n-1]+fib[n-2]
fib[_]:=Print["fib must be called with integer argument"]

fib[4]

DownValues[fib]

fib[5]

fib[2/3]

fib[-1]

Let’s get rid of our old (overly-general) rule

fib[n_Integer]=.

We want to add a condition to the pattern: only match non-negative integers

fib[-1]

fib[n_Integer/;n≥0]:=fib[n-1]+fib[n-2]
fib[_]:=Print["fib must be called with non-negative integer argument"]

fib[-1]

Another example

abs[x_/;x≥0]:=x
abs[x_/;x<0]:=-x

abs[10]

abs[-10]

Example: defining our own derivative function

We want to differentiate arbitrary expressions

This can be achieved using the basic rules of differentiation:

Linearity

Product rule

Chain rule

We first need the concept of an expression being independent of another expression.
This is achieved using FreeQ

FreeQ[(y^2),x]

FreeQ[y*x,x]

Some common base cases:

diff[x_,x_]:=1
diff[c_,x_]/;FreeQ[c,x]:=0

Addition:

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

Scalar multiplication:

diff[c_*f_,x_]/;FreeQ[c,x]:=c*diff[f,x]

Product rule

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

Chain rule

diff[f_[g_],x_]/;!MatchQ[g,x]:=diff[f[g],g]diff[g,x]

(Have to be careful in the above not to match f[x] (so make sure g doesn’t match x), otherwise have infinite recursion)

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

diff[f[x]g[x],x]

poly=Expand[(((x+1))^10)]

diff[poly,x]

D[poly,x]

If we encounter functions whose derivatives we don’t know, it will leave them unevaluated:

diff[Sin[x]Cos[x],x]

diff[Sin[(x^2)],x]

We can easily add rules for trigonometric functions:

diff[Sin[x_],x_]:=Cos[x]
diff[Cos[x_],x_]:=-Sin[x]

diff[Sin[(x^2)],x]

diff[Sin[x]Cos[x],x]

Our definitions will automatically work with function composition:

diff[Sin[(x^2)]Cos[(x^3)],x]

Trace[diff[Sin[(x^2)]Cos[(x^3)],x]]

D[Sin[(x^2)]Cos[(x^3)],x]

Speed: memoization

Recall our Fibonacci sequence

fib[0]:=0
fib[1]:=1
fib[n_]:=fib[n-1]+fib[n-2]

How fast (slow?) is it?

Timing[fib[20]]

Timing[fib[25]]

Timing[fib[30]]

timings=Table[Timing[fib[i]][[1]],{i,1,30}]

timingplt=ListPlot[timings,Joined->True,PlotRange->All]

In fact, the runtime for the recursive algorithm satisfies

({{({{(Subscript[T, n])=(Subscript[T, n-1])+(Subscript[T, n-2])}})}})

and so the time to compute fib[n] is actually given by fib[n]

We can easily find a closed form representation of the Fibonacci sequence (involving the golden ratio) using e.g. the characteristic polynomial

FullSimplify[FunctionExpand[Fibonacci[n]],{n∈Integers}]

factor=(timings[[30]]/(GoldenRatio^30));

estimateplt=Plot[(GoldenRatio^n)*factor,{n,0,30},PlotRange->All,PlotStyle->{Black,Dashed}]

Show[{timingplt,estimateplt}]

We can try to optimize our recursive implementation using memoization

Trace[fib[10]]

Look how many times fib[1] and fib[0] are evaluated! We can speed things up by storing these values

fib[n_]:=fib[n]=fib[n-1]+fib[n-2]

Trace[fib[15]]

Trace[fib[10]]

bench[n_]:=Module[{myfib},
myfib[0]:=0;
myfib[1]:=1;
myfib[m_]:=myfib[m]=myfib[m-1]+myfib[m-2];
myfib[n]]

newtimings=Table[Timing[bench[i]][[1]],{i,1,30}]

newtimingplt=ListPlot[newtimings,Joined->True,PlotRange->All,PlotStyle->Red]

Show[timingplt,newtimingplt]

benchresults=Table[Timing[bench[i]][[1]],{i,1,500}];
benchplot=ListPlot[benchresults,Joined->True,PlotRange->All,PlotStyle->Red]

linearestimateplot=Plot[benchresults[[400]]/400*n,{n,0,500},PlotRange->All,PlotStyle->{Black,Dashed}]

Show[benchplot,linearestimateplot]

We can see that the timings are linear, which is clear given the recurrence

Code compilation

Mathematica code is by default interpreted by the Mathematica kernel

The kernel can also compile a subset of Mathematica code for better performance

e.g. approximate π numerically

Sum[((((-1))^k)/(2k+1)),{k,0,Infinity}]

Truncate this sum at n terms

piapprox[n_]:=4Sum[((((-1))^k)/(2k+1.0)),{k,0,n}]

piapprox[100000]//N

RepeatedTiming[piapprox[100000]]

cf=Compile[{n},4Sum[((((-1))^k)/(2k+1.0)),{k,0,n}],CompilationTarget->"C",RuntimeOptions->"Speed"]

RepeatedTiming[cf[100000]]

Needs["CompiledFunctionTools`"]

CompilePrint[cf]

Generating code

Needs["CCodeGenerator`"]

CCodeStringGenerate[cf,"approximatePi","CodeTarget"->"Standalone"]

SymbolicCGenerate[cf,"approximatePi","CodeTarget"->"Standalone"]

JIT Compilation

Suppose we want to roll our own JIT compilation

We want to speed up Select, which returns those elements of a list that match a pattern

mylist=RandomInteger[{1,100000},1000000];

Timing[Select[mylist,Mod[#,2]==0&]][[1]]

selectJIT[list_,pred_]:=(jitted[pred]=Compile[{{lst,_Integer,1}},Select[lst,pred],CompilationTarget->"C",RuntimeOptions->"Speed"];
selectJIT[list2_,pred]:=jitted[pred][list2];
jitted[pred][list])

Timing[selectJIT[mylist,Mod[#,2]==0&]][[1]]

Timing[selectJIT[mylist,Mod[#,2]==0&]][[1]]

Timing[selectJIT[mylist,Mod[#,17]==0&]][[1]]

Timing[selectJIT[mylist,Mod[#,17]==0&]][[1]]