Chapter 9
Numerics

from
Introduction to Programming with Mathematica
2nd Ed.

by
Richard J. Gaylord
Samuel N. Kamin
Paul R. Wellin

Copyright 1995 TELOS/Springer-Verlag

9.1 Types of Numbers

Code from book

3/27

{Head[1/9], FullForm[1/9]}

{Precision[7], Precision[1/9]}

Precision[92.2226669991111111111122]

2^256 + 2^1024

3/4 /. x_/y_ -> {x, y}

3/4 /. Rational[x_, y_] -> {x, y}

3/4 /. p_Rational -> {Numerator[p], Denominator[p]}

{Head[1.61803], Precision[1.61803]}

$MachinePrecision

Precision[1.23]

MachineNumberQ[1.23]

2.0^100

Precision[%]

pi = N[Pi]

pi - %

pi - 3.14159

{$MaxMachineNumber, $MinMachineNumber}

{Head[Pi], NumberQ[Pi]}

Random[Real, {-Pi, Pi}]

(* Version 3.0 *)
Random[Real, {E, Pi}]

(* Version 3.0 *)
Rationalize[Pi, .0001]

(* Version 3.0 *)
Map[NumericQ, {Pi, E, EulerGamma}]

(* Version 3.0 *)
E^Pi > Pi^E

1/Sin[N[Pi]]

1/Sin[Pi]

z = 3 + 4 I

Head[z]

z + (2 - I)

Re[z]

Im[z]

Conjugate[z]

Abs[z]

Arg[4 I]

34/21 + 2/11

3/4 + 9/4

Precision[10^100 + 1.3]

1 + 2.1 + Pi

Clear[a,b,c,d]
Expand[(a + b - 7.1)(c - d + EulerGamma)]

Precision[N[Sqrt[2], 50] * N[Sqrt[3], 80]]

a = N[2];

b = N[2^99, 30];

{Precision[a], Precision[b], Precision[a b]}

{Accuracy[1.23], Accuracy[12.5]}

Accuracy[1.23 + 12.3]

1.0 + 10^(-25)

Accuracy[%]

x = IntegerDigits[1293]

x^2

Apply[Plus, x^2]

BaseForm[18, 2]

2^^10010

Table[BaseForm[j, 16], {j, 1, 20}]

BaseForm[N[Pi, 5], 2]

RealDigits[N[Pi], 2]

RealDigits[N[Pi], N[GoldenRatio]]

Exercises

Ex.1

complexToPolar[z_] := {Abs[z], Arg[z]}

complexToPolar[3 + 3 I]

z = E^(Pi I/3);

complexToPolar[z]

Ex.3

