2振動子で蔵本モデル

某所で蔵本モデルをルンゲクッタ法を用いてシミュレーションしていたのだが、そんなに難しくなさそうなのでPythonで実装してみることにした。

困ったのは、ルンゲクッタ法にも色々種類があるようだということ。とっかかりなので一番簡単そうな古典的ルンゲクッタ法を用いることにする。
とりあえず振動子2つ版を実装。順繰りにマッチョにkを求める。

今後は実習の合間を使って、ベクトル形式にしてきれいにしたり、複雑化させていこうと思う。
----------------------------------------------------------
# simulating oscillation by Kuramoto model
# 2 unit

import numpy as np
import matplotlib.pyplot as plt

#setting
theta=np.array([0, np.pi]).reshape(1, 2) # start point
omega=np.array([1, 1]) # anglular velocity
time=100 # measurement time
h=0.1 # micro time
K=0.4 # coupling strength

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

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, 2)]

t+=1

axis=np.arange(theta.shape[0])*h
plt.plot(axis, theta[:, 0], label='theta1')
plt.plot(axis, theta[:, 1], label='theta2')
plt.title('oscilation model by kuramoto model with runge kutta method')
plt.xlabel('time')
plt.ylabel('position')
plt.legend()
plt.show()
----------------------------------------------------------

すぐ終わりそう。(理解できる範囲が狭いことによる)

これ終わったら二重振り子やります。