[PYTHON] [Einführung in die Simulation] Sign Wave Mass Game ♬

Ich war auch süchtig nach dieser Zeit. Das folgende Pendel, das auf Twitter veröffentlicht wurde. Wonder of Science @ Wonderofscience

Da jeder unabhängig ist, sieht es aus wie eine Gruppenbewegung. Darüber hinaus erinnert mich die Bewegung an den rekursiven Satz von Poancare. "Das mechanische System kehrt innerhalb einer begrenzten Zeit fast in seinen willkürlichen Ausgangszustand zurück, wenn bestimmte Bedingungen erfüllt sind." Nun, diesmal ist es ein einfaches Pendel, also ist es natürlich, wenn Sie den Zyklus steuern. .. .. Ich hatte jedoch intuitiv das Gefühl, dass es nicht so einfach war.

Ich habe auch die folgenden Schnappschüsse gemacht. Es ist ersichtlich, dass der kürzeste Pendelzyklus ungefähr 1 Sekunde beträgt und der längste Pendelzyklus ungefähr das 4/3-fache des Zyklus zu sein scheint. Allein mit diesen Informationen war jedoch unklar, wie der Pendelzyklus zwischen ihnen so sauber ausgerichtet werden sollte. Nachdem ich gekämpft und verschiedene Dinge getan hatte, las ich das Originalmaterial als Referenz. Es gab von Anfang an einen Link zur ursprünglichen Geschichte, aber ich habe sie erst zu diesem Zeitpunkt nicht gelesen. [Referenz] Originalmaterial Pendulum Waves@Harvard Natural Sciences Lecture Demonstrations Spoiler. "Die Dauer eines vollständigen Tanzzyklus beträgt 60 Sekunden. Die längste Pendellänge wird so eingestellt, dass sie in diesen 60 Sekunden 51-mal vibriert. Die ununterbrochene kurze Pendellänge ist in dieser Zeit weiter. Sorgfältig abgestimmt, um eine Vibration auszuführen, vibriert das 15. Pendel (kürzeste) 65-mal. Sobald alle 15 Pendel zusammen gestartet wurden, sind sie bald nicht mehr synchron. Aufgrund der unterschiedlichen Zyklen der relativen Vibration ändert sich die relative Phase kontinuierlich. Nach 60 Sekunden führen sie jedoch alle ganzzahlige Vielfache der Vibration aus und sind in diesem Moment bereit, den Tanz neu zu synchronisieren und zu wiederholen. . "

Natürlich wäre es ziemlich schwierig, dies mit einem echten Pendel als Referenz zu realisieren. Ich werde es hier jedoch simulieren, also werde ich mehr versuchen. 【Referenz】 [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? b490-4864-923d-51639f6a935f & pf_rd_r = VVJHE2FGVTDYCM3J5C6K & psc = 1 & refRID = VVJHE2FGVTDYCM3J5C6K)

Überprüfung des einfachen Pendels

【Referenz】 Kommentar: Einfache Pendelbewegung Eine numerische Lösung der Differentialgleichung ist in Ordnung, aber ich habe mich entschieden, hier die folgende analytische Lösung zu verwenden.

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

Die obige Tabelle zeigt Folgendes

Codekommentar

Ich habe den ganzen Code als Bonus gesetzt. Es ist fast das gleiche wie die Axes 3D-Animation neulich. Die wesentlichen Teile sind wie folgt.

def genxy(n,L=20,z0=0):
    phi = 0  #Anfangszeitwert
    g = 980  #Schwerkraftbeschleunigung
    x0=2.5  #Anfangsamplitude
    omega = np.sqrt(g/L)
    theta0 = np.arcsin(x0/L)  #Anfangswinkel
    while phi < 600:  #Pi zum Starten ab maximaler Amplitude/Addiere 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  #Teilen Sie 1 Sek. In n

Rufen Sie dies mit dem folgenden Code auf und verwenden Sie ihn Wenn N = 60, ist L die Länge jedes Pendels. z0 ist die Höhe des Drehpunkts. Berechnen Sie die Daten und zeichnen Sie sie in Zeilen.

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

Dies wird angezeigt, indem Sie es mit der folgenden Aktualisierungsfunktion nacheinander mit num angeben.

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)

Dieses Mal kann das gleiche Bild in zwei Dimensionen erhalten werden, aber da das Pendel mit der Maus realistischer gedreht werden kann (indem es in y-Richtung angeordnet wird), zeichne ich es mit einem 3D-Diagramm. Und die wichtigste Berechnung von L wurde unten unter Bezugnahme auf die obige Tabelle durchgeführt. Mit Fc, das von nun an berechnet werden kann, kann die Länge jedes Pendels mit $ L = L_0 * F_c [i] ** 2 $ berechnet werden.

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)

Fortpflanzung, 15 Pendel

0.2-30 sec 30.2-60 sec

Vergrößern, 30 Pendel

Ich habe doppelt so viele versucht. 0.2-30 sec 30.2-60 sec

Eine kleine Diskussion

Wenn es möglich ist, fürchte ich, aber ich denke, dass es schwierig sein wird, Code von Grund auf neu zu schreiben. Es ist schwer, das Wunder loszuwerden, dass es so schön arrangiert werden kann. Und tatsächlich habe ich versucht, die Anzahl in der Reihenfolge von 3 zu erhöhen. Es ist interessant, sich auf die gleiche Weise zu bewegen, wenn es nur wenige gibt. Natürlich können Sie mit dieser Methode auf die gleiche Weise bis zu dem Punkt berechnen, an dem Sie Fc berechnen können, um das obige L zu finden. Und wenn Sie es weiter erhöhen möchten, können Sie eine andere Variation vornehmen, indem Sie das Teil ändern, das in 60 Sekunden ausgerichtet wird. Also habe ich 50 Stück ausprobiert. Zuerst habe ich es einfach gemacht, dataz und linez mit einer for-Anweisung unter Verwendung einer Liste zu berechnen.

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)

Der anzuzeigende Teil wird ebenfalls wie folgt vereinfacht.

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

Bei 50 Stück sind es 3 m, wie unten gezeigt. Das bringt mich dazu, es zu schaffen.

Zusammenfassung

・ Ich habe versucht, ein Sinuswellen-Massenspiel zu simulieren ・ Einfach aber schön

・ Ich freue mich auf verschiedene Erweiterungen ・ Ich möchte die reale Sache bewegen

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  #Anfangszeitwert
    g = 980  #Schwerkraftbeschleunigung
    x0=2.5  #Anfangsamplitude
    omega = np.sqrt(g/L)
    theta0 = np.arcsin(x0/L)  #Anfangswinkel
    while phi < 600:  #Pi zum Starten ab maximaler Amplitude/Addiere 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  #Teilen Sie 1 Sek. In 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)