计算机编程/物理/物体在空间中的运动（分段近似）

<维基文库：源代码

关于 加速物体位置函数 的问题是，尽管它在理论上是正确的，但它在描述物体在空间中的运动方面并不具有很强的适应性。为了纠正这一点，我们采用分段计算的方式。

我们回顾一下 泰勒级数形式的位置函数

$\mathbf {s} (t+t_{0})=\sum {\frac {t^{n}}{n!}}\mathbf {s} ^{(n)}(t_{0})$.

如果我们将时钟同步，使 $t_{0}$ 为零，则可以将其简化为 麦克劳林级数

$\mathbf {s} (t)=\sum {\frac {t^{n}}{n!}}\mathbf {s} _{0}^{(n)}$.

问题的解决方案就在这里。该方程本身描述了基于初始常数因子 $\mathbf {s} _{0}^{(n)}$ 的完整路径。但是，如果我们分段处理它，也就是说，我们计算 $t$ 等于某个小的增量时间 $\Delta t$ 的数值解，并在每次计算循环后将 $t_{0}$ 重新同步到 $t_{previous}$，那么我们就可以使用以下函数的增量形式，通过连续计算来动态构建物体路径。

$\mathbf {s} (t)=\mathbf {s} (\Delta t+t_{previous})=\sum {\frac {\Delta t^{n}}{n!}}\mathbf {s} _{previous}^{(n)}$.

注意，随着 $\Delta t\rightarrow 0$ 以及我们使用的项数越多，精度就越高。

根据初始常数因子$\mathbf {s} _{0}^{(n)}$计算出物体的初始位置变化后，我们只需提供影响物体的任何外部加速度，$\mathbf {s} ^{(2)}$，然后我们可以计算出物体的下一个位置$\mathbf {s} ^{(0)}$和速度$\mathbf {s} ^{(1)}$，以及下一组常数因子$\mathbf {s} ^{(n)}$，其中$\mathbf {s} ^{(0)}$，$\mathbf {s} ^{(1)}$和$\mathbf {s} ^{(2)}$都是已知的。

速度由以下公式计算：

$\mathbf {v} (t)=\sum {\frac {\Delta t^{n}}{n!}}\mathbf {v} _{previous}^{(n)}$

或

$\mathbf {s} ^{(1)}(t)=\sum {\frac {\Delta t^{n}}{n!}}\mathbf {s} _{previous}^{(n+1)}$.

其他$\mathbf {s} ^{(n)}(t)$值可以近似为：

$\mathbf {s} ^{(n)}(t)={\frac {\mathbf {s} ^{(n-1)}(t)-\mathbf {s} _{previous}^{(n-1)}}{\Delta t}}$ 当 $n>2$ 时，

前提是 $\Delta t$ 非常小。

也就是说，

$\mathbf {s} ^{(n)}(t)={\frac {ds^{(n-1)}}{dt}}=\lim _{\Delta t\rightarrow 0}{\frac {\Delta s^{(n-1)}}{\Delta t}}$.

C++

template<class Vector,class Number>
void Motion(Vector *s,Vector *prev_s,Vector a,Number dt,size_t Accuracy)
//s[] is the array of "current derivative of s values" that we are trying to calculate
//prev_s[] is the array of "previous derivative of s values" that is used as the constant factors
{
     size_t n(0);
     Number factor(1);
     
     //Note: the following code for position and velocity can be optimized for speed
     //      but isn't in order to explicitly show the algorithm
     for(s[0]=0;n<Accuracy;n++)
     //Position
     {
          if(n)factor*=(t/n);
          s[0]+=(factor*prev_s[n]);
     }

     for(s[1]=0,n=0,factor=1;n<(Accuracy-1);n++)
     //Velocity
     {
          if(n)factor*=(t/n);
          s[1]+=(factor*prev_s[n+1]);
     }

     s[2]=  a; //Acceleration
            //+PositionDependentAcceleration(s[0])                    //e.g. Newtonian Gravity
            //+VelocityDependentAcceleration(s[1])                    //e.g. Infinite Electro-magnetic field
            //+PositionVelocityDependentAcceleration(s[0],s[1]);      //e.g. Finite Electro-magnetic field

     for(n=3;n<Accuracy;n++)
     //Everything else
     //We're going to assume that dt is so small that s-prev_s is also very small
     //     i.e. approximates the differential
     {
          s[n]=(s[n-1]-prev_s[n-1])/dt;
     }
}

//Calling the function
//...
     Vector *s,*prev_s,*temp;
//...
     //The previous "current s" becomes the current "previous s"
     temp=prev_s;     //Speed optimization (instead of copying all of the s array
     prev_s=s;        //     into the prev_s array, we just swap the pointers)
     s=temp;
     Motion(s,prev_s,a,dt,Accuracy);
//...