1.4 Lesson 4

# Lesson 4. Iteration and Fixed Points
# ==========================
# Newton's method is a particular case of an iteration method.  In general, iteration deals with a sequence 
# defined by x[n+1] = g(x[n]) for some function g.
# Suppose the sequence x[n] converges to some value c as n -> infinity, and g is continuous at c.  Then we 
# can take limits on both sides of the equation and get c = g(c).  So c will be a "fixed point" of the function 
# g.
# You might use an iteration scheme to find a fixed point of g.  Or you could use it to solve an equation 
# that's equivalent to g(c)=c.  Newton's method is an example of this, with g(x) = x - f(x)/D(f)(x), because 
# g(c)=c is equivalent to f(c)=0 (assuming D(f)(c) is not 0).  Now for iteration to work, it must happen that 
# if
# x[0] is close to the fixed point c, the sequence x[n] will converge to c.  This may or may not happen, 
# depending on g. 
# Here is a very simple function, depending on a parameter "a".  We'll use several different values of the 
# parameter and starting points, starting with  a=2.8 and x[0]=1.5.
> g:= x -> a*x - x^2;
> a:= 2.8; x[0]:= 1.5;
# We'll calculate the next 40 terms in the sequence.
> for nn from 1 to 40 do 
>       x[nn]:= g(x[nn-1]) od:
# I used ":" above so we don't see the individual x[n]'s.  Instead we can plot them.
# We want to make a list of the points [n, x[n]], but we don't want to have to type these in individually.  
# The "seq" command makes a sequence of objects, one for each value of an index variable.  Its syntax is
#   seq(expression, index_variable = 
#                                      start_value .. end_value)
# Just as with a do loop, the seq command leaves its index variable with a value, so I use a double letter for 
# the index variable.  
# Enclosing the seq command in square brackets makes a list.
> pts:= [ seq( [nn, x[nn]], nn = 0 .. 40) ]:
> plot(pts);
# It looks very much like the sequence converges to 1.8. This makes sense since 1.8 is a fixed point.  
> g(1.8);
# If you try other starting points x[0], you will find that the sequence converges to 1.8 as long as  0 < x[0] < 
# 2.8.  If x[0] < 0 or x[0] > 2.8, the limit is -infinity.  There is another fixed point at x = 0, but
#  the only ways to get x[n] to converge to 0 are  x[0] = 0 (in which case all x[n] are 0)  or x[0] = 2.8 (in 
# which case x[1] and all later x[n] are 0). If x[0] is close to, but not exactly, 0, x[n] will be close to 0 for a 
# while, but eventually moves away toward either 1.8 or -infinity.
# Of the two fixed points, 1.8 is said to be an "attractor", while 0 is a "repeller".
# A fixed point r of g is stable if we can guarantee that x[n] always stays close to r by taking the starting 
# point x[0] sufficiently close to r.  Otherwise it is unstable.
# A stable fixed point r is an  attractor if x[n] -> r as n -> infinity whenever x[0] is sufficiently close to r.
# An unstable fixed point r is a  repeller if whenever the starting point x[0] is sufficiently close to, but not 
# equal to, r, there will be some x[n] which are at least a certain distance away from r.
#  For us the main classifications of fixed points are attractors and repellers.  We won't consider any 
# examples that are neither one nor the other.
# The main way to tell whether a fixed point r  is an attractor or a repeller is by looking at g'(r).
# A fixed point r of a differentiable function g is an attractor if  |g'(r)| < 1.
# It is a repeller if |g'(r)| > 1. 
# If |g'(r)| = 1, it could be an attractor or a repeller or neither; we need more information to decide which.
# This can be seen as follows: consider the slope of the secant line to the graph of g through  [r,r] and [x[0], 
# g(x[0])] = [x[0], x[1]], where  x[0] is close to (but  not exactly) r.  This slope is (x[1] - r)/(x[0] - r).  It will 
# be arbitrarily  close to g'(r) if x[0] is close enough to r, since the derivative is the limit of the slopes of 
# secant lines.  In particular, if |g'(r)| < 1 we will also have  |(x[1] - r)/(x[0] - r)| < 1, i.e. 
#  |x[1] - r| < |x[0] - r|, if  x[0] is close enough to r.  Thus x[1] will be even closer to r in this case.  That is 
# in itself not enough to make the sequence converge to r, but we can improve the result: 
# if |g'(r)| < 1, take some number p with |g'(r)| < p < 1.  If  x[0] is close enough to r, we will have
#  |(x[1] - r)/(x[0] - r)| < p, i.e.  |x[1] - r| < p |x[0] - r|.
# Since  x[1] will also be "close enough'' to r (if not exactly equal to r),  
#     |x[2] - r| < p |x[1] - r| < p^2 |x[0] - r|.  
# This will go on forever: at each iteration, the distance to r decreases at least by a factor of p, and as 
# n -> infinity  the distance goes to 0 and x[n] converges to r. Thus in this case the fixed point is an
# attractor.
#  On the other hand, if |g'(r)| > 1 a similar argument shows that when x[0] is close enough to (but not equal 
# to) r the distance from r increases at each iteration until x[n] is no longer ``close enough'' to r.  In this 
# case, the fixed point is a repeller.
#  One way to visualize what goes on is to plot a "cobweb" or "staircase" diagram.  We start with the 
# graphs of y = g(x) and y = x.  The intersections of these are the fixed points.  Now, given any x[0], we 
# plot the following sequence of points: [x[0], x[0]], [x[0], x[1]], [x[1], x[1]], [x[1],x[2]], ...
# This starts on the diagonal y = x, then alternately moves vertically to the curve y = g(x) and horizontally 
# to the diagonal.  The basic building block is a function I'll call "stair" that takes x and plots the lines 
# from [x,x] to [x, g(x)] to [g(x),g(x)].
> stair:= x -> plot([[x,x], [x,g(x)], [g(x), g(x)]]);\

