(Sethna, "Entropy, Order Parameters, and Complexity", ex. 4.3)
© 2017, James Sethna, all rights reserved. This exercise was developed in collaboration with Christopher Myers.
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.
Import packages
%matplotlib inline
from matplotlib.pyplot import plot, figure, show, axes, legend, hist
from matplotlib.pyplot import xlim, scatter, setp, xlabel, ylabel, xticks, yticks
from numpy import *
We consider the behavior of the `logistic' mapping from the unit interval $(0,1)$ into itself: \begin{equation} f(x) = 4 \mu x (1-x). \end{equation} (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 $x_0$ as the sequence of points $x_0$, $f(x_0)$, $f(f(x_0))$, \dots, $f^{[n]}(x_0)$, ... Iteration can be thought of as a time step (one iteration of a Poincare return map of the exercise "Jupiter" or one step $\Delta t$ in a time-step algorithm.
def f(x,mu):
"""
Logistic map f(x) = 4 mu x (1-x), which folds the unit interval (0,1)
into itself.
"""
return 4*...
Attracting fixed point. For small $\mu$, 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: \begin{equation} |f(x^* + \epsilon) - x^*| \approx |f'(x^*)| \epsilon < \epsilon, \end{equation} 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 $d{y}/d{t} = F(y)$ will be stable if $d{F}/d{y}$ 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 $\mu = 0.2$; iterate $f$ for some initial points $0<x_0<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 $\mu = 0.4$, and $0.6$. What happens?
Analytics. Find the non-zero fixed-point $x^*(\mu)$ of the map $f(x)$, and show that it exists and is stable for $1/4 < \mu < 3/4$. If you are ambitious or have a computer algebra program, show that there is a stable, period-two cycle for $3/4 < \mu < (1+\sqrt{6})/4$.
x = 0.77
for t in range(100):
x = ...(...,0.2)
print("x goes to", x, "when f(x,0.2) is iterated")
figure()
xs = arange(0,1,0.01)
plot(xs, xs)
for mu in (0.2, 0.4, 0.6):
plot(xs, f(...), label = "mu=%g" %mu)
legend(loc=2)
axes().set_aspect('equal')
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 $\langle O\rangle = O(x^*)$.
For larger values of $\mu$, more complicated things happen. At $\mu=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 $\rho(x)$. The invariant density $\rho(x) \,d{x}$ is the probability that a point on a long trajectory will lie between $x$ and $x+d{x}$. To find it numerically, we iterate a typical point (not an unstable fixed-point or unstable periodic orbit!) $x_0$ a thousand or so times ($N_{\mathrm{transient}}$) to find a point $x_a$ on the attractor, and then collect a long trajectory of perhaps a million points ($N_{\mathrm{cycles}}$). A histogram of this trajectory gives $\rho(x)$. Averaging over this density is manifestly the same as a time average over the trajectory of a million points. We call $\rho(x)$ invariant because it is left the same under the mapping $f$; iterating our million-point approximation for $\rho$ once under $f$ only removes the first point $x_a$ and adds one extra point to the end.
(b) Invariant density. Set $\mu = 1$; iterate $f$ many times, and form a histogram of values giving the density $\rho(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 $\rho(x)$ must be independent of time, verify for $\mu=1$ that the density of points is \begin{equation} \rho(x) = \frac{1}{\pi \sqrt{x(1-x)}}. \end{equation} (You need not derive the factor $1/\pi$, which normalizes the probability density to one.) Plot this theoretical curve with your numerical histogram.
(Hint: The points in a range $d{x}$ around a point $x$ map under $f$ to a range $d{y} = f'(x)\,d{x}$ around the image $y=f(x)$. Each iteration maps two points $x_a$ and $x_b = 1-x_a$ to $y$, and thus maps all the density $\rho(x_a) |d{x}_a|$ and $\rho(x_b) |d{x}_b|$ into $d{y}$. Hence the probability $\rho(y)\,d{y}$ must equal $\rho(x_a) |d{x}_a| + \rho(x_b) |d{x}_b|$, so \begin{equation} \rho(f(x_a)) = \rho(x_a)/|f'(x_a)| + \rho(x_b)/|f'(x_b)|. \end{equation} Substitute the second-to-latest equation for $\rho(x)$ above into this equation. You will need to factor a polynomial.)def IterateList(x,mu,Niter=10,Nskip=0):
"""
Iterate the function f(x,mu) Niter-1 times, starting at x
(or at x iterated Nskip times), so that the full trajectory
contains N points.
Returns the entire list
(x, f(x), f(f(x)), ... f(f(...(f(x))...))).
Can be used to explore the dynamics starting from an arbitrary point
x0, or to explore the attractor starting from a point x0 on the
attractor (say, initialized using Nskip).
For example, you can use Iterate to find a point xAttractor on the
attractor and IterateList to create a long series of points attractorXs
(thousands, or even millions long, if you're in the chaotic region),
and then use
hist(attractorXs, bins=2000, density=True)
to see the density of points.
"""
for i in range(Nskip):
x = f(x,mu)
fiter = [x]
for i in range(Niter-1):
x = ...
fiter.append(x)
return fiter
hist(IterateList(..., 1.0, Niter=1000000, Nskip=1000), bins = 2000, density = True);
xs = arange(0.001,0.999,0.001)
plot(xs, ..., "r-");
Mathematicians call this probability density $\rho(x)\,d{x}$ 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 \begin{equation} \langle O \rangle = \int O(x) \rho(x)\,d{x}. \end{equation} To a mathematician, a measure is a way of weighting different regions when calculating integrals---precisely our $\rho(x)\,d{x}$. Notice that, for the case of an attracting fixed-point, we would have $\rho(x) = \delta(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 $\mu$ slightly smaller than one, our mapping has a rather complex invariant density.
(c) Find the invariant density (as described above) for $\mu=0.9$. Make your trajectory length $N_{\mathrm{cycles}}$ 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 $\rho(x)$ at the edges of the attractor, and also at places inside the attractor (called 'boundaries'). Locate some of the more prominent cusps.
hist(IterateList(..., 0.9, Niter=1000000, Nskip=1000), bins = 2000, density = True);
Analytics of cusps. Notice that $f'(1/2)=0$, so by the equation for $\rho(f(x_a))$ above we know that $\rho(f(x)) \ge \rho(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 $\mu = 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.)
hist(...);
boundaries = IterateList(...,0.9,8,1)
plot(boundaries, 1.+arange(len(boundaries)), 'ro-')
Bifurcation diagram. The evolution of the attractor and its invariant density as $\mu$ 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 $\mu$, for $0.9<\mu<1$. (Pick regularly-spaced $\delta \mu$, run $n_{\mathrm{transient}}$ steps, record $n_{\mathrm{cycles}}$ steps, and plot. After the routine is working, you should be able to push $n_{\mathrm{transient}}$ and $n_{\mathrm{cycles}}$ both larger than 100, and $\delta\mu < 0.01$.)
On the same plot, for the same $\mu$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?
def BifurcationDiagram(muMin=0., muMax=1., deltaMu=0.0002, x0=0.49, Nskip=100,
Niter=512, marker = '.', size=1, color='b'):
xlim((muMin,muMax))
muArray = arange(muMin, muMax, deltaMu)
mus = []
trajs = []
for mu in muArray:
mus.extend([mu]*Niter)
trajs.extend(IterateList(x0, mu, Niter, Nskip))
scatter(mus, trajs, marker = marker, s=size,c=color)
yticks(fontsize=40)
xticks(fontsize=40)
ylabel('x', fontsize=40)
xlabel('mu', fontsize=40)
figure(figsize=(50,20))
BifurcationDiagram(0.9,1.0)
BifurcationDiagram(0.9,1.0,x0=...,Nskip=...,Niter=8,marker='o',color='red',
size=30);