我是 python 新手,我正在尝试实现一个适用于谐波振荡器的速度 Verlet 积分器。正如您从下面我的笔记本中看到的(取自:http ://hplgit.github.io/prog4comp/doc/pub/._p4c-solarized-Python022.html ),欧拉的方法有效,但verlet 积分器无效。我错过了什么?
## Euler's method
from numpy import zeros, linspace, pi, cos, array
import matplotlib.pyplot as plt
import numpy as np
omega = 1
m=1
N=500
dt=0.08
t = linspace(0, N*dt, N+1)
u = zeros(N+1)
vi = zeros(N+1)
X_0 = 2
u[0] = X_0
vi[0] = 0
for n in range(N):
u[n+1] = u[n] + dt*vi[n]
vi[n+1] = vi[n] - dt*omega**2*u[n]
fig = plt.figure()
l1, l2 = plt.plot(t, u, 'b-', t, X_0*cos(omega*t), 'r--')
fig.legend((l1, l2), ('numerical', 'exact'), 'upper left')
plt.xlabel('t')
plt.show()
## Velocity Verlet
for n in range(N):
v_next = vi[n] - 0.5*omega**2*u[n]*dt
u[n+1]= u[n]+ dt*v_next
vi[n+1] = v_next - 0.5*omega**2*u[n]*dt
np.append(u, x_next)
plt.plot(t, X_0*cos(omega*t), label = 'Analytical')
plt.plot(t, u, 'r--', label = 'Verlet')
plt.title('Verlet', fontweight = 'bold', fontsize = 1)
plt.xlabel('t', fontweight = 'bold', fontsize = 14)
plt.ylabel('X', fontweight = 'bold', fontsize = 14)
plt.legend()
plt.show()