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?



Numerical calculation of Arnold tongue
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*](http://www2.warwick.ac.uk/fac/sci/maths/people/staff/sebastian_van_strien/demelo-strien.pdf), section I.1.


You can find a paper by Bruin ([*Numerical determination of the continued fraction expansion of the rotation number*](http://www.ams.org/mathscinet/search/publdoc.html?arg3=&co4=AND&co5=AND&co6=AND&co7=AND&dr=all&pg4=AUCN&pg5=TI&pg6=PC&pg7=ALLF&pg8=ET&r=1&review_format=html&s4=bruin&s5=numerical%20determination&s6=&s7=&s8=All&vfpref=html&yearRangeFirst=&yearRangeSecond=&yrop=eq)) in which he compares different methods on Arnold tongues.


EDIT[update]:
Recently, I wrote for myself some [sage](http://www.sagemath.org/) 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)

```

