# Calculation of the Feigenbaum Constants¶

by Harald Hofstätter, 25.10.2015

In this Julia notebook we calculate the Feigenbaum constants $\alpha=-2.502907875\dots$ and $\delta=4.669201609\dots$ to high precision. The notebook can be downloaded here.

We use the method given in

Keith Briggs, A precise calculation of the Feigenbaum constants, Math. Comp. 57 (1991), 435-439.

## Subroutines needed for the calculation of $\alpha$¶

The following subroutine evaluates $$g(x)=1+\sum_{j=1}^{n}g_jx^{2j}$$ for given $x$ and $g_j$ using Horner's method:

In :
function eval_g(x, g::Vector, n::Int)
x2 = x^2
y = g[n]
for j = n-1:-1:1
y = y*x2 + g[j]
end
y = y*x2 + 1
end

Out:
eval_g (generic function with 1 method)


The following subroutine evaluates the derivative $$g'(x)=\sum_{j=1}^{n}2jg_jx^{2j-1}$$ for given $x$ and $g_j$ using Horner's method:

In :
function eval_diff_g(x, g::Vector, n::Int)
x2 = x^2
y = 2*n*g[n]
for j = n-1:-1:1
y = y*x2 + 2*j*g[j]
end
y = y*x
end

Out:
eval_diff_g (generic function with 1 method)


For a given set of collocation points $\{x_1,\dots,x_n\}\subset(0,1]$, the following system of $n$ nonlinear equations has to be solved for the coefficients $g_1,\dots,g_n$ of $g(x)=1+\sum_{j=1}^{n}g_jx^{2j}$:

$$F_i(g_1,\dots,g_n)=g(1)g(x_i)-g(g(g(1)g(x_i)))=0,\quad i=1,\dots,n$$

For Newton's method we also need the Jacobian $J$ of $F$:

$$J_{ij}=\frac{\partial F_i}{\partial g_j} = g(x_i)+g(1)x_i^{2j}-g(g(1)x_i)^{2j}-g'(g(g(1)x_i))\cdot\big(g(1)^{2j}x_i^{2j}+g'(g(1)x_i)x_i\big)$$

The following subroutine computes $F$ and $J$ for given $g_j$ and $x_i$:

In :
function generate_F_J(F::Vector, J::Matrix, g::Vector, x::Vector, n::Int)
g1 = eval_g(1, g, n)
gx = [eval_g(x[j], g, n) for j=1:n]
gdax = [eval_diff_g(g1*x[j], g, n) for j=1:n]
gax = [eval_g(g1*x[j], g, n) for j=1:n]
gdgax = [eval_diff_g(gax[j], g, n) for j=1:n]
ggax = [eval_g(gax[j], g, n) for j=1:n]
g12 =  g1*g1
x2 = x.*x
gax2 = gax.*gax
F[:] = g1*gx - ggax
for j=1:n
if j==1
x2j = x2
gax2j = gax2
g12j = g12
else
x2j .*= x2
gax2j .*= gax2
g12j .*= g12
end
J[:,j] =  gx + g1*x2j  - gax2j - gdgax.*(g12j*x2j + gdax.*x)
end
end

Out:
generate_F_J (generic function with 1 method)


At each step of the Newton iteration the approximation for the Feigenbaum constant $\alpha=1/g(1)=1/(1+\sum_{j=1}^{n}g_j)$ together with the residuum $\max_i|F_i(g_1,\dots,g_n)|$ for the current values of $g_j$ is printed. The iteration is terminated if the residuum is less than 100 times the machine epsilon. Some heuristics is used to determine the number of valid digits of $\alpha$.

In :
function newton_iteration_for_alpha_g(g::Vector, x::Vector, n::Int, maxit::Int)
T=typeof(g)
J = zeros(T,n,n)
F = zeros(T,n)
alpha = 1/eval_g(1, g, n)
alphaold = alpha
diff = 1000
diffold = diff
for k=1:maxit
generate_F_J(F, J, g, x, n)
g[:] = g - J\F # Newton iteration

res = maximum(abs(F))
alpha = 1/eval_g(1, g, n)
diff = abs(alpha-alphaold)
println(alpha ," ",@sprintf(" %.3e %.3e",res, diff))
alphaold = alpha
if res<100*eps(T) || diff>.1*diffold
break
end
diffold = diff
end
digits = int(floor(-exponent(diff)*log(2)/log(10)))
println(digits, " digits seem to be correct:")
println("α=", string(alpha)[1:digits+3])
return alpha
end

Out:
newton_iteration_for_alpha_g (generic function with 1 method)


## Calculation of $\alpha$¶

We use $x_i=\sqrt{i/n},\ i=1,\dots,n$ as collocation points and (the coefficients of)

$$g(x) = 1-1.5218x^2+1.0481\cdot 10^{-1}x^4+2.67076\cdot 10^{-2}x^6$$

as the starting approximation for Newton's method:

In :
MyType=BigFloat
set_bigfloat_precision(1024)

n=int(get_bigfloat_precision()/5) # For this relation between bigfloat_presision and n the heuristsics
# in newton_iteration_for_alpha_g and power_method_for_delta_h seems to work.

maxit=30
x = sqrt([1:n]/MyType(n)) # collocation points
g = zeros(x)
g[1:3]=[-1.5218 ,1.0481e-1, 2.67076e-2]

alpha=newton_iteration_for_alpha_g(g, x, n, maxit)

-2.5015096545702020561742105864983246542905164138125082092391139293372914859782429002699498562918184318981302067343323080810359047017377660511144769781504542448925188682833665990492904180095959297629208801558365895549853717453608157195488431596489139687669675494593778695681446919234431339766507620962550146109e+00  1.371e-02 6.074e-02
-2.502907674276637625370806463935613885923965033174913673349829908870565207392040241039591520918831200812987370848644319932734152661936820733621440476541789526618308158990916009826611002404040990445080163003096986071421830718083298085967666762191066544182378350897558336345390078655614643661157376708162433625879e+00  3.152e-04 1.398e-03
-2.502907875096094507779349598398278109139444216631215934021496782714341407374084648333210096659233790123641165774080595842944769802073162495889219510276981245165604935716997191924073933032308683178065751919099126022670858127495880203539476892937951440632235140278167052639586539433123990589383438303035902607513e+00  9.603e-08 2.008e-07
-2.502907875095892822283902868205402414393679862622757902087614696561736549360141438026784433433440487350186191371416901791470946233713876836626529997679319784011900519228247531416627089195433191863970851823854631761783969681296587120193458324214195045824844270217072612787669350184428336296374479630788196727543e+00  2.583e-14 2.017e-13
-2.502907875095892822283902873218215786381271376727149964661667669163516916228267066488555646151880570830955765812827439503532072037089031621686548701795521162182110591761167824495632872802636079924887489016196618749344169278952443359614173791721444192106172739212908943510995730784144983313348031055252894752676e+00  1.124e-27 5.013e-27
-2.502907875095892822283902873218215786381271376727149977336192056779235463179590206703299649746433834129595146855023638090551139244064701863970547273246738912499015799001377346534110888503529098315101602476839442748938016499162981483120839554960027762376882297757882196098169717171900087573140985760925239711159e+00  3.081e-54 1.267e-53
-2.502907875095892822283902873218215786381271376727149977336192056779235463179590206703299649746433834129595231869995854723942182377785445179272863314993372578112163594879503744781260997380598671239711737328927665400522774408407521876772998481324024579756872876500478020026313846522263830252789428794900776471939e+00  1.815e-107 8.501e-107
-2.502907875095892822283902873218215786381271376727149977336192056779235463179590206703299649746433834129595231869995854723942182377785445179272863314993372578112163594879503744781260997380598671239711737328927665404401030669831383460009413932236448681146210744161479082800528950863341960371653597073232137299151e+00  8.093e-214 3.878e-213
-2.502907875095892822283902873218215786381271376727149977336192056779235463179590206703299649746433834129595231869995854723942182377785445179272863314993372578112163594879503744781260997380598671239711737328927665404401030669831383460009413932236448423818861516107626403586585775897765197162388211662348336256556e+00  4.024e-308 2.573e-247
246 digits seem to be correct:
α=-2.502907875095892822283902873218215786381271376727149977336192056779235463179590206703299649746433834129595231869995854723942182377785445179272863314993372578112163594879503744781260997380598671239711737328927665404401030669831383460009413932236448


Out:
-2.502907875095892822283902873218215786381271376727149977336192056779235463179590206703299649746433834129595231869995854723942182377785445179272863314993372578112163594879503744781260997380598671239711737328927665404401030669831383460009413932236448423818861516107626403586585775897765197162388211662348336256556e+00 with 1024 bits of precision


Here the value from http://www.plouffe.fr/simon/constants/feigenbaum.txt for comparision:

-2.50290787509589282228390287321821578638127137672714997733619205677923546317959020670329964974643383412959523186999585472394218237 77854451792728633149933725781121635948795037447812609973805986712397117373289276654044010306698313834600094139322364490657889951220 58431725078733774630878534242853519885875000423582469187408204281700901714823051821621619413199856066129382742649709844084470100805 45496779367608881264464068851815527093240075425064971570470475419932831783645332562415378693957125097066387979492654623137

## Subroutines needed for the calculation of $\delta$¶

The Feigenbaum constant $\delta=4.6692016091\dots$ is calculated as the largest eigenvalue of $V^{-T}L$, where $V$ is the Vandermonde matrix with respect to $x_1^2,\dots,x_n^2$ (the squares of the collocationg points from above), and $L$ is defined as

$$L_{ij}=\alpha\big(g(x_i/\alpha))^{2(j-1)}+g'(g(x_i/\alpha))(x_i/\alpha)^{2(j-1)}\big),$$

with the polynomial $g(x)=1+\sum_{j=1}^{n}g_jx^{2j}$ calculated above and $\alpha=1/g(1)$.

The following subroutine calculates the matrix $L$ for given $g_j$ and $x_i$:

In :
function generate_L(L::Matrix, g::Vector , x::Vector, n::Int)
g1 = eval_g(1, g, n)
ax = g1*x
gax = [eval_g(ax[k], g, n) for k=1:n]
gdgax = [eval_diff_g(gax[k], g, n) for k=1:n]
ax2 = ax.^2
gax2 = gax.^2
alpha = 1/g1
for j=1:n
if j==1
ax2j =  1
gax2j = 1
elseif j==2
ax2j = ax2
gax2j = gax2
else
ax2j .*= ax2
gax2j .*= gax2
end
L[:,j] = alpha*(gax2j + gdgax.*ax2j)
end
end

Out:
generate_L (generic function with 1 method)


For the solution of the Vandermonde system $V^T a=f$ we use Algorithm 4.6.1 from

Gene H. Golub and Charles F. van Loan, Matrix Computations, 3rd edition, John Hopkins University Press.

Note that the input vector $f$ is overwritten by the solution vector $a$.

In :
function solve_vander(x::Vector, f::Vector, n::Int)
# Algorithm 4.6.1 from Golub/van Loan
for k=1:n-1
for i=n:-1:k+1
f[i] = (f[i]-f[i-1])/(x[i]-x[i-k])
end
end
for k=n-1:-1:1
for i=k:n-1
f[i] = f[i]-f[i+1]*x[k]
end
end
end

Out:
solve_vander (generic function with 1 method)


The largest eigenvalue $\delta$ and the corresponding eigenvector $h$ of the eigenvalue problem $$V^{-T}Lh=\delta h$$ are computed using the power method, cf. Golub/van Loan, section 7.3. Some heuristics is used to terminate the iteration and to determine the number of valid digits of $\delta$.

In :
function power_method_for_delta_h(h::Vector, g::Vector, x::Vector, n::Int, maxit::Int)
L=zeros(typeof(g),n,n)
generate_L(L, g, x, n)
x2=x.^2
delta=4.669
deltaold=delta
diff=100
diffold=diff
for k=1:maxit
h = L*h
solve_vander(x2, h, n)
delta = h
if mod(k,20)==0
diff = abs(delta-deltaold)
println(k, " ",delta, @sprintf(" %.3e", diff))
deltaold = delta
if diff>.9*diffold
break
end
diffold=diff
end
h = h/delta
end
digits = int(floor(-exponent(diff)*log(2)/log(10)))
println(digits, " digits seem to be correct:")
println("δ=", string(delta)[1:digits+2])
return delta
end

Out:
power_method_for_delta_h (generic function with 1 method)


## Calculation of $\delta$¶

In :
h=zeros(x)
h[1:3]=[1.0, -3.2565e-1, -5.05539e-2]
maxit=1000
delta = power_method_for_delta_h(h, g, x, n, maxit)

20 4.669201609102999926426746612346816322976798004358427711661775266033202248405129779491070533538575943633752656006028662821261199748985724384748529510626878468162288077528724872312169226986559041791619309843822282958360441836698606915691022758417070998550424637905553881132507045407028243859840359177717434828491e+00 2.016e-04
40 4.669201609102990671853203820847709820449049472100067932428062853939235982738979837367783268806092402419920436486078031856810859889275305791531301519497115202488070341052489217890045563183959278417652025575502118690946967103312321257097931524085183519067993860906561133042915303105011463133590455970913563812374e+00 9.255e-15
60 4.669201609102990671853203820466201617258201304673616153799230285192890332466568665652301194923468984969314331357906336664579208783498224564136563211101446291126081950685799091955108749104242829051601033094986616495821387048879476196871232509016706274740610778366650792754767949618951607890132065804030301127351e+00 3.815e-28
80 4.669201609102990671853203820466201617258185577475768633393985376365213337559206644899033678176455346427933608947668435741775506668575909205234861437888777124694123308629032318952527236332652194967327998193855244373813418535200317599891304026378287837751858032080703956514658719026835995507746792635141885799846e+00 1.573e-41
100 4.669201609102990671853203820466201617258185577475768632745651343004161056969818668524635980210551218914429644052104995161135642627212124667423348504495109567189494156525880474993592260546850478583681782115442716259034840667647145227554040258910558404543087018109161090556694501520954279258162941278648811177398e+00 6.483e-55
120 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473824046688878420164958076300920083411222466345449951477183296818001777791689440509747933260632691637892360636537004990308860832835177264076835155732185022940699034134539843021785910238983330803769154575863399634980379032291745519e+00 2.673e-68
140 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402399343320001173832627483032007861014759284765265974800787304270383496913008859176869362975786867162432497609987828829233791973966564956352459328258363099612474753680096909568246845159095537592622434743568712209701233e+00 1.102e-81
160 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165986727546986364363212672834224911012296211871111263283536457885549249006954537015625107632866767596750443818262903977550147556898616637494227217927113628015439161376166890711236819289670392904550563389082e+00 4.542e-95
180 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402394801381716598485518981513448581272133361893301716717185243561810205583961109875989613586367782008579924049322303710575782712997368634861360390297844404948573589969923607736913943919922768826986661275871148405606348e+00 1.872e-108
200 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027935704191276951344665830943037739320283731250797505916046929854907873429837304381340338726135725905961831993396695334201333337384370843404075063583179458496898099104281546082597e+00 7.719e-122
220 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442989022059829655218946787382924448896409135559317130667196566459049017300658394948223960990812674877569958788802961956170824446761723970367719936857316095020561929087e+00 3.182e-135
240 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402394801381716598485518981513440862714202793252231244298889089085994494087051320402232762381254729200045828140135448316704517094566220838561880560687266328949981574240720104191009054867548933297095110066331969874715082e+00 1.312e-148
260 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411755389915073338192521800339607089183359308457793245796319513409576569410758226888429739106052222423880542700128794234495009713857846773e+00 5.407e-162
280 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411532481714219956644755759327350153832544868089009904734075761337726427598254205117052957433460652988077269262504062534181096688697690804e+00 2.229e-175
300 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411532481714219947455644365824172011617637137504526879245489973358460686341096371959409981370332558869052525740674431975914271972142357185e+00 9.189e-189
320 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411532481714219947455644365823793202009561073946527102451336578758317818555458071823844010891415825518120376221710091931032854708154335229e+00 3.788e-202
340 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411532481714219947455644365823793202009561058330575458618295970118332431369604091203206046119926569201601213249782315567617303786007381624e+00 1.562e-215
360 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411532481714219947455644365823793202009561058330575458617652222070385437176527519468262903909473236379642941204596161167675169757857767961e+00 6.437e-229
380 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411532481714219947455644365823793202009561058330575458617652222070385410638680184519253899662908677158706713939627213923219341705770516754e+00 2.654e-242
400 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402394801381716598485518981513440862714202793252231244298889089085994493546323671341153248171421994745564436582379320200956105833057545861765222207038541064112594159351608651357901840171359689026539162945227072177954464e+00 2.446e-246
420 4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411532481714219947455644365823793202009561058330575458617652222070385410636891694051720946214756662673113179131141684681119728955465577156e+00 4.234e-246
245 digits seem to be correct:
δ=4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402394801381716598485518981513440862714202793252231244298889089085994493546323671341153248171421994745564436582379320200956105833057545861765222207038541063


Out:
4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142027932522312442988890890859944935463236713411532481714219947455644365823793202009561058330575458617652222070385410636891694051720946214756662673113179131141684681119728955465577156e+00 with 1024 bits of precision


Here the value from http://www.plouffe.fr/simon/constants/feigenbaum.txt for comparison:

4.669201609102990671853203820466201617258185577475768632745651343004134330211314737138689744023948013817165984855189815134408627142 02793252231244298889089085994493546323671341153248171421994745564436582379320200956105833057545861765222207038541064674949428498145 33917262005687556659523398756038256372256480040951071283890611844702775854285419801113440175002428585382498335715522052236087250291 6788603626745272133990571316068753450834339344461037063094520191158769724322735898389037949462572512890979489867683346116