NERCOMP • Amherst College • 22 Feb 2000

Hero's method

The problem

• Problem:  To estimate (Sqrt[c]) for  c > 0
• Method:  iteration

Function g to iterate

Define g

g[x_]  :=  (1/2) (x  +  (c/x))

Not really a function definition, but a rewrite rule.

Form x_  is a pattern.

g[1]

c  is just a symbol.

g[g[1]]

Multiple arguments

Above, c is global to function g of 1 argument.  Could use function of 2 arguments:

gg[cc_,  x_]  :=  (1/2) (x  +  (cc/x))
gg[2,  1]

Advanced topic:  Everything is an expression.

FullForm[x  +  c/x]

FullForm[ 
      Hold[ g[x_]  :=  (1/2) (x  +  
                (c/x))] ]

FullForm[  g[
          x_]  :=  (1/2) (x  +  (c/x))]

Always get a result, maybe the symbol Null.

Advanced topic:  Head

FullForm[3/2]

Head[3/2]

Advanced topic:  Pure functions

gg[2,  #] &[1]

Using /. for Rule

g[1]

g[1]  /.  c  ->  2

gg[cc,  x0]  /.  {x0  ->  1,  cc  ->  2}

Iterating g

Fix value of c to estimate (Sqrt[2])

c  =  2

c  =  2;

Set ( = )  vs.  SetDelayed ( := )

Variations on evaluating

?g

g[1]        
g[g[1]]

Prefix input form:

g@g[1]

N[g@g[1]]

Postfix input form:

g@g[1]  //  N

Display more digits:

N[g@g[1],  $MachinePrecision]

NumberForm[N[g@g[1]],  $MachinePrecision]

Using lists:

{ N[g[1]],  N[g@g[1]] }

N has the attribute Listable:

N[ {g[1],  g@g[1]} ]

The initial estimate: a floating-point approximate real number

x0  =  1.  ;

Iterating—at last!

{x0,  g[x0],  g@g[x0],  g@g@g[x0],  
      g@ g@g@g[x0]}

List vs. compound expression:

(x0;  g[x0];  g@g[x0];  
      g@g@g[x0];   g@g@g@g[x0])

Reusing previous values:

x0
g[%]
g[%]

"Functional" method for iteration

Nest[g,  x0,  4]

NestList[g,  x0,  4]

Advanced topic:  Lists as arguments

Nest[g,  {1. ,  9. },   4]

Nest[g,  x0,   {2,  4}]

Map[Nest[g,  x0,  #] &,  {2,  4}]

Nest[g,  x0,  #] &  /@  {2,  4}

Fixed points

FixedPoint[g,  x0]

NumberForm[% ,  $MachinePrecision]

FixedPointList[g,  x0]  //  
      NumberForm[#,  $MachinePrecision] &

How good are the estimates?

NumberForm[N[(Sqrt[2])],  $MachinePrecision]
NumberForm[FixedPoint[g,  x0],  $MachinePrecision]

"Arbitrary precision" approximate numbers

Precision[1. ]

Precision[1. ]  ==  $MachinePrecision

1.`50

Precision[%]

x0  =  1.`50;

FixedPoint[g,  x0]

Precision[%]

You cannot get "something for nothing":

NumberForm[FixedPoint[g,  1. ],  50]

Precision[FixedPoint[g,  1. ]]

Procedural programs

Use machine precision again:

x0  =  1. ;

\(hero[c_, \ n_, \ x0_]\  := \ 
	
      Module[{x\  = \ x0}, 
		Do[
          x\  = \ g[x], \ {n}]; 
		x
	]\)

hero[2,  4,  1. ]

\(herolist[c_, \ n_, \ x0_]\  := \ 
      Module[{x\  = \ {x0}}, 
		Do[
          x\  = \ Append[x, \ g[Last[x]]], \ {n}]; x]\)

herolist[2,  4,  1. ]

Advanced topic:  Type checking

Clear[herolist]
\(herolist[c_, \ n_Integer, \ x0_?Positive]\  := \ 
      Module[{x\  = \ {x0}}, 
		Do[
          x\  = \ Append[x, \ g[Last[x]]], \ {n}]; x]\)

herolist[2,  4.1,  1. ]

herolist[2,  4,  -1. ]

herolist[2,  4,  1. ]

Advanced topic:  Defaults

herolist[c_,  n_]  :=  herolist[c,  n,  1. ]
herolist[2,  4]
herolist[2,  4,  0.5]

No change to original definition!

Advanced topic:  Options

Clear[herolist]
Options[herolist]  =  {Iterations  ->  10,  
          Start  ->  1. };

\(herolist[c_, \ opts___?OptionQ]\  := \ 
      Module[{x, \ x0, \ 
          n}, 
		n\  = \ \(Iterations\  /. \ {opts}\)\  /. \ 
            Options[herolist]; 
		x0\  = \ \(Start /. \ {opts}\)\  /. \ 
            Options[herolist]; 
		x\  = \ {x0}; 
		Do[
          x\  = \ Append[x, \ g[Last[x]]], \ {n}]; x]\)

herolist[2]

herolist[2,  Start  ->  0.5]

herolist[2, 
      Start  ->  0.5,   Iterations  ->  4]

Recursive programs

With If test

\(rhero[n_, \ 
        x0_]\  := \ \ Module[{}, 
		If[n\  == \ 0, \ x0, \ 
          g@rhero[n - 1, \ x0]]]\)

rhero[4,  1. ]
rhero[0,  1. ]

With pattern-matching and rewrite rules

rh[0,  x_]  =  x;
rh[n_,  x_]  :=  g@rh[n - 1,  x]

rh[4,  1. ]
rh[0,  1. ]

?rh

Advanced topic:  "Dynamic programming"

drh[0,  x_]  =  x;
drh[n_,  x_]  :=  drh[n,  x]  =  
        g@drh[n - 1,  x]

?drh

drh[1,  1. ]

?drh

drh[2,  1. ]

?drh

Solve exactly

Solve[g[x]  ==  x,  x]

Visualizing fixed points

Plot[{g[x],  x},  {x,  -2,  2}]

Graphics objects

bothGraphs  =  Plot[{g[x],  x},  {x,  -2,  2}]

cstFuncGraph  =  Plot[(Sqrt[2]),  {x,  -2,  2}]

Show[bothGraphs,  cstFuncGraph]

Enhancing graphs

\(\(myg\  = \ 
        Plot[{g[x], \ x}, \ {x, \ \(-2\), \ 2}, \ 
          PlotPoints\  -> \ 80, \ 
          PlotStyle\  -> \ {RGBColor[1, 0, 0], \ RGBColor[0, 0, 1]}, 
          Frame\  -> \ True, \ 
          FrameLabel\  -> \ \ \ {\*""(1/2)(x + (2/x))"", ""}];\)\)

Cobweb diagram

Based upon an original program by Frank Wattenberg.

Advanced topic:  Usage message

\(\(cobweb::usage\  = 
\ \ \ \*""cobweb[f, a, b, n, x0] creates a cobweb diagram for function f with initial value x0 and final value (f^n)(x0); the graph is plotted on the interval [a, b]."";\)\)

?cobweb

The program

\(bounce[f_, \ x_]\  := \ 
      Module[{fx, \ ffx}, \ fx\  = \ f[x]; \ 
        ffx\  = \ 
          f[fx]; 
	{Line[{{x, \ f[x]}, {f[x], \ f[x]}, {f[x], \ 
                f[f[x]]}}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Point[{f[x], \ f[f[x]]}]}]\)

\(web[f_, \ n_, \ x0_]\  := \ 
      Module[{x\  = \ x0, \ steps\  = \ {}}, \ 
	
        Do[\ steps\  = \ Append[steps, \ bounce[f, \ x]]; \ 
          x\  = \ f[x], \ {n}]; steps]\)

cobweb[f_,  lo_,  hi_,  n_,  x0_]  :=  
      Show[Plot[{f[x],  x},  {x,  lo,  hi},  
          PlotRange  ->  {0,  hi},  
          DisplayFunction  ->  Identity],  
        Graphics[{PointSize[0.02],  RGBColor[1, 0, 0], Dashing[{0.01}],  
            Point[{x0,  f[x0]}],  
            web[f,  n,  x0]}], 
AspectRatio  ->  
          Automatic,  Axes  ->  Automatic,  
        DisplayFunction  ->  $DisplayFunction]

Run it

g[x_]  :=  0.5 (x  +  (2/x))

cobweb[g,  0,  10,  4,  9];

Animate  it

cobwebEvolve[f_,  lo_,  hi_,  n_,  x0_]  :=  
      Module[{i},  
        Table[cobweb[f,  lo,  hi,  i,  x0], {i,  0,  n}]]

cobwebEvolve[g,  0,  10,  7,  9];

Main notebook

Copyright © 2000 by Murray Eisenberg.  All rights reserved.