Hello

I've taken it upon myself to finally start learning how game engines work by making a toy engine myself. I've never done any game programming or anything like that but as an avid gamer and someone with a physics degree I've always wondered how these things are done.

Today I implemented a basic spring force using Hooks Law. I was surprised however when over a short time span, maybe 30s, my default integrator of Velocity Verlet caused the spring's motion to escalate wildly. I tried stepping down to Symplectic Euler and the motion became normal. At first I was thinking there might be a bug in my code, but I can't see one so now I'm wondering if the good people of this forum might help shed some light on this. It does not make any sense to me that Symplectic Euler would be a better integrator than Velocity Verlet. Both are time reversible and symplectic, but Velocity Verlet is second order accurate. 

I've been writing up what I learn on a little blog, feel free to check it out to see my current understanding of these methods

Symplectic Euler

Velocity Verlet

Code

Symplectic Euler

void SymplecticEuler::Solve(
            const NetForceAccumulator& net_force_accumulator,
            const std::vector<std::shared_ptr<PhysicsEntity>> &entity_ptrs,
            const std::shared_ptr<PhysicsEntity> entity_ptr)
{
    Vector3Gf xi = entity_ptr->GetPosition();
    Vector3Gf vi = entity_ptr->GetVelocity();
    GLfloat mass = entity_ptr->GetMass();

    Vector3Gf F;
    F.setZero();
    net_force_accumulator.ComputeNetForce(entity_ptrs,entity_ptr,F);

    Vector3Gf vf = vi + m_dt*(1/mass)*F;
    Vector3Gf xf = xi + m_dt*vf;

    entity_ptr->SetNextPosition(xf);
    entity_ptr->SetNextVelocity(vf);
};

Velocity Verlet

void Verlet::Solve(
            const NetForceAccumulator& net_force_accumulator,
            const std::vector<std::shared_ptr<PhysicsEntity>> &entity_ptrs,
            const std::shared_ptr<PhysicsEntity> entity_ptr)
{
    Vector3Gf xi = entity_ptr->GetPosition();
    Vector3Gf vi = entity_ptr->GetVelocity();
    GLfloat mass = entity_ptr->GetMass();

    Vector3Gf Fi,Ff,ai,xf,af,vf;

    Fi.setZero();
    net_force_accumulator.ComputeNetForce(entity_ptrs,entity_ptr,Fi);

    ai = (1/mass)*Fi;
    xf = xi + m_dt*vi + 0.5f*m_dt*m_dt*ai;

    entity_ptr->SetPosition(xf); // Load xf into position slot for computing F(x(t+h))

    Ff.setZero();
    net_force_accumulator.ComputeNetForce(entity_ptrs,entity_ptr,Ff); // Compute F(x(t+h))

    af = (1/mass)*Ff;
    vf = vi + 0.5f*m_dt*(ai + af);

    entity_ptr->SetPosition(xi); // Load xi back into position slot as to not affect force calulations for other entities

    entity_ptr->SetNextPosition(xf);
    entity_ptr->SetNextVelocity(vf);
};

Spring force

void SpringForceGenerator::AccumulateForce(
    const GLfloat k,
    const GLfloat l0,
    const std::shared_ptr<PhysicsEntity> entity1_ptr,
    const std::shared_ptr<PhysicsEntity> entity2_ptr,
    Vector3Gf &F) const
{
    Vector3Gf x1 = entity1_ptr->GetPosition();
    Vector3Gf x2 = entity2_ptr->GetPosition();

    Vector3Gf n = x2 - x1;
    GLfloat l = n.norm();
    n = n/l;

    F += k*(l - l0)*n;
}

Thanks!

I believe that the reason for the enhanced stability of Sympletic/Semi-implicit Euler is that it tends to have a damping effect (loses energy) and undershoots the integral, rather than overshoots. Vertlet might be more accurate but since it doesn't have that damping effect it is less stable. This is why most game physics engines (e.g. Bullet) use Semi-implicit Euler, even though it's not that accurate (it's also faster).

When I once implemented physics for a cannon projectile, Semi-implicit Euler would always undershoot the analytical trajectory, regular Euler would overshoot, and RK4 would be very close to the analytical path. With springs/constraints, RK4 was actually less stable than Semi-implicit Euler. Accuracy != stability.

Thanks for your reply Aressera. I'm aware of the difference between accuracy and stability (I wrote some articles about it as a slightly high level on the blog I linked). I would have expected the region of stability for Velocity Verlet to be nicer than Symplectic Euler, although to be fair I've never seen it pictured. Since both methods are symplectic they should both nearly conserve energy (see sections 1.1.3 and 1.2.2 here) so I would not have expected the damping of Symplectic Euler to have such a pronounced effect. However I guess the degree to which they (nearly) preserve energy need not be the same. I am new to all of this so I'll have to take your word for it. It's good to know I'm not insane after all!

Thanks!