Invariant Measures

(Sethna, “Entropy, Order Parameters, and Complexity”, ex. 4.3)
© 2018, James P. Sethna, all rights reserved.

Liouville’s theorem tells us that all available points in phase space are equally weighted when a Hamiltonian system is averaged over all times. What happens for systems that evolve according to laws that are not Hamiltonian? Usually, the system does not continue to explore all points in its state space; at long times it is confined to a subset of the original space known as the attractor.

We consider the behavior of the `logistic’ mapping from the unit interval (0,1) into itself:
f(x)=4μx(1−x).
(We also study this map in the exercises “Chaos, Lyapunov, and entropy increase”, “Fractal Dimensions”, and “Period doubling”). We talk of the trajectory of an initial point x0 as the sequence of points x0, f(x0), f(f(x0)), \dots, f[n](x0), ... Iteration can be thought of as a time step (one iteration of a Poincare return map of the exercise “Jupiter” or one step Δt in a time-step algorithm.

f[x_,μ_] := 4 ...

Attracting fixed point. For small μ, our mapping has an attracting fixed-point. A fixed-point of a mapping is a value x*=f(x*); a fixed-point is stable if small perturbations shrink after iterating:
|f(x*+ϵ)−x*|≈|f ‘(x*)|ϵ<ϵ,
which happens if the derivative |f ’(x∗)|<1. (For many-dimensional mappings, a sufficient criterion for stability is that all the eigenvalues of the Jacobian have magnitude smaller than one. A continuous time evolution dy/dt=F(y) will be stable if dF/dy is smaller than zero, or (for multidimensional systems) if the Jacobian DF has eigenvalues whose real parts are all less than zero.)

(a) Iteration. Set μ=0.2; iterate f for some initial points 0<x0<1 of your choosing, and convince yourself that they are all attracted to zero. Plot f and the diagonal y=x on the same plot. Are there any fixed-points other than x=0? Repeat for μ=0.4, and 0.6. What happens? than zero.)

(* NestList iterates f(#,μ)& and returns a list of the values 
   f(#,0.2)& is the one-variable 'pure function' corresponding to f(x,μ) where μ=0.2 *)
NestList[f[#,...&,...]

Plot[{f[x,0.2], ..., x},{x,0,1}, AspectRatio->1, PlotLegends->LineLegend["Expressions"]]

(b) Analytics. Find the non-zero fixed-point x*(μ) of the map f(x), and show that it exists and is stable for 1/4<μ<3/4. If you are ambitious or have a computer algebra program, show that there is a stable, period-two cycle for 3/4 < μ < (1+(Sqrt[6]))/4.

xStarSols = Solve[...,x]
xStar = x/.xStarSols[...] // Expand   
(* Stability if |f'[x]|<1 *)
dfAtxStar = D[...] /. x->xStar // Expand
Solve[dfAtxStar==-1, μ]

x2StarSols = Solve[f[f[...]...]==x,x] // Expand
x2Star =...
df2Atx2Star = D[...,x] /. x->x2Star // Expand
Solve[df2Atx2Star==-1, μ]

An attracting fixed-point is the antithesis of Liouville’s theorem; all initial conditions are transient except one, and all systems lead eventually to the same, time-independent state. (On the other hand, this is precisely the behavior we expect in statistical mechanics on the macroscopic scale; the system settles down into a time-independent equilibrium state! All microstates are equivalent, but the vast majority of accessible microstates have the same macroscopic behavior in most large systems.) We could define a rather trivial `equilibrium ensemble’ for this system, which consists of the single point x*; any property O(x) will have the long-time average ⟨O⟩=O(x*).

For larger values of μ, more complicated things happen. At μ=1, the dynamics can be shown to fill the entire interval; the dynamics is ergodic, and the attractor fills the entire set of available states. However, unlike the case of Hamiltonian systems, not all states are weighted equally (i.e., Liouville’s theorem does not hold).

We can find time averages for functions of x in two ways: by averaging over time (many iterates of the map) or by weighting an integral over x by the invariant density ρ(x). The invariant density ρ(x)dx is the probability that a point on a long trajectory will lie between x and x+dx. To find it numerically, we iterate a typical point (not an unstable fixed-point or unstable periodic orbit!) x0 a thousand or so times (Ntransient) to find a point xa on the attractor, and then collect a long trajectory of perhaps a million points (Ncycles). A histogram of this trajectory gives ρ(x). Averaging over this density is manifestly the same as a time average over the trajectory of a million points. We call ρ(x) invariant because it is left the same under the mapping f; iterating our million-point approximation for ρ once under f only removes the first point xa and adds one extra point to the end.

