Friday March 6, 2026

Home | Contact | Support | Programming.. More than just code .... | Physics.... Squishy, Bouncy, Springy, Rigid, ...Bodies...

Physics....

Squishy, Bouncy, Springy, Rigid, ...Bodies...

Looking at the whole picture...calculating the forces to apply each frame for the constraints... everything from objects sittin... — Looking at the whole picture...calculating the forces to apply each frame for the constraints... everything from objects sitting on the floor, friction for rolling, joints and limbs on ragdolls to ropes and wheels on cars.

Project Code
• 2D Vanilla JavaScript and Canvas - 2D/Stacks,Cube,Cars,Motors [LINK]
• 3D Vanilla JavaScript and WebGL - 3D/Cubes/Constraints/Hinges/Ragdolls [LINK]

What are global constraint solvers?

You know that moment when you drop a bunch of cubes into a scene and they all start bouncing around like excited puppies?

And then you stack them... and suddenly everything goes wobbly, like the whole tower is trying to remember how gravity works.

That's usually the moment you realise:

"Oh... maybe solving collisions one at a time isn't going to cut it."

Because at first, the "obvious" approach feels fine. You detect collisions, loop through them, fix each one, move on, and it feels like you're being super clean and systematic about it. You’re basically telling yourself, "Look at me, I’m doing proper physics," because each little pair gets handled, and the scene looks stable for a moment.

Tidy. Logical. Satisfying.

And it works... right up until you try to build anything more interesting than a single cube.

A tiny detour: what a physics step is actually doing

Before we go full matrix-mode, it helps to picture what a physics engine does every frame (or every timestep).

A super simplified loop looks like:

