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

<維基文庫：原始碼

關於 加速物體位置函式 的問題是，雖然理論上是正確的，但它在描述物體在空間中的運動方面並不十分靈活。為了解決這個問題，我們分段進行計算。

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

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