Programming
with
Mathematica

Mathematica Days

University of Michigan
Ann Arbor, Michigan
11 Sept 93

and

Boston University
Boston, Massachusetts
18 Sept 93

Paul Wellin
Department of Mathematics
Sonoma State University
Rohnert Park, CA  94928  USA
wellin@sonoma.edu

1  Examples

1.1  Mean and Variance

mean[p_] := Apply[Plus, p] / Length[p]

variance[p_] := mean[ (p - mean[p])^2 ]

data = Table[Random[Integer, {1, 100}], {1000}];

mean[data] //N

variance[data] //N

1.2  Harmonic Numbers

harmonic[n_] := Sum[1/j, {j, 1, n}]

harmonic[100]

N[%]

harmonic[1000]

N[%]

harmonic[10^4]

nHarmonic[n_] := NSum[1/j, {j, 1, n}]

nHarmonic[10^4]

nHarmonic[10^12]

nHarmonic[10^15]

1.3  Perfect Numbers

Perfect numbers are defined as those numbers equal to the sum of their "proper" divisors.

Divisors[6]

Drop[Divisors[6], -1]

Apply[Plus, %]

perfectQ[n_] := Apply[Plus, Divisors[n]] == 2n

perfectQ[6]

perfectQ[7]

Trace[perfectQ[7]]

TracePrint[perfectQ[7]]

perfect[n_] := Select[Range[n], perfectQ]

perfect[500]

The following function makes an assumption about perfect numbers which is as yet, unproven: that there are no odd perfect numbers. It is known there are none below
10^300.

perfect2[n_?EvenQ, m_] := Select[Range[n,m,2], perfectQ]

perfect2[3, 500]

perfect2[2, 500]

1.4  Root Plotting

plotRoots[poly_, z_, opts___] := 
	ListPlot[{Re[z], Im[z]} /. NSolve[poly == 0, z],
		PlotStyle -> {PointSize[.02], RGBColor[1,0,0]},
		opts]

p5 = x^5 - 1

plotRoots[p5, x]

p[n_] := Apply[Plus, Table[z^j, {j, 0, n}]]

p[2]

Do[
	plotRoots[p[i], z,
		Axes->False,
		AspectRatio->Automatic],
{i,20}];

1.5  CardDeck

deck = Flatten[Outer[List, {c,d,h,s},
			Join[Range[2, 10], {J,Q,K,A}]],
		1]

removeRand[lis_] := 
	Delete[lis, Random[Integer, {1, Length[lis]}]]

deal[n_] := Complement[deck, Nest[removeRand, deck, n]]

deal[5]

deal[n_] := 
(deck = Flatten[Outer[List, {c,d,h,s},
			Join[Range[2, 10], {J,Q,K,A}]],1];
 removeRand[lis_] := 
		Delete[lis, Random[Integer, {1, Length[lis]}]];
 Complement[deck, Nest[removeRand, deck, n]]
)

deal[5]

deal[n_] := Module[{deck, removeRand},
	deck = Flatten[Outer[List, {c,d,h,s},
			Join[Range[2, 10], {J,Q,K,A}]],1];
 	removeRand[lis_] := 
		Delete[lis, Random[Integer, {1, Length[lis]}]];
 	Complement[deck, Nest[removeRand, deck, n]]
	]

deal[5]

1.6  Koch Snowflake  (Michael Trott)

startdreieck = { {{-1,0},{1,0}} , 
                 {{1,0},{0,Sqrt[3]}},
                 {{0,Sqrt[3]},{-1,0}}  } // N