# Here we plot y=x and y=g(x) together.  Instead of plotting a single function, we can plot a set of
# functions (enclosed in braces {}).  On a colour terminal, the curves will be in different colours.  On some 
# of the black-and-white terminals, the second colour may be hard to see, so we add the third input
# colour=white  to the plot command.  You can do this conveniently by removing the "#".
> p:= plot({x, g(x)}, x= 1 .. 2.5
> # , colour=white
> ):
> with(plots):
> display({ p, stair(x[0]), 
>           stair(x[1]), stair(x[2]),           stair(x[3])});
# Clearly the fixed point 1.8 is an attractor.
> D(g)(1.8);
# Each x[n] is closer to 1.8 than the previous one, and the "cobweb" spirals inwards toward [1.8, 1.8].
# On the other hand, what should happen when x[n] starts near the other fixed point, 0?
> D(g)(0);
> p:= plot({x, g(x)}, x=-1 .. 1 
> #, colour=white
> ):
> x[0]:= 0.05:
> for nn from 1 to 3 do x[nn]:= g(x[nn-1]) od:
> display({p, stair(x[0]), stair(x[1]), stair(x[2]), stair(x[3])});
# Now let's change the value of the parameter "a".
> a:= 3.1;
# What are the fixed points, and are they attractors or repellers?
> solve(g(x)=x, x);
> D(g)(0), D(g)(2.1);
# This time both fixed points are repellers.  So what will happen?  We'll start near 2.1 and draw a picture 
# with 50 steps, using "seq" to avoid typing in each individually.
> x[0]:= 2.05:
> for nn from 1 to 49 do 
>   x[nn]:= g(x[nn-1]) od:
> p:= plot({g(x), x}, x = 1.5 .. 2.5):
> display({p, 
>  seq(stair(x[nn]), nn = 0 .. 49) });
# The "cobweb" seems to be approaching a rectangle (actually a square).  The x values at the two sides of 
# the  square, say b and c, have the property that g(b)=c and g(c)=b.  They constitute a "cycle" of period 2.  
# It's almost true for x[48] and x[49].
> g(x[48])=x[49], g(x[49])=x[48]; 
# In that case b and c are roots of g(g(x))=x.  What are those roots?
> s:= solve(g(g(x))=x, x);
# Two of those roots are the fixed points, and the others form the cycle (in this case the last two).
# A point on a cycle of period 2 is a fixed point of x -> g(g(x)).  Are these attractors?
# First of all, note that "@" is the composition sign in Maple, so the function x -> g(g(x) can be written as 
# g@g.
> D(g@g)(s[3]), D(g@g)(s[4]);
# They are stable.  Is it a coincidence that these two derivatives are the same?  Think of the chain rule:
# we should have 
> D(g@g)(s[3]) = D(g)(s[3])*D(g)(g(s[3]));
# And of course, g(s[3]) = s[4], so this is
> D(g)(s[3]) * D(g)(s[4]);
# And the same holds true for D(g@g)(s[4]) since g(s[4]) = s[3].
# If x[0] is close enough to one of the points on the cycle (say b), then as n -> infinity the even-numbered 
# x[n] approach b and the odd-numbered ones approach c. We say the cycle is an attractor for g.
# In general, a periodic point p of period m is a point such that g@...@g(p) = p (with m g's, but not with 
# any smaller positive number of g's).  It is a fixed point of the function g@...@g (which can also be 
# written in Maple as "g@@m").  Together with the points (g@@k)(p) for k = 1 to m-1, which are also 
# periodic points of period m, it forms a cycle of period m.  To test stability of the cycle, we look at 
# D(g@@m)(p), which is the product of the values of D(g) at the m points of the cycle.  If this has absolute 
# value less than 1 the cycle is an attractor, while if its absolute value is greater than 1the cycle is a repeller.
# Next we'll try a = 3.8.
> a:= 3.8;
> solve(g(x)=x, x);
> s:=solve(g(g(x))=x, x);
> D(g)(s[3])*D(g)(s[4]);
> x[0]:= 3.38:
> for nn from 1 to 89 do 
>   x[nn]:= g(x[nn-1]) od:
> plot([seq([nn, x[nn]], nn=0..89)]);
> p:= plot({g(x), x}, x = 0 .. 3.8):
> display({p, 
>  seq(stair(x[nn]), nn = 0 .. 10) });
> display({p, 
>  seq(stair(x[nn]), nn = 10 .. 89)});