# The onset of chaos: Full renormalization-group calculation

(Sethna, "Entropy, Order Parameters, and Complexity", ex. 12.XXX)

In this exercise, we implement Feigenbaum's numerical scheme for finding high-precision values of the universal constants \begin{aligned} \alpha &= -2.50290787509589282228390287322 \\ \delta &= 4.66920160910299067185320382158, \end{aligned} that quantify the scaling properties of the period-doubling route to chaos (Fig. 12.17, Exercise 'Period doubling'). This extends the lowest-order calculation of the companion Exercise 'The onset of chaos: Lowest order renormalization-group for period doubling'}.

Import packages

In [ ]:
%pylab inline
from scipy import *
from scipy.optimize import root
from scipy.linalg import eig

alphaFeigenbaum = -2.502907875095892822283902873218
deltaFeigenbaum = 4.669201609102990671853203821578


Our renormalization group operation (Exercises 'Period doubling and the renormalization group' and the companion Exercise) coarse-grains in time taking $g \to g\circ g$, and then rescales distance $x$ by a factor of $\alpha$. Centering our functions at $x=0$, this leads to $T[g](x) = \alpha g\left(g(x/\alpha)\right)$.

We shall solve for the properties at the onset of chaos by analyzing our function-space renormalization-group by expanding our functions in a power series $$g(x) \approx 1 + \sum_{n=1}^N G_n x^{2 n}.$$ Notice that we only keep even powers of $x$; the fixed point is known to be symmetric about the maximum, and the unstable mode responsible for the exponent $\delta$ will also be symmetric.

In [ ]:
def g(G,x):
"""
Returns 1 + G[0] x^2 + G[1] x^4 + ..., where G_n = G[n-1]
We will sometimes call g with a whole array of x-values.
"""
# enumerate(G) = [[0,G[0]], [1,G[1]], ...], conveniently giving n-1 and Gn in the formula.
# enumerate(G,1) starts the numbering at one
# 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
# array [x1,x2,...] we want an array of values [g(x1),g(x2),...]. sum(M,axis=0) sums up the rows of the matrix.
return 1.+sum([... for n,Gn in enumerate(G,1)],axis=0)

def T(g,G,x,alpha=None):
"""
Returns renormalization-group transform T[g](x).
If alpha is not known, calculate it from g using your result from (a) below.
"""
if alpha is None:
alpha = ...
return ...

def Dg(G,x):
"""
Returns g'(x)
"""
return sum(...)

# Test your functions by plotting them. G = [-1.5, 0, 0, ...] should give T[g] close to g for x<1
x = arange(0,3,0.01)
plot(x,g([-1.5,0.],x))
plot(x,T(g,[-1.5,0],x))


First, we must approximate the fixed point $g^*(x)$ and the corresponding value of the universal constant $\alpha$. At order $N$, we must solve for $\alpha$ and the $N$ polynomial coefficients $G^*_n$. We can use the $N+1$ equations fixing the function at equally spaced points in the positive unit interval: $$T[g^*](x_m) = g^*(x_m), ~~~~~~~~ x_m = m/N,~m = \{0,\dots,N\}.$$ We can use the first of these equations to solve for $\alpha$.

(a) Show that the equation for $m=0$ sets $\alpha = 1/g^*(1)$.

We can use a root-finding routine to solve for $G_n^*$.

(b) Implement the other $N$ constraint equations above in a form appropriate for your method of finding roots of nonlinear equations, substituting your value for $\alpha$ from part (a). Check that your routine at $N=1$ gives values for $\alpha\approx -2.5$ and $G^*_1 \approx -1.5$. (These should reproduce the values from the companion Exercise part (c).)

In [ ]:
def toZero(G):
"""Returns T[g](x) - g(x) for N points [1/N,2/N,...,1], given N terms in its power series Gn"""
N = len(G)
x = linspace(...)
return ...

# Check that your return gives a sensible value for the difference of Tg and g at x=1, for G1 = -1.5
print(toZero([-1.5]))

# Use root to find the best solution for N=1. The values giving zero is returned as root(toZero,[initial values]).x
G1 = root(...,[-1.5]).x

# What do we get for alpha[1]?
1/...


