【ポアンカレ再帰定理】蛇振り子を100個・長周期に拡張して遊んでみた♬ #Python

マスゲームだから、意味ない（当たり前）かなと思いつつ、でも惑星直列などを考えると、こういう計算も実は意味のあるものかなと思い、まとめておくこととします。
なんといってもポアンカレの再帰定理。
「力学系は、ある種の条件が満たされれば、その任意の初期状態に有限時間内にほぼ回帰する」の具体的な事例としての意味があるということだと思いなおしました。
なお、惑星直列の計算とシミュレーションは今回はやりません。
ということで、今回は前回の式を個数と周期を自由に変更したいと思います。
※コード全体はおまけに載せました。
　また、gifアニメーションがエラーが出て貼れないのでYoutubeに載せました

100個を1ダンス120秒のシミュレーション

そして、100個の場合に拡張です。
これだと、前回の式では無理なので、以下のように変更します。

N = 60
z0 =500
L0=980*(120/130)**2/(2*np.pi)**2
Fc0=1.0
Fc=[]
fc=1
Fc.append(fc)
for i in range(1,100,1):
    fc=fc*((131-i)/(130-i))
    Fc.append(fc)

あと、120secで元に戻るようにしたので、結果は以下のようになります。

蛇の振り子100個をシミュレーションして遊んでみた♬

1５個を1ダンス120秒のシミュレーション

今度は、１５個の振り子を120秒という長さの周期のダンスをさせてみる。
これで、実質的にあらゆる振り子の個数、長さ、周期のものが設計できることが分かる。
【蛇の振り子】1５個を1ダンス120秒のシミュレーション作成して遊んでみた♬

主要なコードは、以下で実現できた。

N = 60
L0=980*(120/65)**2/(2*np.pi)**2  #ここの120/65の120が重要
print(L0, 2*np.pi/np.sqrt(980/L0))
Fc0=1.0
z0 =150
Fc=[]
fc=1
Fc.append(fc)
for i in range(1,15,1):
    fc=fc*((66-i)/(65-i))  #65/64で1個1個のずらしを調整
    Fc.append(fc)

dataz=[]
linez=[]
y=0
for j in range(15):
    y += 2  #振り子のy軸方向の間隔
    dataz0 = np.array(list(genxy(N,L=L0*Fc[j]**2,z0=z0,y0=y))).T
    linez0, = ax.plot(dataz0[0, 0:1], dataz0[1, 0:1], dataz0[2, 0:1],'o')
    dataz.append(dataz0)
    linez.append(linez0)

generatorは前回のものと比べると、以下のようにy軸方向に配置できるように拡張している。

def genxy(n,L=20,z0=0,y0=0):
    phi = 0
    g = 980
    x0=5  #2.5
    omega = np.sqrt(g/L)
    theta0 = np.arcsin(x0/L)
    while phi < 600:
        yield np.array([L*np.sin(np.sin(phi*omega+np.pi/2)*theta0), y0,z0-L+L*(1-np.cos(np.sin(phi*omega+np.pi/2)*theta0))])
        phi += 1/n

まとめ

・蛇振り子の個数や長さを自由を変更できるように拡張した
・100個の圧倒的な迫力と、120秒周期という長周期は鑑賞という意味で面白い

・惑星直列など自然界にあるポアンカレの再帰定理の具現化をシミュレーションしたいと思う
・現実の蛇振り子を自由な周期で実現したいと思う

おまけ

from matplotlib import pyplot as plt
import numpy as np
import mpl_toolkits.mplot3d.axes3d as p3
from matplotlib import animation

fig, ax = plt.subplots(1,1,figsize=(1.6180 * 4, 4*1),dpi=200)
ax = p3.Axes3D(fig)

def genxy(n,L=20,z0=0,y0=0):
    phi = 0
    g = 980
    x0=5
    omega = np.sqrt(g/L)
    theta0 = np.arcsin(x0/L)
    while phi < 600:
        yield np.array([L*np.sin(np.sin(phi*omega+np.pi/2)*theta0), y0,z0-L+L*(1-np.cos(np.sin(phi*omega+np.pi/2)*theta0))])
        phi += 1/n

def update(num, data, line,s):
    line.set_data(data[:2,num-1 :num])
    line.set_3d_properties(data[2,num-1 :num])
    ax.set_xticklabels([])
    ax.grid(False)

N = 60
L0=980*(120/65)**2/(2*np.pi)**2
print(L0, 2*np.pi/np.sqrt(980/L0)) #84.6cm 1.846sec
Fc0=1.0   #omega=2pi/T
z0 =150
Fc=[]
fc=1
Fc.append(fc)
for i in range(1,15,1):
    fc=fc*((66-i)/(65-i))
    Fc.append(fc)

dataz=[]
linez=[]
y=0
for j in range(15):
    y += 2
    dataz0 = np.array(list(genxy(N,L=L0*Fc[j]**2,z0=z0,y0=y))).T
    linez0, = ax.plot(dataz0[0, 0:1], dataz0[1, 0:1], dataz0[2, 0:1],'o')
    dataz.append(dataz0)
    linez.append(linez0)

# Setting the axes properties
ax.set_xlim3d([-10., 10])
ax.set_xlabel('X')

ax.set_ylim3d([0.0, 30.0])
ax.set_ylabel('Y')

ax.set_zlim3d([0.0, z0-L0])
ax.set_zlabel('Z')
elev=0 #20. 
azim=90 #35.
ax.view_init(elev, azim)

frames =60*120
fr0=60
s0=0
s=s0

while 1:
    s+=1
    num = s
    for j in range(15):        
        update(num,dataz[j],linez[j],s0)
    """
    if s%fr0==0:
        print(s/fr0, s)
    plt.pause(0.001)    
    """

    if s%(fr0/10)==0:
        ax.set_title("s={}_sec".format(int(10*s/fr0)/10),loc='center')
        plt.pause(0.001)
        print(s/fr0)
        plt.savefig('./pendulum/'+str(int(10*s/fr0))+'.png')
        if s>=s0+frames:
            break

s0=30*120
s=s0+30*120
fr0=60
from PIL import Image,ImageFilter
images = []
for n in range(1,1201,1):
    exec('a'+str(n)+'=Image.open("./pendulum/'+str(n)+'.png")')
    images.append(eval('a'+str(n)))
images[0].save('./pendulum/pendulum_{}_.gif'.format(100),
               save_all=True,
               append_images=images[1:],
               duration=100,
               loop=1)

動画保存

gifアニメーションが貼れないので以下のコードでmp4ファイルを保存し、youtubeにアップします。

from matplotlib import pyplot as plt
import numpy as np
import mpl_toolkits.mplot3d.axes3d as p3
import cv2

def cv_fourcc(c1, c2, c3, c4):
        return (ord(c1) & 255) + ((ord(c2) & 255) << 8) + \
            ((ord(c3) & 255) << 16) + ((ord(c4) & 255) << 24)

OUT_FILE_NAME = "output_video.mp4"
fourcc = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
dst = cv2.imread('./pendulum/1.png')
rows,cols,channels = dst.shape
out = cv2.VideoWriter(OUT_FILE_NAME, int(fourcc), int(10), (int(cols), int(rows)))

from PIL import Image,ImageFilter
for n in range(1,1201,1):
    dst = cv2.imread('./pendulum/'+str(n)+'.png')
    out.write(dst) #mp4やaviに出力します