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