蔵本モデルで遊ぶ

前に作った蔵本モデル(Nコ振動子のコードでN=2)で同期現象をちゃんと見ようということで初期の数値を変えて遊んでいた。
-----------------------------------------------------
#-*- coding:utf-8 -*-
# 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, 0]).reshape(1, -1)  # start point after 'phase reset'
omega=np.array([0.3, 1.1])     # anglular velocity
N=theta.size
time=100       # measurement time
h=0.1        # micro time
K=np.array([[0, 0.8],[0.8, 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()
-----------------------------------------------------

同期してそう。この時の2振動子の位相差をプロットする。
-----------------------------------------------------
dif=np.abs(theta[:, 0]-theta[:, 1])
for i in range(dif.shape[0]):
    if dif[i]>np.pi:
        dif[i]=2*np.pi-dif[i]

plt.plot(axis, dif, label='dif')
plt.show()
-----------------------------------------------------

差がある大きさに収束していく様子で同期していそうである。

次に、角速度の差を大きくしてみる。
omega=[0.3, 1.2]で実行。
-----------------------------------------------------

捕捉しきれず同期できていない。

ここで振動子2つの場合の蔵本モデルの式を確認する。
dΘ1/dt=ω1+(K/2)*sin(Θ2-Θ1)
dΘ2/dt=ω2+(K/2)*sin(Θ1-Θ2)

同期した場合2振動子の速度が等しくなるので

dΘ1/dt-dΘ2/dt=ω1-ω2-K*sin(Θ1-Θ2)=0
から
ω1-ω2=K*sin(Θ1-Θ2)
正弦波は-1～1の値なので
|ω1-ω2|≦K
であれば同期できるということになる。

とはいえルンゲクッタ法の誤差もあるし、いっぱいやれば同期しなくなるンでは？と思い、omega=[0.3, 1.1], time=10000で実行。

おお、同期している。すごーい！上の図はもうただの長方形にしか見えんな。

振動子の動きを円周上でプロットするgifを作ろうとしていたのだけれど、matplotlibがあまりに苦手で実装できなかったでござる。蔵本モデルならこんな感じでもいいんだけど、二重振り子やろうとするならやはりanimationをできるようにならんとなァといったところです。

(追記)
誤差少しできたところでまた収束の方向に向かうんだからそらそうか