1. You start with velocities (how fast things are moving/spinning).
2. You add things like gravity (so velocities change).
3. You detect collisions / contacts (who's touching who).
4. You solve constraints (stop overlaps, enforce joints, limits, etc).
5. You integrate positions (move objects forward using the new velocities).

That "solve constraints" bit is where stacks live or die.

Because stacks aren't just "two cubes colliding."
They're a whole network of "everyone is touching everyone."

Why stacks break the "one collision at a time" idea

Stacks are different!

And the reason is actually really simple once you say it out loud... when you solve contacts "one at a time", you're usually solving them as pairs. Just two bodies. One contact. One little local problem.

Cube A is overlapping Cube B... cool... push them apart.
Cube B is overlapping Cube C... cool... push them apart.
And so on.

Sounds sensible... right?

But here's the catch: a stack isn't a bunch of independent pairs. It's a chain. And if you only look at one link at a time, you don't see the weight of the whole chain.

Picture a stack of ten cubes.

That bottom cube has ten cubes worth of weight going through it.
Not just "itself"... not just "the one above"... all of them.

In real life, the ground contact has to push up with the combined weight of the whole stack. If it only pushes back with "one cube worth", the bottom cube doesn't magically stay put... it gets squashed down, then the next cube settles, then the next... and pretty soon your tower is doing that slow-motion "melt into the floor" thing.

And this is exactly what happens with a naive pair-by-pair solver in a single pass.

Because when you solve the ground contact for the bottom cube, what does it "know"?

It knows:
• bottom cube vs ground
• the current relative velocity between those two
• the current penetration at that one contact

What it doesn't directly know (yet) is:
• "Oh by the way, there are nine more cubes above me and I'm about to have a very bad day."

So you end up with the bottom contact applying a push that makes sense for the pair... but not for the stack.

It's like trying to support a pile of books with your hand while only paying attention to the top book.
You can keep the top one from falling, sure... but your wrist is still going to sink because the load underneath is the real story.

Newton's cradle is the classic desk-toy version of this. You tap one end... and the response shows up on the other end. The middle balls aren't "extra"... they're the path the impulse travels through. A stack is the same idea, just with gravity constantly tapping the whole thing, every frame.

So if you do:

"contact 1... solved!"
"contact 2... solved!"
"contact 3... solved!"

...only once... in a single sweep...

You haven't actually found a consistent answer. You've just nudged a bunch of connected things without letting the information (the load) propagate properly through the chain.

That's why stacks wobble, jitter, or sink with simple solvers: the bottom contact is acting like it's only supporting the cube directly above... because that's all it can see in that moment.

And that's also why iterative/global solvers feel like magic when they work.

They go around again... and again... and again... passing the "load" through the stack until the numbers settle.

Top contact ends up supporting one cube.
Next supports two.
Next supports three.
And the bottom contact finally goes: "Ahhh... right... I'm holding the whole tower."
...and pushes up with the full amount.

That's the difference.

Not "more accurate collision detection"...
but "the solver actually understands the whole stack is one connected problem."

What a global solver does instead

So here's where you stop doing the "pair by pair" thing... and start thinking like the engine.

A global solver basically treats the whole frame like one big situation. Not "contact 12 needs fixing", but:

"Ok... what rules are being broken right now, across the entire world... and what set of nudges would make everything behave at the same time?"

That little shift in mindset is everything.

Instead of chasing problems one-by-one (and accidentally undoing yourself), you first gather all the rules you care about — every contact, every joint, every limit — and you solve them together as one connected system.

And those rules have a name: constraints.

Constraints are just small promises your simulation is trying to keep:

• don't overlap
• don't sink into the ground
• don't stretch this joint
• don't rotate past this angle

Basically: "after the physics step, I want this statement to be true."

Once you've collected those promises, you can finally write them down as math. And it turns out you already know the shape of the math... even if you haven't met it in a physics context before.

The equation you've probably seen before

You've probably come across this in school:

A x = b

That's the classic "linear system."
It's just simultaneous equations wearing a nicer outfit.

Instead of writing:

• equation 1
• equation 2
• equation 3

...with the unknowns tangled together, you bundle the whole thing into one sentence:

"Find the vector \(x\) that makes all these equations true at once."

Same idea. Cleaner packaging. And honestly? Much easier to hand to a computer.

In physics we use the same structure, but we change the letter because it means something very specific:

A λ = b

Same shape. Different letters.
And those letters are basically the whole story.

The pieces: \(b\), \(\lambda\), and \(A\)

\(b\): "What's wrong right now?"

Think of \(\,b\,\) as your solver's to-do list for this frame.

It's the stuff you want to fix:

• "This contact is overlapping by a bit."
• "This joint drifted apart."
• "This hinge rotated past its limit."

If it helps, \(b\) is like the list of complaints your simulation is making right now.

It's not the solution. It's the problem report.
"Here's how wrong things currently are."

\(\lambda\): "What should I do about it?"

\(\,\lambda\,\) (lambda) is what you solve for.

It's the list of impulses the solver decides to apply.

Not forces — impulses.

And if that word feels a bit slippery, here's the friendly version:

• A force is something you apply over time (like pushing a shopping trolley steadily).
• An impulse is like a quick tap that instantly changes velocity (like kicking a ball).

Constraint solvers love impulses because they operate "per frame."
They're basically saying:

"Given the mess you're in right now, what velocity changes should I apply so these rules stop being broken?"

Each constraint gets one or more \(\lambda\) values depending on what it's enforcing:

• contact normal: one \(\lambda\)
• friction: often two more (sideways directions)
• more complex joints: more components

So \(\lambda\) is the solver's "action plan."
A big list of "push here by this much."

\(A\): "How connected is everything?"

And then there's \(\,A\,\).

This is the part that makes it global.

\(A\) describes how constraints influence each other.

So if you push on cube #5, cube #4 reacts, which changes cube #3, and so on.
That ripple is exactly what \(A\) captures.

A nice way to think about it is:

• \(b\) is the problem list
• \(\lambda\) is the fix list
• \(A\) is the awkward truth that fixes aren't independent

It's the dependency graph between fixes.

Like:

"If I tweak this constraint, which other constraints does it accidentally mess with?"

Where \(A\) comes from (no magic involved)

You don't invent \(A\) out of thin air.
You build it from two very physical ideas:

A = J M⁻¹ J^T

This looks intense, but the meaning is actually pretty intuitive once you translate it into normal human language.

\(\,M^{-1}\,\) is the inverse mass matrix (i.e., How stubborn is each body?)

It tells you how much an impulse will change a body's velocity:

• Heavy body: tiny velocity change
• Light body: big velocity change
• Infinite mass / pinned body: no change

And because impulses can also spin objects, \(M^{-1}\) also includes rotational inertia.

A really concrete picture:

• Push a bowling ball → it barely moves.
• Push a ping-pong ball → it zooms away.
• Push a nailed-down object → nothing happens.
• Push at the edge of a box → it spins.

That "spin part" is inertia showing up.

\(\,J\,\) is the Jacobian (i.e., What does this constraint care about?).

I like to think of it as the constraint's "measuring tool."

A constraint is basically asking:

"Along my direction, how fast are things trying to move right now?"

For a contact constraint, the direction is the contact normal \(\mathbf{n}\).
So it cares about velocity along the normal.

If two bodies are moving into each other along the normal → that's a problem.
If they're separating → no problem.

For a joint, it's measuring relative velocity between anchor points.
For a hinge limit, it's measuring relative angular velocity past the limit direction.

So \(J\) is the thing that turns "raw body velocities" into "constraint-specific motion."

That's all it is: the wiring between "what bodies are doing" and "what the constraint complains about."

What a body looks like to the solver (it's simpler than you think)

This part is super helpful for beginners, because it strips away the fluff.

To the solver, a 3D rigid body is basically a little slot of velocities:

• 3 numbers for linear velocity \(v\)
• 3 numbers for angular velocity \(\omega\)

So in 3D you can think of it as:

\[
\begin{bmatrix}
v_x \\ v_y \\ v_z \\ \omega_x \\ \omega_y \\ \omega_z
\end{bmatrix}
\]

That's the stuff impulses change.

The solver does not care about your mesh.
It does not care about your texture.
It does not care that you named it "Cube_Anthony_Final_FINAL2".

It cares about:

"How fast are you moving?"
"How fast are you spinning?"
"Cool. Here's an impulse. Now you move differently."

Once you see bodies like that — as velocity packets — the rest of the solver suddenly feels less mystical.

A really grounded Jacobian example (contact at center vs edge)

If a cube is sitting flat on the ground and you push it straight up at its center, it doesn't spin.
It just goes up.

If you push it straight up at a corner, it spins.

Same direction, different point → different rotational effect.

That's why Jacobian rows usually have two parts:

• linear part: "this direction affects translation"
• angular part: "this point causes rotation too"

For one body, you'll often see something like:

\[
J = [\, \mathbf{n},\; (r \times \mathbf{n}) \,]
\]

You don't need to worship the cross product here.
Just read it as:

• \(\mathbf{n}\): "what direction am I pushing in?"
• \(r \times \mathbf{n}\): "am I pushing off-center, causing spin?"

That's it. That's the whole vibe.

So what does \(A = J M^{-1} J^T\) feel like?

Here's the plain-English version, with the "shape" intact:

1. \(J^T\) takes a constraint impulse and turns it into body-space impulses.
(Like: "this constraint wants to push — what does that mean for the bodies?")

2. \(M^{-1}\) turns impulses into velocity changes.
(Heavy bodies change less, light bodies change more.)

3. \(J\) measures those new velocities back in constraint-space.
(Like: "after pushing, how did the constraint situation change?")

So \(A\) is basically:

"If I push on constraint \(i\), how does that affect constraint \(j\)?"

And because most constraints only touch a couple of bodies, most constraints don't affect most other constraints.

That's why \(A\) is usually sparse (mostly zeros).
Your tower might be big, but each contact only connects nearby pieces.

The practical twist: you don't solve it directly

Now here's the part that sneaks up on you.

Even though the math looks like a clean "solve it once" system, you don't usually solve \(A \lambda = b\) exactly.

That would be expensive.

Like: "your laptop becomes a space heater" expensive.

So real-time physics engines usually use an iterative method, most commonly Gauss–Seidel (or something in that family).

And the vibe is very "tighten it gradually."

Gauss–Seidel, explained like you're adjusting a wobbly chair

Imagine your constraints are screws on a wobbly chair.

You don't tighten one screw once and declare victory.

You tighten screw 1 a bit.
Then screw 2.
Then screw 3.

But tightening screw 3 changes the load, so screw 1 might be a little loose again.

So you go around again.

That's iterative solving.

In solver terms, each iteration looks like:

1. pick one constraint
2. compute a little correction to \(\lambda\)
3. apply that impulse (so velocities change)
4. move to the next constraint
5. loop again

Each pass makes the system more consistent.
More passes = less wobble.

And yes: the number of iterations is one of the big "feel" knobs in a physics engine.

Why contacts force you to upgrade the math (enter the LCP)

So far, \(A\lambda = b\) feels like a nice clean system.

But then you hit contacts.

Contacts have a rule that joints don't:

They can push, but they cannot pull.

A contact constraint is like a polite bouncer at a club:

• If you're trying to walk into the wall: "Nope, push back."
• If you're already separating: "Cool, do nothing."
• It never grabs you and drags you into the wall.

So if your solver produces a negative impulse at a contact, that would mean "pull them together," which is physically wrong.

That's the moment you switch from "just a linear system" to a Linear Complementarity Problem (LCP).

Instead of only solving:

\[
A \lambda = b
\]

you add the contact rules:

\[
A \lambda + b \ge 0,\quad
\lambda \ge 0,\quad
\lambda^{T}(A \lambda + b) = 0
\]

If you want the beginner translation:

• \(\lambda \ge 0\): "only push"
• \(A\lambda + b \ge 0\): "don't push through the contact / don't create negative gap"
• \(\lambda^{T}(A \lambda + b) = 0\): "either you're pushing, or you're separated — not both"

That last line is the magic "either-or" rule.

It's saying:

• If the bodies are not touching (gap > 0), the contact force should be 0.
• If the bodies are touching (gap = 0), the contact force can be > 0.

Never "gap > 0 and force > 0" at the same time.
Never "gap < 0 and force < 0" weirdness.

It's the solver enforcing: push exactly when needed, otherwise chill.

A super relatable example: a book on a table

Put a book on a table.

Gravity pulls down.

The table pushes up.

That upward push is the normal force — a contact constraint doing its job.

Now lift the book slightly.

The moment it's not touching the table, the normal force becomes zero.

The table doesn't chase the book upward with a magical pulling force.

That's complementarity.

Either:

• contact is active and pushing
or
• contact is inactive and doing nothing

Why we talk about 3D (and how 2D fits in)

Quick reassurance: if you're thinking "this seems very 3D," you're right — but it's the same structure in 2D.

The only real difference is how many velocity components each body has.

In 3D:

• 3 linear + 3 angular = 6

In 2D:

• 2 linear + 1 angular = 3

Same idea, fewer knobs.

If you understand the 3D version, 2D is basically the same thing with half the furniture removed.

The whole picture (in plain language)

So here's the flow you can keep in your head:

• Your simulation has a bunch of rules (constraints).
• You gather them all up each frame.
• You figure out "what's wrong right now" (\(b\)).
• You solve for "what impulses fix it" (\(\lambda\)).
• You account for how everything is connected (\(A\)).
• For contacts, you add "only push" rules (LCP).
• You solve it iteratively so your computer doesn't melt.

And once this clicks, constraint solving stops feeling like mysterious math...

...then starts feeling like a big, slightly fussy, but very logical machine you can actually reason about.

Which is exactly when your stacks stop wobbling, your joints stop stretching, and your simulation finally feels... right.

Third Gear: a tiny 1D stack you can solve by hand (with actual numbers)

Alright — now let's make this real.

We're going to do a full toy example in 1D (just up/down), with a two-block stack on the ground, and we'll write down the matrices with the values inside.

But this time we're not just dumping numbers and hoping they "look right."
We're going to keep doing the sanity checks that make solvers finally click:

"Do these values make sense?"
"Is this what I'd expect physically?"
"What would I change if I wanted different behaviour?"

That mindset is the difference between "math soup" and "ohhhh... ok I get it."

The scene (1D, vertical only)

Picture this:

• Block 1 (bottom) is sitting on the ground.
• Block 2 (top) is sitting on Block 1.
• Everything is perfectly aligned, so we only care about vertical motion.

We'll use easy numbers:

• \(m_1 = 1\)
• \(m_2 = 1\)
• timestep \(dt = 1\) (just to keep arithmetic simple)
• gravity gives a velocity change of \(-10\) each step

Start at rest:

\[
v =
\begin{bmatrix}
v_1 \\ v_2
\end{bmatrix}
=
\begin{bmatrix}
0 \\ 0
\end{bmatrix}
\]

Step 1: apply gravity (the "free" velocities)

Before constraints do anything, gravity tries to pull both blocks down:

\[
v_{\text{free}} =
\begin{bmatrix}
-10 \\ -10
\end{bmatrix}
\]

Meaning: both bodies are now moving downwards.
If you did nothing else, they'd happily fall straight through the ground and each other.

So constraints exist to stop that.

Step 2: define the constraints (what rules are we enforcing?)

In 1D, a contact constraint is basically:

"You are not allowed to move into the thing you're touching."

We have two contacts, so two constraints.

Constraint 0: Block 1 vs ground

Ground is "fixed" (infinite mass), and block 1 is the moving thing.

We'll use a simple sign convention:

• positive \(v\) means "up"
• negative \(v\) means "down"

So the constraint velocity is simply:

\[
c_0 = v_1
\]

If \(c_0 < 0\), block 1 is closing into the ground → constraint should push back.

Constraint 1: Block 2 vs Block 1

This one cares about the relative motion between the two blocks.

If block 2 is moving down faster than block 1, then block 2 is "catching up" and pushing into it.

So:

\[
c_1 = v_2 - v_1
\]

If \(c_1 < 0\), top block is closing into the bottom block → constraint should push back.

So far this is just "what we measure" for each constraint.

Step 3: build \(J\) (the Jacobian) — and check the signs make sense

The Jacobian maps body velocities into constraint velocities:

\[
c = J v
\]

We have 2 constraints and 2 bodies, so \(J\) is \(2 \times 2\).

From our definitions:

• \(c_0 = v_1\)
So constraint 0 depends on body 1 with coefficient \(+1\), and not at all on body 2.

Row 0:

\[
[\,1,\;0\,]
\]

Sanity check: yep. The ground contact should only "see" block 1 directly.

• \(c_1 = v_2 - v_1\)
So it depends on \(v_1\) with coefficient \(-1\) and on \(v_2\) with coefficient \(+1\).

Row 1:

\[
[\,-1,\;1\,]
\]

Sanity check: perfect.
- If both move together, \(v_2=v_1\) → \(c_1=0\).
- If top falls faster, \(v_2 - If top moves up relative to bottom, \(c_1>0\) (separating).

So the Jacobian is:

\[
J =
\begin{bmatrix}
1 & 0 \\
-1 & 1
\end{bmatrix}
\]

Those \(1\), \(0\), \(-1\) aren't random. They're literally encoding:

• who participates in the constraint
• and whether their motion increases or decreases the measured "closing speed"

Quick rule of thumb:

• \(+1\) → "this body moving up increases the constraint velocity"
• \(-1\) → "this body moving up decreases it"
• \(0\) → "this body doesn't matter here"

Step 4: build \(M^{-1}\) (inverse mass) — and what it physically means

In 1D with no rotation, inverse mass is just a diagonal matrix:

\[
M^{-1} =
\begin{bmatrix}
1/m_1 & 0 \\
0 & 1/m_2
\end{bmatrix}
\]

With \(m_1=m_2=1\):

\[
M^{-1} =
\begin{bmatrix}
1 & 0 \\
0 & 1
\end{bmatrix}
\]

Meaning: an impulse of \(+10\) changes velocity by \(+10\) for either body.

If you made the top block heavier, say \(m_2=2\), this becomes:

\[
M^{-1} =
\begin{bmatrix}
1 & 0 \\
0 & 1/2
\end{bmatrix}
\]

Meaning: same impulse gives half the velocity change on the heavier block.
Exactly what you'd expect.

Step 5: build the global matrix \(A\) — and read what its values are telling you

We use:

\[
A = J M^{-1} J^{T}
\]

Here \(M^{-1}\) is identity, so:

\[
A = J J^T
\]

Compute:

\[
J =
\begin{bmatrix}
1 & 0 \\
-1 & 1
\end{bmatrix},
\quad
J^T =
\begin{bmatrix}
1 & -1 \\
0 & 1
\end{bmatrix}
\]

\[
A =
\begin{bmatrix}
1 & 0 \\
-1 & 1
\end{bmatrix}
\begin{bmatrix}
1 & -1 \\
0 & 1
\end{bmatrix}
=
\begin{bmatrix}
1 & -1 \\
-1 & 2
\end{bmatrix}
\]

Now pause. This matrix is insanely informative.

What do the entries mean?

• \(A_{00} = 1\)
Ground constraint affects itself with strength 1.
Makes sense: it touches one body of mass 1.

• \(A_{11} = 2\)
Block-block constraint affects itself with strength 2.
Why bigger? Because it involves two bodies. When you push at the interface, you can change both velocities, so you've got more "leverage."

• \(A_{01} = A_{10} = -1\)
This is the stack coupling.
It says the two constraints are linked through block 1. If you apply impulse at the interface, it affects the ground contact too.

So even before solving anything, \(A\) is already whispering:

"These constraints aren't independent. They share a body. They're coupled."

Step 6: compute \(b\) — and check the story it tells

A common choice for the RHS (in velocity solving form) is:

\[
b = J v_{\text{free}}
\]

Compute:

\[
v_{\text{free}}=
\begin{bmatrix}
-10 \\ -10
\end{bmatrix}
\]

\[
b =
\begin{bmatrix}
1 & 0 \\
-1 & 1
\end{bmatrix}
\begin{bmatrix}
-10 \\ -10
\end{bmatrix}
=
\begin{bmatrix}
-10 \\
0
\end{bmatrix}
\]

Interpretation:

• \(b_0 = -10\) → the bottom block is closing into the ground at speed 10. Bad.
• \(b_1 = 0\) → the two blocks are falling together, so there's no closing between them (yet).

That's exactly what you'd expect. Gravity pulls both equally, so the top isn't "catching up" to the bottom.

Step 7: solve the contact rules (LCP) — and why it gives the numbers you'd expect

For contact constraints, we enforce:

\[
w = A\lambda + b
\]
\[
w \ge 0,\quad \lambda \ge 0,\quad \lambda^T w = 0
\]

Write it out:

\[
A=
\begin{bmatrix}
1 & -1 \\
-1 & 2
\end{bmatrix},
\quad
\lambda=
\begin{bmatrix}
\lambda_0 \\ \lambda_1
\end{bmatrix},
\quad
b=
\begin{bmatrix}
-10 \\ 0
\end{bmatrix}
\]

\[
w=
\begin{bmatrix}
\lambda_0-\lambda_1-10 \\
-\lambda_0+2\lambda_1
\end{bmatrix}
\]

If the stack is resting, you expect both contacts to be active and end up with zero closing velocity:

\[
w_0=0,\quad w_1=0
\]

Solve:

\[
\lambda_0-\lambda_1-10=0
\]
\[
-\lambda_0+2\lambda_1=0
\Rightarrow \lambda_0=2\lambda_1
\]

Plug in:

\[
2\lambda_1-\lambda_1-10=0
\Rightarrow \lambda_1=10
\Rightarrow \lambda_0=20
\]

So:

\[
\lambda=
\begin{bmatrix}
20 \\ 10
\end{bmatrix}
\]

What do these \(\lambda\) values mean?

This is the part where it should feel obvious, not magical.

• \(\lambda_1 = 10\) is the impulse needed at the interface to support the top block.
The top block wants to fall at \(-10\), so you need +10 impulse there to cancel it.

• \(\lambda_0 = 20\) is the impulse needed at the ground to support both blocks.
It has to cancel the bottom block's \(-10\) and also transmit the load from the top.

So ground contact gets double because it's carrying the whole stack.

That is exactly what you expect:
• top contact supports only the top block
• ground supports everything

Step 8: apply impulses back to bodies — and verify we get what we wanted

Velocity correction is:

\[
\Delta v = M^{-1} J^T \lambda
\]

\[
J^T=
\begin{bmatrix}
1 & -1 \\
0 & 1
\end{bmatrix}
\]

Compute:

\[
J^T \lambda =
\begin{bmatrix}
1 & -1 \\
0 & 1
\end{bmatrix}
\begin{bmatrix}
20 \\ 10
\end{bmatrix}
=
\begin{bmatrix}
10 \\ 10
\end{bmatrix}
\]

Since \(M^{-1}\) is identity:

\[
\Delta v=
\begin{bmatrix}
10 \\ 10
\end{bmatrix}
\]

Apply:

\[
v_{\text{post}} = v_{\text{free}}+\Delta v=
\begin{bmatrix}
-10 \\ -10
\end{bmatrix}
+
\begin{bmatrix}
10 \\ 10
\end{bmatrix}
=
\begin{bmatrix}
0 \\ 0
\end{bmatrix}
\]

So after solving, both blocks have zero velocity: they're resting.

Again: exactly what you expect.

Other cases: changing the rule changes the Jacobian (and the solver behaviour)

This is the fun part, because it shows how flexible this framework is.

The solver isn't "a collision solver."
It's a "constraint satisfier."

So if you want different behaviour, you don't rewrite the solver — you rewrite the constraint.

Case A: "Lock these two blocks together"

If you want block 2 and block 1 to move together always (like they're glued), you're enforcing:

\[
v_2 - v_1 = 0
\]

That's the same measurement we used for contact, but the difference is:

• contact only pushes when needed (\(\lambda \ge 0\))
• a lock constraint is allowed to push and pull (so \(\lambda\) can be positive or negative)

So you treat it like an equality constraint (not complementarity).

The Jacobian row is the same:

\[
J_{\text{lock}}=[-1,\;1]
\]

But you don't clamp \(\lambda\ge 0\).

That single change is huge:
• contact = one-way bouncer (only push)
• lock = rigid bar (push or pull)

Case B: "Limit how fast the top block can move relative to the bottom"

Maybe you don't want a full lock. Maybe you want:

• the top block can slide, but only up to some max relative speed
• or you want a damper-like "don't let it slam"

In a velocity solver, you can express that as:

\[
v_2 - v_1 \le v_{\max}
\]

That turns into the same complementarity pattern you saw with contacts, but with a different target.

Rearrange:

\[
(v_2 - v_1) - v_{\max} \le 0
\]

Now the constraint becomes active only when you exceed the limit.

Same Jacobian row \([-1,\;1]\), different RHS / bias.

Case C: "Only prevent downward relative motion (like a one-way clutch)"

This is basically contact behaviour but between the blocks, even if they're not physically colliding.

You enforce:

• if \(v_2 - v_1\) tries to go negative (top closing down), push back
• but if it goes positive (separating), do nothing

Same Jacobian. Same LCP rule. Different intention.

The big lesson

The numbers in \(J\) are not just math. They're you saying:

• which bodies participate
• how their motion combines (plus/minus)
• what direction counts as "violating" the rule

And whether you clamp \(\lambda\) or not decides if the constraint is:

• one-way (contact, inequality, limit)
• two-way (joint/lock, equality)

Once you see that, you can design behaviours on purpose instead of hoping.

If you want a nice next experiment: change \(m_2 = 2\) (heavier top block) and watch how \(\lambda_1\) and \(\lambda_0\) shift. It's a great intuition builder.

Fourth Gear: two 3D rigid bodies with a ball joint (a real "joint constraint" example)

Alright — now we step up from the cute 1D stack into something that feels like "real engine stuff."

Two rigid bodies in 3D, connected by a ball joint.

A ball joint is the simplest "proper" 3D joint to learn, because the rule is so clean:

"These two points — one on body A and one on body B — must stay glued together in space."

Like a shoulder, a hip, or any socket-type joint.

No limits yet. No motors. No friction.
Just: "these anchor points should match."

And the nicest part is you're not learning new machinery here. It's still the same pieces:

• bodies have velocities
• constraints measure motion with a Jacobian
• the solver computes impulses \(\lambda\)
• those impulses change velocities so the rule stops being broken

1) What the solver thinks a 3D body is (quick reminder)

Each rigid body carries 6 velocity numbers:

\[
\mathbf{v} =
\begin{bmatrix}
v_x \\ v_y \\ v_z \\ \omega_x \\ \omega_y \\ \omega_z
\end{bmatrix}
\]

So with two bodies, the solver's velocity vector is 12 numbers:

\[
\mathbf{V} =
\begin{bmatrix}
\mathbf{v}_A \\
\mathbf{v}_B
\end{bmatrix}
\in \mathbb{R}^{12}
\]

That's the thing the solver nudges.

It doesn't directly change positions — it changes velocities, then integration moves positions later.

2) The ball joint rule (what are we enforcing?)

Pick an anchor point on each body, expressed in local space:

• \(\mathbf{r}_A\): anchor offset from A's center of mass
• \(\mathbf{r}_B\): anchor offset from B's center of mass

In world space, those points are:

\[
\mathbf{p}_A = \mathbf{x}_A + R_A \mathbf{r}_A
\]
\[
\mathbf{p}_B = \mathbf{x}_B + R_B \mathbf{r}_B
\]

And the ball joint constraint says:

\[
\mathbf{p}_A - \mathbf{p}_B = \mathbf{0}
\]

That's 3 equations (x, y, z).
So a ball joint contributes 3 constraint rows.

The useful mental picture:

• It prevents separation at the anchor point.
• It does not prevent rotation around that point.

So body B can spin like mad, as long as the anchor stays attached.

3) The constraint error and the "to-do list" \(b\)

Let the position error be:

\[
\mathbf{C} = \mathbf{p}_A - \mathbf{p}_B
\]

If \(\mathbf{C}\) isn't zero, the joint is "stretched" — the anchors don't line up.

In practice, solvers don't just use position error. They use velocity form + a bias term (to correct drift).

Velocity form:

\[
\dot{\mathbf{C}} = \mathbf{v}_{pA} - \mathbf{v}_{pB}
\]

Anchor point velocities:

\[
\mathbf{v}_{pA} = \mathbf{v}_A + \boldsymbol{\omega}_A \times (R_A \mathbf{r}_A)
\]
\[
\mathbf{v}_{pB} = \mathbf{v}_B + \boldsymbol{\omega}_B \times (R_B \mathbf{r}_B)
\]

And a typical "goal" is:

\[
\dot{\mathbf{C}} + \beta \frac{\mathbf{C}}{dt} = \mathbf{0}
\]

That \(\beta\) term is the classic "don't drift forever" knob:
• if the joint drifts a bit, you gently pull it back
• so it doesn't slowly stretch apart over time

So your \(b\) ends up being "current relative anchor velocity + bias."

4) Building the Jacobian \(J\) (and why the signs look exactly how you'd expect)

We want a matrix \(J\) such that:

\[
J \mathbf{V} = \dot{\mathbf{C}}
\]

Remember:
• \(\dot{\mathbf{C}}\) is 3 numbers (x, y, z)
• \(\mathbf{V}\) is 12 numbers (6 for A, 6 for B)

So \(J\) is a \(3 \times 12\) matrix.

Let:
• \(\mathbf{r}_A^w = R_A \mathbf{r}_A\)
• \(\mathbf{r}_B^w = R_B \mathbf{r}_B\)

Define \([r]_\times\) such that:

\[
[r]_\times \omega = r \times \omega
\]

Then the Jacobian is:

\[
J =
\begin{bmatrix}
I_3 & -[\mathbf{r}_A^w]_\times & -I_3 & [\mathbf{r}_B^w]_\times
\end{bmatrix}
\]

The sanity checks you should always do:

• The linear blocks are \(+I_3\) and \(-I_3\).
That's literally "A minus B."

• If \(r=0\) (anchor at center), the angular blocks should disappear.
And they do.

That's the matrix behaving like physics.

5) The inverse mass matrix \(M^{-1}\)

