The onset of chaos: Full renormalization-group calculation

(Sethna, “Entropy, Order Parameters, and Complexity”, ex. 12.XXX)
© 2017, James P. Sethna, all rights reserved.

In this exercise, we implement Feigenbaum’s numerical scheme for finding high-precision values of the universal constants
       α = −2.50290787509589282228390287322
       δ = 4.66920160910299067185320382158,
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’}.

αFeigenbaum = -2.502907875095892822283902873218;
δFeigenbaum = 4.669201609102990671853203821578;

Our renormalization group operation (Exercises ‘Period doubling and the renormalization group’ and the companion Exercise) coarse-grains in time taking g→g∘g, and then rescales distance x by a factor of α. Centering our functions at x=0, this leads to 
            T[g](x) = α g(g(x/α)).
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) ≈ 1 + (Subscript[(Σ^N), n=1])(Subscript[G, n])(x^(2n)).
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 δ will also be symmetric.

(* N us a reserved variable; use Nn instead *)
g[Nn_][x_] := 1+Sum[...,{n,1,Nn}]
T[g_][x_] := α ...
(* We'll also want the derivative of g later *)
Dg[Nn_][x_] := Sum[...,{...}]
(* Test your functions by plotting them.G=[-1.5,0,0,...] should give T[g] close to g *)
G0[1] := -3/2
G0[n_] := 0 /; n ≠ 1
Plot[{g[2][x],T[g[2]][x]/.{α->1/g[2][1]}}/.G->G0,{x,0,2}]

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

(a) Show that the equation for m=0 sets α = 1/g*(1).

We can use a root-finding routine to solve for (Subscript[(SuperStar[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 α from part (a). Check that your routine at N=1 gives values for α≈−2.5 and G∗1≈−1.5. (These should reproduce the values from the companion Exercise part (c).)

Nmax=20;
For[Nn=1,Nn<=Nmax, Nn=Nn+1,
xm = Range[...];
vars = Table[{G[n],G0[n]},{...}];eqns = Table[...==.../. α->1/g[Nn][1],{x,xm}];GStar[Nn] = FindRoot[..., WorkingPrecision->50];
α[Nn] = ...]/.GStar[Nn]]
Table[{Nn, ...},{Nn,1,Nmax}]//MatrixForm
Table[{Nn, αFeigenbaum - ...

Now we need to solve for the renormalization group flows T[g], linearized about the fixed point g(x)=g*(x)+ϵψ(x). Feigenbaum tells us that T[g*+ϵψ]=T[g∗]+ϵℒ[ψ], where ℒ  is the linear operator taking ψ(x) to
                 ℒ[ψ](x) = α ψ( g*(x/α)) + α g*’(g(x/α)) ψ(x/α).
                 
(d) Derive the equation above.

ANSWER HERE

We want to find eigenfunctions that satisfy ℒ[ψ]=λψ. Again, we can expand ψ(x) in a polynomial
                        ψ(x)=(Subscript[(∑^(N-1)), n=0]) (Subscript[ψ, n])(x^(2n))      ((Subscript[ψ, 0]) ≡ 1).
We then approximate the action of ℒ on ψ by its action at N points (Subscript[x̃, i]), that need not be the same as the N points (Subscript[x, m]) we used to find g∗. We shall use  (Subscript[x̃, i]) = (i−1)/(N−1), i=1,…,N. (For N=1, we use (Subscript[x̃, 1])=0.) This leads us to a linear system of N equations for the coefficients (Subscript[ψ, n]), using the previous two equations.
(Subscript[(∑^(N-1)), n=0]) [α g((((Subscript[x̃, i])/α))^(2n))+α g’(g((Subscript[x̃, i])/α)) ((((Subscript[x̃, i])/α))^(2n))] (Subscript[ψ, n]) = λ (Subscript[(∑^(N-1)), n=0]) (Subscript[((x̃)^(2n)), i]) (Subscript[ψ, n])
These equations for the coefficients (Subscript[ψ, n]) of the eigenfunctions of  is in the form of a generalized eigenvalue problem
(Subscript[(∑^(N-1)), n=0]) (Subscript[L, in]) (Subscript[ψ, n]) = λ (Subscript[(∑^(N-1)), n=0]) (Subscript[X, in ])(Subscript[ψ, 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×1 matrices. For N=1, what is your prediction for δ? (These should reproduce the values from the companion Exercise part (d).)

(* Make sure your matrix hasn't transposed rows (i) and columns (n).Each row should give powers of one x_i. *)
(* Avoid 0^0 for N=1 by using 'Evaluate' to set up n=0 column in X *)
xtildes[Nn_] := Range[0,1,1/(Nn-1)]
xtildes[1] := {0.}
X[Nn_] := Table[Evaluate[Table[...,{n,0,Nn-1}]], {xtilde,xtildes[Nn]}]
X[1]
X[3]

L[Nn_] :=Table[Evaluate[Table[α[Nn] g[...]^(...)+α[Nn] Dg[...] (...)^(...) /.GStar[Nn],{n,0,Nn-1}]],{xtilde,xtildes[Nn]}]
L[1]
L[3]//MatrixForm
Eigenvalues[{L[3],X[3]}]

Nmax=20;
For[Nn=1,Nn<=Nmax, Nn=Nn+1,
eigvals = ...[{L[Nn],X[Nn]}];
δ[Nn] = eigvals[[1]]]
Table[...]// ...
Table[...]// ...