N振動子で蔵本モデル

マッチョにkを求めるスタイルは変わらず。そこでループ使って表したとてあまりすっきりしない気がした故。

--------------------------------------------------------
# simulating oscillation by Kuramoto model
# N units

import numpy as np
import matplotlib.pyplot as plt

def theta_mat(theta):
theta=theta.reshape(1, -1)
theta_mat=theta.copy()
while theta_mat.shape[0]<theta_mat.shape[1]:
theta_mat=np.r_[theta_mat, theta]
return theta_mat

#setting
theta=np.array([0, np.pi/2, np.pi]).reshape(1, -1) # start point
omega=np.array([0.3, 1, 0.6]) # anglular velocity
N=theta.size
time=400 # measurement time
h=0.1 # micro time
K=np.array([[0, 0.4, 0.6],[0.4, 0, 0.7],[0.6, 0.7, 0]]) # coupling strength(N*N)

t=0
while t*h<time:
theta_now=theta[t, :].reshape(-1, 1)
k1=(omega+np.sum(np.sin(theta_mat(theta_now)-theta_now)*K, 1)/N).reshape(-1, 1)
theta4k2=theta_now+(h/2)*k1
k2=(omega+np.sum(np.sin(theta_mat(theta4k2)-theta4k2)*K, 1)/N).reshape(-1, 1)
theta4k3=theta_now+(h/2)*k2
k3=(omega+np.sum(np.sin(theta_mat(theta4k3)-theta4k3)*K, 1)/N).reshape(-1, 1)
theta4k4=theta_now+h*k3
k4=(omega+np.sum(np.sin(theta_mat(theta4k4)-theta4k4)*K, 1)/N).reshape(-1, 1)

theta_next=theta_now+(h/6)*(k1+2*k2+2*k3+k4)
theta_next=np.mod(theta_next, 2*np.pi)
theta=np.r_[theta, theta_next.reshape(1, N)]

t+=1

axis=np.arange(theta.shape[0])*h
for i in range(N):
la='theta'+np.str(i)
plt.plot(axis, theta[:, i], label=la)

plt.title('oscilation model by kuramoto model with runge kutta method')
plt.xlabel('time')
plt.ylabel('position')
plt.legend()
plt.show()
--------------------------------------------------------

今までのイメージは「小惑星に別の小惑星が近づいてきてスイングバイする」ようなもの。そのため、青線がゆるやかになるタイミングが若干早くない？と思ったのだがそれは間違い。蔵本モデルでは角速度に影響を与えるのはsinの位相差なので、位相差が垂直つまり±pi/2で大きいのであった。

(11/9 こっそり修正した)
(11/11 またも修正した)