Each body contributes a \(6 \times 6\) block:

\[
M^{-1}_A =
\begin{bmatrix}
m_A^{-1} I_3 & 0 \\
0 & I_A^{-1}
\end{bmatrix}
\quad
M^{-1}_B =
\begin{bmatrix}
m_B^{-1} I_3 & 0 \\
0 & I_B^{-1}
\end{bmatrix}
\]

Global matrix is block diagonal:

\[
M^{-1} =
\begin{bmatrix}
M_A^{-1} & 0 \\
0 & M_B^{-1}
\end{bmatrix}
\]

Same intuition as before:
• heavy = stubborn
• light = twitchy
• pinned = doesn't move

6) From \(J\) and \(M^{-1}\) to impulses \(\lambda\)

Ball joint has 3 rows, so:

\[
\lambda =
\begin{bmatrix}
\lambda_x \\ \lambda_y \\ \lambda_z
\end{bmatrix}
\]

This is a key beginner intuition:

A ball joint doesn't produce "one impulse."
It produces a 3D impulse vector at the joint.

And unlike contacts, it can push and pull, because it's an equality constraint.

7) How the solver applies the joint impulse to bodies

Once the solver picks \(\lambda\), it applies it back through \(J^T\):

\[
\Delta \mathbf{V} = M^{-1} J^T \lambda
\]

