[PYTHON] [Introduction à la simulation] Jeu de masse d'onde de signe ♬

J'étais aussi accro à cette époque. Le pendule suivant qui est sorti sur Twitter. Merveille de la science @ wonderofscience

Puisque chacun est indépendant, cela ressemble à un mouvement de groupe. De plus, le mouvement me rappelle ce théorème récursif de Poancare. "Le système mécanique revient presque à son état initial arbitraire dans un temps fini si certaines conditions sont remplies." Eh bien, cette fois, c'est un simple pendule, donc c'est naturel si vous contrôlez le cycle. .. .. Cependant, j'ai intuitivement senti que ce n'était pas si facile.

J'ai également pris les clichés suivants. On peut voir que le cycle de pendule le plus court est d'environ 1 seconde et que le cycle de pendule le plus long semble être d'environ 4/3 fois le cycle. Cependant, avec cette seule information, il n'était pas clair comment aligner le cycle du pendule entre eux si proprement. Après avoir lutté et fait diverses choses, j'ai lu le matériel original pour référence. Il y avait un lien vers l'histoire originale depuis le début, mais je ne l'ai pas lu seulement pour le moment. [Référence] Matériel original Pendulum Waves@Harvard Natural Sciences Lecture Demonstrations divulgacher. "La durée d'un cycle complet de danse est de 60 secondes. La plus longue longueur du pendule est ajustée pour vibrer 51 fois pendant ces 60 secondes. La courte longueur continue du pendule est encore plus pendant cette période. Soigneusement réglé pour effectuer une vibration, le 15e pendule (le plus court) vibre 65 fois.Une fois que les 15 pendules sont démarrés ensemble, ils seront bientôt désynchronisés. En raison des différents cycles de vibration relative, la phase relative change continuellement, mais après 60 secondes, ils exécutent tous des multiples entiers de la vibration et sont prêts à se resynchroniser et à répéter la danse à ce moment. . "

Bien sûr, il serait assez difficile de réaliser cela avec un vrai pendule comme référence. Cependant, je vais le simuler ici, donc je vais essayer plus. 【référence】 [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-02b-4c b490-4864-923d-51639f6a935f & pf_rd_r = VVJHE2FGVTDYCM3J5C6K & psc = 1 & refRID = VVJHE2FGVTDYCM3J5C6K)

Examen du pendule simple

【référence】 Commentaire: Mouvement à pendule simple Une solution numérique de l'équation différentielle est bien, mais j'ai décidé d'utiliser ici la solution analytique suivante.

\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.À 82 cm\\
T=1.000sec

Le tableau ci-dessus montre ce qui suit

Commentaire de code

Je mets tout le code en bonus. C'est presque la même chose que l'animation Axes 3D de l'autre jour. Les parties essentielles sont les suivantes.

def genxy(n,L=20,z0=0):
    phi = 0  #Valeur de temps initiale
    g = 980  #Accélération gravitationnelle
    x0=2.5  #Amplitude initiale
    omega = np.sqrt(g/L)
    theta0 = np.arcsin(x0/L)  #Angle initial
    while phi < 600:  #Pi pour partir de l'amplitude maximale/Ajouter 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  #Diviser 1sec en n

Appelez-le avec le code suivant et utilisez-le Où N = 60, L est la longueur de chaque pendule. z0 est la hauteur du point d'appui. Calculez dataz et tracez-le en 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')

Ceci est affiché en le spécifiant avec num de manière séquentielle avec la fonction de mise à jour suivante.

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)

Cette fois, la même image peut être obtenue en deux dimensions, mais comme le pendule peut être tourné de manière plus réaliste (en le disposant dans la direction y) avec la souris, je le dessine avec un tracé 3D. Et le calcul le plus important de L a été effectué ci-dessous en référence au tableau ci-dessus. En utilisant Fc qui peut être calculé à partir de maintenant, la longueur de chaque pendule peut être calculée par $ 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 pendules

0.2-30 sec 30.2-60 sec

Agrandir; 30 pendules

J'ai essayé deux fois plus. 0.2-30 sec 30.2-60 sec

Une petite discussion

Si cela peut être fait, j'en ai peur, mais je pense qu'il sera difficile d'écrire du code à partir de zéro. Il est difficile de se débarrasser de la merveille que cela puisse être arrangé si bien. Et, en fait, j'ai essayé d'augmenter le nombre de 3. Il est intéressant de se déplacer de la même manière quand il y en a peu. Bien sûr, avec cette méthode, vous pouvez calculer de la même manière jusqu'au point où vous pouvez calculer Fc pour trouver le L. ci-dessus. Et si vous souhaitez l'augmenter davantage, vous pouvez faire une variation différente en changeant la partie qui s'aligne en 60 secondes. Alors, j'ai essayé 50 pièces. Tout d'abord, j'ai facilité le calcul de dataz et linez avec une instruction for en utilisant une liste.

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)

La partie à afficher est également simplifiée comme suit.

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

Avec cela, dans le cas de 50 pièces, ce sera 3 m comme indiqué ci-dessous. Cela me donne envie de réussir.

Résumé

・ J'ai essayé de simuler un jeu de masse sinusoïdale ・ Simple mais beau

・ J'attends avec impatience diverses extensions ・ Je veux déplacer la vraie chose

prime

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  #Valeur de temps initiale
    g = 980  #Accélération gravitationnelle
    x0=2.5  #Amplitude initiale
    omega = np.sqrt(g/L)
    theta0 = np.arcsin(x0/L)  #Angle initial
    while phi < 600:  #Pi pour partir de l'amplitude maximale/Ajouter 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  #Diviser 1sec en 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)