(b) Invariant density. Set μ=1; iterate f many times, and form a histogram of values giving the density ρ(x) of points along the trajectory. You should find that points x near the boundaries are approached more often than points near the center.

Analytics. Using the fact that the long-time average ρ(x) must be independent of time, verify for μ=1 that the density of points is
                ρ(x)=(1/(π (Sqrt[(1-x) x]))).
(You need not derive the factor 1/π, which normalizes the probability density to one.) Plot this theoretical curve with your numerical histogram. (Hint: The points in a range dx around a point x map under f to a range dy=f ’(x)dx around the image y=f(x). Each iteration maps two points xa and xb = 1−xa to y, and thus maps all the density ρ(xa)|dxa| and ρ(xb)|dxb| into dy. Hence the probability ρ(y)dy must equal ρ(xa)|dxa|+ρ(xb)|dxb|, so
                  ρ(f(xa))=ρ(xa)/|f ’(xa)|+ρ(xb)/|f ‘(xb)|.
Substitute the second-to-latest equation for ρ(x) above into this equation. You will need to factor a polynomial.)

xa = Nest[...];
xs = NestList[...];
hist = Histogram[xs, 2000, "PDF"]
plot = Plot[...,{x,0,1}];
Show[hist,plot]

Mathematicians call this probability density ρ(x)dx the invariant measure on the attractor. (There are actually many possible invariant measures on some attractors; this one is the SRB measure (John Guckenheimer, private communication).) To get the long-term average of any function O(x), one can use
                  ⟨O⟩=∫O(x)ρ(x)dx.
To a mathematician, a measure is a way of weighting different regions when calculating integrals---precisely our ρ(x)dx. Notice that, for the case of an attracting fixed-point, we would have ρ(x)=δ(x-x*). (The case of a fixed-point then becomes mathematically a measure with a point mass at x*.)

Cusps in the invariant density. At values of μ slightly smaller than one, our mapping has a rather complex invariant density.

(c) Find the invariant density (as described above) for μ=0.9. Make your trajectory length Ncycles big enough and the bin size small enough to see the interesting structures. Notice that the attractor no longer fills the whole range (0,1); locate roughly where the edges are. Notice also the cusps in ρ(x) at the edges of the attractor, and also at places inside the attractor (called ‘boundaries’). Locate some of the more prominent cusps.

xa = Nest[...];
xs = NestList[...];
hist09 = Histogram[xs, 2000, "PDF"]

Analytics of cusps. Notice that f ’(1/2)=0, so by the equation for ρ(f(xa)) above we know that ρ(f(x))≥ρ(x)/|f ‘(x)| must have a singularity near x=1/2; all the points near x=1/2 are squeezed together and folded to one side by f. Further iterates of this singularity produce more cusps; the crease after one fold stays a crease after being further stretched and kneaded.

(d) Set μ=0.9. Calculate f(1/2), f(f(1/2)), ... and compare these iterates to the locations of the edges and cusps from part~(c). (You may wish to include them both on the same plot.)

cusps = ListPlot[Table[{Nest[...,...,n],n},{n,1,8}], Joined->True]
Show[hist09, cusps]

Bifurcation diagram. The evolution of the attractor and its invariant density as μ varies are plotted in the bifurcation diagram. One of the striking features in this plot are the sharp boundaries formed by the cusps.

(e) Bifurcation diagram. Plot the attractor as a function of μ, for 0.9<μ<1. (Pick regularly-spaced δμ, run ntransient steps, record ncycles steps, and plot. After the routine is working, you should be able to push ntransient and ncycles both larger than 100, and δμ<0.01.)

On the same plot, for the same μs, plot the first eight images of x=1/2, that is, f(1/2),f(f(1/2)),.... Are the boundaries you see just the cusps? What happens in the bifurcation diagram when two boundaries touch?

points[μ_,ntransient_:100, ncycles_:512] := Block[{xa = ..., mus = Table[μ,{n,1,ncycles}]},
Transpose[{mus,NestList[...]}]]
BifurcationDiagram = ListPlot[Flatten[Table[points[μ],{μ,...}],1]];
cusps = Flatten[Table[{μ,Nest[f[#,μ]&,0.5,n]},{n,1,8},{μ,0...}],1];
cuspDiagram = ListPlot[cusps, PlotStyle->Red];
Show[BifurcationDiagram, cuspDiagram]