Physically:
• body A gets \(+\lambda\)
• body B gets \(-\lambda\)
• and if the anchor is off-center you also get torques

That's exactly "equal and opposite forces at a point."

8) What changes if you want different behaviour?

• If you want to lock orientation too, you add 3 angular constraints → fixed joint.
• If you want softness, you add compliance/damping terms → springy joint.
• If you want just a distance constraint, you reduce it to 1 row along the line between anchors.

Same machinery, different rule.

Fifth Gear: the same ball joint, but with actual numbers you can crunch

Now we do the exact same ball joint, but with real numbers you can work through.

To keep it beginner-friendly (and not turn into a 12-page algebra ambush), we'll pick "nice" values:

• both bodies have mass 1
• both bodies have identity inertia
• both bodies have identity orientation (so local = world for anchors)

That way the numbers you see are there because the physics demands them — not because a rotation matrix snuck in when you weren't looking.

0) What we're building

Two bodies, A and B, connected by a ball joint.

The rule:

"The anchor point on A and the anchor point on B must move together."

Velocity form:

\[
\dot{\mathbf{C}} = \mathbf{v}_{pA} - \mathbf{v}_{pB}
\]

We'll compute \(\lambda\) (the impulse vector) that fixes the joint motion.

1) Pick the body properties

Masses:

