[seqfan] Re: Feigenbaum Function Coefficients

[Richard Mathar]

Followup on http://list.seqfan.eu/pipermail/seqfan/2010-August/005639.html :

There is an option of a brute-force standard approach that expands
everything into a finite Taylor series with some (preferentially large) upper
cutoff in the orders, nmax. If we insert
  a = 1+C[1]+C[2]+C[3]+.. +C[nmax]
  g(x) = 1+C[1]*x^2+C[2]*x^4+...+C[nmax]*x^(2*nmax)
into
a*g = g(g(a*x))
and expand both sides in series of x, then equalize the coefficients
in front of all the powers of x, we get nonlinear system of equations
where essentially some multinomials in the C(.) are set to constants.
[These can be written down analytically-ugly, because insertion of g(x)
into the g(g(a*x)) gives a binomal-type of double sum.]

I've tested this with nmax=4, and the output of Maple in this case
is the estimate:

{C[3] = -0.004610361600, C[4] = 0.00006788410799, C[5] = -0.4195271136 e-6  ,
    C[1] = -1.376969275, C[2] = 0.1244085255}

One would probably use some more dedicated vector-Newton approach
to solve the system of equations iteratively, because all the derivatives with
respect to the unknowns are "easy."


# Maple source code:
xlim := 5 ; # some cut-off value for the maximumpower in x
g := add(1+C[n]*x^(2*n),n=1..xlim) ;   # g=1+C(n)*x^(2n), finite ansatz
a := 1+add(C[n],n=1..xlim) ;  # a=1+sum C(n) as required by g(1)=a

gax := subs(x=a*x,g) ;   # expansion of of g(a*x)
gax := expand(gax) ;

ggax := subs(x=gax,g) ; # expansion of g(g(a*x))

# need to solve a*g=g(g(a*x)), which is a*g-ggax = 0
zer := a*g-ggax ;

# Set the coefficients of x^(2n),n=0,1,2,.. in a*g-g(g(a*x)) to zero,
# and build the set of the taylor coefficients for that purpose,
sset := { coeftayl(zer,x=0,0) } ;
for i from 2 to xlim do
	sset := sset union { coeftayl(zer,x=0,2*i)} ;
	print("tayl",i) ;
end do;
# attempt to solve numerically given estimates
fsolve(sset,{C[1]=-1.52763299703630145,C[2]=+0.104815194787,C[3]=+0.0267056705,C[4]=-0.0035274096,C[5]=+0.00008160096654753174}) ;