As promised, a (mostly) beginner-friendly thread on how to code simple physical simulations like those I've been posting recently.

Disclaimer: I am by no means an expert, I still have a lot to learn, and the methods I describe are among the simplest possible.

Don't be afraid of terminology, it is here so that you could google this stuff easily.

Of course, questions are welcome.

All units are SI.

I'm having 2D in mind but everything equally applies to 3D.

Use any programming language / engine / environment you like. You'll want to visualize the results, though.

I'd rather post it as a blog post, but I don't have a blog, so let's do a twitter thread.

Start with freely moving balls. The Newton equations of motion for a single particle are F=ma, or if we separate the velocity

v' = F/m
q' = v

where F is the force (a vector), m is mass, v is velocity (a vector), q is position (a vector as well), and ' denotes time derivative.

For many particles, we have an equation for each particle's position & velocity. We also need starting positions & velocities. This is called an initial-value problem for second-order ordinary differential equation. We are going to solve ("integrate") it numerically.

Popular methods for doing this are explicit Euler, Verlet, Runge-Kutta, etc. There are also implicit methods that (as far as I know) are hard to use for this kind of simulations, and semi-implicit methods that will come up a bit later.

I'll stick to explicit Euler for now. This method consists in using just the first two terms in the Taylor expansion: x(t+Δt) = x(t) + x'(t)⋅Δt.

So, if we have an equation x' = f(x), we do

x(t+Δt) = x(t) + f(x(t))⋅Δt

then

x(t+2Δt) = x(t+Δt) + f(x(t+Δt))⋅Δt

and so on.

