The Quadratic Map

Code for Chapter Opening Art

.  Iterating Functions

NestList[Cos,1.5,20]

NestList[Cos,Nest[Cos,1.5,100],10]

Cos[0.739085]

NestList[Cos,1,5]

Table[Nest[Cos,c,100],{c,RandomReal[{-100,100},{10}]}]

Plot[Evaluate[NestList[Cos,x,5]],{x,0,π},PlotRange->{-0.2,1.5},Frame->True,
Axes->False,PlotStyle->Thick,GridLines->{{},{{0.739,Dashed}}},FrameTicks->{Range[0,π,(π/2)],Range[-0.5,2,0.5],None,{0.739}}]

Options[CobwebPlot] = {WorkingPrecision->MachinePrecision};

CobwebPlot[f_,{x_,a_,b_},x0_,n_Integer, opts___]:=CobwebPlot[f,{x,a,b},x0,{0,n}, opts];

CobwebPlot[f_,{x_,a_,b_},x0_,{n0_,n_}, opts___]:=Module[{fcn,data,start,wp},
wp=WorkingPrecision/.{opts}/.Options[CobwebPlot];
fcn=If[wp===MachinePrecision,Compile[x,f],Function[{x}, N[f, wp]]];
start=Nest[fcn,x0,n0];
data=NestList[fcn,start,n-1];
Show[Plot[f,{x,a,b},Evaluate[Sequence @@ FilterRules[{opts},Options[Plot]]],PlotStyle->{Red,Thickness[0.004]}],
Graphics[{{Thickness[0.004],Black,
Line[Partition[Riffle[data,data],2,1]]},
{PointSize[Large],Point[{start,start}]},
{Dashing[{0.013,0.023}],Thickness[0.004],Line[{{a,a},{b,b}}]}}],
Sequence @@ FilterRules[{opts},Options[Graphics]],PlotRange->All, Frame->True, Axes->False, FrameTicks->{Automatic,Automatic,False,False}]]

(Subscript[f, r_])[x_]:=r x (1-x)

CobwebPlot[(Subscript[f, 2])[x],{x,0,1},0.1,10]

Grid[Partition[{
CobwebPlot[(Subscript[f, 2])[x],{x,0,1},0.1,10],
       CobwebPlot[(Subscript[f, 3.45])[x],{x,0,0.9},0.1,100],
       CobwebPlot[(Subscript[f, 3.45])[x],{x,0.3,0.9},0.1,{300,50}],
       CobwebPlot[(Subscript[f, 3.839])[x],{x,0,1},0.1,{150,50}],
       CobwebPlot[(Subscript[f, 4])[x],{x,0,1},0.1,{200,75}],
       CobwebPlot[(Subscript[f, 4])[x],{x,0,1},0.1,600]},2]]

GraphicsColumn[Table[GraphicsRow[{
CobwebPlot[(Subscript[f, r])[x],{x,0,1},0.1,100, PlotRange->{{0.2, 0.85},{0.3,0.85}}, FrameTicks->{{.3,.8},{.3,.8},None,None}],CobwebPlot[Nest[(Subscript[f, r]),x,2],{x,0,1},0.1,100, PlotRange->{{0.2, 0.85},{0.5,0.8}}, FrameTicks->{{.3,.8},{.3,.8},None,None}]},  AspectRatio-> 1/2,Frame->True,PlotRegion->{{0,1},{0,1}}, ImagePadding->Automatic, ImageSize->300,PlotLabel->StringForm["r = ``",PaddedForm[r,{3,2}]]],{r,2.9,3.1,0.05}]]

Manipulate[Grid[{
{CobwebPlot[Nest[(Subscript[f, r]),x,2],{x,0,1},0.1,100, PlotRange->{{0.2, 1},{0.5,0.8}}],
CobwebPlot[(Subscript[f, r])[x],{x,0,1},0.1,100, PlotRange->{{0.2, 0.8},{0.3,0.75}}]}}],{r,2.9,3.3,0.01}]

.  The Four Flavors of Real Numbers

(20/49)+(31/17)

1.3+1.12345678955555 (10^(-7))//InputForm

Precision[%]

1.+(1/(10^17))-1

N[1+(1/(10^17))-1]

