{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Invariant measures

\n",
"\n",
"(Sethna, \"Entropy, Order Parameters, and Complexity\", ex. 4.3)\n",
"\n",
"© 2017, James Sethna, all rights reserved. This exercise was developed in collaboration with Christopher Myers."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Liouville's theorem tells us that all available points in phase space are\n",
"equally weighted when a Hamiltonian system is averaged over all times. What\n",
"happens for systems that evolve according to laws that are not Hamiltonian?\n",
"Usually, the system does not continue to explore all points in its state\n",
"space; at long times it is confined to a subset of the original space known\n",
"as the *attractor*.\n",
"\n",
"Import packages"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%pylab inline\n",
"from scipy import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We consider the behavior of the `logistic' mapping from the unit\n",
"interval $(0,1)$ into itself:\n",
"\\begin{equation}\n",
"f(x) = 4 \\mu x (1-x).\n",
"\\end{equation}\n",
"(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\n",
"$x_0$, $f(x_0)$, $f(f(x_0))$, \\dots, $f^{[n]}(x_0)$, ... Iteration\n",
"can be thought of as a time step (one iteration of a Poincare return map\n",
"of the exercise \"Jupiter\" or one step $\\Delta t$ in a time-step algorithm."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def f(x,mu):\n",
" \"\"\"\n",
" Logistic map f(x) = 4 mu x (1-x), which folds the unit interval (0,1)\n",
" into itself.\n",
" \"\"\"\n",
" return 4*..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*Attracting fixed point.*\n",
"For small $\\mu$, our mapping has an attracting fixed-point. A fixed-point\n",
"of a mapping is a value $x^* = f(x^*)$; a fixed-point is stable if\n",
"small perturbations shrink after iterating:\n",
"\\begin{equation}\n",
"|f(x^* + \\epsilon) - x^*| \\approx |f'(x^*)| \\epsilon < \\epsilon,\n",
"\\end{equation}\n",
"which happens if the derivative $|f'(x^*)|<1$.\n",
"(For many-dimensional mappings, a sufficient criterion for\n",
" stability is that all the eigenvalues of the Jacobian have magnitude\n",
" smaller than one. A continuous time evolution $d{y}/d{t} = F(y)$\n",
" will be stable if $d{F}/d{y}$ is smaller than zero, or (for multidimensional\n",
" systems) if the Jacobian $DF$ has eigenvalues whose real parts are all\n",
" less than zero.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(a) **Iteration.** * Set $\\mu = 0.2$; iterate $f$ for some\n",
"initial points $0\n",
"\n",
"***Analytics.** * Find the non-zero fixed-point $x^*(\\mu)$ of the\n",
"map $f(x)$, and show that it exists and is stable for\n",
"$1/4 < \\mu < 3/4$. If you are ambitious or have a computer\n",
"algebra program, show that there is a stable, period-two cycle for\n",
"$3/4 < \\mu < (1+\\sqrt{6})/4$.*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = 0.77\n",
"for t in range(100):\n",
" x = ...(...,0.2)\n",
"print(\"x goes to\", x, \"when f(x,0.2) is iterated\")\n",
"\n",
"figure()\n",
"xs = arange(0,1,0.01)\n",
"plot(xs, xs)\n",
"for mu in (0.2, 0.4, 0.6):\n",
" plot(xs, f(...), label = \"mu=%g\" %mu)\n",
"legend(loc=2)\n",
"axes().set_aspect('equal')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An attracting fixed-point is the antithesis of Liouville's theorem; all\n",
"initial conditions are transient except one, and all systems lead eventually\n",
"to the same, time-independent state. (On the other hand, this is precisely\n",
"the behavior we expect in statistical mechanics *on the macroscopic scale*;\n",
"the system settles down into a time-independent equilibrium state! All\n",
"microstates are equivalent, but the vast majority of accessible microstates\n",
"have the same macroscopic behavior in most large systems.)\n",
"We could\n",
"define a rather trivial `equilibrium ensemble' for this system, which\n",
"consists of the single point $x^*$; any property $O(x)$ will have the\n",
"long-time average $\\langle O\\rangle = O(x^*)$.\n",
"\n",
"For larger values of $\\mu$, more complicated things happen. At $\\mu=1$, the\n",
"dynamics can be shown to fill the entire interval; the dynamics is ergodic,\n",
"and the attractor fills the entire set of available states. However, unlike the\n",
"case of Hamiltonian systems, not all states are weighted equally (i.e.,\n",
"Liouville's theorem does not hold).\n",
"\n",
"We can find time averages for functions of $x$ in two ways: by averaging\n",
"over time (many iterates of the map) or by weighting an integral over $x$\n",
"by the *invariant density* $\\rho(x)$. The invariant density $\\rho(x) \\,d{x}$\n",
"is the probability that a point on a long trajectory will\n",
"lie between $x$ and $x+d{x}$. To find it numerically, we iterate a typical point\n",
" (not an unstable fixed-point or unstable periodic orbit!)\n",
"$x_0$ a thousand or so times ($N_{\\mathrm{transient}}$) to find a point\n",
"$x_a$ on the attractor, and then collect a long trajectory of perhaps a\n",
"million points ($N_{\\mathrm{cycles}}$). A histogram of this trajectory\n",
"gives $\\rho(x)$. Averaging over this density is manifestly the same as\n",
"a time average over the trajectory of a million points. We call $\\rho(x)$\n",
"invariant because it is left the same under the mapping $f$;\n",
"iterating our million-point approximation for $\\rho$ once under $f$\n",
"only removes the first point $x_a$ and adds one extra point to the end.\n",
"\n",
"(b) **Invariant density.** *Set $\\mu = 1$; iterate $f$ many times,\n",
"and form a histogram of values giving the density $\\rho(x)$ of points\n",
"along the trajectory. You should find that points $x$ near the boundaries are\n",
"approached more often than points near the center.\n",
"\n",
"***Analytics.** Using the fact that the long-time average $\\rho(x)$ must\n",
"be independent of time, verify for $\\mu=1$ that the density of points is\n",
"\\begin{equation}\n",
"\\rho(x) = \\frac{1}{\\pi \\sqrt{x(1-x)}}.\n",
"\\end{equation}\n",
" (You need not derive the factor $1/\\pi$, which normalizes the\n",
" probability density to one.)\n",
"Plot this theoretical curve with your numerical histogram.\n",
"\n",
"(Hint: The points in a range $d{x}$ around a point $x$ map under $f$ to a range\n",
"$d{y} = f'(x)\\,d{x}$ around the image $y=f(x)$. Each iteration maps two\n",
"points $x_a$ and $x_b = 1-x_a$ to $y$, and thus maps all the density\n",
"$\\rho(x_a) |d{x}_a|$ and $\\rho(x_b) |d{x}_b|$ into $d{y}$.\n",
"Hence the probability $\\rho(y)\\,d{y}$ must equal\n",
"$\\rho(x_a) |d{x}_a| + \\rho(x_b) |d{x}_b|$, so\n",
"\\begin{equation}\n",
"\\rho(f(x_a)) = \\rho(x_a)/|f'(x_a)| + \\rho(x_b)/|f'(x_b)|.\n",
"\\end{equation}\n",
"Substitute the second-to-latest equation for $\\rho(x)$ above into this equation.\n",
"You will need to factor a polynomial.)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def IterateList(x,mu,Niter=10,Nskip=0):\n",
" \"\"\"\n",
" Iterate the function f(x,mu) Niter-1 times, starting at x \n",
" (or at x iterated Nskip times), so that the full trajectory \n",
" contains N points.\n",
" Returns the entire list\n",
" (x, f(x), f(f(x)), ... f(f(...(f(x))...))).\n",
"\n",
" Can be used to explore the dynamics starting from an arbitrary point\n",
" x0, or to explore the attractor starting from a point x0 on the\n",
" attractor (say, initialized using Nskip).\n",
"\n",
" For example, you can use Iterate to find a point xAttractor on the\n",
" attractor and IterateList to create a long series of points attractorXs\n",
" (thousands, or even millions long, if you're in the chaotic region),\n",
" and then use\n",
" hist(attractorXs, bins=2000, normed=1)\n",
" to see the density of points.\n",
" \"\"\"\n",
" for i in range(Nskip):\n",
" x = f(x,mu)\n",
" fiter = [x]\n",
" for i in range(Niter-1):\n",
" x = ...\n",
" fiter.append(x)\n",
" return fiter\n",
"\n",
"hist(IterateList(..., 1.0, Niter=1000000, Nskip=1000), bins = 2000, normed = 1);\n",
"\n",
"xs = arange(0.001,0.999,0.001)\n",
"plot(xs, ..., \"r-\");\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Mathematicians call this probability density $\\rho(x)\\,d{x}$ the\n",
"*invariant measure* on the attractor.\n",
" (There are actually many possible invariant measures on\n",
" some attractors; this one is the SRB measure (John Guckenheimer,\n",
" private communication).)\n",
"To get the long-term average of any function $O(x)$, one can use\n",
"\\begin{equation}\n",
"\\langle O \\rangle = \\int O(x) \\rho(x)\\,d{x}.\n",
"\\end{equation}\n",
"To a mathematician, a measure is a way of weighting different\n",
"regions when calculating integrals---precisely our $\\rho(x)\\,d{x}$.\n",
"Notice that, for the case of an attracting fixed-point, we would\n",
"have $\\rho(x) = \\delta(x^*)$.\n",
" (The case of a fixed-point then becomes mathematically a\n",
" measure with a point mass at $x^*$.)\n",
"\n",
"**Cusps in the invariant density.**\n",
"At values of $\\mu$ slightly smaller than one, our mapping has a rather complex\n",
"invariant density.\n",
"\n",
"(c) * Find the invariant density (as described above) for $\\mu=0.9$.\n",
"Make your trajectory length $N_{\\mathrm{cycles}}$ big enough and the\n",
"bin size small enough to see the interesting structures. Notice that\n",
"the attractor no longer fills the whole range\n",
"$(0,1)$; locate roughly where the edges are. Notice also the\n",
"cusps in $\\rho(x)$ at the edges of the attractor, and also at places\n",
"inside the attractor (called 'boundaries'). Locate\n",
"some of the more prominent cusps. *"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"hist(IterateList(..., 0.9, Niter=1000000, Nskip=1000), bins = 2000, normed = 1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Analytics of cusps.**\n",
"Notice that $f'(1/2)=0$, so by the equation for $\\rho(f(x_a))$ above\n",
"we know that $\\rho(f(x)) \\ge \\rho(x) / |f'(x)|$ must have a singularity\n",
"near $x=1/2$; all the points near $x=1/2$ are squeezed together and folded\n",
"to one side by $f$. Further iterates of this singularity produce more cusps;\n",
"the crease after one fold stays a crease after being further stretched and\n",
"kneaded.\n",
"\n",
"(d) *Set $\\mu = 0.9$.\n",
"Calculate $f(1/2)$, $f(f(1/2))$, ... and compare these iterates\n",
"to the locations of the edges and cusps from part~(c). (You may wish\n",
"to include them both on the same plot.)*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"hist(...);\n",
"\n",
"boundaries = IterateList(...,0.9,8,1)\n",
"plot(boundaries, 1.+arange(len(boundaries)), 'ro-')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Bifurcation diagram.**\n",
"The evolution of the attractor and its invariant density as $\\mu$ varies are\n",
"plotted in the *bifurcation diagram*. One of the striking features\n",
"in this plot are the sharp boundaries formed by the cusps.\n",
"\n",
"(e) **Bifurcation diagram.** *Plot the attractor as a function of $\\mu$,\n",
"for $0.9<\\mu<1$. (Pick regularly-spaced $\\delta \\mu$, run\n",
"$n_{\\mathrm{transient}}$ steps, record $n_{\\mathrm{cycles}}$ steps, and plot.\n",
"After the routine is working, you should be able to push\n",
"$n_{\\mathrm{transient}}$ and $n_{\\mathrm{cycles}}$ both larger than 100, and\n",
"$\\delta\\mu < 0.01$.)\n",
"\n",
"On the same plot, for the same $\\mu$s, plot the first eight images\n",
"of $x=1/2$, that is, $f(1/2), f(f(1/2)), ...$. Are the boundaries\n",
"you see just the cusps? What happens in the bifurcation diagram\n",
"when two boundaries touch? *"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def BifurcationDiagram(muMin=0., muMax=1., deltaMu=0.0002, x0=0.49, Nskip=100,\n",
" Niter=512, marker = '.', size=1, color='b'):\n",
" xlim((muMin,muMax))\n",
" muArray = arange(muMin, muMax, deltaMu)\n",
" mus = []\n",
" trajs = []\n",
" for mu in muArray:\n",
" mus.extend([mu]*Niter)\n",
" trajs.extend(IterateList(x0, mu, Niter, Nskip))\n",
" scatter(mus, trajs, marker = marker, s=size,c=color)\n",
" setp(axes().get_xticklabels(),fontsize=40)\n",
" setp(axes().get_yticklabels(),fontsize=40)\n",
" ylabel('x', fontsize=40)\n",
" xlabel('mu', fontsize=40)\n",
" return fig\n",
" \n",
"figure(figsize=(50,20))\n",
"BifurcationDiagram(0.9,1.0)\n",
"BifurcationDiagram(0.9,1.0,x0=...,Nskip=...,Niter=8,marker='o',color='red',\n",
" size=30);"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.1"
}
},
"nbformat": 4,
"nbformat_minor": 0
}