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

Feigenbaum constant calculation.

For more info:

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) 11 years ago  # | flag

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)
Python recipes (4591)
FB36's recipes (148)

Required Modules

  • (none specified)

Other Information and Tasks