Midpoint Algorithm - Error C++

Question

I am currently finding an issue with my midpoint algorithm, I have been given a formula to work out the midpoint algorithm, but it seems to not correctly orbit. The aim is to get a more accurate version of the EulerMethod to create an more approximate orbit across all frame times. As the EulerMethod will cause a bit of an inaccurate orbit based on different frame-rates.

1st use Euler’s method to approximate half - way position / velocity :
position after half frame = current position +
half frame time * current velocity
velocity after half frame = current velocity +
half frame time * current acceleration

Calculate forces at half - way position to get half - way acceleration
Then use Euler’s method again with updated info :
position next frame = current position + frame time * halfway velocity
velocity next frame = current velocity + frame time * halfway acceleration

The method causes the 2nd sphere to move diagonally infinitely to the bottom right.

The Algorithm

void midpointmethod(CVector3& position, CVector3& velocity, CVector3& centre, float updateTime)
{
    CVector3 accel = OrbitAcceleration(position, Length(velocity), centre);
    //position = position + updateTime / 2 * velocity * updateTime;
    //velocity = velocity * updateTime / 2 * accel.Length() * updateTime;

    CVector3 halfpos = position + (updateTime / 2) * velocity;
    CVector3 halfveloc = velocity + (updateTime / 2) * accel; 

    // calculate forces ? 
    CVector3 halfwayaccel = accel / 2;


    CVector3 posnexf = position + updateTime * halfveloc;
    CVector3 velosnexf = velocity + updateTime * halfwayaccel;

    position = posnexf;
    velocity = velosnexf;

}

Answer 1

I'm glad you could fix your issue. The midpoint algorithm is more accurate than standard Euler for orbits, being a second-order algorithm, but it still has the problem of gaining energy over time (meaning that your orbiting object will move closer or farther from the planet due to numerical innacuracy, though it may take a while depending on your time step). I think a better algorithm to look at would be "semi-implicit" Euler, which will preserve the energy in your system and allow for a better orbit. Here is a great resource: https://www.youtube.com/watch?v=kxWBXd7ujx0. Or, if you can afford to ignore velocity, the Verlet algorithm also works great for orbits. Cheers :)

I wanted this to be a comment, but I don't have enough rep ...

Answer 2

Figured the issue, It was rather trivial, I misunderstood the halfwayacceleration value, and wasn't actually getting the halfwayacceleration between the two points, but rather a different point which was not drawing it back into the sphere.

CVector3 halfpos = position + (updateTime / 2) * velocity;
CVector3 halfveloc = velocity + (updateTime / 2) * accel; 

// calculate forces ? 
CVector3 halfwayaccel = OrbitAcceleration(halfpos, Length(halfveloc), centre);


CVector3 posnexf = position + updateTime * halfveloc;
CVector3 velosnexf = velocity + updateTime * halfwayaccel;

position = posnexf;
velocity = velosnexf;