Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

Hello. I am working on investigation of family of dynamical systems on the torus $$\dot{x}=\cos(x)+b\cos(t)+a$$ $$\dot{t}=1$$ and it's Poincare map $$P:(x,0) \rightarrow (P(x),2\pi=0)$$ I need to find Arnold tongues of map $P$. I tried simple calculation of solution using Runge-Kutta formulas, then iterating and checking rotation number, but it's not working effectively. Arnold tongues was first calculated in 1970s so maybe there is effective algorithm of doing it?

share|improve this question
    
I'm working on the same family now; Rudolf, do we work in one seminar?;) –  Olga Sep 27 '13 at 10:03
add comment

1 Answer

up vote 8 down vote accepted

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.

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)
share|improve this answer
    
@Michele Thanks a lot! Send this algorithm to me please. –  Rudolf Dec 14 '12 at 15:30
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.