{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Period doubling

\n",
"\n",
"(Sethna, \"Entropy, Order Parameters, and Complexity\", ex. 12.9)\n",
"\n",
"© 2017, James Sethna, all rights reserved."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this exercise, we use renormalization-group and scaling methods to\n",
"study the *onset of chaos*. There are several routes by which a\n",
"dynamical system can start exhibiting chaotic motion; this exercise\n",
"studies the *period-doubling cascade*, first extensively investigated by Feigenbaum.\n",
"\n",
"Import packages"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%pylab inline\n",
"from scipy import *\n",
"from scipy.optimize import brentq"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Chaos is often associated with dynamics which stretch and fold; when\n",
"a batch of taffy is being pulled, the motion of a speck in the taffy depends\n",
"sensitively on the initial conditions. A simple representation of this\n",
"physics is provided by the map\n",
"\\begin{equation}\n",
"f(x) = 4 \\mu x (1-x)\n",
"\\end{equation}\n",
"restricted to the domain $(0,1)$. It takes $f(0) = f(1) = 0$, and\n",
"$f(1/2)=\\mu$.\n",
"Thus, for $\\mu=1$ it precisely folds the unit interval in half, and stretches\n",
"it to cover the original domain."
]
},
{
"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 ..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The study of dynamical systems (*e.g.*, differential equations and\n",
"maps like the logistic map $f(x)$ above often focuses on the behavior after\n",
"long times, where the trajectory moves along the\n",
"*attractor*.\n",
"We can study the onset and behavior of chaos in our system by observing the\n",
"evolution of the attractor as we change $\\mu$. For small enough $\\mu$,\n",
"all points shrink to the origin; the origin\n",
"is a stable fixed-point which attracts the entire interval $x\\in(0,1)$. For\n",
"larger $\\mu$, we first get a stable fixed-point inside the interval, and\n",
"then *period doubling*."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(a) Iteration *Set $\\mu = 0.2$; iterate $f$ for some\n",
"initial points $x_0$ of your choosing, and convince yourself that they all\n",
"are attracted to zero. Plot $f$ and the diagonal $y=x$ on the same plot.\n",
"Are there any fixed-points other than $x=0$? Repeat for\n",
"$\\mu = 0.3$, $\\mu = 0.7$, and $0.8$. What happens?*\n",
"\n",
"Note: We write functions like Iterate abstractly, in terms of a function g(x,eta), so that we can also examine $f_{\\rm sin}(x,B)$ where $\\eta=B$ instead of $\\eta=\\mu$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def Iterate(x0, N, eta, g=f):\n",
" \"\"\"\n",
" Iterate the function g(x,eta) N times, starting at x=x0.\n",
" Return g(g(...(g(x))...)). Used to find a point on the attractor\n",
" starting from some arbitrary point x0.\n",
"\n",
" Calling Iterate for the Feigenbaum map f^[1000] (x,0.9) at mu=0.9 would look like\n",
" Iterate(x, 1000, 0.9, f)\n",
"\n",
" We'll later be using Iterate to study the sine map\n",
" fSin(x,B) = B sin(pi x)\n",
" so passing in the function and arguments will be necessary for\n",
" comparing the logistic map f to fSin.\n",
"\n",
" Inside Iterate you'll want to apply g(x0, eta).\n",
" \"\"\"\n",
" for i in range(N):\n",
" x0 = ...\n",
" return x0\n",
"\n",
"[print(mu, \"->\", Iterate(random.rand(4),100,mu,f)) for mu in [0.2,0.3,0.8,0.9]];"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"figure()\n",
"xs = arange(0,1,0.01)\n",
"plot(xs, xs)\n",
"for mu in (0.2, 0.3, 0.8, 0.9):\n",
" plot(xs, f(...), label = \"mu=%g\" %mu)\n",
"legend(loc=2)\n",
"axes().set_aspect('equal')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*On the same graph, plot $f$, the diagonal $y=x$, and the segments\n",
"$\\{x_0,x_0\\}$, $\\{x_0,f(x_0)\\}$, $\\{f(x_0),f(x_0)\\}$,\n",
"$\\{f(x_0),f(f(x_0))\\}$, ...\n",
"(representing the convergence of the trajectory to the attractor).\n",
"See how $\\mu = 0.7$ and $0.8$ differ. Try other values of $\\mu$.*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def IterateList(x,eta,Niter=10,Nskip=0,g=f):\n",
" \"\"\"\n",
" Iterate the function g(x, eta) 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, g(x), g(g(x)), ... g(g(...(g(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",
" pylab.hist(attractorXs, bins=500, normed=1)\n",
" pylab.show()\n",
" to see the density of points.\n",
" \"\"\"\n",
" x = Iterate(x,Nskip,eta,g)\n",
" xs = [x]\n",
" for i in range(Niter-1):\n",
" x = ...\n",
" xs.append(x)\n",
" return xs\n",
"\n",
"def PlotIterate(mu,Niter=100,Nskip=0,x0=0.49):\n",
" \"\"\"\n",
" Plots g, the diagonal y=x, and the boxes made of the segments\n",
" [[x0,x0], [x0, g(x0)], [g(x0), g(x0)], [g(x0), g(g(x0))], ...\n",
" \n",
" Notice the xs and the ys are just the trajectory with each point \n",
" repeated twice, where the xs drop the final point and the ys drop\n",
" the initial point\n",
" \"\"\"\n",
" figure()\n",
" title('mu = '+ str (mu))\n",
" xarray = arange(0.,1.,0.01)\n",
" plot(xarray,xarray,'r-') # Plot diagonal\n",
" plot(xarray,...,'g-',linewidth=3) # Plot function\n",
" traj = IterateList(x0,mu,Niter,Nskip)\n",
" doubletraj = array([traj,traj]).transpose().flatten()\n",
" xs = doubletraj[...] # Drops last point\n",
" ys = doubletraj[...] # Drops first point\n",
" plot(xs, ys, 'b-', linewidth=1, antialiased=True)\n",
" axes().set_aspect('equal')\n",
" \n",
"for mu in (0.7, 0.75, 0.8, 0.9):\n",
" PlotIterate(mu)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*By iterating the map many times, find a point $a_0$ on the\n",
"attractor. As above, then plot the successive iterates of $a_0$ for\n",
"$\\mu = 0.7$, $0.8$, $0.88$, $0.89$, $0.9$, and $1.0$.*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"for mu in (0.7, 0.75, 0.8, 0.9):\n",
" PlotIterate(mu, Nskip=...)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see at higher $\\mu$ that the system no longer settles into a\n",
"stationary state at long times. The fixed-point where $f(x)=x$ exists for\n",
"all $\\mu>1/4$, but for larger $\\mu$ it is no longer *stable*.\n",
"If $x^*$ is a fixed-point (so $f(x^*)=x^*$) we can add a small perturbation\n",
"$f(x^*+\\epsilon) \\approx f(x^*) + f'(x^*) \\epsilon = x^* + f'(x^*) \\epsilon$;\n",
"the fixed-point is stable (perturbations die away) if $|f'(x^*)|<1$.\n",
" (In a continuous evolution, perturbations die away if\n",
" the Jacobian of the derivative at the fixed-point has all negative\n",
" eigenvalues. For mappings, perturbations die away if all eigenvalues\n",
" of the Jacobian have magnitude less than one.)\n",
"\n",
"In this particular case, once the fixed-point goes unstable the motion\n",
"after many iterations becomes periodic, repeating itself after *two*\n",
"iterations of the map---so $f(f(x))$ has two new fixed-points. This\n",
"is called *period doubling*. Notice that\n",
"by the chain rule $d{\\,f(f(x))}/d{x} = f'(x) f'(f(x))$, and indeed\n",
"\\begin{align}\n",
"\\frac{d{\\,f^{[N]}}}{d{x}} &= \\frac{d{\\,f(f(\\dots f(x)\\dots))}}{d{x}}\\\\\n",
" &= f'(x)f'(f(x))\\dots f'(f(\\dots f(x)\\dots)), \\nonumber\n",
"\\end{align}\n",
"so the stability of a period-$N$ orbit is determined by the product of the\n",
"derivatives of $f$ at each point along the orbit."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(b) Analytics: *Find the 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 the period-two cycle is stable for\n",
"$3/4 < \\mu < (1+\\sqrt{6})/4$.*\n",
"\n",
"\n",
"(c) Bifurcation diagram: *Plot the attractor as a function of $\\mu$,\n",
"for $0<\\mu<1$..\n",
"(Pick regularly-spaced $\\delta \\mu$, run $N_{\\mathrm{transient}}$ steps,\n",
"record $N_{\\mathrm{cycles}}$ steps, and plot. After the routine is working,\n",
"you should be able to push $N_{\\mathrm{transient}}$ and $N_{\\mathrm{cycles}}$\n",
"both larger than 100, and $\\delta\\mu < 0.01$.) \n",
"Also on the bifurcation\n",
"diagram, plot the line $x=1/2$ where $f(x)$ reaches its maximum.\n",
"*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def BifurcationDiagram(etaMin=0., etaMax=1., deltaEta=0.001, g=f, x0=0.49, Nskip=100, Niter=128):\n",
" figure(figsize=(50,20))\n",
" xlim((etaMin,etaMax))\n",
" etaArray = arange(etaMin, etaMax, deltaEta)\n",
" etas = []\n",
" trajs = []\n",
" for eta in etaArray:\n",
" etas.extend([eta]*Niter)\n",
" trajs.extend(IterateList(...))\n",
" scatter(etas, trajs, marker = '.', s=0.2)\n",
" plot([0.,1.],[0.5,0.5],'g-')\n",
" setp(axes().get_xticklabels(),fontsize=40)\n",
" setp(axes().get_yticklabels(),fontsize=40)\n",
" ylabel('x', fontsize=40)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"BifurcationDiagram()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"BifurcationDiagram(...[Zoom in; decrease deltaEta and increase Nskip until it looks nice])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*\n",
"Also plot the attractor for another one-humped map\n",
"\\begin{equation}\n",
"f_{\\mathrm{sin}}(x) = B \\sin(\\pi x),\n",
"\\end{equation}\n",
"for $0***"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def fSin(x, B):\n",
" return ...\n",
"\n",
"BifurcationDiagram(g=fSin)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice the complex, structured, chaotic region for large $\\mu$. How do we get from a stable\n",
"fixed-point $\\mu<3/4$ to chaos? The onset of chaos in this system\n",
"occurs through a cascade of ***period doublings*.\n",
"There is the sequence of bifurcations as $\\mu$ increases---the\n",
"period-two cycle starting at $\\mu_1=3/4$, followed by a\n",
"period-four cycle starting at $\\mu_2$, period-eight at $\\mu_3$---a whole\n",
"period-doubling cascade. The convergence appears geometrical, to a fixed-point\n",
"$\\mu_\\infty$:\n",
"\\begin{equation}\n",
"\\mu_n \\approx \\mu_\\infty - A \\delta^{-n},\n",
"\\end{equation}\n",
"so\n",
"\\begin{equation}\n",
"\\delta = \\lim_{n\\to\\infty} (\\mu_{n-1} - \\mu_{n-2})/(\\mu_n - \\mu_{n-1})\n",
"\\end{equation}\n",
"and there is a similar geometrical self-similarity along the $x$ axis,\n",
"with a (negative) scale factor $\\alpha$ relating each generation of the\n",
"tree.\n",
"\n",
"In the exercise 'Invariant Measures', we explained the boundaries in\n",
"the chaotic region as images of $x=1/2$. These special points are also\n",
"convenient for studying period-doubling.\n",
"Since $x=1/2$ is the maximum in the curve, $f'(1/2)=0$.\n",
"If it were a fixed-point (as it is for $\\mu=1/2$), it would not only\n",
"be stable, but unusually so: a shift by $\\epsilon$ away from the fixed\n",
"point converges after one step of the map to a distance\n",
"$\\epsilon f'(1/2) + \\epsilon^2/2 f''(1/2) = O(\\epsilon^2)$.\n",
"We say that such a fixed-point is *superstable*. (The superstable points \n",
"are the values of $\\mu$ in the figure above which intersect the green line x=1/2.) \n",
"If we have a period-$N$ orbit that passes through $x=1/2$, so that\n",
"the $N$th iterate\n",
"$f^N(1/2) \\equiv f(\\dots f(1/2) \\dots ) = 1/2$, then the orbit is also\n",
"superstable, since the derivative\n",
"of the iterated map is the product of the derivatives along the orbit,\n",
"and hence is also zero.\n",
"\n",
"These superstable points happen roughly half-way between the period-doubling\n",
"bifurcations, and are easier to locate, since we know that\n",
"$x=1/2$ is on the orbit. Let us use them to investigate the geometrical\n",
"convergence and self-similarity of the period-doubling bifurcation diagram\n",
"from part~(d). We will measure both the superstable values of $\\mu$ and the \n",
"size of the centermost 'leaf' in the bifurcation diagram (crossed by the line\n",
"$x=1/2$ where $f(x)$ takes its maximum).\n",
"For this part and part~(h), you will need a routine that\n",
"finds the roots $G(y)=0$ for functions $G$ of one variable $y$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(d) The Feigenbaum numbers and universality: * Numerically,\n",
"find the values of $\\mu^s_n$ at which the $2^n$-cycle is superstable (the intersections \n",
"of the attractor with the green line $x=1/2$), for\n",
"the first few values of $n$.* (Hint: Define a function\n",
"$G(\\mu) = f_\\mu^{[2^n]}(1/2)-1/2$, and find the root as a function\n",
"of $\\mu$. In searching for $\\mu^s_{n+1}$, you will want to search in a range\n",
"$(\\mu^s_n+\\epsilon, \\mu^s_{n}+(\\mu^s_n-\\mu^s_{n-1})/A)$ where\n",
"$A\\sim3$ works pretty well. Calculate $\\mu^s_0$ and $\\mu^s_1$ by hand.)\n",
"Also, find the separation $1/2 - f^{[n-1]}(1/2,\\mu^s_{n})$, the opening between the two edges of\n",
"the leaf crossed by the maximum of $f$ (green line above).\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def FindSuperstable(g, Niter, etaMin, etaMax, xMax=0.5):\n",
" \"\"\"\n",
" Finds superstable orbit g^[nIterated](xMax, eta) = xMax\n",
" in range (etaMin, etaMax).\n",
" Must be started with g-xMax of different sign at etaMin, etaMax.\n",
"\n",
" Notice that nIterated will usually be 2^n. (Sorry for using the\n",
" variable n both places!)\n",
"\n",
" Uses optimize.brentq, defining a temporary function, G(eta)\n",
" which is zero for the superstable value of eta.\n",
"\n",
" G iterates g nIterated times and subtracts xMax\n",
" (basically Iterate - xMax, except with only one argument eta)\n",
" \"\"\"\n",
" def G(eta):\n",
" return Iterate(...)-xMax\n",
" eta_root = brentq(G, etaMin, etaMax)\n",
" return eta_root\n",
"\n",
"\n",
"def GetSuperstablePointsAndIntervals(g, eta0, eta1, nMax = 11, xMax=0.5):\n",
" \"\"\"\n",
" Given the parameters for the first two superstable parameters eta_ss[0]\n",
" and eta_ss[1], finds the next nMax-1 of them up to eta_ss[nMax].\n",
" Returns dictionary eta_ss\n",
"\n",
" Usage:\n",
" Find the value of the parameter eta_ss[0] = eta0 for which the fixed\n",
" point is xMax and g is superstable (g(xMax) = xMax), and the\n",
" value eta_ss[1] = eta1 for which g(g(xMax)) = xMax, either\n",
" analytically or using FindSuperstable by hand.\n",
" mus = GetSuperstablePoints(f, 9, eta0, eta1)\n",
"\n",
" Searches for eta_ss[n] in the range (etaMin, etaMax),\n",
" with etaMin = eta_ss[n-1] + epsilon and\n",
" etaMax = eta_ss[n-1] + (eta_ss[n-1]-eta_ss[n-2])/A\n",
" where A=3 works fine for the maps I've tried.\n",
" (Asymptotically, A should be smaller than but comparable to delta.)\n",
" \"\"\"\n",
" eps = 1.0e-10\n",
" A = 3.\n",
" eta_ss = {}\n",
" eta_ss[0] = eta0\n",
" eta_ss[1] = eta1\n",
" leafSize = {}\n",
" leafSize[1] = ...\n",
" for n in arange(2, nMax+1):\n",
" etaMin = ...\n",
" etaMax = ...\n",
" nIterated = 2**n\n",
" eta_ss[n] = FindSuperstable(...)\n",
" leafSize[n] = Iterate(xMax, 2**(...), eta_ss[n], g)-xMax\n",
" return eta_ss, leafSize\n",
"\n",
"mu0 = FindSuperstable(f,1,0.3,1.0,0.5)\n",
"mu1 = FindSuperstable(f,2,mu0 + 1.e-6,1.0,0.5)\n",
"mus, leafSizes = GetSuperstablePointsAndIntervals(f,mu0,mu1)\n",
"\n",
"print(mus)\n",
"print(leafSizes)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*\n",
"Calculate the ratios $(\\mu^s_{n-1}-\\mu^s_{n-2})/(\\mu^s_n-\\mu^s_{n-1})$;\n",
"do they appear to converge to the Feigenbaum number\n",
"$\\delta = 4.6692016091029909\\dots$? Estimate $\\mu_\\infty$ by using your last two\n",
"values of $\\mu^s_n$, your last ratio estimate of $\\delta$,\n",
"and the equations \n",
"$\\mu_n \\approx \\mu_\\infty - A \\delta^{-n}$ and\n",
"$\\delta = \\lim_{n\\to\\infty} (\\mu_{n-1} - \\mu_{n-2})/(\\mu_n - \\mu_{n-1})$\n",
"above. *In the superstable orbit with $2^n$ points, the nearest point to $x=1/2$\n",
"is $f^{[2^{n-1}]}(1/2)$. (This is true because, at the previous superstable orbit,\n",
" $2^{n-1}$ iterates returned us to the original point $x=1/2$.)\n",
"*\n",
"Calculate the ratios of the amplitudes\n",
"$f^{[2^{n-1}]}(1/2)-1/2$ at successive\n",
"values of $n$; do they appear to converge to the universal value\n",
"$\\alpha = -2.50290787509589284\\dots$? \n",
"*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nMax = 11; # Note: Value of mu for n>9 not reliable\n",
"\n",
"def ExponentEstimates(g, eta_ss, leafSizes, xMax=0.5, nMax=nMax):\n",
" \"\"\"\n",
" Given superstable 2^n cycle values eta_ss[n], calculates\n",
" delta[n] = (eta_{n-1}-eta_{n-2})/(eta_{n}-eta_{n-1})\n",
" and alpha[n] = leafSizes[n-1]/leafSizes[n]\n",
"\n",
" Also extrapolates eta to etaInfinity using definition of delta and\n",
" most reliable value for delta:\n",
" delta = lim{n->infinity} (eta_{n-1}-eta_{n-2})/(eta_n - eta_{n-1})\n",
"\n",
" Returns delta and alpha dictionaries for n>=2, and etaInfinity\n",
" \"\"\"\n",
" delta = {}\n",
" alpha = {}\n",
" for n in range(2, nMax+1):\n",
" delta[n] = ...\n",
" alpha[n] = ...\n",
" etaInfinity = eta_ss[nMax] + ...\n",
" return delta, alpha, etaInfinity\n",
"\n",
"deltas, alphas, muInfinity = ExponentEstimates(f, mus, leafSizes, nMax=11)\n",
"\n",
"deltas[nMax], alphas[nMax], muInfinity\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*\n",
"Calculate the same ratios for\n",
"the map $f_2(x) = B \\sin(\\pi x)$; do $\\alpha$ and $\\delta$ appear\n",
"to be universal (independent of the mapping)?\n",
"*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"B0 = FindSuperstable(fSin,1,...)\n",
"B1 = FindSuperstable(fSin,2,B0+1.e-6,...)\n",
"Bs, leafSizesSin = GetSuperstablePointsAndIntervals(...)\n",
"deltaSin, alphaSin, BInfinity = ExponentEstimates(...)\n",
"\n",
"deltaSin[nMax], alphaSin[nMax], BInfinity"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The limits $\\alpha$ and $\\delta$ are independent of the map, so long as\n",
"it folds (one hump) with a quadratic maximum. They are the same, also,\n",
"for experimental systems with many degrees of freedom which undergo the\n",
"period-doubling cascade. This self-similarity and universality suggests\n",
"that we should look for a renormalization-group explanation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(e) Coarse-graining in time. \n",
"*\n",
"Plot $f(f(x))$ vs.\\ $x$ for $\\mu = 0.8$, together with the line $y=x$.\n",
"Notice that the period-two\n",
"cycle of $f$ becomes a pair of stable fixed-points for $f^{[2]}$.<\\em>\n",
"(We are coarse-graining in time---removing every other point in the time\n",
"series, by studying $f(f(x))$ rather than $f$.)\n",
"** Compare the plot with that for $f(x)$ vs.\\ $x$ for $\\mu = 0.5$. Notice\n",
"that the region zoomed in around $x=1/2$ for $f^{[2]}=f(f(x))$ looks quite\n",
"a bit like the entire map $f$ at the smaller value $\\mu=0.5$.\n",
"Plot $f^{[4]}(x)$ at $\\mu = 0.875$; notice again the small one-humped map\n",
"near $x=1/2$.\n",
"*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"x = arange(0.,1.,0.01)\n",
"plot(...))\n",
"plot(...)\n",
"title('mu='+str(0.8))\n",
"axes().set_aspect('equal')\n",
"\n",
"figure()\n",
"plot(x,f(x,0.5))\n",
"plot(x,x)\n",
"title('mu='+str(0.5))\n",
"axes().set_aspect('equal')\n",
"\n",
"figure()\n",
"plot(...)\n",
"plot(...)\n",
"title('mu='+str(0.875))\n",
"axes().set_aspect('equal')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The fact that the one-humped map reappears in smaller form just after\n",
"the period-doubling bifurcation is the basic reason that succeeding\n",
"bifurcations so often follow one another. The fact that many things are\n",
"universal is due to the fact that the little one-humped maps have a\n",
"shape which becomes *independent of the original map* after several\n",
"period-doublings.\n",
"\n",
"Let us define this renormalization-group transformation $T$, taking\n",
"function space into itself. Roughly speaking, $T$ will take the\n",
"small upside-down hump in $f(f(x)),$\n",
"invert it, and stretch it to cover the interval from $(0,1)$. Notice\n",
"in your graphs for part~(g) that the\n",
"line $y=x$ crosses the plot $f(f(x))$ not only at the two points on the\n",
"period-two attractor, but also (naturally) at the old fixed-point\n",
"$x^*[f]$\n",
"for $f(x)$. This unstable fixed-point plays the role for $f^{[2]}$ that\n",
"the origin played for $f$; our renormalization-group rescaling must map\n",
"$(x^*[f], f(x^*)) = (x^*,x^*)$ to the origin. The corner of the window\n",
"that maps to $(1,0)$ is conveniently located at $1-x^*$, since our map\n",
"happens to be symmetric about $x=1/2$. \n",
"(For asymmetric maps, we would need to locate this other corner\n",
" $f(f(x_c))=x^*$ numerically. As it happens, breaking this symmetry is\n",
" irrelevant at the fixed-point.)\n",
"For a general one-humped map $g(x)$ with fixed-point\n",
"$x^*[g]$ the side of the window is thus of length\n",
"$2(x^*[g]-1/2)$. To invert and stretch, we must thus rescale by a factor\n",
"$\\alpha[g] = -1/(2(x^*[g]-1/2))$. Our renormalization-group transformation\n",
"is thus a mapping $T[g]$ taking function space into itself, where\n",
"\\begin{equation}\n",
"T[g](x) = \\alpha[g]\n",
" \\left(g\\left(g(x/\\alpha[g]+x^*[g])\\right) - x^*[g]\\right).\n",
"\\end{equation}\n",
"(This is just rescaling $x$ to squeeze into the window, applying $g$ twice,\n",
"shifting the corner of the window to the origin, and then rescaling by\n",
"$\\alpha$ to fill the original range $(0,1) \\times (0,1)$.)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"(f) Scaling and the renormalization group:\n",
"*\n",
"Write routines that calculate $x^*[g]$ and $\\alpha[g]$, and define\n",
"the renormalization-group transformation $T[g]$. Plot $T[f]$,\n",
"$T[T[f]]$,... and compare them. Are we approaching a fixed-point $f^*$\n",
"in function space?\n",
"*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def XStar(g, mu):\n",
" \"\"\"\n",
" Finds fixed point of one-humped map g, which is assumed to be\n",
" between xMax and 1.0.\n",
" \"\"\"\n",
" def gXStar(x,mu): return ...\n",
" return brentq(gXStar, 1/2, 1.0, args=(mu,))\n",
"\n",
"def Alpha(g, mu):\n",
" \"\"\"\n",
" Finds the (negative) scale factor alpha which inverts and rescales\n",
" the small inverted region of g(g(x)) running from (1-x*) to x*.\n",
" \"\"\"\n",
" gWindowMax = XStar(g, mu)\n",
" gWindowMin = ...\n",
" return 1.0/(gWindowMin-gWindowMax)\n",
"\n",
"class T:\n",
" \"\"\"\n",
" Creates a new function T[g] from g, implementing Feigenbaum's\n",
" renormalization-group transformation of function space into itself.\n",
"\n",
" We define it as a class so that we can initialize alpha and xStar,\n",
" which otherwise would need to be recalculated each time T[g] was\n",
" evaluated at a point x.\n",
"\n",
" Usage:\n",
" Tg = T(g, args)\n",
" Tg(x) evaluates the function at x\n",
" \"\"\"\n",
" def __init__(self, g, mu):\n",
" \"\"\"\n",
" Stores g and args.\n",
" Calculates and stores xStar and alpha.\n",
" \"\"\"\n",
" self.mu = mu\n",
" self.xStar = XStar(g, mu)\n",
" self.alpha = Alpha(g, mu)\n",
" self.g = g\n",
"\n",
" def __call__(self, x, mu):\n",
" \"\"\"\n",
" Defines xShrunk to be x/alpha + x*\n",
" Evaluates g2 = g(g(xShrunk))\n",
" Returns expanded alpha*(g2-xStar)\n",
" \"\"\"\n",
" xShrunk = .../self.... + self....\n",
" g2 = self.g(self.g(...), ...)\n",
" return ... * (g2 - ...)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def PlotTIterates(g, mu, nMax = 2):\n",
" \"\"\"\n",
" Plots g(x), T[g](x), T[T[g]](x) ...\n",
" on the same plot, for x from (0,1)\n",
" \"\"\"\n",
" figure(figsize=(10,10))\n",
" xs = arange(0.01, 1., 0.01)\n",
" plot(xs,xs)\n",
" axes().set_aspect('equal')\n",
" Tg = {}\n",
" Tg[0] = lambda x, mu: g(x, mu)\n",
" for n in range(1,nMax+1):\n",
" Tg[n] = T(Tg[n-1], mu=mu)\n",
"\n",
" for n in range(nMax+1):\n",
" plot(xs, Tg[n](xs,mu))\n",
"\n",
"PlotTIterates(f,muInfinity,nMax=3)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"This explains the self-similarity; in particular, the value\n",
"of $\\alpha[g]$ as $g$ iterates to $f^*$ becomes the Feigenbaum number\n",
"$\\alpha = -2.5029\\dots$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(g) Universality and the renormalization group:\n",
"*\n",
"Using the sine function above, compare\n",
"$T[T[f_{\\mathrm{sin}}]]$ to $T[T[f]]$ at their onsets of chaos.\n",
"Are they approaching the same fixed-point?\n",
"*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"figure(figsize=(10,10))\n",
"xs = arange(0.001, 1., 0.01)\n",
"plot(xs,xs)\n",
"axes().set_aspect('equal')\n",
"T2F = T(T(f, mu=muInfinity),mu=muInfinity)\n",
"plot(xs,T2F(xs,muInfinity),\"r-\")\n",
"T2Fsin = T(T(fSin, mu=BInfinity),mu=muInfinity)\n",
"plot(xs,T2Fsin(xs,BInfinity),\"g-\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By using this rapid convergence in function space, one can prove both\n",
"that there will (often) be an infinite geometrical series of\n",
"period-doubling bifurcations leading to chaos, and that this series will share\n",
"universal features (exponents $\\alpha$ and $\\delta$ and features) that\n",
"are independent of the original dynamics."
]
}
],
"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
}