I'm trying to implement a physics engine based upon Catto's 2005 paper. It's working quite nicely, and I have collision caching with warmstarting working. However, I still have too much jitter to properly stack boxes - they'll shake themselves apart eventually.

I have read that the problem is Baumgarte stabilisation (beta=0.1) introducing extra momentum which adds up over time. A solution which is mentioned in one of Catto's lecture slides is to use split impulses and store extra pseudo-velocities to separate the velocity constraint from the correction of the positional drift. I get the general idea - but I'm finding it difficult to dig up any concrete details (only 5ish sparse slides - don't think it is in the original paper) of this method or code to learn from. In particular I'm a bit confused how to clamp the lambdas when they are separated. Can anyone point me in the right direction?

Btw, this is for a pico-8 project - so a more complicated solver is probably out and I can't use another library (and the fun is in the learning ?). I've had a look in Box2D and Box2D-lite, but I don't think this is the approach taken (I could be mistaken).

Thanks

If you google some combination of “Erin Catto split impulse pseudo velocity” you should find a few threads about this, here are a couple to get you started:

https://pybullet.org/Bullet/phpBB3/viewtopic.php?t=1644

https://www.gamedev.net/forums/topic/613945-sequential-impulse-bias-velocities/

IMO for pico-8 you might want to consider a Position Based Dynamics solver (there are several papers about it here https://matthias-research.github.io/pages/publications/publications.html ) since IME this involves a bit less code.

(Note that most people use PBD in a cloth/deformable context but it's possible to make a rigid body simulation using the same ideas – it's a bit lengthy to explain here, but essentially each rigid body is modeled as 2 particles + 1 distance constraint, which defines a rigid 2d frame that can translate+rotate. You can then define constraint geometry (ie the “r” vector in Box2D terms (IIRC) -- the local-space position of the constrained points) in this frame so that “impulses” (I don't know what the term is for position-level deltas) applied to the embedded geometry are converted into impulses applied to the 2 particles.)