\[
m_A = 1,\quad m_B = 1
\Rightarrow m_A^{-1}=1,\quad m_B^{-1}=1
\]

Inverse inertia:

\[
I_A^{-1} = I_3,\quad I_B^{-1} = I_3
\]

So each body's inverse mass block is:

\[
M_A^{-1} =
\begin{bmatrix}
I_3 & 0 \\
0 & I_3
\end{bmatrix},
\quad
M_B^{-1} =
\begin{bmatrix}
I_3 & 0 \\
0 & I_3
\end{bmatrix}
\]

Global inverse mass matrix:

\[
M^{-1} =
\begin{bmatrix}
M_A^{-1} & 0 \\
0 & M_B^{-1}
\end{bmatrix}
\]

2) Choose the anchors

Anchors offset along the x-axis:

\[
\mathbf{r}_A = (0.5,\ 0,\ 0),\quad
\mathbf{r}_B = (-0.5,\ 0,\ 0)
\]

This is a great choice because it gives clean cross products, and you'll actually see rotation matter in y/z.

3) Choose a current state (velocities + a tiny position error)

Body A is spinning around z:

\[
\mathbf{v}_A=(0,0,0),\quad \boldsymbol{\omega}_A=(0,0,2)
\]

Body B is not moving:

\[
\mathbf{v}_B=(0,0,0),\quad \boldsymbol{\omega}_B=(0,0,0)
\]

