シミュレーションのアニメーション

二重振り子やるにしても可視化できなけりゃ面白くないので、とりあえずアニメーションのgifを作れるようになろうとしたもの。また蔵本モデル使う。

update関数内でtからt+1の計算を行えば、今の様にデータをすべて保持する必要がなさそうだということはわかったのだが、昨日までのコードをそのまま使った方が楽なので、一旦全部計算してます。
plt.axes().set_aspect('equal', 'datalim')をネットから見つけ出すのに時間がかかった。
そしてそれよりもimagemagickの設定の方がさらに時間かかった。
--------------------------------------------------------
#-*- coding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

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

#making model data
#setting

theta=np.array([0, 0.5*np.pi, np.pi]).reshape(1, -1)
omega=np.array([np.pi, 1.5*np.pi, 2.0*np.pi])     # anglular velocity
N=theta.size
time=30       # measurement time
h=0.05        # micro time
K=np.array([[0, 1.0, 1.5],[1.0, 0, 0.5], [1.5, 0.5, 0]])  # coupling strength(N*N)

t=1
while t*h<time:
    theta_now=theta[t-1, :].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

def update(i):
    plt.clf()
    circle=plt.Circle((0, 0), radius=1.0, fc='y')
    plt.gca().add_patch(circle)
    circle1=plt.Circle((np.cos(theta[i, 0]), np.sin(theta[i, 0])), radius=0.1, fc='b')
    plt.gca().add_patch(circle1)
    circle2=plt.Circle((np.cos(theta[i, 1]), np.sin(theta[i, 1])), radius=0.1, fc='r')
    plt.gca().add_patch(circle2)
    circle3=plt.Circle((np.cos(theta[i, 2]), np.sin(theta[i, 2])), radius=0.1, fc='g')
    plt.gca().add_patch(circle3)
 
    plt.xlim(-1.5, 1.5)
    plt.ylim(-1.5, 1.5)
    plt.axes().set_aspect('equal', 'datalim')

fig=plt.figure()
ani=animation.FuncAnimation(fig, update, frames=theta.shape[0])
ani.save('oscillation.gif', writer='imagemagick', fps=np.int(1/h))
--------------------------------------------------------
Exception ignored in: <bound method TimerQT.__del__ of <matplotlib.backends.backend_qt5.TimerQT object at 0x000001D6221C4BA8>>
Traceback (most recent call last):
  File "C:\Users\_____\Anaconda3\lib\site-packages\matplotlib\backends\backend_qt5.py", line 201, in __del__
TypeError: 'method' object is not connected

なんかエラー吐かれたけどgif自体は出力されてるのでよしとする。
(分かる人通りかかったら教えてください)

何はともあれアニメーションできましたとさ。

好奇心から振動子を増やしてみる。
いい感じで同期しないかなーと思い、Kを若干大きめに。
振動子をグラデーションにしようかとも思ったが面倒なのでパス。
--------------------------------------------------------
#-*- coding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

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

#making model data
#setting

theta=np.random.rand(100).reshape(1, -1)
theta=theta*2*np.pi
omega=np.random.rand(100)*np.pi
N=theta.size
time=30       # measurement time
h=0.05        # micro time
K=np.random.rand(10000).reshape(100, 100)+1
K=K+K.T
for i in range(K.shape[0]):
    K[i, i]=0

t=1
while t*h<time:
    theta_now=theta[t-1, :].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

def update(i):
    plt.clf()
 
    circle=plt.Circle((0, 0), radius=1.0, fc='y')
    plt.gca().add_patch(circle)
 
    for j in range(theta.shape[1]):
        circles=plt.Circle((np.cos(theta[i, j]), np.sin(theta[i, j])), radius=0.05, fc='b')
        plt.gca().add_patch(circles)
     
    plt.xlim(-1.5, 1.5)
    plt.ylim(-1.5, 1.5)
    plt.axes().set_aspect('equal', 'datalim')

fig=plt.figure()
ani=animation.FuncAnimation(fig, update, frames=theta.shape[0])
ani.save('manyoscillation.gif', writer='imagemagick', fps=np.int(1/h))
--------------------------------------------------------

即固まった。

ところで蔵本モデルは弱結合を仮定しているのに逆行する様子がみられますが、振動子が多い場合にも逆行はタブーなんでしょうか。

蔵本モデルで遊ぶ

前に作った蔵本モデル(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をできるようにならんとなァといったところです。

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

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 またも修正した) 

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()
----------------------------------------------------------

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

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