Revisions to Numerical calculation of Arnold tongue

added 1509 characters in body

Source Link

edited Mar 19, 2013 at 12:18

1.3k
15
25

EDITEDIT[update]: Recently, I wrote for myself some sage lines implementing the algorithm I described you. This is my second version, now working for rational numbers too : I was originally interested only in irrational rotation numbers, so these lines do not work properly for rational numbers, but I think you can adapt it very easily (comments are welcome to improve it!).

L=9L=8 #length for cf-expansion, depending on your computer, 8 or 9 suggested for a try run
A=100000 #maximum size of single element of the sequence

def partfrac(x):
    return x-floor(x)

##### computing rational approximations given continued fraction expansion
# input b=a continued fraction expansion
# input l=L length of computed expansion
def rational_approximation(b,l):
    p=[0,1]
    q=[1,b[1]]
    for i in range(1,l+1):
        p.append(b[i+1]*p[i]+p[i-1])
        q.append(b[i+1]*q[i]+q[i-1])
    return simplify(p[l+1]/q[l+1])


#computing rotation number of a given circle map f
def rotation(f):
    a=[0]
    orbit=[]
    orbit.append(partfrac(f(0)))
    if orbit[0]==0 :
        print 'map with a fixed point'
        return 0
    
    def shift(x):  #set f(0) as the origin + 1
        if partfrac(x)>orbit[0]:
            return partfrac(x)-1
        return partfrac(x)
        
    def first_return(p,pre_p,y):
        x=shift(f(y))
        while x<pre_p or x>p:
            x=shift(f(x))
        return x

    a.append(1)
    x=orbit[0]
    
    if shift(f(orbit[0]))==0:
        returnprint 'map with fixedperiodic point'point of order 2'
        return 1/2
    
    if shift(f(orbit[0]))<0:
        while shift(f(x))<0:
            a[1]=a[1]+1
            x = shift(f(x))
            if a[1]>10000a[1]>A:
                returnprint 'approximatively 0'
                return 0
            if shift(f(x))==0:
                print 'periodic point'
                a[1]=a[1]+1
                return 1/a[1]
        orbit.append(shift(x))
        z = shift(f(x))
        a.append(0)
        while z>0:
            y = z
            z = first_return(shift(orbit[0]),shift(orbit[1]),z)
            a[2]=a[2]+1
            if a[2]>10000a[2]>A:
                returnprint 'approximatively 0'rational'
                return 1/a[1]
            if z==0:
                print 'periodic point'
                a[2]=a[2]+1
                return rational_approximation(a,1)
        orbit.append(y)
    
    if shift(f(orbit[0]))>0:
        def shift(y):  #set f(0) as the origin
            if partfrac(y)>=orbit[0]:
                return partfrac(y)-1
                return partfrac(y)
        orbit.append(orbit[0]-1)
        a.append(0)
        while shift(f(x))>0:
            a[2] = a[2] + 1
            x = shift(f(x))
            if a[2]>10000a[2]>A:
                returnprint 'approximatively rational'
                return 1/a[1]
            if shift(f(x))==0:
                print 'periodic point'
                a[2]=a[2]+1
                return rational_approximation(a,1)
        orbit.append(shift(x))    
        z = shift(f(x))
    
    for i in range(1,L):
        a.append(0)
        if shift(orbit[i+1])<shift(orbit[i]):
            while z>0:
                y = z
                z = first_return(shift(orbit[i]),shift(orbit[i+1]),z)
                a[i+2]=a[i+2]+1
                if a[i+2]>10000a[i+2]>A:
                    returnprint 'approximatively rational'
                    return rational_approximation(a,i)
                if z==0:
                    print 'periodic point'
                    a[i+2]=a[i+2]+1
                    return rational_approximation(a,i+1)
        if shift(orbit[i+1])>shift(orbit[i]):
             while z<0:
                 y = z
                 z = first_return(shift(orbit[i+1]),shift(orbit[i]),z)
                 a[i+2]=a[i+2]+1
                 if a[i+2]>10000a[i+2]>A:
                    return print 'approximatively rational'
        orbit.append             return rational_approximation(ya,i)
 
    print a  #shows the cf expansion

    #computing rational approximationsif z==0:
    p=[0,1]                 print 'periodic point'
    q=[1,a[1]]                 a[i+2]=a[i+2]+1
    for i in range              return rational_approximation(1a,L+1i+1):
        p orbit.append(a[i+1]*p[i]+p[i-1]y)
     
    q.append(a[i+1]*q[i]+q[i-1])print a
    return simplifyrational_approximation(p[L+1]/q[L+1]a,L) #returns the rational approximation

EDIT: Recently, I wrote for myself some sage lines implementing the algorithm I described you. I was interested in irrational rotation numbers, so these lines do not work properly for rational numbers, but I think you can adapt it very easily (comments are welcome to improve it!).

L=9 #length for cf-expansion

def partfrac(x):
    return x-floor(x)

#computing rotation number
def rotation(f):
    a=[0]
    orbit=[]
    orbit.append(partfrac(f(0)))

    def shift(x):  #set f(0) as the origin + 1
        if partfrac(x)>orbit[0]:
            return partfrac(x)-1
        return partfrac(x)

    def first_return(p,pre_p,y):
        x=shift(f(y))
        while x<pre_p or x>p:
            x=shift(f(x))
        return x

    a.append(1)
    x=orbit[0]

    if shift(f(orbit[0]))==0:
        return 'map with fixed point'

    if shift(f(orbit[0]))<0:
        while shift(f(x))<0:
            a[1]=a[1]+1
            x = shift(f(x))
            if a[1]>10000:
                return 'approximatively 0'
        orbit.append(shift(x))
        z = shift(f(x))
        a.append(0)
        while z>0:
            y = z
            z = first_return(shift(orbit[0]),shift(orbit[1]),z)
            a[2]=a[2]+1
            if a[2]>10000:
                return 'approximatively 0'
        orbit.append(y)

    if shift(f(orbit[0]))>0:
        def shift(y):  #set f(0) as the origin
            if partfrac(y)>=orbit[0]:
                return partfrac(y)-1
            return partfrac(y)
        orbit.append(orbit[0]-1)
        a.append(0)
        while shift(f(x))>0:
            a[2] = a[2] + 1
            x = shift(f(x))
            if a[2]>10000:
                return 'approximatively rational'
        orbit.append(shift(x))
        z = shift(f(x))

    for i in range(1,L):
        a.append(0)
        if shift(orbit[i+1])<shift(orbit[i]):
            while z>0:
                y = z
                z = first_return(shift(orbit[i]),shift(orbit[i+1]),z)
                a[i+2]=a[i+2]+1
                if a[i+2]>10000:
                    return 'approximatively rational'
        if shift(orbit[i+1])>shift(orbit[i]):
            while z<0:
                y = z
                z = first_return(shift(orbit[i+1]),shift(orbit[i]),z)
                a[i+2]=a[i+2]+1
                if a[i+2]>10000:
                    return 'approximatively rational'
        orbit.append(y)
 
    print a  #shows the cf expansion

    #computing rational approximations
    p=[0,1]
    q=[1,a[1]]
    for i in range(1,L+1):
        p.append(a[i+1]*p[i]+p[i-1])
        q.append(a[i+1]*q[i]+q[i-1])
    return simplify(p[L+1]/q[L+1]) #returns the rational approximation

EDIT[update]: Recently, I wrote for myself some sage lines implementing the algorithm I described you. This is my second version, now working for rational numbers too : I was originally interested only in irrational rotation numbers (comments are welcome to improve it!).

L=8 #length for cf-expansion, depending on your computer, 8 or 9 suggested for a try run
A=100000 #maximum size of single element of the sequence

def partfrac(x):
    return x-floor(x)

##### computing rational approximations given continued fraction expansion
# input b=a continued fraction expansion
# input l=L length of computed expansion
def rational_approximation(b,l):
    p=[0,1]
    q=[1,b[1]]
    for i in range(1,l+1):
        p.append(b[i+1]*p[i]+p[i-1])
        q.append(b[i+1]*q[i]+q[i-1])
    return simplify(p[l+1]/q[l+1])