Sanity check before we compute anything:
A is rotating, and its anchor is offset in +x, so that anchor point should have upward motion in +y. So you expect the joint to complain about y.

Now add a little drift (anchors slightly separated in y):

\[
\mathbf{C} = (0,\ 0.1,\ 0)
\]

Bias term:

\[
\text{bias} = \beta \frac{\mathbf{C}}{dt}
\]

Pick:

\[
dt = 1,\quad \beta = 0.2
\Rightarrow \text{bias} = (0,\ 0.02,\ 0)
\]

So we want a tiny "closing" behaviour in y to correct drift gently.

4) Write the velocity vector \(V\) (12 numbers)

Order:

\[
\mathbf{V}=
\begin{bmatrix}
\mathbf{v}_A \\ \boldsymbol{\omega}_A \\ \mathbf{v}_B \\ \boldsymbol{\omega}_B
\end{bmatrix}
=
\begin{bmatrix}
v_{Ax}\\v_{Ay}\\v_{Az}\\\omega_{Ax}\\\omega_{Ay}\\\omega_{Az}\\
v_{Bx}\\v_{By}\\v_{Bz}\\\omega_{Bx}\\\omega_{By}\\\omega_{Bz}
\end{bmatrix}
\]

With our values:

\[
\mathbf{V}=
\begin{bmatrix}
0\\0\\0\\0\\0\\2\\
0\\0\\0\\0\\0\\0
\end{bmatrix}
\]