(c) Use a root-finding routine to calculate $\alpha$ for $N=1, \dots, 9$. Start the search at $G^*_1 = -1.5$, $G^*_n = 0$ ($n>1$) to avoid landing at the wrong fixed point. (If it is convenient for you to use high-precision arithmetic, continue to higher $N$.) To how many decimal places can you reproduce the correct value for $\alpha$ at the beginning of this exercise?

In [ ]:
# Fill dictionary with your values of alpha[N] for N = 1...9
# Also keep your values for the fixed point function G[N] for use in calculating delta
alpha = {}
G = ...
Nmax = 9
for N in range(1,Nmax):
G0 = zeros(N)
G0[0]=-1.5
G[N] = root(...
alpha[N] = ...

print(array([(N,...) for N in range(1,Nmax)]))

# Calculate how far they deviate from alphaFeigenbaum
[(N,alphaFeigenbaum-...


Now we need to solve for the renormalization group flows $T[g]$, linearized about the fixed point $g(x)=g^*(x)+\epsilon \psi(x)$. Feigenbaum tells us that $T[g^*+\epsilon \psi] = T[g^*] + \epsilon \mathcal{L}[\psi]$, where $\mathcal{L}$ is the linear operator taking $\psi(x)$ to $$\mathcal{L}[\psi](x) = \alpha \psi(g^*(x/\alpha)) + \alpha {g^*}'(g(x/\alpha) \psi(x/\alpha).$$

(d) Derive the equation above.

We want to find eigenfunctions that satisfy $\mathcal{L}[\psi] = \lambda \psi$. Again, we can expand $\psi(x)$ in a polynomial $$\psi(x) = \sum_{n=0}^{N-1} \psi_n x^{2n} ~~~~~~ (\psi_0 \equiv 1).$$ We then approximate the action of $\mathcal{L}$ on $\psi$ by its action at $N$ points $\tilde{x}_i$, that need not be the same as the $N$ points $x_m$ we used to find $g^*$. We shall use $\tilde{x}_i = (i-1)/(N-1)$, $i=1,\dots,N$. (For $N=1$, we use $\tilde{x}_1 = 0$.) This leads us to a linear system of $N$ equations for the coefficients $\psi_n$, using the previous two equations. $$\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$$ These equations for the coefficients $\psi_n$ of the eigenfunctions of $\mathcal{L}$ is in the form of a generalized eigenvalue problem $$\sum_n L_{in} \psi_n = \lambda X_{in} \psi_n.$$ The solution to the generalized eigenvalue problem can be found from the eigenvalues of $X^{-1} L$, but most eigenvalue routines provide a more efficient and accurate option for directly solving the generalized equation given $L$ and $X$.

(e) Write a routine that calculates the matrices $L$ and $X$ implicitly defined by the previous two equations. For $N=1$ you should generate $1\times1$ matrices. For $N=1$, what is your prediction for $\delta$? (These should reproduce the values from the companion Exercise part (d).)

In [ ]:
def X(N):
"""Returns X_{in} = xtilde_i**(2n)"""
# Make sure your matrix hasn't transposed rows (i) and columns (n). Each row should have
xtildes = linspace(0.,1.,N)
return array([[... for n in range(N)] for xtilde in xtildes])

print(X(1))
print(X(3))

def Ln(xtildes,n,alpha,G):
"""Returns one column of L, given the array of xtilde values"""
return alpha*g(...)**(...)+alpha*Dg(...)*(...)**(...)

# Test Ln on the one-element column for N=1: does it give a reasonable value for delta?
print('delta[1] should be the entry in ', Ln(array([0.]),0,alpha[1],G[1]))

def L(N):
"""Builds an array Lin from the columns Ln"""
# Again, make sure your matrix has rows (i) and columns (n). You may need to use M.transpose.
xtildes = ...
return array([Ln(...) for n in range(N)]).transpose()

print(L(1))

print(L(3))

eig(L(3),X(3))


(f) Solve the generalized eigenvalue problem for $L$ and $X$ for $N=1,\dots,9$. To how many decimal places can you reproduce the correct value for $\delta$ at the beginning of this exercise?

In [ ]:
# Fill dictionary with your values of alpha[N] for N = 1...9
delta = {}
Nmax = 9
for N in range(1,Nmax):
eigvals, eigvecs = ...
delta[N] = real(eigvals[0])