[PYTHON] [Introduction to simulation] Sine wave mass game ♬

I was addicted to this time as well. The following pendulum that came out on Twitter. Wonder of Science @ wonderofscience

Since each is independent, it looks like a group movement. Moreover, the movement reminds me of that Poincare recurrence theorem. "A dynamical system almost returns to its arbitrary initial state within a finite time if certain conditions are met." Well, this time it's a simple simple pendulum, so it's natural if you control the cycle. .. .. However, I intuitively felt that it was not so easy.

I also took the following snaps. You can see that the shortest pendulum cycle is about 1 second, and the longest pendulum cycle seems to be about 4/3 times as long. However, with this information alone, it was unclear how to align the pendulum cycle between them so neatly. After struggling and doing various things, I read the original material for reference. There was a link to the original story from the beginning, but I didn't read it only at this time. [Reference] Original material Pendulum Waves@Harvard Natural Sciences Lecture Demonstrations spoiler. "The duration of one complete cycle of dance is 60 seconds. The longest pendulum length is adjusted to vibrate 51 times in this 60 seconds. The length of consecutive short pendulums is further during this period. It is carefully adjusted to perform one vibration, so the 15th pendulum (shortest) vibrates 65 times. When all 15 pendulums are started together, they are immediately out of sync. Due to the different cycles of relative vibration, the relative phase changes continuously. However, after 60 seconds, they all perform an integral multiple of the vibration, and at that moment they are ready to resynchronize and repeat the dance. . "

Of course, it would be quite difficult to realize this with a real pendulum as a reference. However, I'm going to simulate it here, so I'll try more. 【reference】 [Jaburiko @ Amazon](https://www.amazon.co.jp/ScienceGeek-DIY-%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3% 83% B3-% E3% 83% 90% E3% 83% A9% E3% 83% B3% E3% 82% B9% E3% 83% 89% E3% 83% BC% E3% 83% AB% E3% 83 % 9C% E3% 83% BC% E3% 83% AB-% E6% 8C% AF% E5% AD% 90% E6% B3% A2-% E3% 82% A4% E3% 83% B3% E3% 83 % 86% E3% 83% AA% E3% 82% A2-% E7% 89% A9% E7% 90% 86% E5% AD% A6% E3% 82% A8% E3% 83% 8D% E3% 83% AB% E3% 82% AE% E3% 83% BC / dp / B085T26F1G / ref = pd_sbs_21_1 / 355-3562875-4369549? _Encoding = UTF8 & pd_rd_i = B085T26F1G & pd_rd_r = e634f512-037b-4112-b7c2-b2e8 b490-4864-923d-51639f6a935f & pf_rd_r = VVJHE2FGVTDYCM3J5C6K & psc = 1 & refRID = VVJHE2FGVTDYCM3J5C6K)

Review of the simple pendulum

【reference】 Commentary: Simple pendulum movement Numerical solutions of differential equations are fine, but I decided to use the following analytical solutions here.

\theta = \theta_0\sin(\omega t+\delta)  \\
\omega = \sqrt{\frac{g}{L}}\\
T=\frac{2\pi}{\omega}\\
\theta = \theta_0\sin(\frac{2\pi t}{T}+\delta)  \\
T;sec \\ t; sec\\
L=\frac{g}{\omega ^2}=g(\frac{T}{2\pi})^2\\
g=980 cm/s^2\\
L=24.At 82 cm\\
T=1.000sec

The table above spoilers is as follows

Code commentary

I put the whole code as a bonus. It's almost the same as the Axes 3D animation the other day. The essential parts are as follows.

def genxy(n,L=20,z0=0):
    phi = 0  #Initial time value
    g = 980  #Gravitational acceleration
    x0=2.5  #Initial amplitude
    omega = np.sqrt(g/L)
    theta0 = np.arcsin(x0/L)  #Initial angle
    while phi < 600:  #Pi for starting from maximum amplitude/Add 2
        yield np.array([L*np.sin(np.sin(phi*omega+np.pi/2)*theta0), 0,z0-L+L*(1-np.cos(np.sin(phi*omega+np.pi/2)*theta0))])
        phi += 1/n  #Divide 1sec into n

Call this with the following code and use it Here, N = 60, L is the length of each pendulum. z0 is the height of the fulcrum. Calculate dataz and plot it in linez.

dataz = np.array(list(genxy(N,L=L0*Fc[0]**2,z0=z0))).T
linez, = ax.plot(dataz[0, 0:1], dataz[1, 0:1], dataz[2, 0:1],'o')

This is displayed by specifying it with num sequentially with the following update function.

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)

This time, I can get the same picture in 2D, but I draw it in 3d plot because I can rotate the pendulum more realistically (by arranging it in the y direction) with a mouse. And the most important calculation of L was carried out below with reference to the above table. Using the Fc that can be calculated from now on, the length of each pendulum can be calculated by $ L = L_0 * F_c [i] ** 2 $.

L0=980*(60/65)**2/(2*np.pi)**2
Fc=[]
fc=1
Fc.append(fc)
for i in range(1,15,1):
    fc=fc*((66-i)/(65-i))
    Fc.append(fc)

Reproduction; 15 pendulums

0.2-30 sec 30.2-60 sec

Enlarged; 30 pendulums

I tried twice as many. 0.2-30 sec 30.2-60 sec

A little discussion

If it can be done, I think it will be a problem if you actually write the code from scratch. It's hard to get rid of the wonder that it can be arranged so beautifully. And, in fact, I tried to increase the number in order from 3. It is interesting to move in the same way when there are few. Of course, with this method, the same calculation can be performed up to the point where Fc can be calculated to obtain the above L. And if you want to increase it further, you can make a different variation by changing the part that aligns in 60 seconds. So, I tried 50 pieces. First, I made it easy to calculate dataz and linez with a for statement using a list.

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

The part to be displayed is also simplified as follows.

    for j in range(50):        
        update(num,dataz[j],linez[j],s0)

With this, in the case of 50 pieces, it will be 3 m as shown below. This makes me want to make it.

Summary

・ I tried to simulate a sine wave mass game ・ Simple but beautiful

・ I'm looking forward to various expansions ・ I want to move the real thing

bonus

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):
    phi = 0  #Initial time value
    g = 980  #Gravitational acceleration
    x0=2.5  #Initial amplitude
    omega = np.sqrt(g/L)
    theta0 = np.arcsin(x0/L)  #Initial angle
    while phi < 600:  #Pi for starting from maximum amplitude/Add 2
        yield np.array([L*np.sin(np.sin(phi*omega+np.pi/2)*theta0), 0,z0-L+L*(1-np.cos(np.sin(phi*omega+np.pi/2)*theta0))])
        phi += 1/n  #Divide 1sec into 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 = 360
z0 =50
L0=980*(60/65)**2/(2*np.pi)**2
Fc=[]
fc=1
Fc.append(fc)
for i in range(1,15,1):
    fc=fc*((66-i)/(65-i))
    Fc.append(fc)

dataz=[]
linez=[]
for j in range(50):
    dataz0 = np.array(list(genxy(N,L=L0*Fc[j]**2,z0=z0))).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([-10.0, 10.0])
ax.set_ylabel('Y')

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

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