""" File: Kuramoto.py Author: Martin Sauer Email: sauer@math.tu-berlin.de Date: June 24, 2014 Simulation of Kuramoto oscillators with disorder in both frequencies and connection matrix """ import numpy as np from matplotlib import pyplot as plt from matplotlib import animation from matplotlib import patches """ Class definitions """ class KuramotoNetwork: """KuramotoNetwork class This class defines a network of Kuramoto oscillators with possible disorder in the individual frequencies and also the connection matrix. init_state is [N, dt, theta0_j, omega_j, K_ij] where N : number of oscillators, dt : stepsize of numerical integration of the ODE, theta0_j : vector of length N of initial phases, omega_j : vector of length N with individual characteristic frequencies, K_ij : NxN matrix of the coupling strength. """ def __init__(self, init_state = [64, 0.05, np.random.uniform(0, 2 * np.pi, 64), np.random.normal(0, 1, 64), 1. / 64 * np.ones([64, 64])] ): self.size = init_state[0] self.state = np.asarray(init_state[2], dtype='float') self.omega = np.asarray(init_state[3], dtype='float') self.coupling = np.asarray(init_state[4], dtype='float') self.stepsize = init_state[1] self.time_elapsed = 0 def timestep(self): """ Time evolution of the network""" # Numerical scheme is the fourth order Runge-Kutta method """ k1 = np.zeros(self.size) for i in xrange(0, self.size): k1[i] = np.dot(self.coupling[i,], np.sin(self.state - self.state[i])) k2 = np.zeros(self.size) for i in xrange(0, self.size): k2[i] = np.dot(self.coupling[i,], np.sin(self.state - self.state[i] + 0.5 * self.stepsize * (k1 - k1[i]))) k3 = np.zeros(self.size) for i in xrange(0, self.size): k3[i] = np.dot(self.coupling[i,], np.sin(self.state - self.state[i] + 0.5 * self.stepsize * (k2 - k2[i]))) k4 = np.zeros(self.size) for i in xrange(0, self.size): k4[i] = np.dot(self.coupling[i,], np.sin(self.state - self.state[i] + self.stepsize * (k3 - k3[i]))) self.state += self.stepsize * self.omega + self.stepsize / 6 * (k1 + 2 * k2 + 2 * k3 + k4) """ # Numerical scheme is the Heun method (second order) k1 = np.zeros(self.size) for i in xrange(0, self.size): k1[i] = np.dot(self.coupling[i,], np.sin(self.state - self.state[i])) k2 = np.zeros(self.size) for i in xrange(0, self.size): k2[i] = np.dot(self.coupling[i,], np.sin(self.state - self.state[i] + self.stepsize * (k1 - k1[i]))) self.state += self.stepsize * self.omega + self.stepsize / 2 * (k1 + k2) # Numerical scheme is the Euler method (first order) """ k1 = np.zeros(self.size) for i in xrange(0, self.size): k1[i] = np.dot(self.coupling[i,], np.sin(self.state - self.state[i])) self.state += self.stepsize * (self.omega + k1) """ self.state = np.mod(self.state, 2 * np.pi) self.time_elapsed += self.stepsize def position(self): pos = np.exp(self.state * complex(0,1.)) xy = [np.real(pos), np.imag(pos)] return xy def order_param(self): mean = 1./self.size * sum(np.exp(self.state * complex(0,1.))) r = np.abs(mean) psi = np.angle(mean) return r, psi def position_cosine(self): return np.cos(self.state) class KuramotoNetworkConnectionDisorder(KuramotoNetwork): """KuramotoNetworkConnectionDisorder class Introducing a special form of disorder in the connections between oscillators using K_ij = K/N + sigma/sqrt(N) * N(0,1) iid random numbers. init_state is [N, dt, K, sigma] where N : number of oscillators, dt : stepsize of numerical integration of the ODE, K : coupling strength, sigma : strength of coupling disorder. """ def __init__(self, init_state = [64, 0.05, 1, 1] ): self.size = init_state[0] self.state = np.random.uniform(0, 2 * np.pi, self.size) self.omega = np.random.normal(0, 1, self.size) self.K = np.asarray(init_state[2], dtype='float') self.sigma = np.asarray(init_state[3], dtype='float') self.coupling = (self.K / self.size * np.ones([self.size, self.size]) + self.sigma / np.sqrt(self.size) * np.random.normal(0, 1, [self.size, self.size])) self.stepsize = init_state[1] self.time_elapsed = 0 def local_field(self): field = (1 / self.sigma * np.dot(self.coupling, np.exp(self.state * complex(0,1.)))) return [np.real(field), np.imag(field)] def mean_local_field(self): field = (1 / self.sigma * np.dot(self.coupling, np.exp(self.state * complex(0,1.)))) return np.mean(np.abs(field)) """ Definitions of functions for animation """ def animation_kuramoto(N, K): """ Animation of Kuramoto oscillators on the circle with plots of the order parameters evolving in time. N : network size K : coupling constant. """ network = KuramotoNetwork([N, 0.05, np.random.uniform(0, 2 * np.pi, N), np.random.normal(0, 1, N), 1. * K / N * np.ones([N, N])]) fig = plt.figure() # First axes object for plot on circle ax_circ = plt.axes([0.05,0,0.4,1], aspect='equal', xlim=(-1.5, 1.5), ylim=(-1.5, 1.5)) # empty plot of dots for the oscillators dots, = ax_circ.plot([], [], 'o') # plot of unit circle unit_circle = patches.Circle((0,0), 1, ls='dashed', lw=1, fill=False, color='k') ax_circ.add_patch(unit_circle) # empty plot of line representing the order parameters line, = ax_circ.plot([], [], color='r', lw=2) # Second axes object for evolution of r ax_r = plt.axes([0.55, 0.55, 0.4, 0.4], xlim=(0, 1), ylim=(-0.1, 1.1)) line_r, = ax_r.plot([], [], lw=2) # Third axes object for evolution of psi ax_psi = plt.axes([0.55, 0.05, 0.4, 0.4], xlim=(0, 1), ylim=(-np.pi - 0.1, np.pi + 0.1)) line_psi, = ax_psi.plot([], [], lw=2) # Create the actual animation object def init(): dots.set_data([], []) line.set_data([], []) line_r.set_data([], []) line_psi.set_data([], []) return dots, line, line_r, line_psi, def animate_kuramoto(t): # Do the time evolution network.timestep() # Update plot data dots.set_data(network.position()) [r, psi] = network.order_param() line.set_data([0, r * np.cos(psi)], [0, r * np.sin(psi)]) line_r.set_data(np.append(line_r.get_xdata(), network.time_elapsed), np.append(line_r.get_ydata(), r)) ax_r.set_xlim(right=1+network.time_elapsed) line_psi.set_data(np.append(line_psi.get_xdata(), network.time_elapsed), np.append(line_psi.get_ydata(), psi)) ax_psi.set_xlim(right=1+network.time_elapsed) return dots, line, line_r, line_psi, ani = animation.FuncAnimation(fig, animate_kuramoto, frames=100, interval=1000*network.stepsize, blit=False, init_func=init) plt.show() def asymptotics(N, K1, K2, T, dt): """ Long time analysis of the order parameter for two different coupling strengths N : network size K1 : coupling constant network 1 K2 : coupling constant network 2 T : time horizon dt : stepsize """ network1 = KuramotoNetwork([N, 0.05, np.random.uniform(0, 2 * np.pi, N), np.random.normal(0, 1, N), 1. * K1 / N * np.ones([N, N])]) network2 = KuramotoNetwork([N, 0.05, np.random.uniform(0, 2 * np.pi, N), np.random.normal(0, 1, N), 1. * K2 / N * np.ones([N, N])]) [r1, r2] = np.zeros([2, T]) for t in xrange(T): network1.timestep() network2.timestep() tmp = network1.order_param() r1[t] = tmp[0] tmp = network2.order_param() r2[t] = tmp[0] plt.plot(xrange(T), r1, lw=2) plt.plot(xrange(T), r2, lw=2) plt.show() #animation_kuramoto(200, 2.) asymptotics(200, 0.5, 2., 1000, 0.05)