MATH332 — Assignment 3, 2015

Exercise 1: Logistic Map, λ = 2 (super-attracting)

Define the map; evaluate it's derivative:

(Subscript[Q, λ_])[x_]:=λ x(1-x);
D[(Subscript[Q, λ])[x],x]//Expand//Factor

Plot a graph, showing the most salient features:

Plot[{(Subscript[Q, 2])[x],x},{x,0,1},Frame->True,AspectRatio->1,PlotRange->All,PlotStyle->{{Hue[0],Thickness[0.008`]},{Hue[(2/3)]}},GridLines->{{0,(1/2),1},{0,(1/2),1}},Background->White]

a

Locate the fixed-points:

(Subscript[fixedpoints, 2])=Solve[(Subscript[Q, 2])[x]==x]

Check their stability by calculating the derivative (slope of the graph) at each fixed-point:

(D[(Subscript[Q, 2])[x],x]//Expand)/.(Subscript[fixedpoints, 2])

The fixed point at x=0 is unstable, while the one at x=(1/2) is stable and attracting; indeed super-attracting, since the slope is 0 there.

b

Change coordinates:

(Subscript[x, n_]):=(1/2)(1-(Subscript[y, n]))

Solve[(Subscript[x, n])==x,(Subscript[y, n])]

Rewrite the action of (Subscript[Q, 2])in terms of (Subscript[y, n]):

1-2(Subscript[Q, 2])[(Subscript[x, n])]//Expand

Hence (Subscript[y, n+1])=  ((Subscript[y, n])^2), for every n∈ℕ, so that   (Subscript[y, n])=  ((((Subscript[y, 0])))^(2^n)).

c

With |(Subscript[y, 0])|<1, this gives convergence to 0 at a rate that is clearly faster than exponential.
That is,  ((((Subscript[y, 0])))^(2^n))<  ((Subscript[y, 0])^(2n)) for all  n>1.

This gives a closed-form expression in terms of  x :

(Subscript[x, n])/.{(Subscript[y, n])-> (((1-2x))^(2^n))}

?ListPlot

Manipulate[ListPlot[{Table[(A^(2^n)),{n,numterms}],Table[(B^(2n)),{n,numterms}]},PlotStyle->{Hue[0],Hue[2/3]},Joined->True,PlotRange->{0,(10^upper)}],{{A,.5},0,1},{{B,.5},0,1},{{numterms,20},Range[10,50,5]},{upper,Range[-5,5]}]

d

If |(Subscript[y, 0])|=1, then...

(Subscript[x, 0])/.{{(Subscript[y, 0])->1},{(Subscript[y, 0])->-1}}

1-2(Subscript[Q, 2])[(Subscript[x, 0])]/.{{(Subscript[y, 0])->1},{(Subscript[y, 0])->-1}}

(Subscript[x, 1])/.{(Subscript[y, 1])->1}

In either case we drop immediately into the fixed-point at  0 , if not there already.

Exercise 2: Logistic Map, λ = 4 (chaotic)

(Subscript[Q, 4])[x]

Plot a graph, showing the most salient features:

Plot[{(Subscript[Q, 4])[x],x},{x,0,1},Frame->True,AspectRatio->1,PlotRange->All,PlotStyle->{{Hue[0],Thickness[0.008`]},{Hue[(2/3)]}},GridLines->{{0,(1/2),(3/4),1},{0,(3/4),1}},Background->White]

a:  fixed-points, stability

Locate the fixed-points:

(Subscript[fixedpoints, 4])=Solve[(Subscript[Q, 4])[x]==x]

Check their stability:

(D[(Subscript[Q, 4])[x],x]//Expand)/.(Subscript[fixedpoints, 4])

Both fixed-points are unstable.

b:  rewrite as angles

(Subscript[x, n_]):=((Sin[(Subscript[θ, n])])^2)

Notice that up to 4 different angles, modulo 2π, correspond to the same value of  (sin^2)θ .
How are they related?

c:  closed-form expression

Notice that:

Sin[2(Subscript[θ, n])]//TrigExpand

So we deduce that (Subscript[θ, n+1])=2(Subscript[θ, n]).  (More precisely, (Subscript[θ, n+1])=±2(Subscript[θ, n]) modulo π .)

This leads to a closed-form expression:  (Subscript[x, n])= (sin^2)((2^n)(Subscript[θ, 0])).

fixed-point for occurs for θ = 1/3 ℼ = 1/6 × 2π  as well as  0 × 2π ;
2-cycle  occurs for θ ∈ {1/5 × 2π, 2/5 × 2π } 
3-cycle occurs for  θ ∈ {1/9 × 2π, 2/9 × 2π, 4/9 × 2π } and  θ ∈ {1/7 × 2π, 2/7 × 2π, 4/7 × 2π }

d:  fixed-points, as angles

Angles corresponding to the fixed-points are easily identified:

fixedangles=Solve/@{2(Subscript[θ, n])==(Subscript[θ, n]),2(Subscript[θ, n])==π-(Subscript[θ, n]),2(Subscript[θ, n])==π+(Subscript[θ, n]),2(Subscript[θ, n])==2π-(Subscript[θ, n])}

(Subscript[x, n])/.%

D[(Subscript[Q, 4])[x],x]//Expand

(%/.{x->(Subscript[x, n])})//TrigReduce

Clearly these are unstable.

%/.fixedangles

e:  2-cycles

(Subscript[cycles, 2])=Solve/@{4(Subscript[θ, n])==(Subscript[θ, n]),4(Subscript[θ, n])==π-(Subscript[θ, n]),4(Subscript[θ, n])==π+(Subscript[θ, n]),4(Subscript[θ, n])==2π-(Subscript[θ, n]),4(Subscript[θ, n])==2π+(Subscript[θ, n]),4(Subscript[θ, n])==3π-(Subscript[θ, n]),4(Subscript[θ, n])==3π+(Subscript[θ, n])}

(Subscript[x, n])/.%

Eliminating the known fixed-points, there is a 2-cycle:  (Subscript[x, ±])= (1/8) (5±(Sqrt[5])).

(((8x-5-(Sqrt[5]))(8x-5+(Sqrt[5])))/4)//Expand//Simplify

Check its stability:

D[Nest[(Subscript[Q, 4]),x,2],x]//Expand

(%/.{x->(Subscript[x, n])})//TrigReduce

%/.{{(Subscript[θ, n])->(π/5)},{(Subscript[θ, n])->((2π)/5)}}

%//Expand

Alternatively:

(D[(Subscript[Q, 4])[x],x]//Expand)/.{x->(Subscript[x, n])}

%/.{{(Subscript[θ, n])->(π/5)},{(Subscript[θ, n])->((2π)/5)}}

Multiply these, to get the slope of the doubly-iterated map:

(Times@@%)//Expand

Clearly this 2-cycle is unstable.

f :  3-cycles

(Subscript[cycles, 3])=Solve/@{8(Subscript[θ, n])==(Subscript[θ, n]),8(Subscript[θ, n])==π-(Subscript[θ, n]),8(Subscript[θ, n])==π+(Subscript[θ, n]),8(Subscript[θ, n])==2π-(Subscript[θ, n]),8(Subscript[θ, n])==2π+(Subscript[θ, n]),8(Subscript[θ, n])==3π-(Subscript[θ, n]),8(Subscript[θ, n])==3π+(Subscript[θ, n]),8(Subscript[θ, n])==4π-(Subscript[θ, n]),8(Subscript[θ, n])==4π+(Subscript[θ, n]),8(Subscript[θ, n])==5π-(Subscript[θ, n])}

(Subscript[x, n])/.%

Eliminating the known fixed-points, there are two 3-cycles:

{{((Sin[(π/9)])^2),((Sin[((2 π)/9)])^2),((Sin[((4 π)/9)])^2)},{((Sin[(π/7)])^2),((Sin[((2 π)/7)])^2),((Sin[((3 π)/7)])^2)}}//N//TableForm

(D[(Subscript[Q, 4])[x],x]//Expand)/.{x->(Subscript[x, n])}

%/.{{{(Subscript[θ, n])->(π/7)},{(Subscript[θ, n])->((2π)/7)},{(Subscript[θ, n])->((4π)/7)}},{{(Subscript[θ, n])->(π/9)},{(Subscript[θ, n])->((2π)/9)},{(Subscript[θ, n])->((4π)/9)}}}

%//N

Multiply these, to get the slope of the triply-iterated map:

Times@@#&/@%

Clearly both 3-cycles are unstable.

See this in a Manipulate widget.

Evaluate[{Cos[#],Sin[#]}&/@( θ((2^(Range[1,6]))))]

Set the number of iterations; move the slider to change the initial angle.
How do you know when you have a cycle?

Manipulate[ListPlot[Evaluate[{Cos[#],Sin[#]}&/@ (θ((2^(Range[0,num])))2π)],PlotStyle->{{Hue[0],PointSize[0.02]}},PlotRange->1.1{{-1,1},{-1,1}},AspectRatio->1,Epilog->{Circle[{0,0},1]}],{{num,3,"number of iterates"},Range[0,20]},{{θ,1/9,"angle, as a fraction of π"},0,(1/4)}]

This next widget shows all 4 equivalent angles.
Now how do you tell when you have a cycle?

Manipulate[ListPlot[Evaluate[{{Cos[#],Sin[#]}&/@ (θ((2^(Range[0,num])))2π),{Cos[-#],Sin[-#]}&/@ (θ((2^(Range[0,num])))2π),{-Cos[#],Sin[#]}&/@ (θ((2^(Range[0,num])))2π),{-Cos[#],-Sin[#]}&/@ (θ((2^(Range[0,num])))2π)}],PlotStyle->{{Hue[0],PointSize[0.02]},{Hue[2/3],PointSize[0.02]},{Hue[1/12],PointSize[0.02]},{Hue[3/4],PointSize[0.02]}},PlotRange->1.1{{-1,1},{-1,1}},AspectRatio->1,Epilog->{Circle[{0,0},1]}]
,{{num,3,"number of iterates"},Range[0,20]},{{θ,1/9,"angle, as a fraction of π"},0,(1/4)}]

g:  k-cycles

There are always solutions to:

{{(2^k)(Subscript[θ, n])==j π -(Subscript[θ, n])},{(2^k)(Subscript[θ, n])==j π +(Subscript[θ, n])}}

(Subscript[cycles, k_]):=Map[Solve[#,(Subscript[θ, n])]&,Table[{{(2^k)(Subscript[θ, n])==j π -(Subscript[θ, n])},{(2^k)(Subscript[θ, n])==j π +(Subscript[θ, n])}},{j,1,(2^(k-1))}],{2}]

For example:

(Subscript[cycles, 4])

Indeed all rational multiples of π , with odd denominators, will arise this way, as angles corresponding to periodic points.

Exercise 3: Exploring period-2 and period-3

Set the working directory as needed, to find the  GraphicalAnalysis.m  package — adjust the cell below to suit your local environment

SetDirectory["~/Courses/MATH332/2015/"]

Needs["GraphicalAnalysis`"]

define the function to be iterated, showing the 1st 3 iterates

g[x_] := 3.84x(1-x) ;Plot[{x,g[x],g[g[g[x]]]},{x,0,1},Background->White]

Exercise 4: Obtain an estimate of the Feigenbaum number

Here is the working used by one of the students in 2007. He used intersections of Q-curves.
Note how approximations to nearly 13-digit accuracy were obtained.
Execute the cells to regenerate any missing images. (Make sure the package has loaded — see above!)

f[x_,a_] := a x (1-x)

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.0`,4`},Iterations->4,Transients->3,Slices->500],PlotRange->{0.0,0.99`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.569944`,3.569948`},Iterations->2000,Transients->0,Slices->1000],PlotRange->{0.83097,0.83107`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.6785735104280`,3.6785735104285`},Iterations->5,Transients->4,Slices->1000],PlotRange->{0.728155493652`,0.728155493655`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.5748049387590`,3.5748049387595`},Iterations->15,Transients->11,Slices->500],PlotRange->{0.816204957858,0.81620495786`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.5709859403412`,3.570985940342`},Iterations->31,Transients->23,Slices->500],PlotRange->{0.8377011469285,0.837701146929`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.5701684724960`,3.5701684724965`},Iterations->63,Transients->47,Slices->500],PlotRange->{0.8291323435781,0.8291323435782`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.569993388559`,3.5699933885596`},Iterations->127,Transients->95,Slices->500],PlotRange->{0.83255376369,0.83255376371`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.5699558913250`,3.5699558913255`},Iterations->255,Transients->191,Slices->500],PlotRange->{0.83118458453,0.83118458455`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.5699478605640`,3.5699478605648`},Iterations->447,Transients->319,Slices->500],PlotRange->{0.8309431113,0.830943111310`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.569945772245`,3.56994577252`},Iterations->2000,Transients->1000,Slices->500],PlotRange->{0.831059835,0.83105984`},Background->White]

This one shows a bifurcation and the indistinctness of where it actually occurs — it takes lots of bytes:

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.5699444`,3.5699446`},Iterations->2000,Transients->1000,Slices->500],PlotRange->{0.831054,0.831056`},Background->White]

Show[Bifurcation[f[x,a],{x,0.5`},{a,3.0`,4`},Iterations->100,Transients->7,Slices->500],PlotRange->{0.0,0.99`},Background->White]

values for μ

(Subscript[μ, 1])	3.6785735104283      using curves  4,5

(Subscript[μ, 2])        3.5925721841070         using curves  7,9

(Subscript[μ, 3])	3.5748049387592      using curves 11,15

(Subscript[μ, 4])	3.5709859403415      using curves 23,31

(Subscript[μ, 5])	3.5701684724963      using curves 47,63

(Subscript[μ, 6])	3.5699933885592      using curves 95,127

(Subscript[μ, 7])	3.5699558913252      using curves 191,255

(Subscript[μ, 8])	3.5699478605646      using curves 319,447

(Subscript[μ, 9])	3.5699461406220      using curves >600,<1000

(Subscript[μ, 10])	3.5699457722451     using curves >1000,<2000

Feigenbaum estimates

?SetPrecision

using  (Subscript[μ, 1]) to  (Subscript[μ, 3])

SetPrecision[(3.6785735104283-3.5925721841070)/(3.5925721841070 -3.5748049387592),13]

using  (Subscript[μ, 2]) to  (Subscript[μ, 4])

SetPrecision[(3.5925721841070-3.5748049387592)/(3.5748049387592 -3.5709859403415),13]

using  (Subscript[μ, 3]) to  (Subscript[μ, 5])

SetPrecision[(3.5748049387592-3.5709859403415)/(3.5709859403415 -3.5701684724963),13]

using  (Subscript[μ, 4]) to  (Subscript[μ, 6])

SetPrecision[(3.5709859403415-3.5701684724963)/(3.5701684724963 -3.5699933885592),13]

using  (Subscript[μ, 5]) to  (Subscript[μ, 7])

SetPrecision[(3.5701684724963-3.5699933885592)/(3.5699933885592 -3.5699558913252),13]

using  (Subscript[μ, 6]) to  (Subscript[μ, 8])

SetPrecision[(3.5699933885592-3.5699558913252)/(3.5699558913252 -3.5699478605646),13]

using  (Subscript[μ, 7]) to  (Subscript[μ, 9])

SetPrecision[(3.5699558913252-3.5699478605646)/(3.5699478605646-3.569946140622),13]

using  (Subscript[μ, 8]) to  (Subscript[μ, 10])

SetPrecision[(3.5699478605646-3.569946140622)/(3.569946140622 -3.5699457722451),17]

Even though the calculations have been carried out to 13-digit accuracy, the value is consistent only to 3–4 decimal places, as lying between say  4.6689  and 4.6693 .

The actual value, computed by other means, is roughly  4.669201609...  (Wikipedia) which is indeed within the estimated bounds.

Example error - calculation

SetPrecision[{((1.71994269-0.00000002)/(0.3683769+0.00000002)),(1.71994269/0.3683769),((1.71994269+0.00000002)/(0.3683769-0.00000002))},13]//ColumnForm

estimates for the end of period-doubling (Subscript[μ, ∞]) based on the second and third last estimates for the Feigenbaum constant using the last two  identified "leaves" for μ using a GP.

using 2nd-last estimate, sum a geometric series:

SetPrecision[(3.569946140622 -3.5699457722451)/(1-1/4.669202681851495),13]

subtract it from the right-hand end (larger) value:

SetPrecision[3.569946140622-% ,13]

using 3rd-last estimate, sum a geometric series:

SetPrecision[(3.569946140622 -3.5699457722451)/(1-1/4.6692008226080395),13]

SetPrecision[3.569946140622-%,13]

and now the same but back one step

SetPrecision[(3.5699478605646-3.569946140622)/(1-1/4.669202681851495),13]

SetPrecision[3.5699478605646-%,13]

SetPrecision[(3.5699478605646-3.569946140622)/(1-1/4.6692008226080395),13]

SetPrecision[3.5699478605646 -%,13]

These are consistent with values in the range  3.569945671848 – 3.569945671871 .

The actual value is :  3.5699456718709...  (OEIS) which indeed is just within this range; thus estimates made this way are accurate to 11+ decimal digits.