EckedranDeterministisch[line_List]:=
Module[{d,norm,senk},
 d = line[[2]]-line[[1]];
 norm = Sqrt[d.d];
 senk = {d[[2]],-d[[1]]}; 
 { {line[[1]], line[[1]]+d/3},
   {line[[1]] + d/3,
    line[[1]] + d/2 + senk/norm Sqrt[3] norm/6},
   {line[[1]] + d/2 + senk/norm Sqrt[3] norm/6,
    line[[1]] + 2d/3},
   {line[[1]] + 2d/3,line[[2]]}} //N ]

st[1] = startdreieck;
Do[st[i] = Flatten[ EckedranDeterministisch /@ 
                                 st[i-1], 1], 
  {i,2,5}]

Show[Graphics[
	Table[{Hue[0.8 (i-1)/5],
           Polygon[(i (#-{0,1/Sqrt[3]}) )&/@ 
                                (#[[1]]& /@ 
                       Partition[Drop[Append[
                       hg=Flatten[st[i],1],
                                 hg[[1]] ],1],2]) ]},
    {i,4,1,-1}]
    ],
    AspectRatio->Automatic]

2  Recursive Programming

2.1  Simple sum function

sum[n_] := sum[n-1] + n
sum[0] = 0

sum[2]

Table[sum[j], {j,20}]

Table[Sum[j, {j,1,i}], {i,20}]

Map[sum, Range[20]]

2.2  Fibonacci numbers and dynamic programming

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

fib[1] = fib[2] = 1;

?fib

Table[fib[j], {j, 15}]

?fib

fib2[n_] := fib2[n] = fib2[n-1] + fib2[n-2]
fib2[1] = fib2[2] = 1;

?fib2

Table[fib2[j], {j,12}]

?fib2

fib2[500]

$RecursionLimit

$RecursionLimit = 1000

Clear[fib2]
fib2[n_] := fib2[n] = fib2[n-1] + fib2[n-2]
fib2[1] = fib2[2] = 1;

fib2[500]

2.3  Maxima

maxima[{}] := {}
maxima[{x_, r___}] := 
	Join[{x}, Select[maxima[{r}], (# > x)&]]

maxima[{4, 2, 7, 3, 4, 9, 14, 11, 17}]

3  Functional Programming

3.1  Simple sum functions

Timing[ Sum[k, {k, 10000}] ]

?FoldList

?Fold

FoldList[Plus, 0, Range[10]]

Timing[ Fold[Plus, 0, Range[10000]]]

Timing[ Apply[Plus, Range[10000]] ]

Timing[ Plus @@ Range[10000] ]

3.2  Base conversion

l = IntegerDigits[12509, 2]

r = Reverse[l];
Sum[r[[j]] 2^(j-1), {j, 1, Length[l]}]

baseConvert[list_List, base_:2]:=
	Fold[(base #1 + #2)&, 0, list]

baseConvert[l]

baseConvert[{a, b, c, d, e}, x]

Expand[%]

3.3  Maxima

maxima[x_] := Union[Rest[FoldList[Max, 0, x]]]

maxima[{4, 2, 7, 3, 4, 9, 14, 11, 17}]

4  Iterative Programming

4.1  Iteration by "Do"

sum = 0;
Do[ sum += i, {i,10000} ];
sum

4.2  Iteration by "For"

sum = 0;
For[ sum = 0; i = 1, 
     	i <= 10000, 
     	i++,
     	sum += i];
sum

4.3  Iteration by "While"

sum = 0;
i = 1;
While[ i <= 10000, sum += (i++)];
sum

5  Evaluation of Expressions

5.1  Order of evaluation

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

fib[1] = fib[2] = 1;

?fib

Trace[fib[3]]

Specific rules before general rules

- Evaluate head of expression
- Evaluate each element in turn
- Apply transformations associated with Attributes (Orderless, Listable, and Flat)
- Apply any user-supplied definitions
- Apply any built-in defintions
- Evaluate the result

x = 2;
Trace[ a x^2 + b x + c ]

Clear[x]
Log[x y]

(
Unprotect[Log];
Log[a_ b_] := Log[a] + Log[b];
Protect[Log];
)

Log[x y]

??Log

5.2  Immediate and delayed evaluation

x = Random[]

y := Random[]

?x

?y

Table[x,{10}]

Table[y,{10}]

6  Graphics Programming

6.1  Simple closed paths

coords = Table[{Random[],Random[]},{10}]

points = Map[Point, coords]

Show[Graphics[{PointSize[.02], points}]]

Show[Graphics[
	{PointSize[.025], RGBColor[1,0,0], Map[Point,coords]}],
	Axes->Automatic];

lines = Line[coords]

Show[Graphics[{PointSize[.02], 
	{RGBColor[1,0,0], points}, lines}]];

path = Line[AppendTo[coords, First[coords]]];
Show[Graphics[{PointSize[.02], 
	{RGBColor[1,0,0],points}, path}]];

base = coords[[Random[Integer,{1,Length[coords]}]]]

The following function finds the angle between the horizontal axis and point "a", with point "b" being taken as the vertex of the angle.

angle[a_List, b_List] := Apply[ArcTan, (b-a)]

remain = DeleteCases[coords, base];

s = Sort[remain, (angle[base, #1] <= angle[base, #2])&]

p = Join[{base}, s, {base}]

path = Line[p];

Show[Graphics[{PointSize[.02], 
	{RGBColor[1,0,0],points}, path}]]

SimpleClosedPath

simpleClosedPath[l_] := Module[{points, base, sorted, path},
	points = {PointSize[.02], RGBColor[1,0,0], Map[Point,l]};
	(* base = l[[ Random[Integer, {1, Length[l]}] ]]; *)
	base = Last[Sort[l,(#2[[2]] > #1[[2]])&]];
	remain = DeleteCases[l,base];
	angle[a_,b_]:= Apply[ArcTan, (b - a)];
	sorted = Sort[remain, (angle[base, #1] <= angle[base, #2])&];
	path = Line[Join[{base},sorted,{base}]];
	
	Show[Graphics[{{RGBColor[1,0,0],points}, path}]]
]

data = Table[{Random[],Random[]},{25}];
simpleClosedPath[data];

data = Table[{Random[],Random[]},{100}];
simpleClosedPath[data];

6.2  Random audio walk

(* Equal tempered C major scale *)
cmajor = Table[N[261.62558 2^(j/12)],{j,0,12}]

We want to do a random walk across our list of tones. We will allow steps of 0, 1, 2, -1,-2, steps.

step20 = Table[Random[Integer, {-2,2}], {2}];

Need to add 1, so that we get no Part[expr, 0] which would return the Head, not a position.

FoldList[Plus, 5, step20]

pos = Mod[FoldList[Plus, 5, step20], 11] + 1

brown = cmajor[[pos]]

tones[freq_] := Play[Sin[2 Pi t freq], {t, 0, .5}]

Map[tones, brown];

7  Additional sources of information on programming in Mathematica.

"Programming in Mathematica", R. Maeder, Addison-Wesley Publishing Co., 1992.

"Introduction to Programming with Mathematica",  R. Gaylord, S. Kamin, 
  P. Wellin, TELOS/Springer-Verlag, 1993.

"Mathematica: A Practical Approach", N. Blachman, Prentice-Hall, 1992.

"The Mathematica Journal, Miller Freeman, San Francisco.

"Mathematica in Education", TELOS/Springer-Verlag. To order a subscription or receive a complimentary copy of Mathematica in Education, send a message to:
  
  	Paul Wellin, editor
  	Mathematica in Education
  	Department of Mathematics
  	Sonoma State University
  	Rohnert Park, CA  94928  USA
  	
  	email: MathInEd@groucho.sonoma.edu
  	fax:  707-664-2505