HW 2 Prob 2 - Iterating Maps
This is a sample of attempts to simply obtain some periodic orbits of the logisitic map for periods 3, 4, 6 and 8. In some cases one can be guided by the plots. Every iteration of the logistic map doubles the degree of the function representing the nth iteration. This tends to give a large number of wiggles and makes root finding much harder, requiring also a larger number of digits of precision. This can be seen below. Finally, a simple bisection type method allows one to zoom in on the starting value for higher period orbits. See the use of this method in obtaining a period 8 orbit.
> restart:with(plots): Digits:=10:
Function f(x) and Kth Composition F(x) = (f^K)(x)
>
f:=x->r*x*(1-x):
> C:=proc(k,g) if k=1 then g(x) else g(C(k-1,g)) fi end:
Set K here for period you are seeking and change r until you see the desired intersections in the graph.
> K:=3: F3:=(unapply (C(K,f),x)):
> K:=4: F4:=(unapply (C(K,f),x)):
> K:=6: F6:=(unapply (C(K,f),x)):
> K:=8: F8:=(unapply (C(K,f),x)):
Period 3
> r:=3.9: plot({F3(x),f(x),x},x=0..1,color=[blue,red,green]);
Here is an animation showing when period three orbits are born as the parameter r is changed.
> r0:=3.7: r1:=3.9: N:=10: dr:=(r1-r0)/N:
> for i from 0 to N do r:=evalf(r0+i*dr): P0:=plot({f(x),x},x=0..1,color=red): T:=textplot([.2,1,cat(`r = `,convert(r,string))], align={ABOVE,RIGHT},color=blue,font=[TIMES,ROMAN,14]): P[i]:=display({plot(F3(x),x=0..1,color=blue),T,P0}): od:
> display({seq(P[i],i=1..N)},insequence=true);
Zooming in one fixed point of F3(x) for r near the birth value.
> r:=3.8285: plot({F3(x),x},x=0.158..0.161,color=[blue,red,green],view=[0..0.2,0.158..0.161]);
Candidates for period 3 - less the fixed points plus verification that we have a period three orbit.
> X:={solve(F3(x)-x)} minus {solve(f(x)-x)};
> X[1]; f(X[1]); f(f(X[1]));
> X[2]; f(X[2]); f(f(X[2]));
Period 4 Orbits
Animation sequence as above.
> r0:=3.6: r1:=4.0: N:=20: dr:=(r1-r0)/N:
> for i from 0 to N do r:=evalf(r0+i*dr): P0:=plot({f(f(x)),x},x=0..1,color=red): T:=textplot([.2,1,cat(`r = `,convert(r,string))], align={ABOVE,RIGHT},color=blue,font=[TIMES,ROMAN,14]): P[i]:=display({plot(F4(x),x=0..1,color=blue),T,P0}): od:
> display({seq(P[i],i=1..N)},insequence=true);
> r:=3.5: plot({F4(x),f(x),f(f(x)),x},x=0.0..1,color=[red,blue],view=[0..0.15,0.0..1]);
> r:=3.45: X:={solve(F4(x)-x)} minus {solve(f(f(x))-x)}: for i to nops(X) do if conjugate(X[i])=X[i] then print(i, ` `, X[i]): fi: od:
Using first root yields only period 2
> J:=1: X[J]; f(X[J]); f(f(X[J])); F3(X[J]); F4(X[J]);
Found a Period 4
> J:=2: X[J]; f(X[J]); f(f(X[J])); F3(X[J]); F4(X[J]);
Period 6
This is done in a similar fashion to the period 4 results.
> r0:=3.8: r1:=4.0: N:=20: dr:=(r1-r0)/N:
> for i from 0 to N do r:=evalf(r0+i*dr): P0:=plot({f(f(x)),x},x=0..1,color=red): T:=textplot([.2,1,cat(`r = `,convert(r,string))], align={ABOVE,RIGHT},color=blue,font=[TIMES,ROMAN,14]): P[i]:=display({plot(F6(x),x=0..1,color=blue),T,P0}): od:
> display({seq(P[i],i=1..N)},insequence=true);
> r:=3.999: plot({F6(x),f(x),F3(x),x},x=0.0..1,color=[red,blue],view=[0..0.15,0.0..1]);
> r:=3.999: plot({F6(x),x},x=0.14..0.18,color=[red,blue],view=[0..0.15,0.0..1]);
> Digits=20: r:=3.999: X:={fsolve(F6(x)-x,x=0..1)} minus {fsolve(F3(x)-x)}: for i to nops(X) do print(i, ` `, X[i]): od:
> J:=1: X[J]; f(X[J]); f(%); f(%); f(%); f(%); f(%);
> J:=2: X[J]; f(X[J]); f(%); f(%); f(%); f(%); f(%);
> J:=3: X[J]; f(X[J]); f(%); f(%); f(%); f(%); f(%);
> J:=4: X[J]; f(X[J]); f(%); f(%); f(%); f(%); f(%);
> J:=5: X[J]; f(X[J]); f(%); f(%); f(%); f(%); f(%);
> J:=6: X[J]; f(X[J]); f(%); f(%); f(%); f(%); f(%);
Period 8
> r0:=3.55: r1:=4.0: N:=20: dr:=(r1-r0)/N:
> for i from 0 to N do r:=evalf(r0+i*dr): P0:=plot({f(f(x)),x},x=0..1,color=red): T:=textplot([.2,1,cat(`r = `,convert(r,string))], align={ABOVE,RIGHT},color=blue,font=[TIMES,ROMAN,14]): P[i]:=display({plot(F8(x),x=0..1,color=blue),T,P0}): od:
> display({seq(P[i],i=1..N)},insequence=true);
> r:=3.55: plot({F8(x),F4(x),x},x=0.0..1,color=[red,blue],view=[0..0.15,0.0..1]);
> r:=3.55: plot({F8(x),F4(x),x},x=0.32..0.38,color=[red,blue,green],view=[0..0.15,0.33..0.38]);
Another approch to finding the root.
> Digits:=20: r:=3.55: z:=.35: for i to 50 do z1:=z: z2:=F8(z): z:=.5*(z1+z2): od:
> z; f(z); f(%); f(%); f(%); f(%); f(%); f(%); f(%);
>