#computing rotation number of a given circle map f
def rotation(f):
    a=[0]
    orbit=[]
    orbit.append(partfrac(f(0)))
    if orbit[0]==0 :
        print 'map with a fixed point'
        return 0
    
    def shift(x):  #set f(0) as the origin + 1
        if partfrac(x)>orbit[0]:
            return partfrac(x)-1
        return partfrac(x)
        
    def first_return(p,pre_p,y):
        x=shift(f(y))
        while x<pre_p or x>p:
            x=shift(f(x))
        return x

    a.append(1)
    x=orbit[0]
    
    if shift(f(orbit[0]))==0:
        print 'map with periodic point of order 2'
        return 1/2
    
    if shift(f(orbit[0]))<0:
        while shift(f(x))<0:
            a[1]=a[1]+1
            x = shift(f(x))
            if a[1]>A:
                print 'approximatively 0'
                return 0
            if shift(f(x))==0:
                print 'periodic point'
                a[1]=a[1]+1
                return 1/a[1]
        orbit.append(shift(x))
        z = shift(f(x))
        a.append(0)
        while z>0:
            y = z
            z = first_return(shift(orbit[0]),shift(orbit[1]),z)
            a[2]=a[2]+1
            if a[2]>A:
                print 'approximatively rational'
                return 1/a[1]
            if z==0:
                print 'periodic point'
                a[2]=a[2]+1
                return rational_approximation(a,1)
        orbit.append(y)
    
    if shift(f(orbit[0]))>0:
        def shift(y):  #set f(0) as the origin
            if partfrac(y)>=orbit[0]:
                return partfrac(y)-1
                return partfrac(y)
        orbit.append(orbit[0]-1)
        a.append(0)
        while shift(f(x))>0:
            a[2] = a[2] + 1
            x = shift(f(x))
            if a[2]>A:
                print 'approximatively rational'
                return 1/a[1]
            if shift(f(x))==0:
                print 'periodic point'
                a[2]=a[2]+1
                return rational_approximation(a,1)
        orbit.append(shift(x))    
        z = shift(f(x))
    
    for i in range(1,L):
        a.append(0)
        if shift(orbit[i+1])<shift(orbit[i]):
            while z>0:
                y = z
                z = first_return(shift(orbit[i]),shift(orbit[i+1]),z)
                a[i+2]=a[i+2]+1
                if a[i+2]>A:
                    print 'approximatively rational'
                    return rational_approximation(a,i)
                if z==0:
                    print 'periodic point'
                    a[i+2]=a[i+2]+1
                    return rational_approximation(a,i+1)
        if shift(orbit[i+1])>shift(orbit[i]):
             while z<0:
                 y = z
                 z = first_return(shift(orbit[i+1]),shift(orbit[i]),z)
                 a[i+2]=a[i+2]+1
                 if a[i+2]>A:
                     print 'approximatively rational'
                     return rational_approximation(a,i)
                if z==0:
                     print 'periodic point'
                     a[i+2]=a[i+2]+1
                     return rational_approximation(a,i+1)
         orbit.append(y)
     
    print a
    return rational_approximation(a,L)

added 2952 characters in body; added 7 characters in body

Source Link

edited Dec 14, 2012 at 16:01

user47274

1.3k
15
25

RecentlyEDIT: Recently, I wrote for myself some sage lines implementing the algorithm I described you. I was interested in irrational rotation numbers, so these lines do not work properly for rational numbers, but I think you can sendadapt it to you if youvery easily (comments are interestedwelcome to improve it!).

L=9 #length for cf-expansion

def partfrac(x):
    return x-floor(x)

#computing rotation number
def rotation(f):
    a=[0]
    orbit=[]
    orbit.append(partfrac(f(0)))

    def shift(x):  #set f(0) as the origin + 1
        if partfrac(x)>orbit[0]:
            return partfrac(x)-1
        return partfrac(x)

    def first_return(p,pre_p,y):
        x=shift(f(y))
        while x<pre_p or x>p:
            x=shift(f(x))
        return x

    a.append(1)
    x=orbit[0]

    if shift(f(orbit[0]))==0:
        return 'map with fixed point'

    if shift(f(orbit[0]))<0:
        while shift(f(x))<0:
            a[1]=a[1]+1
            x = shift(f(x))
            if a[1]>10000:
                return 'approximatively 0'
        orbit.append(shift(x))
        z = shift(f(x))
        a.append(0)
        while z>0:
            y = z
            z = first_return(shift(orbit[0]),shift(orbit[1]),z)
            a[2]=a[2]+1
            if a[2]>10000:
                return 'approximatively 0'
        orbit.append(y)

    if shift(f(orbit[0]))>0:
        def shift(y):  #set f(0) as the origin
            if partfrac(y)>=orbit[0]:
                return partfrac(y)-1
            return partfrac(y)
        orbit.append(orbit[0]-1)
        a.append(0)
        while shift(f(x))>0:
            a[2] = a[2] + 1
            x = shift(f(x))
            if a[2]>10000:
                return 'approximatively rational'
        orbit.append(shift(x))
        z = shift(f(x))

    for i in range(1,L):
        a.append(0)
        if shift(orbit[i+1])<shift(orbit[i]):
            while z>0:
                y = z
                z = first_return(shift(orbit[i]),shift(orbit[i+1]),z)
                a[i+2]=a[i+2]+1
                if a[i+2]>10000:
                    return 'approximatively rational'
        if shift(orbit[i+1])>shift(orbit[i]):
            while z<0:
                y = z
                z = first_return(shift(orbit[i+1]),shift(orbit[i]),z)
                a[i+2]=a[i+2]+1
                if a[i+2]>10000:
                    return 'approximatively rational'
        orbit.append(y)

    print a  #shows the cf expansion

    #computing rational approximations
    p=[0,1]
    q=[1,a[1]]
    for i in range(1,L+1):
        p.append(a[i+1]*p[i]+p[i-1])
        q.append(a[i+1]*q[i]+q[i-1])
    return simplify(p[L+1]/q[L+1]) #returns the rational approximation

