The onset of chaos: Full renormalization-group calculation

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

© 2017, James Sethna, all rights reserved.

In this exercise, we implement Feigenbaum's numerical scheme for finding high-precision values of the universal constants \begin{equation} \begin{aligned} \alpha &= -2.50290787509589282228390287322 \\ \delta &= 4.66920160910299067185320382158, \end{aligned} \end{equation} 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 \begin{equation} g(x) \approx 1 + \sum_{n=1}^N G_n x^{2 n}. \end{equation} 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: \begin{equation} T[g^*](x_m) = g^*(x_m), ~~~~~~~~ x_m = m/N,~m = \{0,\dots,N\}. \end{equation} 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 out your alphas
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 \begin{equation} \mathcal{L}[\psi](x) = \alpha \psi(g^*(x/\alpha)) + \alpha {g^*}'(g(x/\alpha) \psi(x/\alpha). \end{equation}

(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 \begin{equation} \psi(x) = \sum_{n=0}^{N-1} \psi_n x^{2n} ~~~~~~ (\psi_0 \equiv 1). \end{equation} 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. \begin{equation} \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 \end{equation} These equations for the coefficients $\psi_n$ of the eigenfunctions of $\mathcal{L}$ is in the form of a generalized eigenvalue problem \begin{equation} \sum_n L_{in} \psi_n = \lambda X_{in} \psi_n. \end{equation} 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])

# Print out your deltas
print(array([...

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