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

二重振り子やるにしても可視化できなけりゃ面白くないので、とりあえずアニメーションの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))
--------------------------------------------------------

即固まった。

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