convert[digits_List, base_:2] := 
	Fold[(base #1 + #2)&, 0, digits]

IntegerDigits[13,2]

convert[%]

IntegerDigits[129,16]

convert[%,16]

This function is essentially an implementation of Horner's method for fast polynomial multiplication.

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

Expand[%]

Ex.4

sumsOfCubes[n_Integer] := Apply[Plus, IntegerDigits[n]^3]

sumsOfCubes[16]

Ex.5

sumsOfSums[n_Integer, iter_]:=
	NestList[SumsOfCubes, n, iter]

Here we see that 4 enters into a cycle.

sumsOfSums[4,12]

In fact, it appears as if powers of 2 enter cycles.

sumsOfSums[8,12]

sumsOfSums[32,12]

Here are some other sums.

sumsOfSums[7,12]

sumsOfSums[372,12]

Ex.6

Ex.7

sumsOfPowers[n_, p_] :=
	Apply[Plus, IntegerDigits[n]^p]

sumsOfPowers[123,5]

Ex.8

Ex.9

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 (using the Convert function from Exercise 2 of Section 9.1).

convert[l, 2]

survivor[n_]:= Module[{p},
	p = RotateLeft[IntegerDigits[n, 2], 1];
	Fold[(2 #1 + #2)&, 0, p]]

survivor[100]

You could of course do the same thing without the symbol p, but it is a bit less readable.

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

survivor2[100]

9.2 Random Numbers

Code from book

Random[]

Random[Integer, {0, 100}]

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

Needs["Graphics`Graphics`"]     (* Load definition of BarChart *)

(* Load definition of Frequencies *)
Needs["Statistics`DataManipulation`"]

BarChart[Frequencies[numbers]]

Sum[i^2, {i, 1, 3}]

Table[Sum[i^2, {i, 1, n}], {n, 5}]

numb = Random[Integer, {1, 1000}]

Sum[i^2, {i, 1, numb}]

(2 numb + 1)(numb + 1)numb / 6

Clear[x]
Table[(Sum[i^2, {i,x}] == (2x + 1)(x + 1)x/6 ) /.
   x -> Random[Integer, {1, 1000}], {10}]

Exercises

Ex.1

rollEm := {Random[Integer,{1,6}], Random[Integer,{1,6}]}

rollEm

Ex.2

Here is a density plot of a 50x50 table with each site taking on a value of 1 to 256.

ListDensityPlot[
	Table[Random[Integer, {1,2^8}], {50}, {50}]]

Ex.3

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[data=NestList[linearCong, 0.0, 5000]]

ListPlot[Abs[Fourier[ data ]]]

Ex.4

Clear[quadraticCong]
quadraticCong[a_,b_,c_,m_,x_]:=
	Mod[a x^2 + b x + c, m]

quadraticCong[1,2,3,2^10,0]

The parameters must be chosen carefully to ensure sequences of reasonably long cycles.

NestList[quadraticCong[1,2,3,2^10,#]&, 0, 15]

Ex.5a

chiSquare[l_List] := 
	Module[{m=Length[l], n=Max[l]},
		Sum[(Count[l,i] - m/n)^2, {i,1,n}] / (m/n)]

Ex.5b

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

chiSquare[data]

N[%]

Notice that the statistic is quite far from 2Sqrt[n] or n. This is a particularly pathological sequence. You can see a cycle of length 63.

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

Here are those positions that contain the number 1.

Position[%,1]

Ex.6

Probably the most straightforward implementation of the middle-square method is the following (for details, see Knuth Vol.2):

middleSquare[n_, b_:10]:= 
	Mod[Floor[n^2/10^(b/2)], 10^b]

middleSquare[1234567890]

You can see that 2-digit starting values often lead to the repeating cycle 0,0,0,0...

NestList[middleSquare[#,2]&, 43, 6]

9.3 Precision and Accuracy

Code from book

{Precision[1.23], Accuracy[1.23]}

n = 12345.6789101112

{Precision[n], Accuracy[n]}

{N[Pi], Precision[N[Pi]]}

x1 = N[Pi, 2]

Sin[3.1]

Sin[x1]

InputForm[x1]

a = SetPrecision[1/3, 30]

b = SetPrecision[a, 50]

RealDigits[N[1/3], 2]

2^^.010101

FullForm[%]

N[Sqrt[2]]^200

% - Sqrt[2]^200

Table[N[Sqrt[2]]^j - Sqrt[2]^j, {j, 100, 1000, 100}]

Exercises

Ex.1

First note the exact result.

Cos[Pi/2]

Table[Cos[N[Pi/2, j]], {j, 1, 10}]

It is instructive to look at the Precision of each result.

Precision[Cos[N[Pi/2, 10]]]

So, each of these numbers is forced to be a machine-precision real. If you increased the precision past machine-precision, you would see different results.

Table[Cos[N[Pi/2, j]], {j, 10, 50, 10}]

Ex.2

(10^13 - 10^13) + 1.0/3.0

Attributes[Plus]

FullForm[Hold[(10^13 - 10^13) + 1.0/3.0]]

Tracing the computation shows what is happening in Mathematica's evaluation procedure. The next to last line below shows that Mathematica reorders the numbers placing 0.333333 in the middle of the two exact numbers. At this stage, Mathematica performs a floating point computation at machine precision using default accuracy on the first two numbers -10000000000000  and 0.333333.

Trace[(10^13 - 10^13) + 1.0/3.0]

-10000000000000 + 0.333333

FullForm[%]

Notice the decreased accuracy due to adding these two numbers (one with many digits to the right of the decimal and the other with many digits to the left of the decimal).

{Precision[%], Accuracy[%]}

9.4 Numerical Computations

Code from book

findRoot[fun_, init_, eps_] :=
   Module[{a = init, funa = fun[init]},
         While[Abs[funa] > eps,
               a = N[a - funa / fun'[a]];
               funa = fun[a]
              ];
         a]

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

secant[f, -3.1, -1.8]

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

m = {{0, 3}, {3, 0}};
b = {5, 6};

m . {2, 5/3}

pivot[S_] := Module[{p, ST1},
   ST1 = Abs[Transpose[S][[1]]];
   p = Position[ST1, Max[ST1]][[1, 1]];
   Join[{S[[p]]}, Delete[S, p]]]

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]

HilbertMatrix[n_] := Table[1/(i + j - 1), {i, n}, {j, n}]

HilbertMatrix[3] //TableForm

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

NHilbertMatrix[3] //TableForm

h50 = NHilbertMatrix[50];

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

xLS = LinearSolve[h50, b50];

xGE = solve[h50, b50];

totalerrorGE = Apply[Plus, Abs[h50.xGE - b50]]

xPP = solvePP[h50, b50];

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

Exercises

Ex.1

Clear[newton]
newton[fun_, init_]:=
   Module[{p=Precision[init], fi=fun[init], n=0},
         x = SetPrecision[init, p + 12];
         While[n++;
               Accuracy[fi] - Precision[fi] < p,
               x = x - fun[x]/fun'[x];
               fi = fun[x]
              ];
         Print["The zero is ", N[x,p], " after ",
              n, " iterations"]
   ]

f[x] := Cos[x]

newton[f, 1.5]

Ex.2

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

f[x_] := x - Sin[x]

secant[f, -.1, .1, 20]

Ex.3

Clear[x]

FindRoot[E^x - x - 1 == 0, {x, .1}]

FindRoot[E^x - x - 1 == 0, {x, .1},
	DampingFactor -> 2]

Ex.4

(# - f[#] / f'[#])&

Clear[newton]
newton[f_, x0_] := 
	FixedPoint[(# - f[#] / f'[#])&, x0]

f[x_]:=x^2-2

newton[f, N[Sqrt[3], 50]]

Ex.6

v = {1, -3, 2};

normInfinity[v_?VectorQ] := Max[Abs[v]]

normInfinity[v]

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

normInfinity[m_?MatrixQ]:=
   normInfinity[Apply[Plus, Abs[Transpose[m]]]]

normInfinity[m]

HilbertMatrix[n_Integer] :=
	Table[1/(i+j-1), {i,n},{j,n}]

normInfinity[ HilbertMatrix[10] ]

Ex.7

conditionNumber[m_?MatrixQ] :=
   normInfinity[m] normInfinity[Inverse[m]]

conditionNumber[HilbertMatrix[3]]

Map[conditionNumber[HilbertMatrix[#]]&, Range[10]]

N[%]