For our case, 'x' is the vector of all particles' positions & velocities, and f contains forces (for v') and velocities (for q'). There's no need to store and handle 'x' explicitly, though. We do

v(t+Δt) = v(t) (since we have no forces yet)
q(t+Δt) = q(t) + q(t)⋅Δt

Thus, we store positions & velocities of all particles and update them using the formulas above over certain time period. What Δt to use and when to update the physics is entirely up to the application. I prefer to keep the virtual simulation time equal to real application time.

Sometimes I do a single physics update before each rendering frame, taking Δt = 1 / (frame rate). As we'll see later, smaller Δt are required sometimes. I've used Δt = 1 ms, thus doing ~16 physics updates before each rendering frame at 60 Hz refresh rate.

At this point, we have something like this:

Let's add some gravity. Gravity is the same for all objects, and points down (or wherever you prefer, really). I'm using a gravity value of 10 m/s². Add it to the update code:

v(t+Δt) = v(t) + g⋅Δt
q(t+Δt) = q(t) + v(t)⋅Δt

(g is the gravity vector)

Note that it is not the same as doing

v ⭠ v + g⋅Δt
q ⭠ q + v⋅Δt

which corresponds to

v(t+Δt) = v(t) + g⋅Δt
q(t+Δt) = q(t) + v(t+Δt)⋅Δt

This is a different method called semi-implicit Euler or symplectic Euler and actually conserves energy better.

From now on I'll stick to the assignment (⭠) notation, and we are using symplectic Euler scheme. With gravity, we get something like

Let's make some walls. For a wall to function properly, we need collision detection & collision response. Test each ball for collision with each wall, i.e. if the distance from the ball center to the wall is less than the ball radius.

Note that we need a signed distance here, i.e. for a wall at x = 0 we need x < radius, not abs(x) < radius, for the collision to be reported. Otherwise, we'll miss objects that happened to tunnel far through the wall somehow (e.g. due to high speed).

Suppose a collision occurs. We want the ball to bounce off the wall. Speaking of velocities, this means inverting the normal component of velocity and keeping the tangential component unchanged.

Let n be the inward normal of the wall. Reflecting the velocity through the wall can be done via the formula

v ⭠ v - 2⋅n⋅(v⋅n)

Here, (v⋅n) is the dot product of vectors. This part is the collision response. If multiple collisions occur, handle all of them sequentially.

Put the collision detection & response after the update code, to get

(notice that they all have the same bouncing frequency, since the acceleration is the same)

This collision is completely elastic, which is not that realistic. With many objects colliding this way & many other forces (e.g. undamped springs) you'll soon have an incomprehensible mess of stuff flying around & vibrating. We want the energy to dissipate (decrease) somehow.

Let's make the collisions non-elastic. This is as simple as changing the collision response to

v ⭠ v - (1+e)⋅n⋅(v⋅n)

Here, 'e' is the elasticity parameter. e=1 means completely elastic collision (normal velocity is reflected), e=0 is inelastic (normal velocity is killed).

With e=0.5 we get

A few more hints on collision. Firstly, it is a good idea to accept the collision only if the object is moving towards the wall, i.e. (v⋅n) < 0. Otherwise, an object "sunk" into the wall but already moving away would be considered colliding, it's velocity will get reflected, ...

... and it will start moving towards the wall. On next update, the velocity will be flipped again, and so on, flipping on each update and never really leaving the wall.

(the thread is in the process of being continued)

Secondly, I usually offset the object from the wall so that no collision occurs anymore. That is:

p ⭠ p + n⋅d

where 'd' is the penetration depth (e.g. d = -(x - radius) for a particle at [x,y] colliding with a wall x = 0). This prevents objects from sinking into the wall.

I forgot to mention that for simplicity I assume all masses to be equal to 1.

Let's do some springs. A spring connects two balls and has some predefined rest length. If its length is different, a force pulls the spring to return it to rest length, governed by Hooke's law.

So, for a pair of balls with positions p₁ and p₂, calculate the relative vector r=p₂-p₁, distance d=|r|, and normalized direction vector n=r/d. If the spring rest length is l, the force is

F = -K⋅n⋅(d - l)

E.g. if d > l, the force will contract the spring.

Since a spring connects two balls, split the force between them. Force affects velocity, so put this to the velocity update rule

v ⭠ v + (g + F/2)⋅Δt

If a single ball is connected to multiple springs, add all the spring forces to its velocity update.

With K=100, we get this (I disabled gravity to make the spring more apparent)

Note that it oscillates forever (unless it hits the wall, in which case it may transform spring oscillation energy into velocity or rotation). As I've said earlier, we want energy to dissipate. For springs, this means we want some damping.

Spring force accelerates the spring towards the rest length. Spring damping eases out this acceleration, gradually decreasing the oscillation. Suppose for a moment that the first ball is fixed p₁ = const. Then, the damping equation looks like

v₂' = -C⋅v₂

That is, the velocity decreases with time, proportional to the velocity magnitude and the damping constant C. With two balls we can't simply apply this equation to both velocities: this would violate impulse & torque conservation laws.

Instead, we need to take the relative velocity Δv = v₂-v₁ and decompose it into normal & tangential components, and apply the equation to normal velocity only.

Let n be the normalized direction vector from the spring equation. Normal velocity is

vₙ = n⋅(Δv⋅n)

Using explicit Euler for the damping equation, we get an impulse that should be added to the spring:

j = - C⋅vₙ⋅Δt

This impulse should be divided between the two balls, as well. The full spring update rule becomes

v ⭠ v + (g + F/2 - C⋅vₙ/2)⋅Δt

Again, for each spring connected to the ball, we get another (F/2 - C⋅vₙ/2) term.

With C=10, we get nice damping:

Notice how it stops oscillating but continues moving & rotating, thanks to conservation laws.

If you try larger values of K or C, you may find that your simulation becomes wildly unstable. Let's look closer at the velocity update, and assume p₁ = const again. Focus on damping:

v₂ ⭠ v₂ - C⋅vₙ/2⋅Δt

Assume vₙ = v₂. The rule rewrites as

v₂ ⭠ v₂⋅(1 - C⋅/2⋅Δt)

If it so happens that C⋅/2⋅Δt > 1, we get a problem: instead of decreasing the velocity, the rule flips it. If C⋅/2⋅Δt > 2, it even increases the magnitude of velocity. So, we really want C⋅Δt to be sufficiently smaller than 1.

In other words, we want Δt sufficiently smaller than 1/C.

The same kind of reasoning applies to the spring force, except that it involves a second order equation (q''=F/m, remember?) so it is more forgiving: we want (Δt)² smaller than 1/K.

With Δt = 0.001 (1 ms), K = 10000 and C = 200, I got this (I'm moving it with a mouse)

Notice that with K and C being that large, the spring appears to be completely stiff until pressed in the corner.

At this point, we can already create almost rigid structures:

And even freely rotating wheels (if we forbid one ball to move):

The balls don't bounce off the ground anymore. This is ok, though: the energy transforms into springs tension and rapidly dissipates away.

They also slide on the ground as if it was ice. Time to change that by adding friction.

Bouncing changes normal velocity. Friction, on the other hand, changes tangential velocity.

Let's revisit our bouncing equations. n is the inward normal to the wall. The normal velocity is

vₙ = n⋅(v⋅n)

and the tangential velocity is all that's left

vₜ = v - vₙ

Friction works similarly to damping: the equation is

vₜ' = -f⋅vₜ

where f is the friction coefficient. Collision response transforms into

v ⭠ v - (1+e)⋅vₙ - f⋅vₜ⋅Δt

With f = 100, we get somewhat realistic ground friction:

It still feels a bit unnatural, though. The reason is that friction is actually tied to bounce impulse through Coulomb's law of friction. For more details on that, see https://gamedevelopment.tutsplus.com/tutorials/how-to-create-a-custom-2d-physics-engine-friction-scene-and-jump-table--gamedev-7756

Now, let's do ball-ball collisions. In fact, the collision response code is almost the same as for ball-wall collisions, only the detection phase is different.

Two balls collide if the distance between their centers is less then the sum of their radii.

If collision occurs, we need to remember about conservation laws again. Take the normal direction between the balls

n = (p₂-p₁)/|p₂-p₁|

and the relative velocity

Δv = v₂-v₁

As before, we have normal & tangential velocity:

vₙ = n⋅(Δv⋅n)
vₜ = v - vₙ

Bounce the balls back & add friction, be careful with the signs, and don't forget to split the impulse between the two balls:

v₁ ⭠ v₁ + vₙ⋅(1+e)/2 + f⋅vₜ⋅Δt/2
v₂ ⭠ v₂ - vₙ⋅(1+e)/2 - f⋅vₜ⋅Δt/2

The bouncing part is more important at the beginning, so feel free to test with f=0.

To get the signs right, assume ball 1 fixed and ball 2 moving directly into ball 1, draw all the vectors on a piece of paper, and think of what should happen to the velocities after the impact.

You can also offset the balls to resolve the collision:

p₁ ⭠ p₁ - n⋅d/2
p₂ ⭠ p₂ + n⋅d/2

where d is penetration depth

d = 2⋅radius - |p₂-p₁|

We should get real collisions now:

We have a ton of different velocity & position updates coming from various collisions, forces, and dampings, applied at some order, altering the effect of each other. Collisions may resolve differently based on balls order. At this point it's a good idea to reorganize everything.

We start the update routine with old values of positions & velocities. For each particle, track the impulse (velocity change) and offset (position change) applied at this update. All calculations use the old positions & velocities, and add to this cumulative impulses & offsets.

At the end of the update routine, actually add the impulses & offsets to the velocities & positions.

Ball-stick collisions are a bit trickier. A decent physics engine would trick the sticks as real solid objects, with their own positions, velocities, and orientations. We only have sticks as virtual connection between balls, so we have to hack a bit.

First, collision detection. To get the distance between a ball 1 & a stick connecting balls 2-3, find the normalized stick direction vector

Δq = (q₃-q₂) / |q₃-q₂|

and rotate it 90° by the transformation

(x,y) ⭢ (y,-x)

This way you get the stick normal vector n. The abs of dot product n⋅(q₁-q₂) is the distance between the stick 2-3 and the ball 1. A collision occurs if this distance is less than ball radius.

Collision response happens the same way as in the ball-ball case, but treating the stick as a virtual ball of mass 2 and interpolated velocity. That is, if the collision happens closer to ball 2, the virtual ball velocity is closer to v₂.

Compute t = Δq⋅(q₁-q₂)/|q₃-q₂| to get the interpolation parameter.

Having mass 2 means that the stick gets only 1/3 of the full collision impulse, and the ball gets the other 2/3. Also, if the collision happens closer to ball 2, the impulse should apply to it more.

Let j be the full collision impulse (including bounce & friction). The update rule becomes

v₁ ⭠ v₁ + j⋅2/3
v₂ ⭠ v₂ - j⋅(1-t)/3
v₃ ⭠ v₃ - j⋅t/3

(the signs may differ based on how you define j)

This is hacky but the result looks plausible:

Now you can make different types of balls & springs, e.g. non-colliding springs, and make a ton of interesting simulations. The end!

For a more in-depth source, see Physics-Based Animation by Erleben et. al, or just google "game physics engine tutorial", there are thousands of them!