Resonant Sequences

(from http://journals.aps.org/prstab/abstract/10.1103/PhysRevSTAB.17.014001)



In [1]:

    
!date
%matplotlib inline









    



Thu Apr  9 15:28:26 PDT 2015



In [21]:

    
from __future__ import division

def fareySequence(N, k=1):
    """
    Generate Farey sequence of order N, less than 1/k
    """
    # assert type(N) == int, "Order (N) must be an integer"
    a, b = 0, 1
    c, d = 1, N
    seq = [(a,b)]
    while c/d <= 1/k:
        seq.append((c,d))
        tmp = int(math.floor((N+b)/d))
        a, b, c, d = c, d, tmp*c-a, tmp*d-b
    return seq


def resonanceSequence(N, k):
    """
    Compute resonance sequence

    Arguments:
        - N (int): Order
        - k (int): denominator of the farey frequency resonances are attached to
    """
    a, b = 0, 1
    c, d = k, N-k
    seq = [(a,b)]
    while d >= 0:
        seq.append((c,d))
        tmp = int(math.floor((N+b+a)/(d+c)))
        a, b, c, d = c, d, tmp*c-a, tmp*d-b
    return seq


def plotResonanceDiagram(N, figsize=(10,10)):
    import matplotlib.pyplot as plt

    ALPHA = 0.5/N

    plt.figure(figsize=figsize)
    ticks = set([])
    for h, k in fareySequence(N, 1):
        ticks.add((h,k))
        for a, b in resonanceSequence(N, k):
            if b == 0:
                x = np.array([h/k, h/k])
                y = np.array([0, 1])
            elif a== 0:
                x = np.array([0, 1])
                y = np.array([h/k, h/k])
            else:
                m = a/b
                cp, cm = m*h/k, -m*h/k
                x = np.array([0, h/k, 1])
                y = np.array([cp, 0, cm+m])
            plt.plot(  x,   y, 'b', alpha=ALPHA) # seqs. attached to horizontal axis
            plt.plot(  y,   x, 'b', alpha=ALPHA) # seqs. attached to vertical axis
            # also draw symetrical lines, to be fair (otherwise lines in the 
            # lower left traingle will be duplicated, but no the others)
            plt.plot(  x, 1-y, 'b', alpha=ALPHA) 
            plt.plot(1-y,   x, 'b', alpha=ALPHA)
    plt.xlim(0, 1)
    plt.ylim(0, 1)
    plt.xticks([h/k for h,k in ticks], [r"$\frac{{{:d}}}{{{:d}}}$".format(h,k) for h,k in ticks], fontsize=15)
    plt.yticks([h/k for h,k in ticks], [r"${:d}/{:d}$".format(h,k) for h,k in ticks])
    plt.title("N = {:d}".format(N))

Generating resonance sequences is fast

Try it!

Note: in the original paper there was a minor mistake. Eq. (8) read

$$ \Bigg( \Big\lfloor \frac{N+b+a}{d} \Big\rfloor c - a , \Big\lfloor \frac{N+b+a}{d} \Big\rfloor d - b \Bigg) $$

but it should read

$$ \Bigg( \Big\lfloor \frac{N+b+a}{d+c} \Big\rfloor c - a , \Big\lfloor \frac{N+b+a}{d+c} \Big\rfloor d - b \Bigg) $$

I've contacted Rogelio Tomás (the author) and he agreed with the correction (an erratum will be sent to the publication).



In [22]:

    
N = 5
for k in set([k for _,k in fareySequence(N, 1)]):
    print "N={}, k={}:".format(N, k)
    print "\t", resonanceSequence(N, k)









    



N=5, k=1:
	[(0, 1), (1, 4), (1, 3), (1, 2), (2, 3), (1, 1), (3, 2), (2, 1), (3, 1), (4, 1), (1, 0)]
N=5, k=2:
	[(0, 1), (2, 3), (2, 2), (2, 1), (4, 1), (2, 0)]
N=5, k=3:
	[(0, 1), (3, 2), (3, 1), (3, 0)]
N=5, k=4:
	[(0, 1), (4, 1), (4, 0)]
N=5, k=5:
	[(0, 1), (5, 0)]

..., but plotting can be slow for large N (N > 10)

Try it, but be patient ... (lots of lines to plot)



In [23]:

    
# from matplotlib2tikz import save as save_tikz
plotResonanceDiagram(5, figsize=(12,12))
# save_tikz('resonanceDiagram_N7.tikz')

def plotSolution(a,b,c,color='r'):
    x = [c/a, 0, (c-b)/a, 1]
    y = [0, c/b, 1, (c-a)/b]
    plt.plot(x, y, color=color, alpha=0.5, linewidth=4)

    
# plot some example solutions
if True:
    # solutions for (x,y) = (0.5, 1)
    plotSolution( 4, -1,  1)
    plotSolution(-2,  2,  1)
    plotSolution(-2,  3,  2)
    plotSolution( 4,  1,  3)
    plotSolution( 2,  1,  2)

    # solutions for (x,y) = (0.5, 0.5)
    plotSolution( 3, -1,  1, 'g')
    plotSolution(-1,  3,  1, 'g')
    plotSolution( 3,  1,  2, 'g')
    plotSolution( 1,  3,  2, 'g')
    plotSolution( 1,  1,  1, 'g')



In [ ]: