[Calcul scientifique / technique par Python] Solution numérique d'un problème d'oscillateur harmonique unidimensionnel par vitesse Méthode de Berle

[Méthode Speed Berlet](https://ja.wikipedia.org/wiki/%E3%83%99%E3%83%AC%E3%81%AE%E6%96%B9%E6%B3% Selon 95), l'équation de mouvement de l'oscillateur harmonique à une dimension, $ M d ^ 2x / dt ^ 2 = F (x, t) = -kx $ sous la condition initiale $ x = 0 $, $ dx / dt = 1.0 $ Résolvez avec. Soit la masse du point de qualité M = 1 kg et la constante du ressort $ k = 1,0 $ N / m.

Dans la méthode de Berle de vitesse, les coordonnées x et la vitesse v ($ = dx / dt $) sont mises à jour comme suit [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

Elle se caractérise par une meilleure conservation de l'énergie totale que la méthode Lunge-Kutta [1,2]. De plus, la méthode Speed Berle est mathématiquement équivalente à la méthode Berle normale, mais elle résiste aux erreurs d'arrondi et est une expression numériquement plus robuste. Pour ces raisons, la méthode de vitesse de Berle est souvent utilisée dans les simulations de dynamique moléculaire.

import numpy as np
import matplotlib.pyplot as plt

"""
Solution par vitesse Méthode Berley
"""



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

M = 1.0 #Masse de masse[Kg]
k=1.0  #Constante du ressort Constante du ressort[N/m]

t0 = 0.0 #Valeur initiale du temps de simulation
t1 = 20.0 #Temps de simulation maximum
N = 200 #Nombre d'incréments de temps de t1 à t0: 
del_t = (t1-t0)/N #Pas de temps pas de temps

tpoints = np.arange(t0, t1, del_t) #del le point de temps de t0 à t1_Généré par incréments de t
xpoints = [] #Liste de stockage des coordonnées x dans le temps
vpoints = []#Liste de stockage des coordonnées de vitesse v par heure
Etot=[] #Liste de stockage horaire de l'énergie mécanique

# initial condition
x0 = 0.0 #Condition initiale de x
v0 = 1.0 #Condition initiale de vitesse 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) #Énergie mécanique(mv^2/2 + k x^2/2)

    xx=x  #Stockage de l'étape précédente de x
    x +=  v*del_t + 0.5 * f(x,t)*(del_t**2) #Mise à jour des coordonnées par la méthode Speed Berley
    v += 0.5*del_t*(f(xx , t) + f(x, t+del_t)) #Mise à jour de la vitesse par la méthode Speed Berley
     
# for plot  #Solution exacte x= sin(t)Comparaison avec
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()

##Solution exacte v= cos(t)Comparaison avec
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()


#Dessin de pleine énergie Etot. La solution exacte est 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()

Addendum: avantages et inconvénients de la méthode simplifiée [4]

** Avantages **: Calcul haute précision de la dynamique du système conservateur, symétrie d'inversion temporelle (de l'équation de Newton), conservabilité de l'énergie mécanique à long terme, etc.

** Inconvénients **: Il n'est pas pratique de changer le pas de temps (le contrôle automatique de la précision n'est pas possible), et les fonctionnalités ne sont pas utilisées dans les systèmes sans conservation (méthode température constante / pression constante, lorsque la force dépendant de la vitesse sort, etc.

Les références

[1] Rihei Endo, «Solving Common Differential Equations»

[2] Satoshi Yinyama, "À propos de la méthode d'intégration symplectique"

[3] Sur le Web, ici est facile à comprendre. Dans le livre, par exemple, de Harvey [Introduction to Computational Physics](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) et Kunin Computer Physics % E7% AE% 97% E6% A9% 9F% E7% 89% A9% E7% 90% 86% E5% AD% A6-Steven-Koonin / dp / 4320032918) etc. ont une explication élémentaire et facile à comprendre.

[4] Shuichi Nose, "Simulation de dynamique moléculaire / méthode d'intégration numérique"