Introduction to Mathematica, Part II

Elements of Mathematica Programs

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

Note above that line breaks without ";" are points where the instructions get executed. This three-line sequence has five executions and produces two lines of output and the instructions get executed sequentially.

Clear[a]

Loops

?Do

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

a

In this case, Do does not produce screen output. The following will and is a simple illustration of how formatted output can be programmed:

a=0.1; Do[a=2a;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=2a;a=a^a; Print["iteration is ", i, " and a is ", a]]

?While

?Table

Clear[a]

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

Note above that we have asked Mathematica to create a list of two-element lists...

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

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

In the instruction above, note that the variable "datatable" is a list. In Mathematica, variables can be numbers, expressions, lists, plots, ...  We will see that this feature is very useful.

Plotting Data from Lists, Using More Mathematica Functions

ListPlot[datatable]

Options[ListPlot]

Note that the options are written as Rules.

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

?*Minimum*

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

By going into the Help Browser, you can see that the output of FindMinimum is a list, the first element of which is the functions minimum value, and the second is a Rule specifing where the minimum occurs.

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

Localized Variables and Switches

Local Variables

In the above, the values of a and  xvalue 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 xvalue without clearing—what will the value be?

xvalue
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

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

xvalue
a
solution

Switches: If, Which

?If

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

?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 and ``Placeholder Variables''

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}

Here's a way to put the last input command into words: Search list AList for any pattern that matches "2 × anything" and replace "2 × anything" with a, then output the new AList. Study the following examples carefully:

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}

The third member of this result may be a bit surprising; it has to do with the way Mathematica internalizes expressions. In this case it represents the fraction as a (b^(-1))and it performs the replacement on the exponent (note that it must treat the "-" sign independently of the "1", as it returns "-One" rather than "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 use 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

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

Functions: Using Patterns in Delayed Assignment (:=)

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; (*This is not a good way to define a function, we will see why later*)

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; (*This is not a good way to define a function, we will see why later*)

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

f[2,3] (*should now be 4^3, which is probably not what the programmer had in mind*)

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]

Using Recursion to Program

Functions that call themselves...

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]

try it out...

factorial[8]

Ooops, This isn't what was expected, but upon reflection it is correct--we forgot to define a part of the rule. (Note also that the message window produced an error about recursion limits) Add the second part of the definition. Here, we don't use delayed evaluation (:=) because we want to assign a value immediately.

factorial[0] = 1;

factorial[120]

factorial[257]

Here is where the recursion limit comes in : our function keeps on calling itself (i.e., recursively). Unless a limit is set the program would keep running forever. Mathematica builds in a limit to how many times a function will call itself:

$RecursionLimit

$RecursionLimit =2^11

Speed versus Memory in Functions

?Timing

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

Using immediate assignment in a function: spending memory to buy time: Each time the function is called, it makes an extra assignment so that previous values can be recalled if needed.

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

This version takes a bit longer the first time, because we are storing data in memory ...

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

But, the next time it is called, the result is much faster.

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

Clear[factorial]

Smarter Functions: Further Restriction on Patterns and Conditions

Restrictions on Patterns

The factorial function is pretty good, but not foolproof as the next few lines will show.

Clear[factorial]

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

The next line will cause an error to appear on the message screen.

factorial[Pi]

The remedy is to restrict the pattern:

Clear[factorial]

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

This time it doesn' t produce an error, and returns a value indicating that it is leaving the function in symbolic form for values it doesn' t know about.

factorial[Pi]

Functions and Patterns with Tests

However, the definition of factorial still needs some improvement--the next line will cause an error.

factorial[-5]

Clear[factorial]

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

factorial[12]

factorial[Pi]

Examples with Graphics and  Conditional Definitions

As a another 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

Here we write our own version, we don' t "name'' the pattern because it is not needed in the function definition.  It is a bit harder to read this way, but I use it here to be instructive.

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

Plot[HeyWhatsYourSign[argument],{argument,-π,E},PlotStyle->Thick]

Plot[{1/x,HeyWhatsYourSign[x]/x},{x,-1,1},
BaseStyle->{FontSize->18,FontFamily->"Helvetica"},PlotStyle->{{Hue[1],Thickness[0.02]},{Hue[0.66],Thickness[0.01]}}]

Functions with Conditional Definitions

In thermodynamics, x ln(x) appears frequently in expressions that involve entropy.  The variable x is restricted to 0 ≤x≤1.

XLogX[x_] := x Log[x]/;(x>0 && x ≤1)
XLogX[0] =XLogX[0.0]= Limit[xsmall Log[xsmall], xsmall->0]

XLogX[1.2]

Plot[XLogX[x]+ XLogX[1-x],{x,-1,2},PlotStyle->Thick]