Introduction to Programming 
with 
Mathematica

By
Richard J. Gaylord 
(gaylord@ux1.cso.uiuc.edu)

Samuel N. Kamin 
(kamin@cs.uiuc.edu)

Paul R. Wellin 
(wellin@sonoma.edu)

Copyright © 1993 by TELOS/Springer-Verlag

Chapter 9 – Numerics

(Mathematica code)

Functions defined in this Notebook

josephus...201
J3 (Josephus problem)...
linearCong...204
ChiSquare...205
newton...209
secant...211
solvePP...213-218
vectorNorm...218
matrixNorm...218

Initialization (evaluate to turn off spell messages)

Off[General::spell]

Off[General::spell1]

Needs["Graphics`Graphics`"]

Needs["Statistics`DataManipulation`"]

Section 9.1  Types of Numbers

r = Table[Random[Integer, {0,9}], {1000}];

BarChart[Frequencies[r]]; (* Make sure you read in the Data
							  Manipulation package first *)

Exercises

Josephus problem

The Josephus number of 100 is 73.

josephus[n_] :=
	Nest[Rest[RotateLeft[#, 1]]&, Range[n], n-1]

josephus[100]

First let's put 100 in base 2,

BaseForm[100,2]

and then get a list of its digits:

IntegerDigits[100,2]

Do a binary shift of one unit,

l = RotateLeft[IntegerDigits[100, 2], 1]

and convert back from base 2 to base 10:

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

J3[n_]:= Module[{p},
	p = Reverse[RotateLeft[IntegerDigits[n, 2], 1]];
	Sum[p[[j]] 2^(j-1),{j, 1, Length[p]}]
	]

J3[100]

J3[10^9] //Timing

Here though is the ultimate functional approach to a binary shift.  If you look carefully, you will see Horner's method for multiplying polynomials lurking behind the scenes.  Fold  is a nice little function!

J4[n_Integer]:= Fold[(2 #1 + #2)&, 0, 
	RotateLeft[IntegerDigits[n, 2], 1]]

J4[10^9] //Timing

Section 9.2  Random Numbers

Linear congruential method

Clear[x,b,m]
linearCong[x_] = Mod[b x + 1, m]

m = 100; b = 15; NestList[linearCong, 5, 10]

m = 381
b = 15;
data = NestList[linearCong, 0, 60]

m = 381;
b = 15;
ListPlot[NestList[linearCong, 0.0, 5000]];

ListPlot[Abs[Fourier[NestList[linearCong, 0.0, 5000]]]];

m = 2^16;
b = 27421;
ListPlot[NestList[linearCong, 0.0, 5000]];

ListPlot[Abs[Fourier[NestList[linearCong,0.,5000]]]];

Exercises

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

ChiSquare[l_List,r_]:= Module[{n=Length[l], i=1, s=0},
	While[i <= r,
		s += ((Count[l,i] - n/r)^2)/(n/r);
		i++;
	];
	{s, N[s]}
]

ChiSquare[d, 100]

Section 9.4  Numerical Computations

Newton's method revisited

newton[fun_, init_]:= 
   Module[{p = Precision[init], fi = fun[init]},
         x = SetPrecision[init, p + 12];
         While[Accuracy[fi] - Precision[fi] < p,
               x = x - fun[x]/fun'[x];
               fi = fun[x]
              ];
         N[x, p]
        ]

f[x_] := x^2 - 2

newton[f, 1.0]

x0 = N[Sqrt[3], 40];

newton[f, x0]

Secant method

secant[f_, a_, b_]:= Module[{x1 = a, x2 = b},
	While[Abs[ f[x2] ] > 10^(-10),
		  df = (f[x2] - f[x1])/(x2 - x1);
		  {x1, x2} = {x2, x2 - f[x2]/df}
	];
	x2]

f[x_] := Abs[x + 3]

FindRoot[f[x]==0, {x, .1}]

Solve[f[x]==0, x]

secant[f, -3.1, -1.8]

Gaussian elimination revisited

solvePP[m_, b_] := 
	Module[{elimx1, pivot, subtractE1, solve},
		elimx1[T_] := Map[subtractE1[T[[1]], #]&, Rest[T]];
		pivot[Q_] := Module[{p, ST1},
			ST1 = Abs[Transpose[Q][[1]]];
			p = Position[ST1, Max[ST1]][[1, 1]];
			Join[{Q[[p]]}, Delete[Q, p]]];

		subtractE1[E1_, Ei_] := 
			Module[{w = Ei[[1]]/E1[[1]]},
				Module[{newE1 = w * Rest[E1]},
					Rest[Ei] - newE1]];

		solve[{{a11_, b1_}}] := {b1/a11};
		
		solve[S_] :=
			Module[{S1 = pivot[S]},
				Module[{E1 = First[S1],
					x2toxn = solve[elimx1[S1]]},
					Module[{b1 = Last[E1],
					a11 = First[E1],
					a12toa1n = Drop[Rest[E1], -1]},
					Join[{(b1 - a12toa1n . x2toxn)/a11},
				x2toxn]]]];
	
		solve[Transpose[Append[Transpose[m], b]]]
	]

m = {{0, 3}, {3, 0}};

b = {5, 6};

solvePP[m, b]

hilbert[n_] := Table[1.0/N[i + j - 1], {i, n}, {j, n}]

hilbert[3] //MatrixForm

h50 = hilbert[50];

b50 = Table[Random[], {50}];

xLS = LinearSolve[h50, b50];

totalerrorLS = Apply[Plus, Abs[h50.xLS - b50]]

xPP = solvePP[h50, b50]

totalerrorPP = Apply[Plus, Abs[h50.xPP - b50]]

Exercises (vector and matrix norms)

v = {1, -3, 2};

Abs[v]

m = {{1,2,3}, {1,0,2}, {2,-3,2}}

vectorNorm[v_VectorQ] := Max[Apply[Plus, Abs[v]]];

vectorNorm[v]

matrixNorm[m_MatrixQ] := Max[Map[vectorNorm, m]]

matrixNorm[hilbert[3]]

matrixNorm[m]