Entertaining Science done at
LASSP.
(* The logistic map is an inverted parabola, depending on a parameter mu. *)
f[x_] := 4 mu x (1-x)
mu := 1.
Plot[f[x],{x,0,1}]
x0 = 0.1
y0 = 0.1000000001
(* NestList makes a list of values {x, f[x], f[f[x]], ... }: *)
NestList[f,x0,35] NestList[f,y0,35] Diff = %-%%(* Notice that the two curves separate exponentially with time. *)
ListPlot[Diff, PlotJoined->True] ListPlot[Log[Abs[Diff]], PlotJoined->True](* The Lyapunov exponent is the initial slope on the log plot. We can find the slope using Fit. We remove the final points using Drop, since they are far separated, bouncing around independently and not still separating exponentially. *)
Fit[Log[Abs[Take[Diff,30]]],{1,t},t]
Plot[%,{t,0,35}]
ListPlot[Log[Abs[Diff]], PlotJoined->True]
Show[%,%%]
mu := 0.2
Plot[{x, f[x]},{x,0,1}]
flist = NestList[f,0.7,30]
(* We can define a fancy routine to plot the iterates. *)
PlotIterate[x0_,n_,options___] := {flist = NestList[f,x0,n];
itplot = ListPlot[
Flatten[Drop[Transpose[
{Transpose[{flist,flist}], Transpose[{flist,RotateLeft[flist]}]}
],-1],1],
PlotJoined->True, DisplayFunction->Identity];
rplot = Plot[f[x],{x,0,1}, DisplayFunction->Identity];
Show[itplot,rplot,DisplayFunction->$DisplayFunction,options]}
PlotIterate[0.7,30]
(* For mu > 1/4, the origin becomes unstable. (Why?) A new stable fixed
point forms. *)
mu := 0.3
Plot[{x,f[x]},{x,0,0.3}]
NestList[f,0.1,30]
PlotIterate[0.7,30,PlotRange->{{0,0.3},{0,0.3}}]
(* Things start happening for larger mu. *)
(* Above mu = 1/2, we get alternation on the way to convergence. *)
mu := 0.7 NestList[f,0.1,30] PlotIterate[0.1,30,PlotRange->All](* Above mu = 3/4, the fixed point becomes unstable, and a period--two cycle starts. *)
mu := 0.8 NestList[f,0.1,30] PlotIterate[0.1,20](* We can see the attractor more clearly if we first iterate the map a few times to get onto the attractor. Nest iterates the map without forming a list. *)
Nest[f,0.1,20] PlotIterate[%,20](* Period Doubling Continues ... *)
mu := 0.89 PlotIterate[Nest[f,0.1,40],20](* And Chaos Begins *)
mu := 0.9 PlotIterate[Nest[f,0.1,100],20] mu := 1.0 PlotIterate[0.1,100]
Clear[mu] f[x_] := 4 mu x (1-x) Solve[f[x]==x,x] // Simplify f'[1-1/(4 mu)] // Simplify Solve[%==-1,mu](* Show that there is a stable period--two cycle for 3/4 < mu < (1 + sqrt[6])/4. *)
TwoPeriods = Solve[f[f[x]]==x,x] // Simplify(* The first solution is zero: the unstable fixed point is also, naturally, an unstable period 2 cycle. The second is the unstable fixed point xstar. The third and fourth are the period--two cycle. The period two cycle itself becomes unstable when it's derivative goes through minus 1. *)
D[f[f[x]],x] /. TwoPeriods // Simplify Solve[%[[3]] == -1,mu] // Simplify(* Show that the amplitude (x2-x1) of this cycle is proportional to sqrt[mu - 3/4] when mu-3/4 << 1. *)
x /. TwoPeriods
%[[3]]-%[[4]] // Simplify
Series[%,{mu,3/4,2}]
(* Show that the period--two cycle is superstable when mu = (1 + sqrt[5])/4:
a "superstable" cycle contains the point x=1/2, so that the Lyapunov
exponent is zero. *)
mu = (1 + Sqrt[5])/4 f[1/2] f[%] // Simplify D[f[f[x]],x] /. x -> 1/2
Clear[mu]
f[x_,mu_] := 4 mu x (1-x)
Attractor[mu_,ntransient_,ncycle_] := Transpose[{Table[mu,{ncycle+1}],
NestList[f[#,mu]&,Nest[f[#,mu]&,0.1,ntransient],ncycle] }]
Attractor[mu_] := Attractor[mu,100,32]
Attractor[0.7]
Attractor[0.8]
Attractor[0.9]
(* 74.5 Second *)
Attractors = Table[Attractor[mu],{mu,0,1,0.01}];
ListPlot[Flatten[Attractors,1]]
Show[%,PlotRange->{{0.8,1.0},{0,1}}]
(* 45.2 Second *)
Attractors = Table[Attractor[mu,100,64],{mu,0.87,0.92,0.001}];
ListPlot[Flatten[Attractors,1]]
f[1/2] == 1/2 mu[0] = 0.5 f[f[1/2]] == 1/2 Solve[%,mu](* There are three roots here. We can use FindRoot to find a numerical solution; if we start near the expected solution, perhaps we'll get it. *)
FindRoot[f[f[1/2]]==1/2,{mu,0.9}]
mu[1] = mu /. %
(* To get the later bifurcation points mu[n], we need a function which
iterates f many times. In Mathematica, the standard method (Nest)
actually expands the composition out algebraically before plugging
in: iterating f 2^12 times involves generating an expression of
order 2^(2^12) -1 long in mu. We use old--fashioned iteration to
speed things up. *)
(* Unfortunately, Mathematica does things algebraically unless we insist
otherwise. We insist that x and mu be real numbers before iterating! *)
f[x_Real,mu_Real,NTimes_] := Block[{y = x},
Do[y = 4 mu y (1-y) ,{NTimes}]; y]
mu[2] = mu /. FindRoot[f[1/2,mu,2^2]==1/2,
{mu,{0.899,0.9}}]
(* As n increases, there are more and more roots to f^(n)[1/2] = 1/2.
To find the right one, we have to be careful. What we can do is
use the fact that the roots are spaced in an approximate geometric
progression: mu[n] ~ mu[n-1] + (mu[n-1] - mu[n-2]) / delta: *)
delta[n_] := (mu[n-1] - mu[n-2]) / (mu[n] - mu[n-1])
delta[2]
mu[n_] := mu[n] = mu /.
FindRoot[f[1/2,mu,2^n]==1/2,
{mu,
{mu[n-1]+(mu[n-1]-mu[n-2])/delta[n-1],
mu[n-1]+(mu[n-1]-mu[n-2])/(delta[n-1]+.01)}
}, AccuracyGoal->14]
mu[3]
delta[3]
mu[4]
delta[4]
mu[5]
mu[6]
mu[7]
Table[mu[n],{n,0,7}]
Table[delta[n],{n,2,7}]
(* delta is converging nicely to the universal value. Try alpha: the
ratio of sizes of successive forks: *)
alpha[n_] := alpha[n] =
(f[0.5,mu[n-1],2^(n-2)]-1/2) / (f[0.5,mu[n],2^(n-1)]-1/2)
Table[alpha[n],{n,1,7}]
(* It, too, is converging to the universal value. It is convenient to
extrapolate the geometrical series for mu[n] to its limit: *)
muInfinity = mu[7] + (mu[7] - mu[6]) / (delta[7]-1)(* What does it mean that these constants are universal? It means that they are unchanged when we change f! Let's now try f2[x] = B sin[pi x] *)
f2[x_] := B Sin[Pi x]
B[0] = B /. FindRoot[f2[1/2] == 1/2,{B,0.9}]
B[1] = B /. FindRoot[f2[f2[1/2]]==1/2,{B,0.9}]
f2[x_Real,B_Real,NTimes_] := Block[{y = x},
Do[y = B Sin[N[Pi] y] ,{NTimes}]; y]
B[2] = B /. FindRoot[f2[1/2,B,2^2]==1/2,{B,{0.899,0.9}}]
delta2[n_] := (B[n-1] - B[n-2]) / (B[n] - B[n-1])
B[n_] := B[n] = B /.
FindRoot[f2[1/2,B,2^n]==1/2,
{B,
{B[n-1]+(B[n-1]-B[n-2])/delta2[n-1],
B[n-1]+(B[n-1]-B[n-2])/(delta2[n-1]+.01)}
}, AccuracyGoal->14]
Table[B[n],{n,0,6}]
BInfinity = B[6] + (B[6] - B[5]) / (delta2[6]-1)
Table[delta2[n],{n,2,6}]
alpha2[n_] := alpha2[n] =
(f2[0.5,B[n-1],2^(n-2)]-1/2) / (f2[0.5,B[n],2^(n-1)]-1/2)
Table[alpha2[n],{n,1,6}]
(* The limits delta and alpha are the same! *)
mu = 0.5
Plot[{x,f[x]},{x,0,1}]
Plot[{x,f[f[x]]},{x,0,1}]
(* At mu = 0.8, there is the (unstable) fixed point, and the period--two
cycle. *)
mu = 0.8
Plot[{x,f[x]},{x,0,1}]
Plot[{x,f[f[x]]},{x,0,1}]
(* f[f[x]] has fixed points where f[x] does, and in addition has fixed points
where f has period--two cycles. *)
(* Notice that f[f[x]] has a downward--pointing hump near x=1/2, quite similar
to the upward--pointing original map. The next period doubling occurs by
analogy to the first one: *)
mu = 0.875
Plot[{x,f[x]},{x,0,1}]
Plot[{x,f[f[x]]},{x,0,1}]
Plot[{x,f[f[f[f[x]]]]},{x,0,1}]
(* We can make this analogy precise, by defining a new map T[g] which is
a blown--up, rescaled version of g[g[x]]. We'll see that, for
mu = muInfinity, that T (as a map from function space into itself)
has a fixed point. *)
mu = muInfinity(* First, the new origin of T[g] will be the fixed point of g: *)
xstar[g_] := (x /. FindRoot[g[x]==x,{x,0.5,0,1}])
(* Next, we must expand both the x and g[x] axes by a scale factor, chosen
so that xstar will stretch to zero. The scale factor will end up
converging to alpha, so we name it appropriately: *)
Clear[alpha] alpha[g_] := - 1 / (2 (xstar[g] - 1/2))(* Test these on the original map: *)
xstar[f] alpha[f](* We now define T[g][x]: *)
T[g_][x_] := alpha[g] ( g[g[(x-1/2)/alpha[g] + 1/2]] - g[xstar[g]] )(* Try it once. *)
T[f][x]
% // Simplify
Plot[f[x],{x,0,1}]
Plot[Release[T[f][x]],{x,0,1}]
Show[%,%%]
(* Try it twice. *)
T[T[f]][x]
Plot[f[x],{x,0,1}]
Plot[Release[T[f][x]],{x,0,1}]
Plot[Release[T[T[f]][x]],{x,0,1}]
Show[%,%%,%%%]
(* Notice that T[T[f]] is almost identical to T[f]. One can also check
that T[T[f2]] is almost identical to T[T[f]]: *)
B = BInfinity
T[f2][x]
Plot[Release[f[x]],{x,0,1}]
Plot[Release[f2[x]],{x,0,1}]
Show[%,%%]
Plot[Release[T[f][x]],{x,0,1}]
Plot[Release[T[f2][x]],{x,0,1}]
Show[%,%%]
Plot[Release[T[f][x]-T[f2][x]],{x,0,1}]
Plot[Release[T[T[f]][x]-T[T[f2]][x]],{x,0,1}]
(* By using this rapid convergence in function space, one can prove both
that there will be an infinite geometrical series of period--doubling
bifurcations leading to the onset of chaos, and that this series
will share universal features (alpha, delta, and others) which are
independent of the original dynamics. *)
Entertaining Science done at LASSP.
Jim Sethna, sethna@lassp.cornell.edu
Statistical Mechanics: Entropy, Order Parameters, and Complexity,
now available at
Oxford University Press
(USA,
Europe).