(a=N[π,30])//InputForm

Precision[(a^100000)]

Nest[(Subscript[f, 4]),N[1/10,30],30]//InputForm

a=N[(π/2)+(1/(10^7)),30];
Precision[Sin[a]]

Sin[Interval[{-π,π}]]

g[x_] := Sin[x]+Cos[x];
g[Interval[{-π,π}]]

Nest[(Subscript[f, 4]),Interval[N[1/10,30]],30]

N[Sin[(10^50)]]

N[Sin[(10^50)],20]

Sin[((Interval[N[10,75]])^50)]

$MaxExtraPrecision=500;
N[Sin[(10^100)],20]

Nest[(Subscript[f, 4]), 0.1, 100]

Nest[(Subscript[f, 4]),N[1/10,70],100]

Log[10.,Numerator[Nest[(Subscript[f, 4]),1/10, 21]]]

Nest[(Subscript[f, 4]), Interval[N[1/10, 100]], 100]

70-Round[Precision/@NestList[(Subscript[f, 4]),N[(1/10),70],100]]

δ=(10^(-150));x0=(1/10);
Table[Δ=Abs[Nest[(Subscript[f, 4]),N[x0-(δ/2),200],n]-Nest[(Subscript[f, 4]),N[x0+(δ/2),200],n]];Round[Log[10,(Δ/δ)]],{n,0,100}]

r=31/8;
Row[{CobwebPlot[(Subscript[f, r])[x],{x,0,1},1/10,{100,300}, PlotStyle->{Red,Thick}, ImageSize->180],Spacer[20],
        CobwebPlot[(Subscript[f, r])[x],{x,0,1},1/10,{100,300},WorkingPrecision->1000,PlotStyle->{Red,Thick}, ImageSize->180]}]

.  Attracting and Repelling Cycles

Clear[r];x/.Solve[(Subscript[f, r])[x]==x,x]

Table[Sort[NestList[(Subscript[f, 3.45]),Nest[(Subscript[f, 3.45]),Random[],500],3]],{10}]

