0

我是 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()
4

1 回答 1

0

您没有正确实现速度 Verlet 步骤,第二次速度更新使用新计算的位置,而不是旧位置,

vi[n+1] = v_next - 0.5*omega**2*u[n+1]*dt

这个小的变化应该会恢复二阶和能量/幅度保存。


还删除了先前更改的剩余部分,没有x_next附加内容。

在此处输入图像描述

于 2020-12-05T09:41:58.753 回答