5) Build the Jacobian \(J\) with numbers (and sanity check it)

We use:

\[
J =
\begin{bmatrix}
I_3 & -[\mathbf{r}_A]_\times & -I_3 & [\mathbf{r}_B]_\times
\end{bmatrix}
\]

Skew matrices.

For \(\mathbf{r}_A=(0.5,0,0)\):

\[
[\mathbf{r}_A]_\times =
\begin{bmatrix}
0 & 0 & 0\\
0 & 0 & -0.5\\
0 & 0.5 & 0
\end{bmatrix}
\]

For \(\mathbf{r}_B=(-0.5,0,0)\):

\[
[\mathbf{r}_B]_\times =
\begin{bmatrix}
0 & 0 & 0\\
0 & 0 & 0.5\\
0 & -0.5 & 0
\end{bmatrix}
\]

So:

\[
J=
\begin{bmatrix}
1&0&0&\;\;0&0&0&\;-1&0&0&\;\;0&0&0\\
0&1&0&\;\;0&0&0.5&\;0&-1&0&\;\;0&0&0.5\\
0&0&1&\;\;0&-0.5&0&\;0&0&-1&\;\;0&-0.5&0
\end{bmatrix}
\]

Sanity checks you should do (seriously — always do these):

• x row has no angular terms → correct, because the lever arm is along x.
• y row involves \(\omega_z\) with 0.5 → correct, z-rotation makes y motion at an x-offset point.
• z row involves \(\omega_y\) with -0.5 → correct, y-rotation makes z motion.

