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:

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

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

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

Herein lies the solution to the problem. The equation, as is, describes the complete path based on the initial constant factors ${𝐬}_0^{(n)}$. However if we treat this piece-meal, as in we calculate the numerical solution for $t$ equal some small incremental time $\mathrm{\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:

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

Notice that we gain accuracy as $\mathrm{\Delta }t\to 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 ${𝐬}_0^{(n)}$, we need only provide any external accelerations that affect the object, ${𝐬}^{(2)}$, from which we can calculate the next position ${𝐬}^{(0)}$ and velocity ${𝐬}^{(1)}$ of the object as well as the next set of constant factors ${𝐬}^{(n)}$, of which ${𝐬}^{(0)}$, ${𝐬}^{(1)}$ and ${𝐬}^{(2)}$ are already known.

Velocity is calculated by

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

or

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

The other ${𝐬}^{(n)}(t)$ values can be approximated by

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

provided that $\mathrm{\Delta }t$ is very small.

That is,

${𝐬}^{(n)}(t)=\frac{d{s}^{(n-1)}}{dt}=\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{\mathrm{\Delta }{s}^{(n-1)}}{\mathrm{\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);
//...