計算機程式設計/物理/物體在空間中的運動（分段近似）

<維基文庫：原始碼

關於 加速物體位置函式 的問題是，儘管它在理論上是正確的，但它在描述物體在空間中的運動方面並不具有很強的適應性。為了糾正這一點，我們採用分段計算的方式。

我們回顧一下 泰勒級數形式的位置函式

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