[Scientific / technical calculation by Python] Numerical solution of one-dimensional harmonic oscillator problem by velocity Verlet method

[Speed Verlet method](https://ja.wikipedia.org/wiki/%E3%83%99%E3%83%AC%E3%81%AE%E6%96%B9%E6%B3% According to 95), the equation of motion of the one-dimensional harmonic oscillator, $ M d ^ 2x / dt ^ 2 = F (x, t) =-kx $, is set under the initial condition $ x = 0 $, $ dx / dt = 1.0 $. Solve with. Let the mass of the mass point be M = 1 kg and the spring constant $ k = 1.0 $ N / m.

In the velocity Verlet method, the coordinates x and velocity v ($ = dx / dt $) are updated as follows [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

It is characterized by higher total energy conservation (having symplectic property) than Runge-Kutta method [1,2]. In addition, the velocity Verlet method is mathematically equivalent to the normal Verlet method, but it is resistant to rounding errors and has a numerically more robust expression. For these reasons, the velocity Verlet method is often used in molecular dynamics simulations.

import numpy as np
import matplotlib.pyplot as plt

"""
Solving by the velocity Berley method
"""



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

M = 1.0 #Mass mass[Kg]
k=1.0  #Spring constant spring constant[N/m]

t0 = 0.0 #Initial value of simulation time
t1 = 20.0 #Maximum simulation time
N = 200 #Number of time increments from t1 to t0: 
del_t = (t1-t0)/N #Time step time step

tpoints = np.arange(t0, t1, del_t) #del the time point from t0 to t1_Generated in t increments
xpoints = [] #List for storing x-coordinates in time
vpoints = []#List for storing velocity v coordinates by hour
Etot=[] #List for storing mechanical energy by hour

# initial condition
x0 = 0.0 #Initial condition of x
v0 = 1.0 #Initial condition of velocity 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) #Mechanical energy(mv^2/2 + k x^2/2)

    xx=x  #Storage of the previous step of x
    x +=  v*del_t + 0.5 * f(x,t)*(del_t**2) #Coordinate update by velocity Berley method
    v += 0.5*del_t*(f(xx , t) + f(x, t+del_t)) #Speed update by the speed Berley method
     
# for plot  #Exact solution x= sin(t)Comparison with
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()

##Exact solution v= cos(t)Comparison with
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()


#Drawing of full energy Etot. The exact solution is 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: Advantages and Disadvantages of Symplectic Method [4]

** Advantages **: High-precision calculation of the dynamics of conservative systems is possible, time-reversal symmetry (as Newton's equation has), long-term mechanical energy conservation, etc.

** Disadvantages **: Inconvenient to change the time step (automatic accuracy control is not possible), features are not utilized in non-preservation systems (constant temperature / constant pressure method, when a force that depends on speed comes out, etc.).

References

[1] Rihei Endo, "Solving Ordinary Differential Equations"

[2] Satoshi Yinyama, "About Symplectic Integral Method"

[3] On the web, here is easy to understand. In the book, for example, by 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) and by Kunin Computer Physics % E7% AE% 97% E6% A9% 9F% E7% 89% A9% E7% 90% 86% E5% AD% A6-Steven-Koonin / dp / 4320032918) etc. have an elementary and easy-to-understand explanation.

[4] Shuichi Nosé, "Molecular dynamics simulation / numerical integration method"