(x/.FindRoot[Nest[(Subscript[f, 3.45]),x,4]==x,{x,#},AccuracyGoal->12]&)/@%[[1]]

CobwebPlot[(Subscript[f, 3.45])[x],{x,0,1},(49/69)+0.01,100]

(Subscript[f, 345/100])[(49/69)]

Solve[(Subscript[f, (345/100)])[x]==(49/69)]

n=80;
orbit=Transpose[{Range[0,n],NestList[(Subscript[f, 4]),N[(1/10),n],n]}];
orbit1=Transpose[{Range[0,n],NestList[(Subscript[f, 4]),N[(1/10)+(1/(10^17)),n],n]}];
ListLinePlot[{orbit,orbit1},PlotStyle->{{Thickness[0.013],Orange},{Thickness[0.004],Black}},AspectRatio->0.25]

TrigIteration[n_, x0_]:=(1/2)(1- Cos[(2^n) ArcCos[1-2x0]])

N[{TrigIteration[5,(1/10)],Nest[(Subscript[f, 4]),(1/10),5]}]

(x/.RSolve[{x[n+1]==4 x[n] (1-x[n]),x[0]==x0},x,n][[1]])[n]//Quiet

TrigExpand[TrigIteration[n+1,x]==4 TrigIteration[n,x] (1-TrigIteration[n,x])]

fval=(Subscript[f, 4])[((Sin[((π v)/2)])^2)];
Simplify[(2/π) ArcSin[(Sqrt[fval])],0≤v≤(1/2)]

PowerExpand[Simplify[(2/π) ArcSin[(Sqrt[fval])],(1/2)<v≤1],Assumptions->(1/2)<v≤1]

stringOp[bits_]:=If[bits[[1]]==0,Rest[bits],1-Rest[bits]]
Nest[stringOp,Flatten[Table[{0,1,0},{10}]],3]

v=FromDigits[{{{0,1,0}},0},2]

ToX[v_]:=((Sin[((π v)/2)])^2);
FromX[x_]:=(2/π)ArcSin[(Sqrt[x])];
x0=ToX[v]

CobwebPlot[(Subscript[f, 4])[x],{x,0,1},x0,10]

v=31/64;
x0=ToX[v]

N[NestList[(Subscript[f, 4]),N[x0,20],7],6]

CobwebPlot[(Subscript[f, 4])[x],{x,0,1},x0+(10^(-10)),80,WorkingPrecision->200]

.  Measuring Instability: The Lyapunov Exponent

Lyapunovλ[r_,x0_,n_,δ_]:=With[{prec=Max[2 n,-Log[10,δ] 2+20]},
Δ=Abs[Nest[(Subscript[f, r]),N[x0-(δ/2),prec],n]-Nest[(Subscript[f, r]),N[x0+(δ/2),prec],n]];N[(1/n)Log[(Δ/δ)],4]]/; Union[Precision /@ {r, x0, n,δ} ]== {∞}

(Subscript[f, 4])[.1]

ListLinePlot[data=Table[{d,Lyapunovλ[4, 1/10,100,(10^(-d))]}, {d,10,100,10}], Frame->True, Axes->False,PlotStyle->{Thick,Black}, PlotRange->{-0.1,1}, Epilog ->{PointSize[0.02], Point[data]}, FrameTicks -> {{#, HoldForm[(10^(#))]} & /@ Range[20,200,20],Automatic, None, None}]

λByDer[r_,x0_,n_]:=(1/n)Log[Abs[(∂_x)Nest[(Subscript[f, r]),x,n]]]/.x->x0
{λByDer[4,0.1,10],Lyapunovλ[4,(1/10),10,(1/(10^10))]}

Manipulate[r1=Rationalize[r];ListLinePlot[Table[{d,Lyapunovλ[r1, 1/10,100,(10^(-d))]}, {d,10,100,10}], Frame->True, Axes->False,PlotStyle->{Thick,Black}, FrameTicks -> {{#, HoldForm[(10^(#))]}& /@ Range[20,200,20],Automatic, None, None},PlotRange->{{5,105},{-3,1}}, PlotLabel-> StringForm["r = ``",N[r,3]]], {{r,4},2,4, 1./100}]

ListLinePlot[Table[Table[{d,Lyapunovλ[r,(1/10),100,(10^(-d))]},{d,10,40,10}],{r,2,4,(1/100)}],Frame->True,Axes->False,PlotStyle->{Thick,Black},FrameTicks->{({#1,HoldForm[(10^(-#1))]}&)/@Range[10,40,10],Automatic,None,{0,{Log[2],"log 2"}}},PlotRange->{{8,42},{-3,1}}]

Monitor[ListLinePlot[Table[{r,Lyapunovλ[r,(1/10),300,(10^(-100))]},{r,1,4,(3/20000)100}], Frame->True, PlotStyle->{Thickness[0.0001],Black},
PlotRange->All, Filling->{1->{0, Red}}], ProgressIndicator[r,{1,4}]]

.  Bifurcations

cf=Compile[{{r,_Real,1},{x,_Real,1}},r+x];
cf[{1.,2.,3.},{2.,3.,4.}]

cf=Compile[{{r,_Real,1},{x,_Real,1}},r x (1-x)];
rVals=Range[2.,4,((4-2)/(10-1))];
Transpose[NestList[cf[rVals,#1]&,Nest[cf[rVals,#1]&,(0.1&)/@rVals,10],4]]

NestList[(Subscript[f, 4]),Nest[(Subscript[f, 4]),0.1,10],4]

Column[Transpose[{rVals, Transpose[NestList[cf[rVals,#]&, Nest[cf[rVals, #]&, 
          Array[ 0.1 &,10] , 10], 4]]}] ]

Options[BifurcationPlot]={PlotPoints->2000,PlotStyle->{Black,PointSize[0.001]}};
BifurcationPlot[f_,{r_,a_,b_},{x_,x0_},{iter0_,iterShow_},opts___]:=Module[{sty,n,makePts,cf,rVals,data},
{sty,n}={PlotStyle,PlotPoints}/.{opts}/.Options[BifurcationPlot];
makePts[{s_,v_}]:=({s,#}&)/@v;
cf=Compile[{{r,_Real,1},{x,_Real,1}},Evaluate[f]];
rVals=Range[N[a],b,((b-a)/(n-1))];
data=Transpose[{rVals,Transpose[NestList[cf[rVals,#]&,Nest[cf[rVals,#1]&,Array[x0&,n],iter0], iterShow]]}];
Graphics[Append[Flatten[{sty}],Point[Flatten[makePts/@data,1]]],Sequence@@FilterRules[{opts},Options[Graphics]], 
AspectRatio->1/3,Frame->True,FrameTicks->{Automatic,Range[0,1,.5], None, None},Axes->None]];

BifurcationPlot[(Subscript[f, r])[x],{r,2.8,4},{x,0.5},{500,100},FrameTicks->{{3,3.45,3.57,3.839,4},Automatic,{{69/20, "69/20"}},None},GridLines->{{3,3.45,3.57,3.839},None}]

Options[BifurcationPlotEfficient]={PlotPoints->2000,PlotStyle->{Black,PointSize[0.001]}};
BifurcationPlotEfficient[f_,{r_,a_,b_},{x_,x0_},{iter0_,iterShow_},opts___]:=Module[{ n,cf,rVals,data  },
{sty,n}={PlotStyle,PlotPoints}/.{opts}/.Options[BifurcationPlot];
makePts[{s_,v_}]:=({s,#}&)/@v;
cf=Compile[{r,x},Evaluate[f]];
com1=Compile[{rr}, ((1+rr+(Sqrt[-3-2 rr+(rr^2)]))/(2 rr))];
com2=Compile[{rr}, ((1+rr-(Sqrt[-3-2 rr+(rr^2)]))/(2 rr))];
iter[rr_]:=  Which[rr<3,{1-1/rr},
rr<3+45/100,{com1[rr], com2[rr]},True,NestList[cf[rr,#]&,Nest[cf[rr,#]&,x0,iter0],iterShow]];
rVals=Range[N[a],b,((b-a)/(n-1))];
data = {#,iter[#]}& /@ rVals;
ListPlot[Flatten[makePts/@data,1],Sequence@@FilterRules[{opts},Options[Graphics]], 
Frame->True,FrameTicks->{Automatic,Range[0,1,0.5], None, None},Axes->None,PerformanceGoal->"Speed",PlotStyle->sty]];
bifPlot=BifurcationPlotEfficient[(Subscript[f, r])[x],{r,2.9,4},{x,0.5},{500,32}];

Manipulate[Column[{CobwebPlot[(Subscript[f, r])[x],{x,0,1},0.1,{300,20},ImageSize->200],Show[bifPlot,GridLines->{{{r,Red}},{}},ImageSize->200]}],{{r,3.84},2.9,4,Appearance->"Labeled"}]

bifPts={3,3.45,3.54409,3.56441,3.56876,3.56969,3.56989,3.569943,3.5699451};
limitr=3.569945557391439;
logMin=Log[10,limitr-2.5];
n=2000; x0=0.1; iter0=1000;iterShow=64;
makePts[{s_,v_}]:=({s,#1}&)/@v;
cf=Compile[{{r,_Real,1},{x,_Real,1}},Evaluate[(Subscript[f, r])[x]]];
rValsLog=Range[logMin,-3,(((-3-logMin))/(n-1))];
rVals=limitr-(10^rValsLog);

data=Transpose[{-rValsLog,Transpose[NestList[cf[rVals,#]&,Nest[cf[rVals,#]&,Array[x0 &,n],iter0],iterShow]]}];

Graphics[{PointSize[0.0013],Point[Flatten[makePts/@data,1]]},Frame->True,AxesOrigin->{0,0.3},GridLines->{-Log[10,limitr-bifPts],None},FrameTicks->{({-#,NumberForm[limitr-(10^(#)),4]} &)/@Log[10,limitr-bifPts],{0.3,0.6,0.9},None,None},PlotRange->{0.3,0.9}]

eFunc[x_]:=(1/2)(E/(E-1))((E^(-(((2 x-1))^2)))-(1/E));
g[x_] := eFunc[Sin[π (x/2)]];
Plot[eFunc[x],{x,0,1}]
Plot[g[x],{x,0,1}]

BifurcationPlot[r g[x],{r,0.95,2},{x,1/10},{1000,128},PlotPoints->1000,PlotRangePadding->{{0,0.03},{0.1,0.1}},FrameTicks->{Automatic,{0,1,2},None,None}]