# 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)});