Those 0s and 0.5s are the physics showing up in the matrix.

6) Compute the constraint velocity and build \(b\)

\[
\dot{\mathbf{C}} = J\mathbf{V}
\]

Result:

\[
\dot{\mathbf{C}}=
\begin{bmatrix}
0\\
1\\
0
\end{bmatrix}
\]

Perfect: the complaint is in y, as expected.

Add bias:

\[
\mathbf{b} = \dot{\mathbf{C}} + \text{bias}
=
\begin{bmatrix}
0\\1\\0
\end{bmatrix}
+
\begin{bmatrix}
0\\0.02\\0
\end{bmatrix}
=
\begin{bmatrix}
0\\1.02\\0
\end{bmatrix}
\]

So the solver wants to cancel the y separation speed and add a tiny corrective pull.

7) Build \(A = J M^{-1} J^T\)

With our identity-ish \(M^{-1}\):

\[
A = J J^T
=
\begin{bmatrix}
2&0&0\\
0&2.5&0\\
0&0&2.5
\end{bmatrix}
\]

Read it like a human:

• x has 2 → both bodies translate, that's \(1+1\).
• y and z have 2.5 → translation plus some rotation leverage from the 0.5 offset.

8) Solve for \(\lambda\)

Equality joint solve:

\[
A\lambda + \mathbf{b} = 0
\Rightarrow A\lambda = -\mathbf{b}
\]

So:

\[
\lambda = -A^{-1}\mathbf{b}
\]

Diagonal solve:

\[
\lambda_x = -0/2 = 0
\]
\[
\lambda_y = -1.02/2.5 = -0.408
\]
\[
\lambda_z = -0/2.5 = 0
\]

So:

\[
\lambda=
\begin{bmatrix}
0\\-0.408\\0
\end{bmatrix}
\]

Meaning: the joint applies an impulse purely in y, and the sign is exactly what you'd want to oppose the upward relative motion.

9) Apply the impulse back to the bodies

\[
\Delta \mathbf{V} = M^{-1} J^T \lambda
\]

With identity-ish \(M^{-1}\):

\[
\Delta \mathbf{V} = J^T \lambda
\]

Result:

\[
\Delta \mathbf{V}=
\begin{bmatrix}
0\\-0.408\\0\\
0\\0\\-0.204\\
0\\+0.408\\0\\
0\\0\\-0.204
\end{bmatrix}
\]

Read it:

• A gets pushed down in y.
• B gets pushed up in y.
• Both get a z angular change because the impulse is applied at an offset (torque).

10) Verify the result

Post-solve relative anchor velocity becomes:

\[
\dot{\mathbf{C}}_{\text{post}} =
\begin{bmatrix}
0\\-0.02\\0
\end{bmatrix}
\]

That's not a mistake — it's the bias doing its job.

You asked for:

\[
\dot{\mathbf{C}} + \text{bias} = 0
\]

So the solver aims for:

\[
\dot{\mathbf{C}}_{\text{post}} = -\text{bias} = (0,-0.02,0)
\]

Meaning: a tiny closing speed to correct the drift gently over time.

11) Quick "other cases"

If body B is pinned (bolted to the world)

Set:

\[
m_B^{-1}=0,\quad I_B^{-1}=0
\]

Then the same joint impulse only affects A.
B becomes "the world," and doesn't budge.

If you only want to constrain one axis

Keep only (say) the y row of the Jacobian.
Now you've built a "slide in a plane/line" style constraint instead of a full ball joint.

If you want to lock rotation too (fixed joint)

Add 3 more constraints on relative angular motion/orientation.

Ball joint (3) + angular lock (3) = 6 constraints total → rigid weld.

That's how you move from "socket joint" to "glued together."

Only the Beginning

At this point, you should be getting an idea of what constraint solvers are and how they work! Of course, the theory is always different from the practice! For example, you'll have times when your matrix will contain invalid values or overlapping constraints (need to add small correction biases)... not to mention all the different types of Jacobian constraints, limits, ... and more!

The following gives you some full implemented prototypes (2D and 3D) for a constraint solver - in vanilla JavaScript that you can look at!

Project Code
• 2D Vanilla JavaScript and Canvas - 2D/Stacks,Cube,Cars,Motors [LINK]
• 3D Vanilla JavaScript and WebGL - 3D/Cubes/Constraints/Hinges/Ragdolls [LINK]

Advert (Support Website)

Visitor: