Welcome, guest | Sign In | My Account | Store | Cart

Feigenbaum constant calculation.

http://en.wikipedia.org/wiki/Feigenbaum_constant

Python, 38 lines
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38``` ```# Feigenbaum constant calculation # FB - 201011151 # http://en.wikipedia.org/wiki/Feigenbaum_constant # Source: # http://keithbriggs.info/documents/how-to-calc.pdf import math def b(a, k): if k == 0: return 0 else: return a - math.pow(b(a, k - 1), 2) def b_prime(a, k): if k == 0: return 0 else: return 1 - 2 * b_prime(a, k - 1) * b(a, k - 1) # main print 'Feigenbaum constant calculation:' maxIt = 10 maxItJ = 10 ai_1 = 1 ai_2 = 0 di_1 = 3.2 print 'i', 'di' for i in range(2, maxIt): ai = ai_1 + (ai_1 - ai_2) / di_1 for j in range(maxItJ): ai = ai - b(ai, math.pow(2, i)) / b_prime(ai, math.pow(2, i)) # di = (ai_1 - ai_2) / (ai - ai_1) print i, di # prepare for the next iteration di_1 = di ai_2 = ai_1 ai_1 = ai ```

#### 1 comment

FB36 (author) 13 years, 3 months ago

This version does not use recursive functions, so it is way faster:

``````# FB - 201011206
# http://en.wikipedia.org/wiki/Feigenbaum_constant
# http://keithbriggs.info/documents/how-to-calc.pdf
print 'Feigenbaum constant calculation:'
maxIt = 13
maxItJ = 10
a_1 = 1.0
a_2 = 0.0
d_1 = 3.2
print 'i ', 'd'
for i in range(2, maxIt):
a = a_1 + (a_1 - a_2) / d_1
for j in range(maxItJ):
x = 0
y = 0
for k in range(2 ** i):
y = 1 - 2 * y * x
x = a - x * x
a = a - x / y
d = (a_1 - a_2) / (a - a_1)
print "%02i %f" % (i, d)
d_1 = d
a_2 = a_1
a_1 = a
``````
 Created by FB36 on Tue, 16 Nov 2010 (MIT)

### Required Modules

• (none specified)