# Lesson 6. Cycles and Chaos # ===================== # # Last time we tried to use fsolve to find a 5-cycle for a=3.8, but the results were rather strange. > g:= x -> a*x-x^2; a:= 3.8: 2 g := x -> a x - x a := 3.8 > s:=fsolve((g@@5)(x)-x, x, 0 .. 3.8); s := 0, .8323303583, .8994824740, 1.691755494, 2.152259899, 2.470081536, 2.608964680, 2.800000000, 3.107369082, 3.285007042, 3.546364943, 3.566634226 # Only three roots show up instead of 12. Are they accurate? We calculate (g@@5)(x)-x for each. > seq((g@@5)(s[jj])-s[jj], jj=1..12); -8 -8 -8 -7 -8 -7 -7 0, .17*10 , .60*10 , .6*10 , -.25*10 , .4*10 , .12*10 , 0, .12*10 , -8 -7 -8 -.4*10 , .13*10 , .2*10 # I don't really know why fsolve has trouble with this (other than the fact that this is a very # complicated polynomial). # ( I later found out that using the fraction a = 38/10 instead of the floating-point a = 3.8 would produce all # 12 roots, with good accuracy, although computation # time would be significant. ) # Here again is the graph of (g@@5)(x)-x. The 5-cycles are the points (besides the fixed points 0 and 2.8) # where it crosses the x axis. > plot((g@@5)(x)-x, x=0..3.8);\ # There should be one point near x=1.7. Here's a closeup of the graph near that point. > plot((g@@5)(x) - x, x = 1.6 .. 1.8); # If fsolve doesn't work, maybe we can use Newton's method directly. Here's the derivative of g. > gp:= x -> a - 2*x; gp := x -> a - 2 x # And here's the Newton formula using f(x) = (g@@5)(x)-x. > newt:= unapply(x - ((g@@5)(x)-x)/(gp(x)*gp(g(x))*gp(g(g(x)))*gp(g(g(g(x))))*gp(g(g(g(g(x)))))-1), x); 2 2 2 newt := x -> x - (791.35168 x - 208.5136 x - 54.872 (3.8 x - x ) 2 2 2 2 - 14.44 (14.44 x - 3.8 x - (3.8 x - x ) ) - 3.8 2 2 2 2 2 2 2 (54.872 x - 14.44 x - 3.8 (3.8 x - x ) - (14.44 x - 3.8 x - (3.8 x - x ) ) ) 2 2 2 ^2 - (208.5136 x - 54.872 x - 14.44 (3.8 x - x ) 2 2 2 2 - 3.8 (14.44 x - 3.8 x - (3.8 x - x ) ) - 2 2 2 2 2 2 2 (54.872 x - 14.44 x - 3.8 (3.8 x - x ) - (14.44 x - 3.8 x - (3.8 x - x ) ) ) / 2 2 2 2 ^2)^2) / (gp(x) gp(3.8 x - x ) gp(14.44 x - 3.8 x - (3.8 x - x ) ) gp( / 2 2 2 2 2 2 2 54.872 x - 14.44 x - 3.8 (3.8 x - x ) - (14.44 x - 3.8 x - (3.8 x - x ) ) ) 2 2 2 gp(208.5136 x - 54.872 x - 14.44 (3.8 x - x ) 2 2 2 2 - 3.8 (14.44 x - 3.8 x - (3.8 x - x ) ) - 2 2 2 2 2 2 2 (54.872 x - 14.44 x - 3.8 (3.8 x - x ) - (14.44 x - 3.8 x - (3.8 x - x ) ) ) ^2) - 1) # We start with 1.6 and do a few iterations of "newt". To ensure accuracy, it's prudent to increase Digits. > x[0]:= 1.6; Digits:= 14: x[0] := 1.6 Digits := 14 > for jj from 1 to 6 do x[jj]:= newt(x[jj-1]) od; x[1] := 1.6873351200278 x[2] := 1.6917045340808 x[3] := 1.6917554864697 x[4] := 1.6917554935152 x[5] := 1.6917554935176 x[6] := 1.6917554935210 # Except for the last couple of digits, which are affected by roundoff error, we have the solution. # We calculate the rest of the 5-cycle. Note that p[6] is (except for the last couple of digits) the same # as p[1]. > p[1]:= x[6]; p[1] := 1.6917554935210 > for jj from 2 to 6 do p[jj]:= g(p[jj-1]) od; p[2] := 3.5666342255213 p[3] := .832330358321 p[4] := 2.4700815362370 p[5] := 3.2850070420417 p[6] := 1.691755493494 # Presumably the 5-cycle is repelling. We check by calculating the derivative of g@@5 at one of its # points. > gp(p[1])*gp(p[2])*gp(p[3])*gp(p[4])*gp(p[5]); -9.3624536117925 # For what "a" will there be an attracting 5-cycle? Certainly if the point a/2 is on a cycle, that cycle is # attracting, because g'(a/2)=0. > a:= 'a'; > (g@@5)(a/2); > plot((g@@5)(a/2)-a/2, a=2 .. 4); > fsolve((g@@5)(a/2)=a/2, a = 2 .. 4); > as:= "; > seq((g@@5)(as[jj]/2) - as[jj]/2, jj=1..4); > a:= as[2]; > x[0]:= a/2; > for jj from 1 to 10 do x[jj]:= g(x[jj-1]) od; # What is really going on? How do things change so much for different values of "a"? # > for ii from 1 to 50 do > a:= 2 + ii*0.04; > x:= a/2; > for jj from 1 to 50 do x:= g(x) od: > x[0]:= x; > for jj from 1 to 50 do x[jj]:= g(x[jj-1]) od: > plist[ii]:= seq([a,x[jj]], jj=0..50); > od: > plot([seq(plist[ii], ii=1..50)], style=POINT, symbol=POINT); > a:= 3.8008; > x[0]:= 2.1; > for jj from 1 to 99 do x[jj]:= g(x[jj-1]); od: > plot([seq([nn,x[nn]], nn=0..99)]); > for jj from 1 to 300 do x[0]:= g(x[0]) od: > for jj from 1 to 99 do x[jj]:= g(x[jj-1 > ]) od: > plot([seq([nn,x[nn]], nn=0..99)]); > stair:= x -> plot([[x,x],[x,g(x)],[g(x),g(x)]]); > with(plots): > display({seq(stair(x[nn]), nn=0..99), > plot({x,g(x)}, x=0..a)}); >