# Lesson 5. Iteration, Fixed and Periodic Points
# ================================
#
# We continue to study iteration of the following
# family of functions:
> g:= x -> a*x - x^2;
> a:= 2.8;
> stair:= x -> plot([[x,x], [x,g(x)], [g(x), g(x)]]);
# As we saw last time, 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:
> with(plots):
> 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];
# Now b and c are roots of g(g(x))=x. What are
# tho
# se 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 for that function?
# 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 attractors. 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 1
# the cycle is a repeller.
#
# Next we'll try a = 3.8.
> a:= 3.8;
# Here are the fixed points, and the periodic points of period 2.
> solve(g(x)=x, x);
> s:=solve(g(g(x))=x, x);
# Is the cycle an attractor?
> D(g)(s[3])*D(g)(s[4]);
# We start near the 2-cycle and calculate some points.
> x[0]:= 3.38:
> for nn from 1 to 89 do
> x[nn]:= g(x[nn-1]) od:
# Here's the graph of these points.
> plot([seq([nn, x[nn]], nn=0..89)]);
# And now a staircase diagram. The first few steps
# show what happens near the 2-cycle. Try it also for nn = 16 .. 89 to see the later
# development.
> p:= plot({g(x), x}, x = 0 .. 3.8):
> display({p, \
seq(stair(x[nn]), nn = 0 .. 15) });
# It's near the 2-cycle for a while, but moves away after
# a while. No apparent pattern develops. This, it seems, is an example of chaos. Actually,
# there are lots of cycles here, but it seems that all of them are unstable.
#
# Is there a 3-cycle?
> plot((g@@3)(x)-x, x=0 .. a);\
# Apparently not: the only fixed points of g@@3 are the fixed points of g.
# Try for a 5-cycle:
#
> plot((g@@5)(x)-x, x=0..a);\
# There are apparently 12 fixed points of g@@5: the 2 fixed points of g, and two 5-cycles.
# "fsolve", by the way, has trouble with this function (which is really quite complicated).
#
> eq:=(g@@5)(x)-x=0;\
# It's even worse if you try to "simplify" it by expanding.
#
> expand(eq);\
# What does "fsolve" do with it?
#
> s:= fsolve(eq, x);\
# It didn't find all the solutions, only a few, and they are not very accurate.
# (Actually, it seems that the PC and X versions get different answers here,
# and the X version is very inaccurate).
#
> subs(x=s[2], eq);
>
