Landau Stuart Oscillator

Equation Used

$$\begin{equation} \dot{x}_i = (\lambda - x_i^2 - y_i^2)x_i -\omega y_i + \kappa A_{ij} \sum_{i,j \in n} (x_j - x_i) \\ \dot{y}_i = (\lambda - x_i^2 - y_i^2)y_i +\omega x_i \end{equation}$$



In [1]:

    
from numpy import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import networkx as nx
import plotly.plotly as py



In [2]:

    
%matplotlib inline



In [55]:

    
state_init = random.normal(2.0, 1.0,[30,2] )
xs,ys = array(zip(*state_init))
state_init = concatenate((xs,ys))
g = nx.powerlaw_cluster_graph(30,4, 0.7)
adjm = nx.adjacency_matrix(g).todense()
gw = random.normal(2.0, 1.0, 30)
gws = empty(30)
gws.fill(2.0)

t = linspace(0.0, 10.0, 1000)

Adaptive Coupling

3 types of coupling have been implemeted

Constant Coupling
Quadratic Coupling $ k = 100(d - d_0)^2$ here $ d = \sqrt{\mid{x_i^2 - x_j^2 \mid}}$
Linear Coupling $ k = 100 d$



In [56]:

    
def kfunc_const(state_init,val =100.0):
    return val

def kfunc_quad(state_init):
    d0 = 2.0
    dim = len(state_init)/2
    state_init = state_init[:dim]
    state_init = state_init.reshape([1,-1])
    d = sqrt(abs((state_init)**2 - (state_init.T)**2))
    k_mat = 100.0*(d - d0)**2
    return k_mat

def kfunc_lin(state_init):
    dim = len(state_init)/2
    state_init = state_init[:dim]
    state_init = state_init.reshape([1,-1])
    d = sqrt(abs((state_init)**2 - (state_init.T)**2))
    return 10.0*d



In [57]:

    
params_c = (gw, kfunc_const, adjm)
params_q = (gw, kfunc_quad, adjm)
params_l = (gw, kfunc_lin, adjm)
params_c_w = (gws, kfunc_const, adjm)
params_l_w = (gws, kfunc_lin, adjm)



In [58]:

    
#d = nx.degree(g)
#colmap = plt.cm.coolwarm
pos = nx.spring_layout(g, iterations=5000)
#nx.draw(g,pos, nodelist=d.keys(), node_color=[v  for v in d.values()], cmap = plt.cm.coolwarm)
#plt.colorbar(sc)
#plt.draw()
d = g.degree().values()
node_color = d
plt.figure(figsize=(15,10), frameon=False)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
plt.colorbar(nx.draw_networkx_nodes(g,pos,node_color=node_color))
nx.draw_networkx_edges(g, pos)
#unique_url = py.plot_mpl(fig, filename="plotly version of an mpl figure")
plt.show()



In [59]:

    
def coupling(state_init):
    state_init = state_init.reshape([1,-1])
    return state_init - state_init.T



In [60]:

    
def landau_stuart(state_init, t, gw, k_fun, adjm, L=4.0, cpl="x"):
    k = k_fun(state_init)
    
    ydim = len(state_init)/2
    state_init = state_init.reshape(2, ydim)
    xs = state_init[0]
    ys = state_init[1]
    if cpl == "x":
        coupling_matrix = coupling(xs)
    if cpl == "y":
        coupling_matrix = coupling(ys)
    coupling_matrix = multiply(adjm, coupling_matrix)
    coupling_matrix = multiply(k, coupling_matrix)
    coupling_terms = sum(coupling_matrix, axis =1).T
    #coupling_terms = multiply(k,coupling_terms)
    coupling_terms = squeeze(asarray(coupling_terms))      
    xsd = multiply((L - xs**2 -ys**2), xs) - multiply(gw, ys) + coupling_terms
    ysd = multiply((L- xs**2 - ys**2), ys) + multiply(gw,xs) 
    zs = array([xsd, ysd])
    zs = zs.reshape(ydim*2)
    
    
    return zs



In [61]:

    
zs = odeint(landau_stuart, state_init, t, params_l_w)



In [62]:

    
zs.shape









    Out[62]:





(1000, 60)



In [356]:

    
x1, x2,x3,x4,x5,x6,x7,x8,x9,x10, y1, y2,y3,y4,y5,y6,y7,y8,y9,y10 = zip(*zs)



In [74]:

    
mpl_fig = plt.figure()
plt.plot(t,zs)
plt.xlabel("Time")
plt.ylabel("Zs")
unique_url = py.plot_mpl(mpl_fig, filename="plotly version of an mpl figure")
plt.show()









    



The draw time for this plot will be slow for clients without much RAM.



In [63]:

    
def xsys(state_final):
    div = len(state_final[0])/2
    xs = []
    ys = []
    for i in state_final:
        xs.append(i[:div])
        ys.append(i[div:])
    return xs,ys



In [64]:

    
x,y = xsys(zs)



In [67]:

    
mpl_fig = plt.figure(figsize=(50,5))
plt.plot(t,x)
plt.title("Time evolution of x coordinate")
plt.xlabel("Time")
plt.ylabel("x")
#unique_url = py.plot_mpl(mpl_fig, filename="plotly version of an mpl figure")
plt.show()



In [71]:

    
mpl_fig = plt.figure(figsize=(50, 5))
plt.plot(t,y)
plt.title("Time evolution of y coordinate")
plt.xlabel("Time")
plt.ylabel("y")
#unique_url = py.plot_mpl(mpl_fig, filename="plotly version of an mpl figure")
plt.show()



In [68]:

    
plt.figure(figsize=(10,8))
plt.plot(x,y)
plt.xlim((-2,2))
plt.ylim((-2.2,2.5))
plt.title("Y vs X")
plt.show()



In [69]:

    
xinds,yinds = array(zip(*x)), array(zip(*y))



In [70]:

    
plt.figure(figsize=(20,20))
colors = plt.cm.gnuplot(linspace(0, 1, 30))
for i in range(30):
    plt.subplot(6,5,i+1)
    plt.plot(xinds[i], yinds[i], color = colors[i])
    plt.title("x%s vs y%s" %(i+1, i+1))



In [ ]: