用迭代法估计弹簧上的速度

Question

问题是：

考虑一个具有质量和弹簧的系统，如图{a1}所示。弹簧的刚度和物体的质量是已知的。因此，如果弹簧被拉伸，弹簧施加的力可以从Hooke`s law计算，瞬时加速度可以从牛顿运动定律估算。将加速度积分两次，得到弹簧移动的距离，并从初始长度中减去该距离，得到一个新位置，以计算加速度并再次开始循环。因此，当加速度线性下降时，速度在某一特定值处趋于平稳(top right)在该点之后，弹簧压缩&；在这种情况下，减速被忽略

我的问题是如何用python编写代码。到目前为止，我已经编写了一些伪代码

instantaneous_acceleration = lambda x: 5*x/10    # a = kx/m
delta_time = 0.01     #10 milliseconds
a[0] = instantaneous_acceleration(12)     #initial acceleration when stretched to 12 m
v[0] = 0        #initial velocity 0 m/s
s[0] = 12       #initial length 12 m
i = 1
while a[i] > 12:
    v[i] = a[i-1]*delta_time + v[i-1]      #calculate the next velocity
    s[i] = v[i]*delta_time + s[i-1]        #calculate the next position
    a[i] = instantaneous_acceleration (s[i])          #use the position to derive the new accleration
    i = i + 1

非常感谢任何帮助或提示

Answer 1

解析法是获得服从胡克定律的简单系统速度的最简单方法

然而，如果你想要一个物理上精确的数值/迭代方法，我强烈建议不要使用标准的Euler或runge-kutta方法（Patrick87建议）。[更正：如果加速度项的符号已更正，则OPs方法为辛一阶方法。]

您可能希望使用哈密顿方法和合适的辛积分器，如二阶蛙跳（Patrick87也建议）

对于胡克定律，可以表示哈密顿量H=T（p）+V（q），其中p是动量（与速度相关），q是位置（与弦离平衡的距离相关）

有动能T和势能V

T（p）=0.5*p^2/m

V（q）=0.5*k*q^2

您只需要这两个表达式的导数来模拟系统

dT/dp=p/m

dV/dq=k*q

我提供了一个详细的示例（尽管是针对另一个二维系统），其中包括一阶和四阶方法的实现： ^方法0和方法1下的{a1}

这些是辛积分器，它具有能量保持特性，这意味着您可以执行长时间模拟，而不会产生耗散误差

Answer 2

如果你要提前积分，这是个好主意，而且绝对是可行的方法，那么你可以把方程写下来，作为t的函数：

x'' = -kx/m
x'' + (k/m)x = 0
r^2 + k/m = 0
r^2 = -(k/m)
r = i*sqrt(k/m)
x(t) = A*e^(i*sqrt(k/m)t)
     = A*cos(sqrt(k/m)t + B) + i*A*sin(sqrt(k/m)t + B)
     = A*cos(sqrt(k/m)t + B)

从初始条件我们知道

x(0) = 12 = A*cos(B)
v(0) = 0 = -sqrt(k/m)*A*sin(B)

只有当我们选择A=0或B=0或B=Pi时，这些方程中的第二个才是真的

• 如果A=0，则第一个方程没有解
• 如果B=0，则第一个方程的解A=12
• 如果B=Pi，第一个方程的解A=-12

我们可能更喜欢B=0和A=12。这给

x(t) =  12*cos(sqrt(k/m)t)
v(t) = -12*sqrt(k/m)*sin(sqrt(k/m)t)
a(t) = -12*(k/m)cos(sqrt(k/m)t)

因此，在任何增量时间t[n+1]=t[n]+dt，我们可以简单地计算t[n]的精确位置、速度和加速度，而不会累积任何漂移或误差

综上所述，如果你对给定一个任意常微分方程的x（t）、v（t）和a（t）的数值计算感兴趣，那么答案就难多了。有很多很好的方法可以被称为数值积分。欧拉方法是最简单的：

// initial conditions

        t[0] = 0
        x[0] = …
       x'[0] = …
         …
  x^(n-1)[0] = …
    x^(n)[0] = 0


// iterative step

  x^(n)[k+1] = f(x^(n-1)[k], …, x'[k], x[k], t[k])
x^(n-1)[k+1] = x^(n-1)[k] + dt * x^(n)[k]
             …
     x'[k+1] = x'[k] + dt * x''[k]
      x[k+1] = x[k] + dt * x'[k]
      t[k+1] = t[k] + dt

选择的dt值越小，运行固定时间所需的时间就越长，但得到的结果越精确。这基本上是对函数及其所有导数进行黎曼和，直到ODE中涉及的最高导数

更准确的版本是辛普森规则，它做的是相同的事情，但取最后一个时间量的平均值（而不是任一端点的值；上面的示例使用间隔的开始）。保证区间内的平均值比任一端点更接近区间内的真值（除非函数在该区间内为常数，在这种情况下，辛普森至少也一样好）

对于常微分方程来说，最好的标准数值积分方法（假设不需要像蛙跳法这样的方法来提高稳定性）可能是龙格-库塔法。一个足够阶数的自适应时间步长龙格库塔方法通常可以做到这一点，并给出准确的答案。不幸的是，解释Runge-Kutta方法的数学可能过于先进和耗时，无法在这里介绍，但您可以在网上或数字配方（Numeric Recipes）中找到有关这些和其他先进技术的信息，这是一系列关于数字方法的书籍，其中包含许多非常有用的代码示例

即使是龙格-库塔方法也基本上是通过在时间量子内对函数值的精确猜测来工作的。他们只是用更复杂的方法来做，这可以证明减少了每一步的误差

Answer 3

力中存在符号错误，对于弹簧或任何其他振动，它应始终与激励方向相反。纠正这一点会立即产生振荡。然而，您的循环条件现在将永远不会得到满足，因此您还必须进行调整

通过将方法从当前的辛欧拉方法提升到蛙跳Verlet，可以立即提高方法的阶数。您只需将v[i]的解释更改为t[i]-dt/2处的速度。然后，第一次更新使用^中的加速度{{CD3}}，用中点公式从^ {< CD5> }的速度计算{{CD4}}的速度。然后在下一行中，位置更新是一个类似的中点公式，使用位置时间之间中间时间的速度。为了获得这一优势，您只需在代码中进行更改，即使用t[0]处的泰勒展开将初始速度设置为t[0]-dt/2时的速度

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

这显然是一种正确的振荡。精确解是12*cos(sqrt(0.5)*t)，使用它及其导数计算数值解中的误差（记住速度的跳跃）

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

下图显示了预期大小delta_time**2中的错误。