# 1.6 Lesson 6

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

```