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