Personal tools

superstable.py

A python module to compute (and load within your workspace) the first few superstable points for the period doubling cascade.

superstable.py — Python Source, 1Kb

File contents

#!
from pylab import *
nmax = 10
#rseed = array( [ 2.0, 3.2360679771751166, 3.498561699315905 ], float128)

global n

def f(r):
    global n 
    x=float128(0.5)
    for i in range(2**n):
        x=r*x*(1.0-x)
    return x-0.5

def dicho( g, x0, x1):
    f0 = g(x0)
    f1 = g(x1)
    if(f0*f1>0):
        print 'no root in interval',x0,x1
        return []
    while (x1-x0) > 1e-19*(abs(x0)+abs(x1)):
        xm = 0.5*(x0+x1)
        fm = g(xm)
        if( fm*f0 > 0): # fm same sign as f0, root is on (xm,x1)
            x0 = xm
            f0 = fm
        else:
            x1 = xm
            f1 = fm
    return xm

def getsuperstablevalues():
    global n 
    rseed = zeros(nmax,float128)
    
    rseed[0] = 2
    n = 1
    rseed[1]=dicho(f,float128(3),float128(3.4))
    n = 2
    rseed[2]=dicho(f,float128(3.45),float128(3.5))
    
    for i in range(3):
        print "%3d %25.20lf" % (i,rseed[i])
    
    for n in range(3,nmax):
        rmin = rseed[n-1] + (rseed[n-1]-rseed[n-2])/6.0
        rmax = rseed[n-1] + (rseed[n-1]-rseed[n-2])/2.0
        newr = dicho(f,rmin,rmax)
        rseed[n]=newr
        delta = (rseed[n-2]-rseed[n-3])/(rseed[n-1]-rseed[n-2])
        print "%3d %25.20lf %20.15lf" % (n,newr,delta)
    return rseed

if __name__ == '__main__':
    rsuper = getsuperstablevalues()
Document Actions