Chapter 8
Iteration

from
Introduction to Programming with Mathematica
2nd Ed.

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

Copyright 1995 TELOS/Springer-Verlag

8.1 Newton's Method

Graphics

Figure 8.1

$DefaultFont = {"Helvetica", 9};

ftmp[x_] := ((x*1000+10)^3.2)/1000 - 30
x1 = 0.05;
y1 = ftmp[x1];
p0 = {x1, y1};
deltax = .01;
p1 = p0-{deltax, deltax * ftmp'[x1]};
p2 = p0+{deltax, deltax * ftmp'[x1]};
x0 = N[x1-ftmp[x1]/ftmp'[x1]];
q0 = {x0, 0};
newton1 =
  Graphics[Join[InputForm[
      Plot[ftmp[x], {x, 0, .06}, Axes -> True, Ticks -> None,
           DisplayFunction -> Identity]][[1,1]],
      {{Dashing[{.01, .01}], Line[{p0, q0}],
                   Line[{p0, {x1, 0}}]}},
      {Line[{p1, p2}],
               Line[{{x0, -10}, {x0, -50}}],
               Line[{{x1, -10}, {x1, -50}}],
               Text["x - dx", {x0, -75}],
               Text["x", {x1, -75}],
               Line[{{x0, -30}, {x1, -30}}],
               Text["dx", {(x1+x0)/2, -55}],
               Text["dx = f[x]/(f[x]/dx)", {.02, 400}],
               Text["= f[x]/f'[x]", {.02, 300}],
               Line[{{x1+.001, y1}, {x1+.002, y1}}],
           Text["f[x]", {x1+.005, y1}]}],
  Axes->True, Ticks->None]

Show[newton1, DisplayFunction ->$DisplayFunction]

Figure 8.2

x0 = 50;
y0 = f[x0];
x1 = N[x0-f[x0]/f'[x0]];
y1 = f[x1];
x2 = N[x1-f[x1]/f'[x1]];
y2 = f[x2];
x3 = N[x2-f[x2]/f'[x2]];
y3 = f[x3];
deltax = 10.0;
p0 = {x0, y0};
p1 = {x1, y1};
p2 = {x2, y2};
p3 = {x3, y3};
p0a = {x1, 0};
p0b = p0+{deltax, deltax * f'[x0]};
p1a = {x2, 0};
p1b = p1+{deltax, deltax * f'[x1]};
p2a = {x3, 0};
p2b = p2+{deltax, deltax * f'[x2]};
newton2 =
   Graphics[
      Join[ InputForm[
             Plot[f[x], {x, 0, 60}, Axes -> True, Ticks -> None,
              DisplayFunction ->Identity]][[1, 1]],
        {Line[{p0a, p0b}],
               Line[{p1a, p1b}],
               Line[{p2a, p2b}],
               {AbsoluteThickness[0.5], Line[{p0, {x0, -100}}]},
               {AbsoluteThickness[0.5], Line[{p1, {x1, -100}}]},
               {AbsoluteThickness[0.5], Line[{p2, {x2, -100}}]},
               {AbsoluteThickness[0.5], Line[{p3, {x3, -100}}]},
               Text["50", {x0, -200}],
               Text["25.5", {x1, -200}],
               Text["13.7", {x2+1, -200}],
               Text["8.7", {x3-.5, -200}]}],
    Axes -> True, Ticks->None];

Show[newton2, DisplayFunction->$DisplayFunction]

Code from book

FindRoot[x^2 - 50 == 0, {x, 50}]

FindRoot[x^(1/3) == 0, {x, .1}]

f[x_] := x^2 - 50

a0 = 50;

a1 = N[a0 - f[a0] / f'[a0]]

a2 = N[a1 - f[a1] / f'[a1]]

a3 = N[a2 - f[a2] / f'[a2]]

a4 = N[a3 - f[a3] / f'[a3]]

a5 = N[a4 - f[a4] / f'[a4]]

a = 50;

a = N[a - f[a] / f'[a]]

a = N[a - f[a] / f'[a]]

a = N[a - f[a] / f'[a]]

a = N[a - f[a] / f'[a]]

a = N[a - f[a] / f'[a]]

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

a = 50;

a = next[a, f]

a = next[a, f]

a = next[a, f]

a = next[a, f]

a = next[a, f]

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

a

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

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

g[0]

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

a

a = 1.0;
Do[a = next[a, g], {25}]

a

epsilon = .001;
a = 50;

While[Abs[f[a]] > epsilon, a = next[a, f]]

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]

Clear[findRoot]
findRoot[fun_, x_, eps_] :=
   (While[Abs[fun[x]] > eps,
          x = N[x - fun[x] / fun'[x]]
         ];
    x)

findRoot[Sin, .1, .01]

TracePrint[findRoot[Sin, .1, .01], findRoot]

Exercise solutions

Ex.1

next[x_, fun_] := N[x - fun[x] / fun'[x]]
a = 50;
b = 60;
Do[a = next[a, (#^2 - 50)&];
   b = next[b, (#^2 - 60)&],
  {10}]
Print["square root of 50 = ", a]
Print["square root of 60 = ", b]

Ex.2

Notice that the function to compute square root of one number
using Nest is very concise, but the one to compute the square
root of two numbers simultaneously using Nest is verbose.
Since Nest applies just one function to the value computed
at the previous step, we need to put the two values into a
list and then write a function that extracts each element in
the list, applies the two functions to the two elements, and
then puts them back into a list.

nestSqrt[a_] := Nest[next[#, (#^2-a)&]&, a, 10]

nestSqrt2[a_, b_] :=
   Module[{},
      nexta[x_] := next[x, (#^2-a)&]; 
      nextb[x_] := next[x, (#^2-b)&]; 
      Nest[{nexta[#[[1]]], nextb[#[[2]]]}&,
           {a, b},
           10]]

Ex.3

We need to place the two expressions that were in the body of
the Do into a list.  Try copying the body of the Do exactly as
above and see what happens.

a = 50;
b = 60;
Table[{a = next[a, (#^2 - 50)&],
       b = next[b, (#^2 - 60)&]},
     {10}]

Ex.4

fib[n_] :=
  Module[{prev = 0,  (* fib[0] *)
          this = 1,  (* fib[1] *)
          next},
    (* At the ith step, prev = fib[i-1], this = fib[i] *)
    Do[ next = prev+this;  (* fib[i+1] *)
        prev = this;       (* prev = fib[i] *)
        this = next,       (* this = fib[i+1] *)
      {n}];
    prev]

Map[fib, Range[20]]

Ex.5

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

Ex.6

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

f[x_]:=x^2-50

findRootList[f,7.0, 10^(-12)]

Ex.7

We go back to the version of findRoot on page~\pageref{page.findroot}, and add multiple initial values.

findRootList[fun_, inits_, eps_] :=
Module[{a = inits},
   While[Min[Abs[Map[fun, a]]] > eps,
   a = Map[N[# - fun[#] / fun'[#]]&, a]
   ];
   Select[a, (Min[Abs[Map[fun, a]]] == Abs[fun[#]])&]]

findRootList[(#^2-50)&, {25, 50, 75, 100}, .001]

Ex.8

bisect[f_, a_, b_, eps_] :=
   Module[{low = Min[a,b],
           high = Max[a,b],
           mid = N[(a+b)/2],
           fofmid = N[f[(a+b)/2]]},
        While[Abs[fofmid] > eps,
              If[fofmid < 0, low = mid, high = mid];
              mid = N[(low+high)/2];
              fofmid = N[f[mid]]];
        mid]

bisect[(#^2-50)&, 0, 50, .001]

8.2 Vectors and Matrices

Code from book

L = {0, 1, 2, 3, 4}

L[[1]] = 10

L

A = {{1, 2, 3}, {4, 5, 6}}

A[[2, 3]] = 20

A

B = A

B[[1, 2]] = 30

B

A

A[[2, 1]] = 40

B

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

Exercises

Ex.1a

reverse[V_] :=
  Module[{rV = Table[0, {Length[V]}]},
     Do[ rV[[i]] = V[[Length[V]-i+1]],
      {i, 1, Length[V]}];
     rV]

Ex.1b

transpose[A_] :=
  Module[{tA = Table[Table[0, {Length[A]}], {Length[A[[1]]]}]},
    Do[ tA[[i, j]] = A[[j, i]],
      {i, 1, Length[A[[1]]]},
      {j, 1, Length[A]}];
    tA]

Ex.1c

The key to this problem is to use the Mod operator to compute
   the target address for any item from V.  That is, the element
   V[i] must move to, *roughly speaking*, position (n+i) mod length[V].
   The "roughly speaking" is due to the fact that the Mod operator
   returns values in the range 0, ..., length[V]-1, whereas vectors
   are indexed by values 1, ..., length[V].  This causes a little
   trickiness in this problem.

rotateRight[V_, n_] :=
  Module[{rV = Table[0, {Length[V]}]},
    Do[ rV[[1+Mod[n+i-1, Length[V]]]] = V[[i]],
      {i, 1, Length[V]}];
    rV]

Ex.1d

rotateRows[A_] :=
  Module[{tA = Table[0, {Length[A]}]},
    Do[ tA[[i]] = rotateRight[A[[i]], i],
      {i, 1, Length[A]}];
    tA]

Ex.1e

rotateRowsByS[A_, S_] :=
  Module[{tA = Table[0, {Length[A]}]},
    Do[ tA[[i]] = rotateRight[A[[i]], S[[i]]],
      {i, 1, Length[A]}];
    tA]

Ex.1f

compress[L_, B_] :=
  Module[{cL = Table[0, {Count[B, True]}],
          last = 1},
    Do[ If[B[[i]],
           cL[[last]] = L[[i]];  last = last+1,
           NULL (* do nothing *)],
      {i, 1, Length[B]}];
    cL]

Ex.1g

This is just a matter of taking the code for primes shown in the book and replacing the last line by code that is very similar to the compress function just given.

primes[n_] := Module[{primes = Table[i, {i, 1, n}], p = 2},
        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];
        ];
        Module[{justprimes = Table[0, {n-Count[primes, 1]}],
                last = 1},
          Do[ If[primes[[i]] != 1,
                 justprimes[[last]] = primes[[i]];  last = last+1,
                 NULL (* do nothing *)],
            {i, 1, Length[primes]}];
          justprimes]]

8.3 Passing Arrays to Functions

Code from book

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

x = {0, 10, 20}

sqrElement[x, 3]

TracePrint[sqrElement[x, 3], sqrElement]

SetAttributes[sqrElement, HoldFirst]

sqrElement[x, 3]

x

sqrElement[{5, 10, 15, 20}, 2]

Exercises

Ex.1

SetAttributes[swap, HoldFirst]
swap[V_, i_, j_] :=
   Module[{temp = V[[i]]},
     V[[i]] = V[[j]];
     V[[j]] = temp;
     Null]

Ex.2

SetAttributes[reverseInPlace, HoldFirst]

reverseInPlace[V_] :=
   Module[{l = Length[V]},
     Do[swap[V, i, l-i+1], {i, 1, l/2}]]

8.4 Gaussian Elimination Revisited

Code from book

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

solveI[S_] :=
   Module[{xs = {}, n = Length[S]},
         Do[elimInPlace[S, i];
           Print[MatrixForm[S]];
	    Print["\n"],
           {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]

Exercises

Ex.1

The basic idea is simple: before calling elimInPlace, call pivot, which will swap rows if necessary.

SetAttributes[{elimInPlace, solveI, pivot}, HoldFirst]

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

elimInPlace and soln do not change

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

\mv{pivot} is similar to \mv{pivot} written for Exercise 1 on page~148,
except that it swaps the rows {\it in place} instead of returning
an entire new array.  (Note that since rows are the elements
of an array, they can be swapped using swap from the previous
exercise set.)  However, what makes it trickier than the \mv{pivot}
from that exercise is that the argument is the entire array,
together with an index $i$ indicating what part of the array we
are interested in.  Thus, \mv{pivotrow} below gives the pivot row,
considering only rows from $i$ to $n$.  So if \mv{pivotrow = 1}, it really means that $i$ itself is the pivot row ({\it i.e.}, no pivoting is needed);
if \mv{pivotrow > 1}, then row $i$ should be swapped with row $i + \mbox{\mv{pivotrow}} - 1$.

pivot[Q_, i_] :=
   Module[{p, pivotrow, Qsuffix, ST1},
       Qsuffix = Drop[Q, i-1];  (* look only in remaining rows *)
       ST1 = Transpose[Qsuffix][[i]]; (* ith elements of rows *)
       p = Position[ST1, x_/; x!= 0];
       If[p == {},
          Print["Matrix is singular"]; Q,
          pivotrow = p[[1]][[1]];
          If[pivotrow != 1, swap[Q, i, i+pivotrow-1]]
       ]]

m = {{1, 1, 1, 1}, {1, 1, 2, 2}, {1, 2, 2, 1}};

solveI[m]

Ex.2

We return to the non-pivoting version of solveI for this problem.
The change in ludecompI from solveI is simply to omit the
final part of computing the solution vector. However, the
changes to elimInPlace are more significant.

SetAttributes[{elimInPlace, ludecompI}, HoldFirst]

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

For elimInPlace, note first that S is now a square matrix
instead of a rectangular one, since it does not include the
vector B.  Also, we had previously used the multipliers only
to eliminate rows, but now we need to need to save them in
the array.  These are the two changes with respect to the
version of elimInPlace from problem 1.

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

Ex.3

SetAttributes[{elimInPlace, ludecompI}, HoldFirst]

ludecompI[S_] := Module[{xs = {}, n = Length[S]},
        Do[pivot[S, i];
           elimInPlace[S, i],
         {i, 1, n-1}];
        S]

For elimInPlace, note first that S is now a square matrix
instead of a rectangular one, since it does not include the
vector B.  Also, we had previously used the multipliers only
to eliminate rows, but now we need to need to save them in
the array.  These are the two changes with respect to the
version of elimInPlace from problem 1.

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

pivot[Q_, i_] :=
   Module[{p, pivotrow, Qsuffix, ST1},
       Qsuffix = Drop[Q, i-1];  (* look only in remaining rows *)
       ST1 = Transpose[Qsuffix][[i]]; (* ith elements of rows *)
       p = Position[ST1, x_/; x!= 0];
       If[p == {},
          Print["Matrix is singular"]; Q,
          pivotrow = p[[1]][[1]];
          If[pivotrow != 1, swap[Q, i, i+pivotrow-1]]
       ]]

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

ludecompI[m]

m = {{1, 1, 1}, {1, 1, 2}, {1, 2, 2}};

ludecompI[m]