{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fractal dimensions

\n",
"\n",
"(Sethna, \"Entropy, Order Parameters, and Complexity\", ex. 5.16)\n",
"\n",
"© 2017, James Sethna, all rights reserved. This exercise was developed in collaboration with Christopher Myers."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are many strange sets that emerge in science. In statistical mechanics,\n",
"such sets often arise at continuous phase transitions, where self-similar\n",
"spatial structures arise (Chapter 12).\n",
"In chaotic dynamical systems, the attractor (the set of points occupied at long\n",
"times after the transients have disappeared) is often a fractal (called\n",
"a *strange attractor*). These sets are often tenuous and jagged, with\n",
"holes on all length scales, as in percolation (Fig. 1.2).\n",
"\n",
"We often try to characterize these strange sets by a dimension.\n",
"The dimensions of two extremely different sets can be the same;\n",
"the path exhibited by a random walk (embedded in three or more dimensions)\n",
"is arguably a two-dimensional set, but does not locally\n",
"look like a surface. However, if two sets have different spatial dimensions\n",
"(measured in the same way) they are certainly qualitatively different.\n",
"\n",
"There is more than one way to define a dimension. Roughly speaking,\n",
"strange sets are often spatially inhomogeneous, and what dimension you measure\n",
"depends upon how you weight different regions of the set. In this exercise,\n",
"we will calculate the *information dimension** (closely connected to\n",
"the non-equilibrium entropy), and the **capacity dimension**\n",
"(originally called the **Hausdorff dimension}, also sometimes called the\n",
"**fractal dimension**).\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": [
"To generate our strange set---along with some more ordinary sets---we\n",
"will use the logistic map\n",
"\\begin{equation}\n",
"f(x) = 4 \\mu x (1-x).\n",
"\\end{equation}\n",
"The attractor\n",
"for the logistic map is a periodic orbit (dimension zero) at $\\mu=0.8$,\n",
"and a chaotic, cusped density filling two intervals (dimension one)\n",
"at $\\mu=0.9$.\n",
" (See the exercise 'Invariant measures'. The chaotic\n",
" region for the logistic map does not have a strange attractor because\n",
" the map is confined to one dimension; period-doubling cascades for\n",
" dynamical systems in higher spatial dimensions have\n",
" fractal, strange attractors in the chaotic region.\n",
"At the onset of chaos at $\\mu=\\mu_\\infty\\approx 0.892486418$\n",
"(see the exercise 'Period doubling') the dimension becomes intermediate\n",
"between zero and one; this strange, self-similar set is called the\n",
"**Feigenbaum attractor**."
]
},
{
"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": [
"Both the information dimension and the capacity dimension are defined\n",
"in terms of the occupation $P_n$ of cells of size $\\epsilon$\n",
"in the limit as $\\epsilon\\to 0$.\n",
"\n",
"(a) ** Write a routine which, given $\\mu$ and a set of bin sizes $\\epsilon$,\n",
"does the following:\n",
"* Iterates $f$ hundreds or thousands of times (to get onto the attractor).\n",
"* Iterates $f$ a large number $N_{\\mathrm{tot}}$ more times, collecting points\n",
"on the attractor.\n",
"(For $\\mu\\le\\mu_\\infty$, you could just integrate $2^n$ times for $n$ fairly\n",
"large.)\n",
"* For each $\\epsilon$, use a histogram to calculate the probability $P_j$\n",
"that the points fall in the $j$th bin.\n",
"* Return the set of vectors $P_j[\\epsilon]$.\n",
"*\n",
"You may wish to test your routine by using it for $\\mu=1$ (where the\n",
"distribution should look like $\\rho(x) = {1}/{\\pi \\sqrt{x(1-x)}}$,\n",
"see the exercise 'Invariant measures') and $\\mu=0.8$ (where the distribution\n",
"should look like two $\\delta$-functions, each with half of the points)."
]
},
{
"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",
"def GetPn(mu, epsilonList, Niter, Nskip=10000):\n",
" \"\"\"\n",
" Generates probability arrays P_n[epsilon].\n",
" Specifically,\n",
" finds a point on the attractor by iterating Nskip times\n",
" collects points on the attractor of size Niter\n",
" for each epsilon in epsilonList,\n",
" generates bins of size epsilon for the range (0,1) of the function\n",
" bins = scipy.arange(0.0,1.0+eps,eps)\n",
" finds the number of points from the sample in each bin, using\n",
" the histogram function\n",
" numbers, bins = pylab.mlab.hist(sample, bins=bins)\n",
" and computes the probability P_n[epsilon] of being in each bin.\n",
" In the period doubling region the sample should of size 2^n so that\n",
" it covers the attractor evenly.\n",
" \"\"\"\n",
" sample = IterateList(0.1, mu, Niter, Nskip)\n",
" P_n = {}\n",
" for eps in epsilonList:\n",
" bins = arange(0.0, 1.0 + eps, eps)\n",
" numbers, bins = histogram(sample, bins=bins)\n",
" P_n[eps] = ... # Probability\n",
" return P_n\n",
"\n",
"Pn = GetPn(0.8,[0.001],10000)\n",
"plot(Pn[0.001])\n",
"figure()\n",
"Pn = GetPn(1.0,[0.001],10000)\n",
"plot(Pn[0.001])"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"**The capacity dimension.**\n",
"The definition of the capacity dimension is motivated by the idea that\n",
"it takes at least\n",
"\\begin{equation}\n",
"N_{\\mathrm{cover}} = V/\\epsilon^D\n",
"\\end{equation}\n",
"bins of size $\\epsilon^D$ to cover a $D$-dimensional set of volume $V$.\n",
" (Imagine covering the surface of a sphere in 3D with tiny\n",
" cubes; the number of cubes will go as the surface area (2D volume)\n",
" divided by $\\epsilon^2$.)\n",
"By taking logs of both sides we find\n",
"$\\log N_{\\mathrm{cover}} \\approx \\log V + D \\log \\epsilon$.\n",
"The capacity dimension is defined as the limit\n",
"\\begin{equation}\n",
"D_{\\mathrm{capacity}} =\n",
" \\lim_{\\epsilon\\to 0} \\frac{\\log N_{\\mathrm{cover}}}{\\log \\epsilon},\n",
"\\end{equation}\n",
"but the convergence is slow (the error goes roughly as $\\log V/\\log \\epsilon$).\n",
"Faster convergence is given by calculating the slope of $\\log N$ versus\n",
"$\\log \\epsilon$:\n",
"\\begin{align}\n",
"D_{\\mathrm{capacity}} &=\n",
"\\lim_{\\epsilon\\to 0} \\frac{d{\\log N_{\\mathrm{cover}}}}{d{\\log \\epsilon}}\\nonumber \\\\\n",
"&= \\lim_{\\epsilon\\to 0}\\frac{\\log N_{j+1}-\\log N_j}\n",
" {\\log \\epsilon_{i+1}-\\log \\epsilon_i}.\n",
"\\end{align}\n",
"\n",
"(b) *Use your routine from part~(a), write a routine to calculate\n",
"$N[\\epsilon]$ by counting non-empty bins. Plot $D_{\\mathrm{capacity}}$\n",
"from the fast convergence versus the midpoint\n",
"$(1/2)(\\log \\epsilon_{i+1}+\\log \\epsilon_i)$.\n",
"Does it appear to extrapolate to $D=1$ for $\\mu = 0.9$?\\,%\n",
" (In the chaotic regions, keep the number of bins\n",
" small compared to the number of iterates in your sample, or you will\n",
" start finding empty bins between points and eventually get a dimension\n",
" of zero.)\n",
"Does it appear to extrapolate to $D=0$ for $\\mu=0.8$? Plot these two curves\n",
"together with the curve for $\\mu_\\infty$. Does the last one appear to\n",
"converge to $D_1 \\approx 0.538$, the capacity dimension for the Feigenbaum\n",
"attractor gleaned from the literature? How small a deviation from $\\mu_\\infty$\n",
"does it take to see the numerical cross-over to integer dimensions?\n",
"*\n",
"\n",
"**Entropy and the information dimension.**\n",
"The probability density\n",
"$\\rho(x_j) \\approx P_j/\\epsilon\n",
" = ({1}/{\\epsilon})({N_j}/{N_{\\mathrm{tot}}})$. Converting\n",
"the non-equilibrium entropy formula to a sum gives\n",
"\\begin{align}\n",
"S &= -k_B \\int \\rho(x) \\log(\\rho(x)) \\,d{x} \\nonumber\\\\\n",
" &\\approx -\\sum_j \\frac{P_j}{\\epsilon}\n",
" \\log\\left(\\frac{P_j}{\\epsilon}\\right)\\epsilon\\nonumber\\\\\n",
" &= -\\sum_j P_j \\log P_j + \\log \\epsilon\n",
"\\end{align}\n",
"(setting the conversion factor $k_B = 1$ for convenience).\n",
"\n",
"You might imagine that the entropy for a fixed-point would be zero, and the\n",
"entropy for a period-$m$ cycle would be $k_B \\log m$. But this is incorrect;\n",
"when there is a fixed-point or a periodic limit cycle, the attractor is on a\n",
"set of dimension zero (a bunch of points) rather than dimension one. The\n",
"entropy must go to minus infinity---since we have precise information\n",
"about where the trajectory sits at long times. To estimate the\n",
"`zero-dimensional' entropy $k_B \\log m$ on the computer, we would\n",
"use the discrete form of the entropy\n",
"summing over bins $P_j$ instead of integrating over $x$:\n",
"\\begin{equation}\n",
"S_{d=0} = -\\sum_j P_j \\log (P_j) = S_{d=1} - \\log(\\epsilon).\n",
"\\end{equation}\n",
"More generally, the `natural' measure of the entropy for a set with $D$\n",
"dimensions might be defined as\n",
"\\begin{equation}\n",
"S_D = -\\sum_j P_j \\log (P_j) + D \\log(\\epsilon).\n",
"\\end{equation}\n",
"Instead of using this formula to define the entropy, mathematicians use\n",
"it to define the information dimension\n",
"\\begin{equation}\n",
"D_{\\mathrm{inf}} =\n",
" \\lim_{\\epsilon\\to0} \\left(\\sum P_j \\log P_j\\right)/\\log(\\epsilon).\n",
"\\end{equation}\n",
"The information dimension agrees with the ordinary dimension for sets that\n",
"locally look like $\\mathbb{R}^D$. It is different from the capacity\n",
"dimension, which counts each occupied\n",
"bin equally; the information dimension counts heavily occupied parts\n",
"(bins) in the attractor more heavily. Again, we can speed up the\n",
"convergence by noting that the equation for the information dimension\n",
"says\n",
"that $\\sum_j P_j \\log P_j$ is a linear function of $\\log \\epsilon$ with\n",
"slope $D$ and intercept $S_D$. Measuring the slope directly, we find\n",
"\\begin{equation}\n",
"D_{\\mathrm{inf}} =\n",
"\\lim_{\\epsilon\\to 0} \\frac{d{\\sum_j P_j(\\epsilon) \\log P_j(\\epsilon)}}\n",
" {d{\\log \\epsilon}}.\n",
"\\end{equation}\n",
"\n",
"(c) *As in part~(b), write a routine that plots $D_{\\mathrm{inf}}$ using\n",
"the fast definition as a function of the midpoint\n",
"$\\log \\epsilon$, as we increase the number of bins.\n",
"Plot the curves for $\\mu=0.9$, $\\mu=0.8$, and $\\mu_\\infty$. Does the\n",
"information dimension agree with the ordinary one for the first two?\n",
"Does the last one appear to converge to $D_1 \\approx 0.517098$, the\n",
"information dimension for the Feigenbaum attractor from the literature?\n",
"*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def DimensionEstimates(mu, epsilonList, Niter):\n",
" \"\"\"\n",
" Estimates the capacity dimension and the information dimension\n",
" for a sample of points on the line.\n",
" The capacity dimension is defined as\n",
" D_capacity = lim {eps->0} (- log(Nboxes) / log(eps))\n",
" but converges faster as\n",
" D_capacity = - (log(Nboxes[i+1])-log(Nboxes[i]))\n",
" / (log(eps[i+1])-log(eps[i]))\n",
" where Nboxes is the number of intervals of size eps needed to\n",
" cover the space. The information dimension is defined as\n",
" D_inf = lim {eps->0} sum(P_n log P_n) / log(eps)\n",
" but converges faster as\n",
" S0 = sum(P_n log P_n)\n",
" D_inf = - (S0[i+1]-S0[i])\n",
" / (log(eps[i+1])-log(eps[i]))\n",
" where P_n is the fraction of the list 'sample' that is in bin n,\n",
" and the bins are of size epsilon. You'll need to add a small\n",
" increment delta to P_n before taking the log: delta = 1.e-100 will\n",
" not change any of the non-zero elements, and P_n log (P_n + delta)\n",
" will be zero if P_n is zero.\n",
"\n",
" Returns three lists, with epsilonBar (geometric mean of neighboring\n",
" epsilonList values), and D_inf, and D_capacity values for each\n",
" epsilonBar\n",
" \"\"\"\n",
" D_inf = []\n",
" D_capacity = []\n",
" epsilonBar = []\n",
" delta = 1.e-100 # Add to make log finite\n",
"\n",
" P_n = GetPn(mu, epsilonList, Niter)\n",
"\n",
" Nboxes = [] # Number of non-zero P_n\n",
" S0 = [] # Zero-dimensional entropy -sum(P_n log(P_n))\n",
"\n",
" for eps in epsilonList:\n",
" Nboxes.append(sum(P_n[eps] > 0))\n",
" S0.append(-sum(... * log(... + delta)))\n",
"\n",
" epsBar = []\n",
" D_capacity = []\n",
" D_inf = []\n",
" for i in range(len(epsilonList) - 1):\n",
" epsi = epsilonList[i]\n",
" epsiP1 = epsilonList[i + 1]\n",
" epsBar.append(sqrt(epsiP1 * epsi))\n",
" D_capacity_estimate = ...\n",
" D_capacity.append(D_capacity_estimate)\n",
" D_inf_estimate = ...\n",
" D_inf.append(D_inf_estimate)\n",
"\n",
" return epsBar, D_capacity, D_inf\n",
"\n",
"muInfinity = 0.892486418\n",
"def PlotDimensionEstimates(mu=muInfinity, nIter=2**18,\n",
" epsilonList=2.0**arange(-4, -16, -1)):\n",
" epsBar, DCapacity, DInformation = DimensionEstimates(mu,epsilonList,Niter)\n",
" figure()\n",
" semilogx(epsBar, DCapacity, label='Capacity dim')\n",
" semilogx(epsBar, ..., label='Info dim')\n",
" title('Fractal dimension estimates, mu = '+ str (mu))\n",
" ylabel('Dimension')\n",
" xlabel('Bin size')\n",
" legend(loc=3)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"PlotDimensionEstimates(0.9)\n",
"PlotDimensionEstimates(0.8)\n",
"PlotDimensionEstimates(muInfinity)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Most 'real world' fractals have a whole spectrum of different characteristic\n",
"spatial dimensions; they are *multifractal**."
]
}
],
"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
}
*