Introduction to Mathematica, Part II

Mathematica can be used in the same fashion as a traditional programming language, such a fortran or C.  For  frequently used programs or computationally intensive applications, Mathematica is not the best programming language choice, but it can be very effective as a way to test out ideas or generate a single numerical result.

Evaluating a sequence of instrutions (;;;)

a  =  1;
 a  =  a + a;  a = a^a
 a =  a + a;   a = a^a

Clear[a]

Loops

?Do

a = 1;  Do[a = 2  a; a = a^a, {i, 1, 2}]

a

a = 0.1;  
    Do[a = 2  a; a = a^a;  
      Print["iteration is ",  i,  " and a is ",  a], {i, 1, 
        4}]

Clear[a]

?For

For[a = 0.1;  i = 1, i  ≤  4, i++, a = 2  a; a = a^a;  
          Print["iteration is ",  i,  " and a is ",  a]]

?While

?Table

Clear[a]

a = 0.25;
Table[{i, a = 2  a; a = a^a}, {i, 1, 4}]

a = 0.75;
Table[{i, a = 2  a; a = a^a}, {i, 1, 4}]

datatable  = 
      Table[{dx, 
          For[a = dx; i = 1,  i  ≤  4,  i++,  a = 2  a;  
            a = a^a]; Log[a]}, {dx, 0.01, 0.5, 0.01}]

ListPlot[datatable]

Options[ListPlot]

ListPlot[datatable, PlotRange -> {250, 500}]

?*Minimum*

FindMinimum[
      For[a = dx; i = 1,  i  ≤  4,  i++,  a = 2  a;  
        a = a^a]; Log[a], {dx, 0.15, 0.25}]

