Constrained Dynamics (I)

I am starting a mini-series of posts that go through the mathematical derivation necessary to build a minimal physics engine with support for constrained dynamics.

This first post is about unconstrained dynamics, using a free-falling particle as an example.

Upcoming posts will describe constrained dynamics, then build a (single-constraint) pendulum, and finally a (multi-constraint) double-pendulum.

We will restrict ourselves to:

2D.
Only one particle mass.
Only one external force (gravity).
Only one type of constraint (distance).

Rigid bodies with a surface area, collisions, contacts, frictions, etc… will be left out of the picture.

Incredibly useful links

This mini-series is nothing but a digest of knowledge that can be found in great detail here:

Erin Catto’s Box2D publications area.
Physically Based Modeling: Principles and Practice by Andrew Witkin and David Baraff.

In particular:

This presentation by Erin Catto.
These lecture notes by Andrew Witkin.

Unconstrained dynamics

By way of newtonian mechanics:

1st Law: The velocity (or lack thereof) of a particle remains unchanged until a force is applied to it.
2nd Law: When a force is applied to a particle, it experiments an acceleration of magnitude $\frac{‖ F ‖}{m}$ in the direction of $F$ .

As a consequence of the 1st Law, we can describe the state $x$ of a particle of mass $m$ by its position and its velocity.

As a consequence of the 2nd Law, we can’t consider acceleration to be part of the state as it only spawns for as long as an external force is applied to the particle. Force (acceleration) is actually the only agent that causes changes of state.

These three magnitudes are 2D vectors. We will use the following notation:

\begin{aligned} p = p o s i t i o n \\ \dot{p} = \frac{d p}{d t} = v e l o c i t y \\ \ddot{p} = \frac{d \dot{p}}{d t} = \frac{d^{2} p}{d t^{2}} = a c c e l e r a t i o n \end{aligned}

Momentum

Mass and velocity together lead to the definition of linear momentum:

ρ = m \dot{p}

Another way to read the 2nd Law is that a force acting on a particle induces a change on its momentum.

F = \frac{d ρ}{d t} = \frac{d (m \dot{p})}{d t} = m \ddot{p}

We can reason that momentum is a measure of the impulse (force over time) that one must exert in order to stop the particle (i.e., cancel out its current velocity).

For a steady force, this is:

\begin{aligned} Imp = \int_{t_{0}}^{t_{1}} F d t = F Δ t \\ Imp = \int_{t_{0}}^{t_{1}} \frac{d ρ}{d t} d t = Δ ρ \end{aligned}

State over time

The evolution of the particle’s position is so ruled by the differential equation:

p = p_{0} + \int_{t_{0}}^{t_{1}} ({\dot{p}}_{0} + \int_{t_{0}}^{t_{1}} \frac{F}{m} d t) d t

Which can be read as: the particle position is a function of the initial state $x_{0}$ (at $t = t_{0}$ ) and the forces applied to the particle, which cause acceleration, which modifies the velocity, which modifies the position.

The classic and most intuitive way to iteratively solve the above diffential equation is the Semi-Implicit Euler method:

\begin{aligned} {\ddot{p}}_{t} = \frac{F}{m} \\ {\dot{p}}_{t + Δ t} = {\dot{p}}_{t} + {\ddot{p}}_{t} Δ t \\ p_{t + Δ t} = p_{t} + {\dot{p}}_{t} Δ t \end{aligned}

Where each iteration is meant to encompass a discrete time slice of duration $Δ t$ .

These expressions can be read as:

Calc the acceleration caused by the total sum of applied forces at time $t$ .
Update the particle’s velocity with said acceleration (times $Δ t$ ).
Update the particle’s position with the updated velocity (times $Δ t$ ).

The intricacies of numerical solvers for differential equations are outside the scope of this post. But the above method does a fine job at producing a discrete approximation to the true state over time. At least, as long as $Δ t$ is small enough. For starters $Δ t = 1 / 60$ (60 steps per second) is a sensible choice.

Note that the smaller the $Δ t$ , the more closely that we will be approximating the exact analytic solution. More accuracy sounds ideal, but the smaller the delta, the more the steps necessary to integrate over the same total time. Also, tiny deltas may drown in the muddy waters of insufficient floating-point precision. So a value that is small, but not too small, is the right choice.

Free-fall (iterative)

We’re only considering gravity in this post ( $F = g$ ). So the only acceleration is downwards and always steady: $\ddot{p} = \frac{g}{m}$ . The simulation steps become:

\begin{aligned} {\dot{p}}_{t + Δ t} = {\dot{p}}_{t} + \frac{g}{m} Δ t \\ p_{t + Δ t} = {\dot{p}}_{t} Δ t \end{aligned}

Here’s a simple implementation for the particle state and the simulation step:

struct particle_t
{
  explicit particle_t() : m( 1 ), p( 0 ), v( 0 ) {}

  void step( const f64_t dt, const vec2d_c& F )
  {
    const vec2d_c a = ( F / m );

    v += ( a * dt );
    p += ( v * dt );
  }

    f64_t m;  // Mass.
  vec2d_c p;  // Position.
  vec2d_c v;  // Velocity.
};

Here’s a humble animation of the particle in free-fall:

Projectile motion

If we give an initial velocity to the particle ( $t = 0$ ), we end up with projectile motion:

Projectile motion

Free-fall (analytic)

Differential equations are outside the scope of this post. But let’s say at least that an analytic solution to this example, where there’s only one steady external force, is possible and easy.

The analytic solution gives us a function that explicitly describes the position of the particle at any given time. This sounds better than an iterative simulation, although analytic solutions are untractable in non-trivial systems.

The vertical force that gravity exerts on a particle is proportional to its mass, which in turn makes gravity accelerate any particle the same regardless of its mass:

g = [\begin{matrix} 0 \\ m g_{y} \end{matrix}]

Plugging $g$ into the differential equation above:

\begin{aligned} x = x_{0} + \int_{t_{0}}^{t_{1}} ({\dot{x}}_{0} + \int_{t_{0}}^{t_{1}} \frac{0}{m} d t) d t \\ y = y_{0} + \int_{t_{0}}^{t_{1}} ({\dot{y}}_{0} + \int_{t_{0}}^{t_{1}} \frac{m g_{y}}{m} d t) d t \end{aligned}

Which analytic solution is the very familiar projectile motion formula:

\begin{aligned} x = x_{0} + {\dot{x}}_{0} t \\ y = y_{0} + {\dot{y}}_{0} t + g_{y} t^{2} \end{aligned}

This type of motion is also called parabolic throw. As can be seen both in the above animation and by looking at the formula, the ballistic trajectory w.r.t. time is horizontally linear, and vertically parabolic.

Things are about to get exciting. I promise!