How to use midpoint method to integrate a particle

Question

I'm using Euler method to numerically integrate particles in a physics engine.

float inverseMass;
Vector3 position, velocity;
Vector3 sumForces, gravity, inputForce;
void Particle::Integrate(float dt) {
   sumForces = (gravity+inputForce)*dt;
   Vector3 a = sumForces*inverseMass;
   velocity += a;
   position += velocity*dt;
}

I would like to implement a better integrator with Midpoint method. Based on this course https://www.cs.cmu.edu/~baraff/pbm/constraints.pdf, with Midpoint method you take a full euler step, then evaluate the force on the particle at the midpoint and then take a step using that midpoint value.

I implemented the first step.

void Particle::Integrate(float dt) {
   //Take full Euler step
   Vector3 sumForces = (gravity+inputForce)*dt;
   Vector3 a = sumForces*inverseMass;
   Vector3 tempVelocity = a+velocity;
   Vector3 deltaX += tempVelocity*dt;

   Vector3 midpoint = deltaX + position;

}

How do I continue from here though? How do I calculate the force at midpoint? Do I just calculate it at half time step, but then what is the purpose of calculating the midpoint position?

Answer 1

You have to consider that your state contains two main components, position x and velocity v. The method has to be formulated uniformly for the full state,

dx1 = v*dt,           dv1 = acc(t,x,v)*dt
dx2 = (v+0.5*dv1)*dt, dv2 = acc(t+0.5*dt, x+0.5*dx1, v+0.5*dv1)*dt

x = x + dx2,          v = v+dv2.

Having written down the formulas from first principles, you can now see that it is relatively easy to eliminate the dx vectors. Thus it remains to compute

dv1 = acc(t,x,v)*dt,
dv2 = acc(t+0.5*dt, x+0.5*v*dt, v+0.5*dv1)*dt,
x = x + (v+0.5*dv1)*dt,
v = v + dv2