Recently, I wrote for myself some sage lines implementing the algorithm I described you. I can send it to you if you are interested.

EDIT: Recently, I wrote for myself some sage lines implementing the algorithm I described you. I was interested in irrational rotation numbers, so these lines do not work properly for rational numbers, but I think you can adapt it very easily (comments are welcome to improve it!).

L=9 #length for cf-expansion

def partfrac(x):
    return x-floor(x)

#computing rotation number
def rotation(f):
    a=[0]
    orbit=[]
    orbit.append(partfrac(f(0)))

    def shift(x):  #set f(0) as the origin + 1
        if partfrac(x)>orbit[0]:
            return partfrac(x)-1
        return partfrac(x)

    def first_return(p,pre_p,y):
        x=shift(f(y))
        while x<pre_p or x>p:
            x=shift(f(x))
        return x

    a.append(1)
    x=orbit[0]

    if shift(f(orbit[0]))==0:
        return 'map with fixed point'

    if shift(f(orbit[0]))<0:
        while shift(f(x))<0:
            a[1]=a[1]+1
            x = shift(f(x))
            if a[1]>10000:
                return 'approximatively 0'
        orbit.append(shift(x))
        z = shift(f(x))
        a.append(0)
        while z>0:
            y = z
            z = first_return(shift(orbit[0]),shift(orbit[1]),z)
            a[2]=a[2]+1
            if a[2]>10000:
                return 'approximatively 0'
        orbit.append(y)

    if shift(f(orbit[0]))>0:
        def shift(y):  #set f(0) as the origin
            if partfrac(y)>=orbit[0]:
                return partfrac(y)-1
            return partfrac(y)
        orbit.append(orbit[0]-1)
        a.append(0)
        while shift(f(x))>0:
            a[2] = a[2] + 1
            x = shift(f(x))
            if a[2]>10000:
                return 'approximatively rational'
        orbit.append(shift(x))
        z = shift(f(x))

    for i in range(1,L):
        a.append(0)
        if shift(orbit[i+1])<shift(orbit[i]):
            while z>0:
                y = z
                z = first_return(shift(orbit[i]),shift(orbit[i+1]),z)
                a[i+2]=a[i+2]+1
                if a[i+2]>10000:
                    return 'approximatively rational'
        if shift(orbit[i+1])>shift(orbit[i]):
            while z<0:
                y = z
                z = first_return(shift(orbit[i+1]),shift(orbit[i]),z)
                a[i+2]=a[i+2]+1
                if a[i+2]>10000:
                    return 'approximatively rational'
        orbit.append(y)

    print a  #shows the cf expansion

    #computing rational approximations
    p=[0,1]
    q=[1,a[1]]
    for i in range(1,L+1):
        p.append(a[i+1]*p[i]+p[i-1])
        q.append(a[i+1]*q[i]+q[i-1])
    return simplify(p[L+1]/q[L+1]) #returns the rational approximation

Source Link

created Dec 14, 2012 at 14:07

user47274

1.3k
15
25

There is a good way to compute rotation number of a circle homeomorphism (this was the way Poincaré thinked of it): you calculate the rotation number buy its continued fraction in a direct way.

You start from a point $x$ and $f(x)$: this gives you a decomposition of the circle into points that are on the right side of $x$ (in $]x,f(x)[$) and points which are on its left side (in $]f(x),x[$). You look at $f^2(x)$ and you write $R$ if it is on the right side of $x$, $L$ otherwise. Iteranting $f$ you find a sequence of $R$'s and $L$'s. If you get $LLLLR$, for example, you record 4 (this is the number of $L$'s) and you approximate the rotation number of $f$ by $1/4$.

Renormalizing $f$, you iterate this process finding $\rho=[0,a_1,a_2,\ldots,a_k]$.

I won't be more precise here.

Every detail is very well explained in de Melo & van Strien's One-Dimensional Dynamics, section I.1.

You can find a paper by Bruin (Numerical determination of the continued fraction expansion of the rotation number) in which he compares different methods on Arnold tongues.

Recently, I wrote for myself some sage lines implementing the algorithm I described you. I can send it to you if you are interested.

Stack Exchange Network

Return to Answer