<p>力中存在符号错误,对于弹簧或任何其他振动,它应始终与激励方向相反。纠正这一点会立即产生振荡。然而,您的循环条件现在将永远不会得到满足,因此您还必须进行调整</p>
<p>通过将方法从当前的辛欧拉方法提升到蛙跳Verlet,可以立即提高方法的阶数。您只需将<code>v[i]</code>的解释更改为<code>t[i]-dt/2</code>处的速度。然后,第一次更新使用^中的加速度{{CD3}},用中点公式从^ {< CD5> }的速度计算{{CD4}}的速度。然后在下一行中,位置更新是一个类似的中点公式,使用位置时间之间中间时间的速度。为了获得这一优势,您只需在代码中进行更改,即使用<code>t[0]</code>处的泰勒展开将初始速度设置为<code>t[0]-dt/2</code>时的速度</p>
<pre class="lang-py prettyprint-override"><code>instantaneous_acceleration = lambda x: -5*x/10 # a = kx/m
delta_time = 0.01 #10 milliseconds
s0, v0 = 12, 0 #initial length 12 m, initial velocity 0 m/s
N=1000
s = np.zeros(N+1); v = s.copy(); a = s.copy()
a[0] = instantaneous_acceleration(s0) #initial acceleration when stretched to 12 m
v[0] = v0-a[0]*delta_time/2
s[0] = s0
for i in range(N):
v[i+1] = a[i]*delta_time + v[i] #calculate the next velocity
s[i+1] = v[i+1]*delta_time + s[i] #calculate the next position
a[i+1] = instantaneous_acceleration (s[i+1]) #use the position to derive the new acceleration
#produce plots of all these functions
t=np.arange(0,N+1)*delta_time;
fig, ax = plt.subplots(3,1,figsize=(5,3*1.5))
for g, y in zip(ax,(s,v,a)):
g.plot(t,y); g.grid();
plt.tight_layout(); plt.show();
</code></pre>
<p><a href="https://i.stack.imgur.com/JKMMn.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/JKMMn.png" alt="enter image description here"/></a></p>
<p>这显然是一种正确的振荡。精确解是<code>12*cos(sqrt(0.5)*t)</code>,使用它及其导数计算数值解中的误差(记住速度的跳跃)</p>
<pre><code>w=0.5**0.5; dt=delta_time;
fig, ax = plt.subplots(3,1,figsize=(5,3*1.5))
for g, y in zip(ax,(s-12*np.cos(w*t),v+12*w*np.sin(w*(t-dt/2)),a+12*w**2*np.cos(w*t))):
g.plot(t,y); g.grid();
plt.tight_layout(); plt.show();
</code></pre>
<p>下图显示了预期大小<code>delta_time**2</code>中的错误。
<a href="https://i.stack.imgur.com/tKZad.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/tKZad.png" alt="enter image description here"/></a></p>