[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()

t1.png

t2.png

t3.png


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"

Recommended Posts

[Scientific / technical calculation by Python] Numerical solution of one-dimensional harmonic oscillator problem by velocity Verlet method
[Scientific / technical calculation by Python] Numerical solution of eigenvalue problem of matrix by power method, numerical linear algebra
[Scientific / technical calculation by Python] Numerical solution of second-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of one-dimensional and two-dimensional wave equations by FTCS method (explicit method), hyperbolic partial differential equations
[Scientific / technical calculation by Python] Sum calculation, numerical calculation
[Scientific / technical calculation by Python] Numerical solution of two-dimensional Laplace-Poisson equation by Jacobi method for electrostatic potential, elliptic partial differential equation, boundary value problem
[Scientific / technical calculation by Python] 2D random walk (drunk walk problem), numerical calculation
[Science / technical calculation by Python] Numerical solution of first-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (2), harmonic oscillator potential, quantum mechanics
[Scientific / technical calculation by Python] Lagrange interpolation, numerical calculation
[Scientific / technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation
[Scientific / technical calculation by Python] Numerical calculation to find the value of derivative (differential)
[Scientific / technical calculation by Python] Analytical solution to find the solution of equation sympy
[Scientific / technical calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method
[Scientific / technical calculation by Python] Basic operation of arrays, numpy
[Scientific / technical calculation by Python] Monte Carlo integration, numerical calculation, numpy
[Scientific / technical calculation by Python] List of matrices that appear in Hinpan in numerical linear algebra
[Scientific / technical calculation by Python] Numerical integration, trapezoidal / Simpson law, numerical calculation, scipy
[Scientific / technical calculation by Python] Solving simultaneous linear equations, numerical calculation, numpy
[Scientific / technical calculation by Python] Solving the Schrodinger equation in the steady state in the 3D isotropic harmonic oscillator potential by the matrix method, boundary value problem, quantum mechanics
[Scientific / technical calculation by Python] Monte Carlo simulation of thermodynamics of 2D Ising spin system by Metropolis method
[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy
[Scientific / technical calculation by Python] Calculation of matrix product by @ operator, python3.5 or later, numpy
[Scientific / technical calculation by Python] Inverse matrix calculation, numpy
[Scientific / technical calculation by Python] Histogram, visualization, matplotlib
[Scientific / technical calculation by Python] Solving (generalized) eigenvalue problem using numpy / scipy, using library
[Scientific / technical calculation by Python] Drawing animation of parabolic motion with locus, matplotlib
[Scientific / technical calculation by Python] Drawing of 3D curved surface, surface, wireframe, visualization, matplotlib
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (1), well-type potential, quantum mechanics
Numerical calculation of compressible fluid by finite volume method
[Scientific / technical calculation by Python] Logarithmic graph, visualization, matplotlib
[Scientific / technical calculation by Python] Polar graph, visualization, matplotlib
[Scientific / technical calculation by Python] Plot, visualization, matplotlib of 2D data read from file
[Scientific / technical calculation by Python] Drawing, visualization, matplotlib of 2D (color) contour lines, etc.
[Python] Numerical calculation of two-dimensional thermal diffusion equation by ADI method (alternate direction implicit method)
[Scientific / technical calculation by Python] 3rd order spline interpolation, scipy
[Numerical calculation method, python] Solving ordinary differential equations by Eular method
[Scientific / technical calculation by Python] List of usage of (special) functions used in physics by using scipy
[Scientific and technical calculation by Python] Drawing of fractal figures [Sierpinski triangle, Bernsley fern, fractal tree]
[Scientific / technical calculation by Python] Wave "beat" and group velocity, wave superposition, visualization, high school physics
Scientific / technical calculation by Python] Drawing and visualization of 3D isosurface and its cross section using mayavi
[Scientific / technical calculation by Python] 1D-3D discrete fast Fourier transform, scipy
[Science / technical calculation by Python] Comparison of convergence speeds of SOR method, Gauss-Seidel method, and Jacobi method for Laplace equation, partial differential equation, boundary value problem
[Scientific / technical calculation by Python] Solving ordinary differential equations, mathematical formulas, sympy
[Scientific / technical calculation by Python] Generation of non-uniform random numbers giving a given probability density function, Monte Carlo simulation
Calculation of technical indicators by TA-Lib and pandas