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: . This can be reduced to the Maclaurin series if we synchronize our clock so that is zero: . Herein lies the solution to the problem. The equation, as is, describes the complete path based on the initial constant factors . However if we treat this piece-meal, as in we calculate the numerical solution for equal some small incremental time and re-synchronize to 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: . Notice that we gain accuracy as and if we use more and more terms. After calculating the initial change in position of the object based on the initial constant factors , we need only provide any external accelerations that affect the object, , from which we can calculate the next position and velocity of the object as well as the next set of constant factors , of which , and are already known. Velocity is calculated by or . The other values can be approximated by for , provided that is very small. That is, . 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); //...