通过迭代方法估计 spring 上的速度

Question

问题：

考虑一个具有质量和 spring 的系统，如图所示 below. The stiffness of the spring and the mass of the object are known. Therefore, if the spring is stretched the force the spring exerts can be calculated from Hooke`s law and the instantaneous acceleration can be estimated from Newton´s laws of motion. Integrating the acceleration twice yields the distance the spring would move and subtracting that from the initial length results in a new position to calculate the acceleration and start the loop again. Therefore as the acceleration decreases linearly the speed levels off at a certain value (top right) 之后的所有内容，spring 在这种情况下忽略压缩和减速。

我的问题是如何在 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

如果您要预先积分 - 这是一个好主意，并且绝对是可行的方法 - 那么您可以将方程式写成 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 中涉及的最高值。

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

可能用于 ODE 的最佳标准数值积分方法（假设您不需要像 leapfrog 方法这样的方法来提高稳定性）是 Runge Kutta 方法。阶数充足的自适应时间步长 Runge Kutta 方法通常应该可以解决问题并为您提供准确的答案。不幸的是，解释 Runge Kutta 方法的数学可能太高级和耗时，无法在此处涵盖，但您可以在线或在例如Numerical Recipes，一系列关于数值方法的书籍，其中包含许多非常有用的代码示例。

不过，即使是 Runge Kutta 方法，基本上也是通过改进对时间量程内函数值的猜测来工作的。他们只是以更复杂的方式做到这一点，这证明可以减少每一步的错误。

Answer 2

分析方法是获得遵守胡克定律的简单系统速度的最简单方法。

但是，如果您想要物理上准确的 numerical/iterative 方法，我强烈建议不要使用标准 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

我提供了一个详细的示例（尽管是针对另一个二维系统），包括此处的一阶和四阶方法的实现： https://zymplectic.com/case3.html 方法 0 和方法 1

这些是辛积分器，它们具有能量守恒属性，这意味着您可以执行长时间模拟而不会出现耗散误差。

Answer 3

你的力有符号错误，对于 spring 或任何其他振荡，它应该始终与激励方向相反。纠正这一点会立即产生振荡。但是，您的循环条件现在将永远不会得到满足，因此您也必须适应它。

您可以通过将方法从当前的辛 Euler 方法提升到 Leapfrog-Verlet 来立即增加方法的阶数。您只需将 v[i] 的解释更改为 t[i]-dt/2 处的速度。然后第一次更新使用 t[i-1] 处的中间加速度，使用中点公式从 t[i-1]-dt/2 处的速度计算 t[i-1]+dt/2=t[i]-dt/2 处的速度。然后在下一行中，位置更新是一个类似的中点公式，使用位置时间中间时间的速度。要获得此优势，您只需更改代码即可使用 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)，使用它和它的导数来计算数值解中的误差（记住速度的跳跃）给出 via

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。

通过迭代方法估计 spring 上的速度

Estimate velocity on a spring by iterative approach

python

math

physics

numerical-methods

differential-equations