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

<维基文库：源代码

关于 加速物体位置函数 的问题是，虽然理论上是正确的，但它在描述物体在空间中的运动方面并不十分灵活。为了解决这个问题，我们分段进行计算。

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

$\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);
//...