{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The onset of chaos: Full renormalization-group calculation

\n",
"\n",
"(Sethna, \"Entropy, Order Parameters, and Complexity\", ex. 12.XXX)\n",
"\n",
"© 2017, James Sethna, all rights reserved."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this exercise, we implement Feigenbaum's numerical\n",
"scheme for finding high-precision values of the\n",
"universal constants\n",
"\\begin{equation}\n",
"\\begin{aligned}\n",
"\\alpha &= -2.50290787509589282228390287322 \\\\\n",
"\\delta &= 4.66920160910299067185320382158,\n",
"\\end{aligned}\n",
"\\end{equation}\n",
"that quantify the scaling properties of the period-doubling route to chaos\n",
"(Fig. 12.17, Exercise 'Period doubling').\n",
"This extends the lowest-order calculation of\n",
"the companion Exercise 'The onset of chaos: Lowest order renormalization-group for\n",
"period doubling'}.\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 root\n",
"from scipy.linalg import eig\n",
"\n",
"alphaFeigenbaum = -2.502907875095892822283902873218\n",
"deltaFeigenbaum = 4.669201609102990671853203821578"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our renormalization group operation (Exercises 'Period doubling and the\n",
"renormalization group'\n",
"and the companion Exercise) coarse-grains in time\n",
"taking $g \\to g\\circ g$, and then rescales distance $x$ by a factor of $\\alpha$.\n",
"Centering our functions at $x=0$, this leads to\n",
"$T[g](x) = \\alpha g\\left(g(x/\\alpha)\\right)$.\n",
"\n",
"We shall solve for the properties at the onset of chaos by\n",
"analyzing our function-space renormalization-group by expanding our functions\n",
"in a power series\n",
"\\begin{equation}\n",
"g(x) \\approx 1 + \\sum_{n=1}^N G_n x^{2 n}.\n",
"\\end{equation}\n",
"Notice that we only keep even powers of $x$; the fixed point is known to\n",
"be symmetric about the maximum, and the unstable mode responsible\n",
"for the exponent $\\delta$ will also be symmetric."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def g(G,x):\n",
" \"\"\"\n",
" Returns 1 + G[0] x^2 + G[1] x^4 + ..., where G_n = G[n-1]\n",
" We will sometimes call g with a whole array of x-values.\n",
" \"\"\"\n",
" # enumerate(G) = [[0,G[0]], [1,G[1]], ...], conveniently giving n-1 and Gn in the formula.\n",
" # enumerate(G,1) starts the numbering at one\n",
" # sum(M) adds up all the entries of a matrix. This is OK if x is a scalar, but if we send in a whole\n",
" # array [x1,x2,...] we want an array of values [g(x1),g(x2),...]. sum(M,axis=0) sums up the rows of the matrix.\n",
" return 1.+sum([... for n,Gn in enumerate(G,1)],axis=0)\n",
" \n",
"def T(g,G,x,alpha=None):\n",
" \"\"\"\n",
" Returns renormalization-group transform T[g](x).\n",
" If alpha is not known, calculate it from g using your result from (a) below.\n",
" \"\"\"\n",
" if alpha is None:\n",
" alpha = ...\n",
" return ...\n",
"\n",
"def Dg(G,x):\n",
" \"\"\"\n",
" Returns g'(x)\n",
" \"\"\"\n",
" return sum(...)\n",
"\n",
"# Test your functions by plotting them. G = [-1.5, 0, 0, ...] should give T[g] close to g for x<1\n",
"x = arange(0,3,0.01)\n",
"plot(x,g([-1.5,0.],x))\n",
"plot(x,T(g,[-1.5,0],x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we must approximate the fixed point $g^*(x)$ and the corresponding\n",
"value of the universal constant $\\alpha$. At order $N$, we must solve\n",
"for $\\alpha$ and the $N$ polynomial coefficients $G^*_n$. We can use\n",
"the $N+1$ equations fixing the function at equally spaced points in\n",
"the positive unit interval:\n",
"\\begin{equation}\n",
"T[g^*](x_m) = g^*(x_m), ~~~~~~~~ x_m = m/N,~m = \\{0,\\dots,N\\}.\n",
"\\end{equation}\n",
"We can use the first of these equations to solve for $\\alpha$.\n",
"\n",
"(a) *Show that the equation for $m=0$ sets $\\alpha = 1/g^*(1)$.*\n",
"\n",
"We can use a root-finding routine to solve for $G_n^*$.\n",
"\n",
"(b) * Implement the other $N$ constraint equations\n",
"above in a form appropriate for your method of finding\n",
"roots of nonlinear equations, substituting your value for $\\alpha$\n",
"from part (a). Check that your routine at $N=1$ gives\n",
"values for $\\alpha\\approx -2.5$ and $G^*_1 \\approx -1.5$.*\n",
"(These should reproduce the values from the companion Exercise part (c).)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def toZero(G):\n",
" \"\"\"Returns T[g](x) - g(x) for N points [1/N,2/N,...,1], given N terms in its power series Gn\"\"\"\n",
" N = len(G)\n",
" x = linspace(...)\n",
" return ...\n",
"\n",
"# Check that your return gives a sensible value for the difference of Tg and g at x=1, for G1 = -1.5\n",
"print(toZero([-1.5]))\n",
"\n",
"# Use root to find the best solution for N=1. The values giving zero is returned as root(toZero,[initial values]).x\n",
"G1 = root(...,[-1.5]).x\n",
"\n",
"# What do we get for alpha[1]?\n",
"1/..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(c) * Use a root-finding routine to calculate $\\alpha$ for $N=1, \\dots, 9$.\n",
"Start the search at $G^*_1 = -1.5$, $G^*_n = 0$ ($n>1$) to avoid\n",
"landing at the wrong fixed point.*\n",
"(If it is convenient for you to use high-precision arithmetic,\n",
"continue to higher $N$.) * To how many decimal places can you reproduce\n",
"the correct value for $\\alpha$ at the beginning of this exercise?*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Fill dictionary with your values of alpha[N] for N = 1...9\n",
"# Also keep your values for the fixed point function G[N] for use in calculating delta\n",
"alpha = {}\n",
"G = ...\n",
"Nmax = 9\n",
"for N in range(1,Nmax):\n",
" G0 = zeros(N)\n",
" G0[0]=-1.5\n",
" G[N] = root(...\n",
" alpha[N] = ...\n",
"\n",
"# Print out your alphas\n",
"print(array([(N,...) for N in range(1,Nmax)]))\n",
"\n",
"# Calculate how far they deviate from alphaFeigenbaum\n",
"[(N,alphaFeigenbaum-...\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we need to solve for the renormalization group flows $T[g]$, linearized\n",
"about the fixed point $g(x)=g^*(x)+\\epsilon \\psi(x)$. Feigenbaum tells\n",
"us that $T[g^*+\\epsilon \\psi] = T[g^*] + \\epsilon \\mathcal{L}[\\psi]$, where\n",
"$\\mathcal{L}$ is the linear operator taking $\\psi(x)$ to\n",
"\\begin{equation}\n",
"\\mathcal{L}[\\psi](x) = \\alpha \\psi(g^*(x/\\alpha)) + \\alpha {g^*}'(g(x/\\alpha) \\psi(x/\\alpha).\n",
"\\end{equation}\n",
"\n",
"(d) * Derive the equation above.*"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We want to find eigenfunctions that satisfy $\\mathcal{L}[\\psi] = \\lambda \\psi$.\n",
"Again, we can expand $\\psi(x)$ in a polynomial\n",
"\\begin{equation}\n",
"\\psi(x) = \\sum_{n=0}^{N-1} \\psi_n x^{2n} ~~~~~~ (\\psi_0 \\equiv 1).\n",
"\\end{equation}\n",
"We then approximate the action of $\\mathcal{L}$ on $\\psi$ by its action\n",
"at $N$ points $\\tilde{x}_i$, that need not be the same as the $N$\n",
"points $x_m$ we used to find $g^*$. We shall use $\\tilde{x}_i = (i-1)/(N-1)$,\n",
"$i=1,\\dots,N$. (For $N=1$, we use $\\tilde{x}_1 = 0$.)\n",
"This leads us to a linear system of\n",
"$N$ equations for the coefficients $\\psi_n$, using the previous two\n",
"equations.\n",
"\\begin{equation}\n",
"\\sum_{n=0}^{N-1} \\left[\\alpha g(\\tilde{x}_i/\\alpha)^{2 n} + \\alpha g'(g(\\tilde{x}_i/\\alpha)) (\\tilde{x}_i/\\alpha)^{2n}\\right]\\psi_n = \\lambda \\sum_{n=0}^{N-1} \\tilde{x}_i^{2 n} \\psi_n\n",
"\\end{equation}\n",
"These equations for the coefficients $\\psi_n$ of the eigenfunctions\n",
"of $\\mathcal{L}$ is in the form of a * generalized eigenvalue problem*\n",
"\\begin{equation}\n",
"\\sum_n L_{in} \\psi_n = \\lambda X_{in} \\psi_n.\n",
"\\end{equation}\n",
"The solution to the generalized eigenvalue problem can be found from the\n",
"eigenvalues of $X^{-1} L$, but most eigenvalue routines provide a more\n",
"efficient and accurate option for directly solving the generalized equation\n",
"given $L$ and $X$.\n",
"\n",
"(e) * Write a routine that calculates the matrices $L$ and $X$\n",
"implicitly defined by the previous two equations.\n",
"For $N=1$ you should generate $1\\times1$ matrices. For $N=1$, what is\n",
"your prediction for $\\delta$?*\n",
"(These should reproduce the values from\n",
"the companion Exercise part (d).)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def X(N):\n",
" \"\"\"Returns X_{in} = xtilde_i**(2n)\"\"\"\n",
" # Make sure your matrix hasn't transposed rows (i) and columns (n). Each row should have \n",
" xtildes = linspace(0.,1.,N)\n",
" return array([[... for n in range(N)] for xtilde in xtildes])\n",
"\n",
"print(X(1))\n",
"print(X(3))\n",
"\n",
"def Ln(xtildes,n,alpha,G):\n",
" \"\"\"Returns one column of L, given the array of xtilde values\"\"\"\n",
" return alpha*g(...)**(...)+alpha*Dg(...)*(...)**(...)\n",
"\n",
"# Test Ln on the one-element column for N=1: does it give a reasonable value for delta?\n",
"print('delta[1] should be the entry in ', Ln(array([0.]),0,alpha[1],G[1]))\n",
"\n",
"def L(N):\n",
" \"\"\"Builds an array Lin from the columns Ln\"\"\"\n",
" # Again, make sure your matrix has rows (i) and columns (n). You may need to use M.transpose.\n",
" xtildes = ...\n",
" return array([Ln(...) for n in range(N)]).transpose()\n",
"\n",
"print(L(1))\n",
"\n",
"print(L(3))\n",
"\n",
"eig(L(3),X(3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(f) * Solve the generalized eigenvalue problem for $L$ and $X$\n",
"for $N=1,\\dots,9$. To how many decimal places can you reproduce\n",
"the correct value for $\\delta$ at the beginning of this exercise?*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Fill dictionary with your values of alpha[N] for N = 1...9\n",
"delta = {}\n",
"Nmax = 9\n",
"for N in range(1,Nmax):\n",
" eigvals, eigvecs = ...\n",
" delta[N] = real(eigvals[0])\n",
"\n",
"# Print out your deltas\n",
"print(array([...\n",
"\n",
"# Calculate how far they deviate from deltaFeigenbaum\n",
"[(N,..."
]
}
],
"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
}