Introduction to Mathematica 
Programming

The Class

Violeta M. Ivanova, Ph.D., Educational Technology Consultant, MIT Academic Computing
Contact: N42-250J, 617-252-1383, violeta@mit . edu
Class materials: http://web.mit.edu/acmath/math/IntroMathematica
Web site: http://web.mit.edu/violeta/www/IAP2006

References

Mathematica 5.2 Documentation, Wolfram Research Inc., http://documents.wolfram.com/mathematica
W. C. Carter, MIT 3.016 Lecture Notebooks, http://pruffle.mit.edu/3.016, ccarter@mit.edu
J. Glynn and T. Gray, "The Beginner's Guide to Mathematica Version 4" http://www.mathware.com
S. Wolfram, "The Mathematica Book", http://documents.wolfram.com/mathematica/book/

Assignments, Patterns and Rules

Fundamental

Set[=] as in lhs = rhs.

SinCos = 
      Sin[x] Cos[2 y] + Sin[3 y] Cos[x] + Sin[x y] Cos[x + y]
Plot3D[
        SinCos, {x, -π, π}, {y, -π, π}];

The following will not work:

Plot3D[
        SinCos, {u, -π, π}, {v, -π, π}];

SetDelayed[:=] as in lhs := rhs.

SinCosFunction[x_, y_] := 
      Sin[x] Cos[2 y] + Sin[3 y] Cos[x] + Sin[x y] Cos[x + y]
Plot3D[
        SinCosFunction[u, 
          v], {u, -π, π}, {v, -π, π}];

Blank[_] or _h stands for "any variable" or, respectively, "any variable with head of type h"

numfraction[x_Integer, y_Integer] := N[(x/y)]
numfraction[x_Real, y_Real] := (x/y)
numfraction[5, 6]
numfraction[5.654, 3.45]
numfraction[6. , 5]
numfraction["ok", a + ι]

Rule[->] and ReplaceAll[/.] as in expression /. rules lhs->rhs.

Plot3D[
          SinCos, {x, -π, π}, {y, -π, π}] /. {x -> u, y -> v};

Also Important

