Computer Programming/Physics/Motion of an object through space (piecemeal approximation)

<Wikisource:Source code

The problem with the Position of an accelerating body function is that although it is theoretically correct, it's not very adaptive in describing the motion of an object through space. To rectify this we attack the calculations piece-meal.

We refer back to the Taylor series form of the position function:

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

This can be reduced to the Maclaurin series if we synchronize our clock so that $t_0$ is zero:

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

Herein lies the solution to the problem. The equation, as is, describes the complete path based on the initial constant factors $\mathbf{s}^{(n)}_0$. However if we treat this piece-meal, as in we calculate the numerical solution for $t$ equal some small incremental time $\Delta t$ and re-synchronize $t_0$ to $t_{previous}$ after each calculation cycle, then we can dynamically build the path of the object through successive calculations using the following incremental form of the function:

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

Notice that we gain accuracy as $\Delta t\rightarrow 0$ and if we use more and more terms.

After calculating the initial change in position of the object based on the initial constant factors $\mathbf{s}^{(n)}_0$, we need only provide any external accelerations that affect the object, $\mathbf{s}^{(2)}$, from which we can calculate the next position $\mathbf{s}^{(0)}$ and velocity $\mathbf{s}^{(1)}$ of the object as well as the next set of constant factors $\mathbf{s}^{(n)}$, of which $\mathbf{s}^{(0)}$, $\mathbf{s}^{(1)}$ and $\mathbf{s}^{(2)}$ are already known.

Velocity is calculated by

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

or

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

The other $\mathbf{s}^{(n)}(t)$ values can be approximated by

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

provided that $\Delta t$ is very small.

That is,

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