FindMinimum[For[a = dx; i = 1,  
        i  ≤  4,  i++,  a = 2  a;  a = a^a]; 
      Log[a], {dx, 0.15, 0.25}

Lets try and do the above the hard way.  I will use Nest to recursively apply the function 4 times (I am just using a shorthand here, we can ignore the use of Nest for this course...). You can see that it works. Don't worry about it, but if you want to know about it, use the Help Browser to get information about Nest and Pure Functions.

Clear[x]

fx  =  Nest[(2  #)^(2  #) &, x, 4]

Take it derivative and set equal to zero...

dfx  = D[fx, x] // Simplify

Finding the zero of this will not be easy.... but FindRoot claims it can do it...

FindRoot[dfx, {x,  .1,  .3}]

Local Variables

In the above, the values of a and  dx are used repeatly—however, if we are only interested in the resulting minimum, their final values have little meaning or importance.  Suppose, we use the value of a or dx without clearing—what will the value be?

dx
a

CurrentValueofA  =  a;

As this could cause problems if a is used in a function call, say Sin[2 π a], it is a good idea to localize variables.

?Module

dx  =  SnickerDoodle;  
    a  =  HappyGoLucky;
Module[
{dx, a, maxiteration = 4, solution, 
        i}, 
solution  =  
      FindMinimum[
          For[a = dx; i = 1,  i  ≤  maxiteration,  i++,  
            a = 2  a;  a = a^a]; Log[a],  {dx, 0.15, 0.25}]; 
      Print[dx /. 
          solution[[2]]]
]

dx
a
solution

Switches: If, Which

?If

a  =  Prime[23]  +  Prime[62]  +  Prime[104]

Factor[a]

?PrimeQ

If[PrimeQ[a], 
Print[
        a ,  " is a Prime Number"], 
Print[
        a,  " is not Prime, its divisors are ",  
        Divisors[
          a]], 
Print["I have no idea what you are asking \
me to do!"]
]

The above program is ok, but not very useful because it only works for  the current value of a.  It would be more useful to have something that worked for any value of a and could use it over again—that is, turn it into a tool.  This involves patterns and function definitions.

Patterns (_)

Patterns and Pattern Replacement are probably one of the most powerful features of Mathematica.  The underscore _ stands for a pattern,  a symbol x with an underscore (x_) is recognized as anypattern that matches x.  Here are some examples.

AList  =  {first,  second,  third  =  2 first,  
        fourth  =  2 second}

AList /. {2 a_  ->  a}

Clear[a]

AList /. {2 a_  ->  a}

AList /. {p_ ,  q_ ,  r_ ,  s_}  ->  {p ,  p q,  
          p q r,  p q r s}

A qualifier on a pattern can be used to restrict which expressions will match the pattern:

{2, 0.667,  a/b,  Pi} /. {p_Integer  ->  p One}

_ all by itself stands for anything.  x_ also  stands for anything, but gives anything a name for later use.

Alist /. _  ->  AppleDumplings

This is not very useful and it really is doing what you asked for—the first thing it found was AList itself and replaced AList with AppleDumplings. 

x^n_ stands for x to any power, and the power can be refered to as n; x_^n_  stands for any expression raised to a power.  We can used these ideas to find rules that look like derivatives (this is not a good way to make a symbolic derivative, but it is illustrative):

PaulieNoMealX  =  Sum[b[i] x^i, {i, 2, 6}]

PaulieNoMealX /. x^n_  ->  n x^(n - 1)

This is ok, but it has least two problems.  The first is that it would be nice to have the rule work for any polynomial...

DerivRule  =  q_^n_  ->  n q^(n - 1);

PaulineOMealY  =  Sum[c[i] z^i, {i, 2, 6}]

PaulineOMealY /. DerivRule
PaulieNoMealX /. DerivRule

Another problem is that it will not work for first-order and zeroeth-order terms...

PaulENoMiel  =   Sum[c[i] HoneyBee^i, {i, 0, 6}]

PaulENoMiel /. DerivRule

This could be fixed, but it would be much easier to do so by defining functions of a pattern

It is also possible to have a pattern apply conditionally using a condition. Here is an example that pulls out cases where the first member of a pair of numbers is less than the second.

Cases[{{1, 2},  {2, 1},  {a, b},  {2,  84},  
        5}, {first_,  second_} /; first < second]

Defining Functions with Patterns

Defining functions with patterns probably combines the most useful aspects of Mathematica.  Define a function that takes patten matching x as its first argument and an argument matching n as its second argument and returns x to the (n^th) power:

f[x_ ,  a_]  =  x^a;

\(What\ happens\ above\ is\ the\ right–hand–side\ is\ evaluated\ immediately\), and its value is assigned to f[a, b].  It was probably intended that f[2, 3] should always return 8 and that f[y, z] should always return (y^z), \(but\ this\ is\ not\ what\ Mathematica\ was\ asked\ to\ \(\(do\)\(:\)\)\)

f[2, 3]
f[y, z]

This works fine, but suppose we had defined  x ahead of time

x  =  4

f[x_ ,  a_]  =  x^a;

RHS evaluated, it is x^(second argument) which is currently 4^(second argument)

f[2, 3]

f[y, z]

Delayed Assignmet (:=)

The problem with the above functions is that the evaluation is made immediately.  For a function, the argument serves as a placeholder for a future calculation.  You do not wish to have the right–hand–side evaluated until the actual values of left–hand–side are specified.  The remedy is to use delayed assignment :=

x = 4
a  =  ScoobyDoo

f[x_ ,  a_]  :=  x^a

f[2, 5]

f[y, z]

f[x, a]

f[a, x]

Clear[f]

Functional Programming with Rules

The canonical programming example is the factorial function n! = (n)×(n-1) ×(n-2)×⃛×(1) where 0! ≡  1; here is a reasonably clever way to use the fact that (n+1)! = (n+1)×n!

factorial[n_]  :=  n factorial[n - 1]

This is part of definition, try it out...

factorial[8]

Ooops, Mathematica does not know when to quit.  Add the second part of the definition.

factorial[0]  =  1;

factorial[120]

factorial[257]

$RecursionLimit

$RecursionLimit  =  2^24

Timing[factorial[2000]][[1]]

Using immediate assignment in a function: spending memory to buy time.

factorial[n_]  :=  factorial[n]  =  
        n*factorial[n - 1]

Timing[factorial[2000]][[1]]

Timing[factorial[2001]][[1]]

Clear[factorial]

Functions and Patterns with Restricted Rules

The factorial function is pretty good, but not foolproof:

Clear[factorial]

factorial[0]  =  1;  
    factorial[n_]  :=  n*factorial[n - 1]

factorial[Pi]

The remedy is to restrict the pattern:

Clear[factorial]

factorial[0]  =  1;  
    factorial[n_Integer]  :=  n*factorial[n - 1]

factorial[Pi]

Still not perfect:

factorial[-5]

Clear[factorial]

factorial[0]  =  1;  
    factorial[n_Integer?Positive]  :=  n*factorial[n - 1]

factorial[12]

factorial[Pi]

As a last example, let's define the Sign function. It should be -1 when its argument is negative, 0 when its argument is zero, and +1 when its argument is positive.  There are lots of ways to write this function, there is no best way. Whatever works is good.

?Sign

HeyWhatsYourSign[0]  =  0;  HeyWhatsYourSign[_?Positive]  :=  1;  
    HeyWhatsYourSign[_?Negative]  :=  -1;

Plot[
      HeyWhatsYourSign[
        argument], {argument, -π, E}]

Plot[{1/x, HeyWhatsYourSign[x]/x}, {x, -1, 1}, 
      PlotStyle -> {{Hue[1], Thickness[0.02]}, {Hue[0.66], 
            Thickness[0.01]}}]