[Wissenschaftlich-technische Berechnung nach Python] Numerische Lösung des Problems des eindimensionalen harmonischen Oszillators nach der Speed-Berle-Methode

[Verlet-Methode](https://ja.wikipedia.org/wiki/%E3%83%99%E3%83%AC%E3%81%AE%E6%96%B9%E6%B3% Gemäß 95) ist die Bewegungsgleichung des eindimensionalen harmonischen Oszillators $ M d ^ 2x / dt ^ 2 = F (x, t) = -kx $ unter der Anfangsbedingung $ x = 0 $, $ dx / dt = 1,0 $ Löse mit. Die Masse des Qualitätspunktes sei M = 1 kg und die Federkonstante $ k = 1,0 $ N / m.

Bei der Geschwindigkeits-Berle-Methode werden die Koordinaten x und die Geschwindigkeit v ($ = dx / dt $) wie folgt aktualisiert [3].

$x_{n+1} = x_n +v_n \ \delta t +0.5 \ F(x_n,t_n)\ (\delta t)^2/M $ v_{n+1} = v_n+ 0.5 \ \delta t\ (\ F(x_n, t_n) + F(x_{n+1}, t_{n+1}))/M

Es zeichnet sich durch eine höhere Gesamtenergieeinsparung (mit Sympraktizität) aus als die Lunge-Kutta-Methode [1,2]. Darüber hinaus entspricht die Geschwindigkeits-Berle-Methode mathematisch der normalen Berle-Methode, ist jedoch resistent gegen Rundungsfehler und ein numerisch robusterer Ausdruck. Aus diesen Gründen wird die Geschwindigkeits-Berle-Methode häufig in molekulardynamischen Simulationen verwendet.

import numpy as np
import matplotlib.pyplot as plt

"""
Lösung nach der Speed-Berley-Methode
"""



def f(x, t): # 
    return -k*x

M = 1.0 #Masse Masse[Kg]
k=1.0  #Federkonstante Federkonstante[N/m]

t0 = 0.0 #Anfangswert der Simulationszeit
t1 = 20.0 #Maximale Simulationszeit
N = 200 #Anzahl der Zeitinkremente von t1 bis t0: 
del_t = (t1-t0)/N #Zeitschritt Zeitschritt

tpoints = np.arange(t0, t1, del_t) #del den Zeitpunkt von t0 bis t1_In t-Schritten generiert
xpoints = [] #Liste zum zeitlichen Speichern von x-Koordinaten
vpoints = []#Liste zum Speichern von Geschwindigkeits-v-Koordinaten pro Stunde
Etot=[] #Liste zur stündlichen Speicherung mechanischer Energie

# initial condition
x0 = 0.0 #Ausgangszustand von x
v0 = 1.0 #Ausgangszustand der Geschwindigkeit v

x, v = x0, v0


for t in tpoints:
    xpoints.append(x)
    vpoints.append(v)
    Etot.append((M*v**2)/2+(k*x**2)/2) #Mechanische Energie(mv^2/2 + k x^2/2)

    xx=x  #Speicherung des vorherigen Schritts von x
    x +=  v*del_t + 0.5 * f(x,t)*(del_t**2) #Koordinatenaktualisierung nach Geschwindigkeit Berley-Methode
    v += 0.5*del_t*(f(xx , t) + f(x, t+del_t)) #Geschwindigkeitsaktualisierung nach der Geschwindigkeits-Berley-Methode
     
# for plot  #Genaue Lösung x= sin(t)Im Vergleich mit
plt.plot (tpoints, xpoints, 'o',label='Velocity Verlet')
plt.xlabel("t",  fontsize=24)
plt.ylabel("x(t)",  fontsize=24)

tt=np.arange(t0, t1, 0.01*del_t)
plt.plot (tt, np.sin(tt), '-',label='Exact: Sin (t)',color = 'Green')
plt.legend(loc='upper left')

plt.show()

##Genaue Lösung v= cos(t)Im Vergleich mit
plt.plot (tpoints,vpoints, 'o',label='Velocity Verlet')
plt.xlabel("t", fontsize=24)
plt.ylabel("v(t)", fontsize=24)

tt=np.arange(t0, t1, 0.01*del_t)
plt.plot (tt, np.cos(tt), '-',label='Exact: Cos (t)', color = 'red')
plt.legend(loc='upper left')

plt.show()


#Zeichnung von voller Energie Etot. Die genaue Lösung ist Etot=0.5
plt.plot (tpoints, Etot, '-',label='Velocity Verlet')
plt.xlabel("t",  fontsize=24)
plt.ylabel("Etot(t)",  fontsize=24)


plt.legend(loc='upper left')

plt.show()

Nachtrag: Vor- und Nachteile der einfachen Methode [4]

** Vorteile **: Hochpräzise Berechnung der Dynamik des konservativen Systems, der Zeitumkehrsymmetrie (der Newton-Gleichung), der langfristigen mechanischen Energieeinsparbarkeit usw.

** Nachteile **: Unbequem zum Ändern des Zeitschritts (automatische Genauigkeitskontrolle ist nicht möglich), Funktionen werden in Nichtkonservierungssystemen nicht verwendet (Methode mit konstanter Temperatur / konstantem Druck, wenn eine von der Geschwindigkeit abhängige Kraft austritt usw.).

Verweise

[1] Rihei Endo, "Lösen gemeinsamer Differentialgleichungen"

[2] Satoshi Yinyama, "Über die Symplectic Integration Method"

[3] Im Internet ist hier leicht zu verstehen. In dem Buch zum Beispiel von Harvey [Einführung in die Computerphysik](https://www.amazon.co.jp/%E8%A8%88%E7%AE%97%E7%89%A9%E7%90%86] % E5% AD% A6% E5% 85% A5% E9% 96% 80-% E3% 83% 8F% E3% 83% BC% E3% 83% 99% E3% 82% A4-% E3% 82% B4 % E3% 83% BC% E3% 83% AB% E3% 83% 89 / dp / 4894713187) und Kunin Computerphysik % E7% AE% 97% E6% A9% 9F% E7% 89% A9% E7% 90% 86% E5% AD% A6-Steven-Koonin / dp / 4320032918) usw. haben eine elementare und leicht verständliche Erklärung.

[4] Shuichi Nose, "Molekulardynamiksimulation / numerische Integrationsmethode"