I've dug up my old JS physics engine and am trying to make some scenes to test its capabilities. One problem I ran across is that restitution is not modeled correctly. In particular, Newton Cradle does not work as expected and it randomly bounces the balls around, as seen here: https://physics.rivsoft.net/Newton (controls are RMB to rotate, wheel to zoom).

As suggested by multiple sources, I'm adding relative normal velocity times restitution to the bias term (for 1D Jacobian), and using velocity at the start of the step for this. However, for obvious reasons this doesn't work if I have multiple velocity resolution steps. How am I supposed to update the bias between steps in order to correctly model restitution?

Your problem is likely that velocity is not propagating across the balls within a single timestep. So what happens is you get one collision over one timestep, then a bounce, and so on and so forth. This will quickly lead to many balls bouncing in many directions.

In real life velocity can propagate near instantaneously through all balls. You will need to do something like: build a graph of all balls in contact, and iterate many times over this graph while enforcing your constraints. Sequential impulses should be able to converge, so long as you build a system of all balls and let the system propagate velocities over the bodies through many iterations.

One trick that most physics engines have used for a while now, is to let rigid bodies penetrate slightly with a slop parameter. This parameter lets the narrow phase detect that two volumes intersect, without forcing the solver to move geometry out of an intersecting configuration. Often times when the configuration is no longer intersecting, internal graph data structure edges are broken, and the system is fragmented, unable to converge globally.

Hope this makes sense!

It actually is propagating immediately, since there are 8 velocity resolution steps currently. I think that's actually *why* it's not working, because if the velocity can propagate between steps, then the initial velocity used to calculate bias is no longer correct. I need a way to update bias inbetween steps, to account for propagated velocities, which should also force balls to bounce apart. If I reduce it to just one step, it starts working as expected, at least for a smaller number of balls (order starts to matter).

I do have that "slop" parameter, except it works the opposite way and creates contacts when objects are within a certain distance.

Hmm. In that case maybe you should try the shock propagation idea, where you lock a rigid body after it has been solved. This way you can ensure energy travels through the system in one direction. Otherwise sequential impulses might just spread around the velocity in an unrealistic way. It would be nice to see a graph or visualization of the propagation in any case.

How would you tell which way it should travel? I looked up some articles and they seem to suggest topological sort in the direction of gravity, so bottom objects resolve first, but in a general case we don't know where the impulses are coming from - for all we know it could be both ways.

I tried tweaking my solver to introduce restitution at different points, I can get it to propagate the impulses all the way through the chain, but then it goes back and tries to reduce the lambdas because it thinks they're too high (since balls are moving apart). Maybe there's some way to do it properly, by dynamically updating a lower limit on lambdas or something?

Would be cool to see a visualization of velocity propagation against different numbers of iterations on the first frame of collision with all spheres. Where has all the velocity gone with one iteration? With ten? With ten thousand?

The non-penetration constraints only say that the spheres can not inter-penetrate. It says nothing about where specifically all the energy in the system goes, and there can indeed be many "valid" solutions.

Another idea, just to get your cradle looking good: You can try separating the stand-still spheres ever so slightly. This way only one collision is solved in isolation at a time, and will prevent the solver from smearing your energy backwards through the system over many iterations. This might be sort of tricky with discrete collision, but with a continuous stepper should make for a really good looking cradle.

I don't have a way to visualize separate collision steps (rendering is completely detached from the physics engine atm), but I went through it in the debugger and it takes several iterations (like 4, depending on which order it decides to resolve everything in) to move everything nicely from the first two balls to the last two as expected, but then it starts going back and messing everything up - basically whenever it runs across a contact in which the balls are already separating (from restitution), it tries to adjust velocities back.

That gap trick worked, but now after the collision the balls aren't vertical anymore, so they start swinging back slowly, and after several swings it doesn't look so good anymore.

Yes this is expected. Sequential impulses will smear the velocities across the system in a manner dependent on your order of constraint solving. If you imagine there are four balls, each penetrating equally to their neighbor, and then think of how many different solutions could be valid here, there are infinite. So long as in the next step energy is conserved and the constraints are satisfied, this means the solver converged. The balls may teeter mostly to the left, mostly to the right, or almost all energy could be transferred to strictly an end ball.

This sort of thing reminds me of the "pool ball problem". Two pool balls are on a pool table touching tangentially. A third ball with some velocity strikes just between both balls. In a sequential solver one constraint must be satisfied before the other. The order of operations will result in two different solutions. Solving for a "real-life" solution can be difficult as well, since the system formed will be overdetermined with no direct solution.

Ah well, so most engines just ignore this and move on? Unless some system can figure out the propagation graph, but idk if that's possible to do for a generic complex scene.

Pool ball worked fine for me though, the two balls roll away at equal angles and the third ball rolls backwards slightly (or even forward, depending on friction).

Well maybe Dirk has an opinion, but I'm pretty sure modern physics engines are just modelling the constraints and calling it "good enough". Since typically in a scene with more than two rigid bodies human intuition isn't actually all that great at noticing situations without a real solution. Sequential impulses, being iterative, can find the minimum even in overdetermined or otherwise not directly solvable situations, making it "look nice" in most situations.

Even in real life a newton's cradle can fairly easily start behaving like your webpage simulation, and does not behave perfectly for long.

One other thing worth experimenting with, if you have the time, is to use soft contact constraints. Since in real life the balls are not perfectly rigid, maybe you can tune the cradle's soft contact settings and get closer to the behavior you like. Of course this isn't going to get you a general "perfect solution", but could be fun to try. I personally have never implemented a soft contact constraint, so I'm not sure exactly what that would look like.