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 8 ± Iteration

(Mathematica code)

Functions defined in this Notebook...(page no. in book)

next...174
findRoot...177
permute...180
primes2...184
rotateRows...185
rotateRowsByS...185
sqrElement...186
elimInPlace...188
soln...189
solveI...187-188

Initialization (evaluate to turn off spell messages)

Off[General::spell]

Off[General::spell1]

Section 8.1  Newton's Method

f[x_] := x^2 - 50

next[x_, fun_] := N[x - fun[x] / fun'[x]]

a = 50;
Do[a = next[a, f], {10}]

a

a = 50;
Do[ a = next[a, f];
    Print["approximation ", i, ": ", a],
  {i, 1, 10}]

Clear[g]

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

g[0]

a = 1.0;

Do[a = next[a, g], {10}]

a

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

findRoot[f, 50, .001]

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

findRoot[f, 50, 0.001]

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

findRoot[f, 50, 0.001]

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

Section 8.2  Vectors and Matrices

permute

permute[S_, P_] := Module[{T = Table[0, {Length[S]}], p},
	Do[ p = P[[i]];
	  	T[[p]] = S[[i]],
	  {i, Length[S]}];
	  T]

Clear[a,b,c,d]

permute[{"a", "b", "c", "d"}, {3, 2, 4, 1}]

primes

primes2[n_] := Module[{primes = Table[i, {i, 1, n}], p = 2},
	Print["Initially: primes = ", primes, ", p = ", p];
	While[ p!= n+1,
 		  Do[ primes[[i]] = 1, {i, 2p, n, p} ];
		  p = p + 1;
		  While[ p != n+1 && primes[[p]]==1, p = p+1];
		  Print["In loop: primes = ", primes, ", p = ", p];
	];
	Select[primes, (# != 1)&]]

primes2[9]

Exercises

(* Note, rotateRight must first be defined for the 
   following to work. This was an exercise in this
   section.
*)

rotateRows[A_] :=
    Map[(rotateRight[A[[#]], #-1])&, Range[1, Length[A]]]

rotateRowsByS[A_, S_] /; Length[A] == Length[S] :=
    Map[(rotateRight[A[[#]], S[[#]]])&, Range[1, Length[A]]]

Section 8.3  Passing Arrays to Functions

sqrElement[L_, i_] := L[[i]] = L[[i]]^2

x = {0, 10, 20}

sqrElement[x, 3]

SetAttributes[sqrElement, HoldFirst]

Section 8.4  Gaussian Elimination Revisited

SetAttributes[{elimInPlace, solveI}, HoldFirst]

elimInPlace[S_, i_] := Module[{m, n = Length[S]},
	Do[m = S[[j, i]] / S[[i, i]];
		Do[ S[[j, k]] -= m*S[[i, k]],
		{k, i+1, n+1}],
	{j, i+1, n}]]

soln[S_, i_, xs_] := Module[{r = 0, n = Length[S]},
	Do[r = r + S[[i, k]] * xs[[k-i]], {k, i+1, n}];
	(S[[i, n+1]] - r) / S[[i, i]] ]

solveI[S_] := Module[{xs = {}, n = Length[S]},
	Do[elimInPlace[S, i],
	{i, 1, n-1}];
	Do[PrependTo[xs, soln[S, i, xs]], {i, n, 1, -1}];
	xs]

m = {{2, 3, 5, 7}, {3, 5, 7, 11}, {5, 7, 11, 13}}

solveI[m]

m1 = {{2, 3, 5}, {3, 5, 7}, {5, 7, 11}};
b1 = {7, 11, 13};
LinearSolve[m1,b1]