python - 通过迭代方法估计 Spring 的速度

Question

问题:

考虑一个具有质量和 Spring 的系统，如图 below 所示。 。 Spring 的刚度和物体的质量是已知的。因此，如果 Spring 被拉伸(stretch)，则 Spring 施加的力可以根据 Hooke`s law 计算出来。瞬时加速度可以根据牛顿运动定律估算。对加速度进行两次积分即可得到 Spring 移动的距离，并从初始长度中减去该距离，得到一个新位置来计算加速度并再次开始循环。因此，随着加速度线性减小，速度会稳定在某个值 (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 中涉及的最高导数。

辛普森规则的更准确版本，执行相同的操作，但取最后一个时间段的平均值(而不是任一端点的值；上面的示例使用间隔的开始)。保证区间内的平均值比任一端点更接近区间内的真实值(除非函数在该区间内保持恒定，在这种情况下，辛普森至少同样好)。

可能 ODE 的最佳标准数值积分方法(假设您不需要诸如蛙跳方法之类的方法来获得更高的稳定性)是龙格库塔方法。具有足够阶数的自适应时间步龙格库塔方法通常应该可以解决问题并为您提供准确的答案。不幸的是，解释龙格库塔方法的数学可能太先进且耗时，无法在此处介绍，但您可以在线或在例如 中找到有关这些和其他高级技术的信息。 Numerical Recipes，一系列关于数值方法的书籍，其中包含许多非常有用的代码示例。

不过，即使是龙格库塔方法也基本上是通过改进对函数在时间段内的值的猜测来工作的。他们只是以更复杂的方式来做这件事，这证明可以减少每一步的错误。