[BlankSequence[__] or __h stands for "any sequence" or "any sequence with head of type h".

sumx2y[x_, y__] := (x^2) + y;
sumx2y[a, b, c, d]

Clear[symbol] and Clear["form"] clear values and definitions.

f1[1] = 30; 
f1[2] = 20; 
f2[x_] := x 10; 
{f1[1], f1[2], f2[1], f2[4], f2[5]}

Clear[f1]
{f1[1], f1[2], f2[1], f2[4], f2[5]}

Clear["f*"]
{f1[1], f1[2], f2[1], f2[4], f2[5]}

Clear["`Global`*"]

PatternTest[?] as in pattern?test.

squareRoot[x_?NumericQ] := (Sqrt[x])
squareRoot[2]
squareRoot[(π/2) + 3 π]
squareRoot[-2]

Condition[/;] as in pattern/; test.

rootSquare[x_Real] := (Sqrt[x]) /; x ≥ 0
rootSquare[-2]

ReplaceRepeated[//.] as in _expression//. rules performs ReplaceAll[/.] until expression stops changing.

{a, b, c, d, X, c, d, X, a, b} /. a -> 
        (X^2)

{a, b, c, d, X, c, d, X, a, b} /. {X -> Y, 
        a -> (X^2)}

{a, b, c, d, X, c, d, X, a, b} //. {X -> Y, 
        a -> (X^2)}

Increment[++], PreIncrement[++], Decrement[--] and PreDecrement[--]:

a = b = 5; 
{a++, ++b, a, b}

AddTo[+=], SubstractFrom[-=], TimesBy[*=] and DivideBy[/=]:

a *= 2; a

RuleDelayed[:>] as in lhs:>rhs is a rule that transforms lhs to rhs after rhs is evaluated.

m = 1; 
    n Cos[n x] + 
        n Sin[n x] /. n :> π m++

Unset[=.] as in lhs =. removes all rules set for lhs.

SinCos = Sin[x] Cos[2 y];
SinCos
SinCos =. 
SinCos

Tests and Logical Operations

Relational Operators

Equal[==] and SameQ[===]:

a==b

a === b

a = 2; b = 2. ; a == b

Unequal[!=] and UnsameQ[=!=]:

Sin[Pi/4] != Cos[Pi/4]

Clear[a, b]
a =!= b

Less[<], Greater[>], LessEqual[<=] and GreaterEqual[>=]:

Sin[x] > 1

Sin[(π/2)] ≥ 1

NumberQ, NumericQ, IntegerQ, PrimeQ, etc.:

NumberQ[π]

Clear[x]; NumberQ[x] /. x -> 5

x = 5; NumberQ[x]

NumberQ[ι]

NumericQ[π]
NumericQ[ι]

Positive, Negative, NonPositive and NonNegative:

Positive[5]

Positive[Abs[-5]]

NonNegative[0]

WhatSign[_?Positive] := 1; 
WhatSign[_?Negative] := -1; 
WhatSign[0] = 0; 
WhatSign /@ {1, -1, π, E, 0}

Logical Operators

Logical operators are almost always used together with relational operators and/or test functions.

And[&&] as in test1 && test2 && test3 ...:

π > 3 && π < 4

Or[||] as in test1 || test2 || test3 ...:

π ≥ 3 || π < 0 || π > 100

Not[!] as in !test:

Not[PrimeQ[π]]

Program Flow Control

Flow Control Functions

If as in If[condition, iftrue, iffalse] or If[condition, iftrue, iffalse, ifneither]:

WhatSignTwo[x_?NumericQ] := If[x > 0, 1, -1, 0]
WhatSignTwo /@ {1, -1, π, E, 0, y}

Do as in Do[expr, {ntimes}] or Do[expr, {i, imax}] or Do[expr, {i, imin, imax}, {j, jmin, jmax}]:

r = 1; Do[r = (1/(1 + r)), {20}]; r

While as in While[testistrue, body]:

n = 0; 
    While[n < 10, Print["for n = ", n, "  n! = ", n!]; 
      n += 2]

Note also above the use of CompoundExpression[;] as in expr1; expr2; expr3; ...

Switch as in Switch[expr, test1, value1, test2, value2, ...]:

\(squareRoot[x_] := 
      Switch[Head[x], 
Integer, N[\@x], 
Real, \@x, 

        Symbol, \@x]\)
squareRoot /@ {3, 3. , π, ι}

Which as in Which[test1, value1, test2, value2, ...]:

Plot[
        Which[x < 0, (1/2) - ((x^2)/2), x ≥ 0. , 
          (1/2) + 0.05 x], {x, - .5, 3}];

The Importance of Being "Listable"

For as in For[start, test, increment, body]:

list1 = {x,  y,  x y,  x + y}
list2 = {};
For[i = 1, i ≤ Length[list1], ++i, 
      AppendTo[list2, 
        Sin[list1[[
            i]]]]]
list2

Programmers beware: since most of Mathematica's functions are Listable, you almost never need to use a For loop!

list3 = Cos[list1]

If a function is not Listable, you usually can make it so:

SetAttributes[WhatSignTwo,  Listable];
WhatSignTwo[{1, -1, π, E, 0, y}]

Variable Scoping

Module[{x, y, z, ...}, expression] specifies that the occurrences of x, y, z, ... in expression should be treated as local variables.

data = Table[Random[], {5}]

\(randomStatistics[n_] := 
      Module[{data}, 
							\ data = 
          Table[Random[], {n}]; 
							\ Print["Maximum = ", 
          Max[data]]; 
							\ Print["Minimum = ", Min[data]]; 
							\ Print["Average = ", 
          Mean[data]]; 
							\ Print["data in function = ", \ data]]\)

randomStatistics[5]
data

Remove[value] or Remove["form"] removes (unprotected) variables completely. All variables in Mathematica's packages are protected.

?SinCos
Remove["Global`*"]
?SinCos

$Packages shows a list of the packages loaded into the current session.

$Packages

Needs["context`package"] loads the appropriate file for the package in the specified context, if it is not already in $Packages.  If symbols from the package are already being used, you need to use Remove before you can load the package.

Helix[a_, b_, tmin_, tmax_] := 
      ParametricPlot3D[{a Sin[t], a Cos[t], b t}, {t, tmin, 
          tmax}]
Helix[3, 2, -π, π]
?Helix

Needs["Graphics`Shapes`"]

Remove[Helix]
?Helix

File Input and Output

Get[<<] and Needs["context`file`"] are mainly used for packages:

Needs["Miscellaneous`CityData`"]

Import["filename.ext"] and Export["filename.ext"] are for data and images in various formats. See supported formats in the Built-In Functions->Input and Output->Import and Export directory in the Help Browser.

Other useful file output  functions to know are OpenWrite, OpenAppend, Write, WriteString and Close:

SetAttributes[CityDistance, Listable]
citylist = {"Sofia", "London", "Sydney", "Rome", "San \
Francisco", "New York"}
bostondistances = Round[CityDistance["Boston", citylist]]
filetowrite = OpenWrite["distances"]
WriteString[
      filetowrite, "Distances from Boston in [km] to: \n"]
n = 0;
Do[n++; 
Write[filetowrite, 
        citylist[[n]], 
        bostondistances[[
          n]]], 
{Length[
          citylist]}]
Close[filetowrite]

!! distances

Programming in Mathematica

Algorithm Development Environment

Computer languages vary from the lowest level machine language where everything is expressed in terms of ones and zeros, to high-level programming languages like Java and C++; scripting languages like Perl; markup languages like HTML; etc.
Mathematica is an algorithm development environment for symbolic and numeric computation, which has built-in abilities of a high-level programming and scripting language. 
Programming in Mathematica means developing algorithms to solve mathematical problems, which includes writing procedures, using built-in functions and creating new functions, and using patterns and rules.

Procedural Programming

Procedural programming in Mathematica is similar to what many programming languages do well. For example, using relational and logical operators and flow control statements (see above) is typical of procedural programming. Useful functions for procedural programming in Mathematica include: Module, Table, If, Do, While, Switch; the logical operators && and ||; the relational operators ==, !=, >, etc.; and functions such as Map[/@], Apply[/@@], and Print[].

a = b; 
    Do[Print["In step ", i, ": ", a *= a + i], {i, 3}]

Map[f,  {π, 2 π, 3 π}]
g /@ {π, (π/2), (π/3)}

Apply[f,  {π, 2 π, 3 π}]
g @@ {π, (π/2), (π/3)}

Functional Programming

Functional programming in Mathematica is also similar to that in other programming languages: it allows you to use libraries of built-in functions as well as write your own functions. The differences are that Mathematica (1) has so many more built-in functions; (2) can work with higher accuracy; (3) can work symbolically as well as numerically; and (4) is so high-level a language that often a one-line command is enough to perform complex operations.

Series[LegendreP[n, z], {z, 0, 3}]

NestList[(((1 + #))^3) &, x + (y^2),  5]

In the expression above, (((1 + #1))^3)& is a pure function and # is its first argument.

Pattern-Based and Rule-Based Programming

Pattern matching and rule-based programming distinguish Mathematica from other symbolic mathematics programs and are fundamental to almost anything you do in Mathematica.

Cos[x] + Cos[x]^2 + Sin[x] + Sin[x]^2

Cos[x] + Cos[x]^2 + Sin[x] + Sin[x]^2  /.  
      Sin[x]^2 + Cos[x]^2 -> 1

Cos[x] + Cos[x]^2 + Sin[x] + Sin[x]^2 /. 
      Sin[n_]^2 + Cos[n_]^2 -> 1

Cos[a]^2 + Cos[x]^2 + Sin[a]^2 + Sin[y]^2 /. 
      Sin[n_]^2 + Cos[n_]^2 -> 1

Clear[a,  b,  c]
g1 /@ {π, (π/2), (π/3)} /. g1 -> 
        Sin
g2 @@ {a Pi, (π/2) b, (π/3) c} /. g2 -> 
        Plus

log[(E^n)_] = n
x + log[(E^Sin)[x + y]]

arrange[{x__, a_, y__}, a_] := {a, x, y, x, y}
arrange[{1, 2, 3, B, x, y, z}, B]

What is a Mathematica program?

A Mathematica program can be as short as a line, using a built-in function ...

FourierTransform[(t^4) (E^(-(t^2))), t, w]

... or several nested functions.

Plot[
      FourierSinTransform[(E^(-t)), t, w], {w, 0, 
        10}]

If the function you want does not exist in Mathematica, you can write it yourself, using procedures, built-in functions, relational and logical operators, flow control statements, patterns, and rules ...

FlatNumericList[alist_List] := 
      Module[{l, n, newlist, templist}, 
        newlist = {}; 
templist = Flatten[alist]; n = 0; 
        Do[n++; 
          If[NumericQ[
              templist[[
                n]]], 
newlist = 
              Join[newlist, {templist[[
                    n]]}]], 
{Length[
              templist]}]; Return[newlist]]

... and then you can use that function in a one-line command.

Blist = 
      FlatNumericList[{{1, 2, 3}, {"ok", π, ϵ}}]

templist
newlist
alist

You can load packages and use your own functions when you are defining more new functions ...

\(ComputeStatistics[alist_List] := 
      Module[{numericlist}, 
Needs["Statistics`DescriptiveStatistics`"]; 
        numericlist = 
          FlatNumericList[alist]; 
Print["Mean: ", 
          Mean[numericlist]]; 
        Print["Standard Deviation: ", StandardDeviation[numericlist]]; 
        Print["Median: ", Median[numericlist]]; 
        Print["Variance: ", Variance[numericlist]]; 
        If[Length[Flatten[alist]] > Length[numericlist], 
          Print["Excluding non-numeric values in input."], 
          Print["Including all input values."]]]\)

... and then write another one-line command that calls your second function that calls your first function ...

ComputeStatistics[{{π, -2 π, (π/3), 2 π}, {1, 
          (3/5), 4, 6}}]

... and you can also use all of the above in longer programs.

Needs["Statistics`ContinuousDistributions`"]
expdist = ExponentialDistribution[0.5]
exparray = RandomArray[expdist, 5000];
ComputeStatistics[exparray];

Mathematica Packages

More Standard Add-On Packages

The Add-Ons & Links directory in the Mathematica  Help Browser  contains several subdirectories with all available packages that come with Mathematica. To use a function in any of these packages, you need to load the package first. For example, to use function ScatterPlot3D, you need to load the package Graphics3D in Standard Packages->Graphics.
Besides the Standard Packages, there are also many advanced packages in the Add-Ons & Links directory. For example, the MathLink package provides functions to create an interface between Mathematica and other computer programs. The GUIKit package provides functions to create Graphical User Interface for your Mathematica programs. The Palettes in Mathematica are essentially such GUIs for math typesetting and other short cuts to Mathematica functions.  
You need to know a lot more programming that what you learned in this session to use some of the advanced packages. For example, open the file slider.nb to see the Mathematica code to create a simple slider. If you have done this kind of programming in, say, Java, this code will not be too difficult for you to understand.

Packages from the MathSource

BasicBiochemData3 package courtesy of Robert Alberty, Professor Emeritus, MIT Department of Chemistry.

Free packages are available from the MathSource web site at http://www.wolfram.com/mathsource/ and from other sources.
For example, here is a package written by an MIT Professor: 
Mathematical Functions for Thermodynamic Properties of Biochemical Reactants, Robert Alberty, MIT
To use it: 
1. Download the notebook file BasicBiochemData3.nb and the Mathematica source file BasicBiochemData3.m to your working directory. In the notebook file, you can read what the package is about. The actual package is the Mathematica source file with extension .m. This is how all other packages that are already in Mathematica look as well. 
2. To use the downloaded package, load the BasicBiochemData3.m file. Note that if you open that file in a text editor, you can find in the header of the file the exact command for loading the package.

Directory[]

SetDirectory["IAP2006/IntroMathematica/Programming\
"]

Needs["BiochemThermo`BasicBiochemData3`"]

$Packages

Creating New Packages

To create a package, use the BeginPackage[] and EndPackage[] commands. See, for example, the file MyFunctions.nb. To save a package, use File->Save As Special->Package Format. The actual package is then generated automatically by Mathematica in an .m file, such as, for example, the MyFunctions.m file. Once you have created a package, you can load it and use it in another Mathematica session.

Needs["MyPackages`MyFunctions`"]

$Packages