P04: Physics Simulation with Constraints

Build a 2D physics engine that simulates rigid body dynamics with constraints (springs, joints, collisions)—requiring you to solve systems of linear equations each frame.

Overview

Attribute	Value
Difficulty	Advanced
Time Estimate	3-4 weeks
Language	C (recommended) or Python
Prerequisites	Basic physics intuition, matrix operations
Primary Book	“Game Physics Engine Development” by Ian Millington

Learning Objectives

By completing this project, you will:

See vectors as physical state - Position, velocity, and acceleration as mathematical objects
Master Ax = b solvers - Gaussian elimination and iterative methods
Understand projections geometrically - Collision response via vector projection
Grasp matrix conditioning - Why some systems are “hard” to solve
Apply null spaces and least squares - Handle over/under-constrained systems
Experience numerical stability - When math works on paper but fails in code

Theoretical Foundation

Part 1: Vectors as Physical State

Position, Velocity, Acceleration

Every physics simulation begins with vectors representing physical quantities:

┌────────────────────────────────────────────────────────────────┐
│                    PHYSICS VECTORS                             │
├────────────────────────────────────────────────────────────────┤
│                                                                │
│    Position (p)     Velocity (v)     Acceleration (a)         │
│    "where"          "how fast"       "how much speeding up"   │
│                                                                │
│    p = [x]          v = dp/dt        a = dv/dt = d²p/dt²      │
│        [y]                                                     │
│                                                                │
│    For a 2D point:                                             │
│    p = (3, 5) means 3 units right, 5 units up                 │
│    v = (1, 0) means moving right at 1 unit/second             │
│    a = (0, -9.8) means accelerating downward (gravity)        │
│                                                                │
└────────────────────────────────────────────────────────────────┘

Physics Vectors

Newton’s Second Law in Vector Form

F = ma

Force vector → acceleration vector → velocity change → position change

For multiple forces:
a = (1/m) × ΣF = (1/m) × (F_gravity + F_spring + F_collision + ...)

Each force is a vector; we sum them.

The State Vector

For a system of n particles, collect all positions and velocities:

State vector S:

    [p₁]     position of particle 1
    [p₂]     position of particle 2
S = [...]
    [v₁]     velocity of particle 1
    [v₂]     velocity of particle 2
    [...]

For n particles in 2D: S has 4n components (2 position + 2 velocity each)

This state vector evolves over time according to physics equations.

Part 2: Numerical Integration

The Euler Method

The simplest integration scheme:

v_new = v_old + a × Δt
p_new = p_old + v_old × Δt

Or more accurately (semi-implicit Euler):
v_new = v_old + a × Δt
p_new = p_old + v_new × Δt  ← use new velocity

Euler is simple but unstable - errors accumulate, springs explode.

Why Integration is Tricky

┌────────────────────────────────────────────────────────────────┐
│                    INTEGRATION ERROR                           │
├────────────────────────────────────────────────────────────────┤
│                                                                │
│    True path                   Euler approximation             │
│                                                                │
│       ╭──────╮                    ╭────────────╮               │
│      ╱        ╲                  ╱              ╲              │
│     ╱          ╲                ╱                ╲             │
│    ╱            ╲              ╱                  ╲            │
│                                      ↑                         │
│                           Each step overshoots                 │
│                           Error accumulates                    │
│                                                                │
│   For springs/pendulums: Energy grows → explosion!            │
│                                                                │
└────────────────────────────────────────────────────────────────┘

Integration Error

Better Methods: Verlet Integration

Position-based method that conserves energy better:

p_new = 2×p_current - p_old + a × Δt²

Velocity is implicit: v ≈ (p_new - p_old) / (2×Δt)

Advantages:
- Second-order accurate (error ∝ Δt²)
- Energy-conserving (symplectic)
- Simple to implement

RK4 (Runge-Kutta 4th Order)

More accurate but expensive:

k₁ = f(t, y)
k₂ = f(t + Δt/2, y + Δt×k₁/2)
k₃ = f(t + Δt/2, y + Δt×k₂/2)
k₄ = f(t + Δt, y + Δt×k₃)

y_new = y + (Δt/6) × (k₁ + 2k₂ + 2k₃ + k₄)

For this project, semi-implicit Euler or Verlet is sufficient.

Part 3: Collision Detection and Response

Detecting Collisions

For circles, collision detection is simple:

Circle A: center pA, radius rA
Circle B: center pB, radius rB

Distance: d = ||pA - pB||
Collision: d < rA + rB

Penetration depth: δ = (rA + rB) - d
Collision normal: n = (pB - pA) / d

For boxes, polygons, etc., it gets more complex (SAT, GJK).

Collision Response via Projection

The key linear algebra insight: collision response is projection!

┌────────────────────────────────────────────────────────────────┐
│                 COLLISION RESPONSE                             │
├────────────────────────────────────────────────────────────────┤
│                                                                │
│    Before collision:                                           │
│                                                                │
│         v (incoming velocity)                                  │
│          ↘                                                     │
│           ↘                                                    │
│    ─────────●───────── surface (normal n)                      │
│                                                                │
│    Decompose v into components:                                │
│                                                                │
│    v_normal = (v · n) × n     (component along normal)         │
│    v_tangent = v - v_normal   (component along surface)        │
│                                                                │
│    For elastic bounce:                                         │
│    v_new = v_tangent - e × v_normal  (e = coefficient of      │
│                                        restitution)            │
│                                                                │
│    For e = 1 (perfectly elastic):                              │
│    v_new = v - 2(v · n)n      (reflection formula)            │
│                                                                │
└────────────────────────────────────────────────────────────────┘

Collision Response - Vector Projection

The math:

Reflection: v_reflected = v - 2(v · n)n

This is the formula for reflecting a vector across a plane!
The dot product measures how much v points into the surface.
Multiplying by n and subtracting reverses that component.

Impulse-Based Collision

For realistic physics, we compute impulses (instantaneous momentum changes):

For two colliding bodies A and B:

Relative velocity: v_rel = vA - vB
Normal component:  v_n = v_rel · n

If v_n > 0: bodies separating, no collision response needed

Impulse magnitude: j = -(1 + e) × v_n / (1/mA + 1/mB)

Apply impulse:
vA_new = vA + (j/mA) × n
vB_new = vB - (j/mB) × n

Part 4: Constraints and Solving Ax = b

What is a Constraint?

A constraint is a restriction on how bodies can move:

┌────────────────────────────────────────────────────────────────┐
│                    TYPES OF CONSTRAINTS                        │
├────────────────────────────────────────────────────────────────┤
│                                                                │
│    Distance constraint (rod/rope):                             │
│    ||pA - pB|| = L   (maintain fixed distance)                 │
│                                                                │
│    Point constraint (pin joint):                               │
│    pA = pB           (two points coincide)                     │
│                                                                │
│    Angle constraint (hinge limit):                             │
│    θ_min ≤ θ ≤ θ_max (rotation bounded)                        │
│                                                                │
│    Contact constraint (no penetration):                        │
│    (pA - pB) · n ≥ 0  (inequality constraint)                  │
│                                                                │
└────────────────────────────────────────────────────────────────┘

Types of Constraints

Constraint Equations

Most constraints can be written as: C(x) = 0

For a distance constraint between particles at p₁ and p₂:

C = ||p₁ - p₂||² - L²  = 0

Taking the derivative (for velocity-level constraint):

Ċ = 2(p₁ - p₂) · (v₁ - v₂) = 0

This says: relative velocity must be perpendicular to the connecting line.

The Jacobian Matrix

For multiple constraints on multiple bodies, we form the Jacobian:

C = [C₁]      Jacobian J = [∂C₁/∂x₁  ∂C₁/∂x₂  ...]
    [C₂]                   [∂C₂/∂x₁  ∂C₂/∂x₂  ...]
    [...]                  [  ...      ...     ...]

The constraint becomes: J × v = 0  (or = -baumgarte × C for stability)

Where v is the velocity vector of all bodies.

Solving for Constraint Forces

To satisfy constraints, we need forces. The constraint force λ:

Total acceleration: a = M⁻¹(F_external + Jᵀλ)

Constraint: J × a = b  (where b comes from the constraint)

Substituting:
J × M⁻¹ × Jᵀ × λ = b - J × M⁻¹ × F_external

This is a system:  A × λ = b

Where A = J × M⁻¹ × Jᵀ (the constraint mass matrix)

This is the heart of constraint solving: solve Ax = b to find λ, then apply forces.

Part 5: Solving Systems of Linear Equations

The General Problem

A × x = b

Given: A (n×n matrix), b (n vector)
Find:  x (n vector)

For physics:

A encodes how constraints couple to each other
b encodes how much the constraints are violated
x (λ) gives the correction forces

Gaussian Elimination

The classic direct method:

Step 1: Forward elimination (create upper triangular form)
Step 2: Back substitution (solve from bottom up)

Example:
[2  1  1 | 8]     [2  1  1 | 8]     [2  1  1 | 8]
[4  3  3 | 20] →  [0  1  1 | 4] →   [0  1  1 | 4]
[8  7  9 | 46]    [0  3  5 | 14]    [0  0  2 | 2]

Back substitute: z = 1, y = 3, x = 2

Complexity: O(n³) - too slow for large systems every frame!

LU Decomposition

Factor A once, solve quickly for different b:

A = L × U  (Lower triangular × Upper triangular)

Solve: L × y = b  (forward substitution)
Then:  U × x = y  (back substitution)

Both steps are O(n²).

Useful when A doesn’t change (same constraints, different forces).

Iterative Methods: Gauss-Seidel

For constraint solving, iterative methods work well:

def gauss_seidel(A, b, x, iterations):
    for _ in range(iterations):
        for i in range(n):
            sigma = sum(A[i][j] * x[j] for j in range(n) if j != i)
            x[i] = (b[i] - sigma) / A[i][i]
    return x

Why iterative works for physics:

Constraints are local (sparse matrix)
Approximate solution is often good enough
We solve many times per second (converges over frames)

Projected Gauss-Seidel (PGS)

For inequality constraints (contacts), clamp solutions:

def pgs(A, b, x, lo, hi, iterations):
    for _ in range(iterations):
        for i in range(n):
            delta = (b[i] - sum(A[i][j] * x[j] for j in range(n))) / A[i][i]
            x[i] = clamp(x[i] + delta, lo[i], hi[i])
    return x

This handles: “bodies can push apart but not pull together”

Part 6: Matrix Conditioning and Numerical Stability

What is Condition Number?

The condition number κ(A) measures how sensitive Ax = b is to perturbations:

If κ(A) is small (close to 1): well-conditioned, stable
If κ(A) is large: ill-conditioned, small errors → large errors

κ(A) = ||A|| × ||A⁻¹|| = σ_max / σ_min  (ratio of singular values)

Why Conditioning Matters for Physics

┌────────────────────────────────────────────────────────────────┐
│                ILL-CONDITIONED SYSTEMS                         │
├────────────────────────────────────────────────────────────────┤
│                                                                │
│    Scenario: Very stiff springs mixed with soft springs        │
│                                                                │
│    Stiff spring: k = 10000                                     │
│    Soft spring:  k = 1                                         │
│                                                                │
│    Constraint matrix eigenvalues: 10000 and 1                  │
│    Condition number: ~10000                                    │
│                                                                │
│    Result: Gauss-Seidel converges slowly                       │
│            Numerical precision issues                          │
│            Simulation may become unstable                      │
│                                                                │
│    Fixes:                                                      │
│    - Scale units to reduce spread                              │
│    - Use regularization (add small diagonal)                   │
│    - Use higher precision (double instead of float)            │
│    - More iterations                                           │
│                                                                │
└────────────────────────────────────────────────────────────────┘

Ill-Conditioned Systems

Numerical Stability Tips

Use double precision for constraint solving
Add bias/regularization: A + εI instead of A
Warm starting: Use previous frame’s solution as initial guess
Scale appropriately: Keep values in reasonable ranges

Part 7: Special Matrix Structures

Sparse Matrices

Physics matrices are typically sparse (mostly zeros):

For a chain of particles connected by springs:

[×  ×  0  0  0  0]    × = non-zero
[×  ×  ×  0  0  0]    0 = zero
[0  ×  ×  ×  0  0]
[0  0  ×  ×  ×  0]    Each particle only interacts with neighbors.
[0  0  0  ×  ×  ×]
[0  0  0  0  ×  ×]

Store only non-zeros → memory efficient
Skip zero multiplications → faster solving

Symmetric Positive Definite (SPD)

Many physics matrices are SPD:

Symmetric: A = Aᵀ
Positive definite: xᵀAx > 0 for all x ≠ 0

SPD matrices have:

All positive eigenvalues
Cholesky factorization: A = LLᵀ (faster than LU)
Guaranteed convergence for iterative methods

Block Structure

For rigid bodies, matrices have block structure:

[M₁  0   0 ]   [v₁]   [F₁]
[0   M₂  0 ] × [v₂] = [F₂]
[0   0   M₃]   [v₃]   [F₃]

Each Mᵢ is a 3×3 (or 2×2 for 2D) mass/inertia matrix.
Block-diagonal structure allows efficient solving.

Part 8: Least Squares and Overdetermined Systems

When Ax = b Has No Exact Solution

If we have more constraints than degrees of freedom:

3 constraints, 2 unknowns:
[1  0]       [3]
[0  1] × x = [4]
[1  1]       [8]   ← inconsistent!

No x satisfies all three equations exactly.

The Least Squares Solution

Minimize

Ax - b

²:

x* = (AᵀA)⁻¹ Aᵀ b  (the normal equations)

This finds x that minimizes the sum of squared residuals.

For physics: “satisfy constraints as well as possible”

Underdetermined Systems (Null Space)

More unknowns than constraints:

1 constraint, 2 unknowns:
[1  1] × x = [5]

Infinite solutions: x = [5-t, t] for any t

The null space of A contains all directions of "free motion"
that don't violate constraints.

For physics: A system with slack has infinite valid configurations.

Part 9: Putting It Together - The Physics Loop

┌─────────────────────────────────────────────────────────────────────────┐
│                    PHYSICS SIMULATION LOOP                             │
├─────────────────────────────────────────────────────────────────────────┤
│                                                                         │
│  For each time step:                                                   │
│  ───────────────────                                                   │
│                                                                         │
│  ┌─────────────────────────────────────────────────────────────────┐   │
│  │ 1. APPLY EXTERNAL FORCES                                         │   │
│  │    F_total = F_gravity + F_user_input + F_air_resistance + ...   │   │
│  └──────────────────────────────────┬──────────────────────────────┘   │
│                                     │                                   │
│                                     ▼                                   │
│  ┌─────────────────────────────────────────────────────────────────┐   │
│  │ 2. INTEGRATE (predict new positions/velocities)                  │   │
│  │    v_predicted = v + (F/m) × Δt                                  │   │
│  │    p_predicted = p + v_predicted × Δt                            │   │
│  └──────────────────────────────────┬──────────────────────────────┘   │
│                                     │                                   │
│                                     ▼                                   │
│  ┌─────────────────────────────────────────────────────────────────┐   │
│  │ 3. DETECT COLLISIONS                                             │   │
│  │    Find all overlapping/intersecting pairs                       │   │
│  │    Compute collision points, normals, depths                     │   │
│  └──────────────────────────────────┬──────────────────────────────┘   │
│                                     │                                   │
│                                     ▼                                   │
│  ┌─────────────────────────────────────────────────────────────────┐   │
│  │ 4. BUILD CONSTRAINT SYSTEM                                       │   │
│  │    Assemble Jacobian J from all constraints                      │   │
│  │    Form A = J M⁻¹ Jᵀ  and  b = -J v - bias × C                  │   │
│  └──────────────────────────────────┬──────────────────────────────┘   │
│                                     │                                   │
│                                     ▼                                   │
│  ┌─────────────────────────────────────────────────────────────────┐   │
│  │ 5. SOLVE FOR CONSTRAINT FORCES                                   │   │
│  │    Solve A λ = b using Gauss-Seidel or PGS                       │   │
│  │    (possibly with warm starting from previous frame)             │   │
│  └──────────────────────────────────┬──────────────────────────────┘   │
│                                     │                                   │
│                                     ▼                                   │
│  ┌─────────────────────────────────────────────────────────────────┐   │
│  │ 6. APPLY CONSTRAINT IMPULSES                                     │   │
│  │    v = v + M⁻¹ Jᵀ λ                                              │   │
│  │    (constraints now satisfied)                                   │   │
│  └──────────────────────────────────┬──────────────────────────────┘   │
│                                     │                                   │
│                                     ▼                                   │
│  ┌─────────────────────────────────────────────────────────────────┐   │
│  │ 7. UPDATE STATE                                                  │   │
│  │    p = p + v × Δt  (with corrected velocity)                     │   │
│  │    Store for next frame                                          │   │
│  └─────────────────────────────────────────────────────────────────┘   │
│                                                                         │
└─────────────────────────────────────────────────────────────────────────┘

Physics Simulation Loop

Project Specification

What You’re Building

A 2D physics simulator that:

Simulates particles/circles under gravity
Handles collisions with walls and other objects
Supports spring/distance constraints
Supports joint/pin constraints
Renders in real-time

Minimum Requirements

Particle system with gravity
Semi-implicit Euler or Verlet integration
Circle-circle collision detection and response
Circle-wall collision (axis-aligned boundaries)
Distance constraints (springs/rods)
Simple constraint solver (Gauss-Seidel)
Real-time rendering (SDL2 or similar)

Stretch Goals

Rigid body rotation (not just particles)
Friction in collisions
Multiple joint types (hinges, sliders)
Stacking and piles of objects
Continuous collision detection
Broad-phase acceleration (spatial hashing)

Solution Architecture

Data Structures

typedef struct {
    Vec2 position;
    Vec2 velocity;
    float mass;
    float inv_mass;  // 1/mass, or 0 for static objects
    float radius;
} Particle;

typedef struct {
    int particle_a;
    int particle_b;
    float rest_length;
    float stiffness;
} DistanceConstraint;

typedef struct {
    int particle_a;
    int particle_b;
    Vec2 normal;
    float penetration;
} Contact;

typedef struct {
    Particle* particles;
    int num_particles;
    DistanceConstraint* constraints;
    int num_constraints;
    Contact* contacts;
    int num_contacts;
    Vec2 gravity;
    float dt;
} PhysicsWorld;

Module Structure

┌─────────────────────────────────────────────────────────────────┐
│                           main.c                                │
│  - Initialize world                                             │
│  - Main loop: update physics, render                           │
│  - Handle input                                                 │
└─────────────────────────────────────────────────────────────────┘
                               │
        ┌──────────────────────┼──────────────────────┐
        ▼                      ▼                      ▼
┌──────────────┐      ┌──────────────┐      ┌──────────────┐
│   vec.c      │      │  physics.c   │      │  render.c    │
│              │      │              │      │              │
│ vec2_add     │      │ apply_forces │      │ draw_circle  │
│ vec2_sub     │      │ integrate    │      │ draw_line    │
│ vec2_scale   │      │ detect_collisions │ │ draw_world   │
│ vec2_dot     │      │ solve_constraints  │ │              │
│ vec2_length  │      │ step_world   │      │              │
│ vec2_normalize│     │              │      │              │
└──────────────┘      └──────────────┘      └──────────────┘
                               │
        ┌──────────────────────┼──────────────────────┐
        ▼                      ▼                      ▼
┌──────────────┐      ┌──────────────┐      ┌──────────────┐
│  collision.c │      │  constraint.c│      │  solver.c    │
│              │      │              │      │              │
│ circle_circle│      │ distance_init│      │ gauss_seidel │
│ circle_wall  │      │ distance_solve│     │ build_matrix │
│ find_contacts│      │ contact_solve│      │ solve_system │
└──────────────┘      └──────────────┘      └──────────────┘

Physics Engine Code Architecture

Simulation Loop Flow

┌─────────────────────────────────────────────────────────────────────────┐
│                      PHYSICS STEP IMPLEMENTATION                       │
├─────────────────────────────────────────────────────────────────────────┤
│                                                                         │
│  void physics_step(PhysicsWorld* world) {                              │
│                                                                         │
│      // 1. Apply external forces (gravity)                             │
│      for each particle p:                                              │
│          p.velocity += world.gravity * world.dt                        │
│                                                                         │
│      // 2. Predict positions (integrate)                               │
│      for each particle p:                                              │
│          p.position += p.velocity * world.dt                           │
│                                                                         │
│      // 3. Detect collisions                                           │
│      world.num_contacts = 0                                            │
│      detect_particle_collisions(world)                                 │
│      detect_wall_collisions(world)                                     │
│                                                                         │
│      // 4. Solve constraints iteratively                               │
│      for (int iter = 0; iter < NUM_ITERATIONS; iter++) {               │
│          // Solve distance constraints                                 │
│          for each constraint c:                                        │
│              solve_distance_constraint(c, world.particles)             │
│                                                                         │
│          // Solve contact constraints                                  │
│          for each contact c:                                           │
│              solve_contact(c, world.particles)                         │
│      }                                                                  │
│                                                                         │
│      // 5. Derive velocities from position changes (for Verlet)        │
│      for each particle p:                                              │
│          p.velocity = (p.position - p.old_position) / world.dt        │
│  }                                                                      │
│                                                                         │
└─────────────────────────────────────────────────────────────────────────┘

Physics Step Implementation

Phased Implementation Guide

Phase 1: Basic Particle Simulation (Days 1-4)

Goal: Particles falling under gravity

Tasks:

Implement Vec2 math library
Create Particle structure
Implement semi-implicit Euler integration
Add gravity
Render particles with SDL2

Implementation:

void integrate_particles(PhysicsWorld* world) {
    float dt = world->dt;
    for (int i = 0; i < world->num_particles; i++) {
        Particle* p = &world->particles[i];
        // Apply gravity
        p->velocity = vec2_add(p->velocity, vec2_scale(world->gravity, dt));
        // Update position
        p->position = vec2_add(p->position, vec2_scale(p->velocity, dt));
    }
}

Verification:

Particles should fall downward
Velocity should increase over time
Position change should match v × Δt

Phase 2: Wall Collisions (Days 5-7)

Goal: Particles bounce off screen boundaries

Tasks:

Detect when particle exits bounds
Compute collision normal
Implement reflection formula
Add coefficient of restitution

Implementation:

void collide_with_walls(Particle* p, float width, float height, float restitution) {
    // Left wall
    if (p->position.x - p->radius < 0) {
        p->position.x = p->radius;
        p->velocity.x = -p->velocity.x * restitution;
    }
    // Right wall
    if (p->position.x + p->radius > width) {
        p->position.x = width - p->radius;
        p->velocity.x = -p->velocity.x * restitution;
    }
    // Similarly for top/bottom
}

Verification:

Particles bounce off walls
Restitution = 1: energy conserved
Restitution = 0.5: each bounce loses energy

Phase 3: Particle-Particle Collisions (Days 8-11)

Goal: Particles collide with each other

Tasks:

Detect circle-circle overlap
Compute collision normal and penetration
Implement impulse-based response
Handle position correction

Implementation:

void detect_circle_collision(Particle* a, Particle* b, Contact* contact) {
    Vec2 delta = vec2_sub(b->position, a->position);
    float dist = vec2_length(delta);
    float min_dist = a->radius + b->radius;

    if (dist < min_dist) {
        contact->normal = vec2_scale(delta, 1.0f / dist);
        contact->penetration = min_dist - dist;
        contact->valid = true;
    } else {
        contact->valid = false;
    }
}

void resolve_collision(Particle* a, Particle* b, Contact* c, float restitution) {
    Vec2 rel_vel = vec2_sub(a->velocity, b->velocity);
    float vel_along_normal = vec2_dot(rel_vel, c->normal);

    // Don't resolve if separating
    if (vel_along_normal > 0) return;

    // Compute impulse
    float j = -(1 + restitution) * vel_along_normal;
    j /= a->inv_mass + b->inv_mass;

    Vec2 impulse = vec2_scale(c->normal, j);
    a->velocity = vec2_add(a->velocity, vec2_scale(impulse, a->inv_mass));
    b->velocity = vec2_sub(b->velocity, vec2_scale(impulse, b->inv_mass));

    // Position correction (prevent sinking)
    float correction = c->penetration / (a->inv_mass + b->inv_mass) * 0.8f;
    a->position = vec2_sub(a->position, vec2_scale(c->normal, a->inv_mass * correction));
    b->position = vec2_add(b->position, vec2_scale(c->normal, b->inv_mass * correction));
}

Verification:

Particles bounce off each other
Momentum is conserved
No permanent interpenetration

Phase 4: Distance Constraints (Days 12-15)

Goal: Connect particles with springs/rods

Tasks:

Define distance constraint structure
Compute constraint violation
Implement iterative position correction
Visualize constraints

Implementation:

void solve_distance_constraint(DistanceConstraint* c, Particle* particles) {
    Particle* a = &particles[c->particle_a];
    Particle* b = &particles[c->particle_b];

    Vec2 delta = vec2_sub(b->position, a->position);
    float current_length = vec2_length(delta);
    float error = current_length - c->rest_length;

    if (current_length < 0.0001f) return;  // Avoid division by zero

    // Correction direction
    Vec2 correction = vec2_scale(delta, error / current_length);

    // Apply correction based on inverse masses
    float total_inv_mass = a->inv_mass + b->inv_mass;
    if (total_inv_mass < 0.0001f) return;  // Both static

    a->position = vec2_add(a->position, vec2_scale(correction, a->inv_mass / total_inv_mass));
    b->position = vec2_sub(b->position, vec2_scale(correction, b->inv_mass / total_inv_mass));
}

Verification:

Connected particles maintain distance
Multiple iterations improve constraint satisfaction
Chain of particles behaves like rope

Phase 5: Building the Constraint Solver (Days 16-19)

Goal: Unified constraint solving system

Tasks:

Abstract constraint interface
Build constraint Jacobian
Implement Gauss-Seidel solver
Handle mixed constraint types

Implementation:

// Generic constraint solver (simplified)
void solve_constraints(PhysicsWorld* world, int iterations) {
    for (int iter = 0; iter < iterations; iter++) {
        // Distance constraints
        for (int i = 0; i < world->num_constraints; i++) {
            solve_distance_constraint(&world->constraints[i], world->particles);
        }

        // Contact constraints
        for (int i = 0; i < world->num_contacts; i++) {
            solve_contact(&world->contacts[i], world->particles);
        }
    }
}

Verification:

All constraints satisfied (within tolerance)
Stable with many constraints
Convergence improves with iterations

Phase 6: Complex Scenes (Days 20-23)

Goal: Build interesting simulations

Tasks:

Create chain/rope scene
Create bridge scene
Create cloth-like grid
Add user interaction (drag, spawn)

Scenes to implement:

void create_rope_scene(PhysicsWorld* world) {
    // Chain of particles with distance constraints
    for (int i = 0; i < 10; i++) {
        add_particle(world, i * 30, 100, 5);
        if (i > 0) {
            add_distance_constraint(world, i-1, i, 30);
        }
    }
    // Pin first particle
    world->particles[0].inv_mass = 0;  // Static
}

void create_soft_body_scene(PhysicsWorld* world) {
    // Grid of particles with structural and shear constraints
    int grid_size = 5;
    for (int y = 0; y < grid_size; y++) {
        for (int x = 0; x < grid_size; x++) {
            add_particle(world, 200 + x * 30, 100 + y * 30, 5);
        }
    }
    // Add constraints (structural, shear, bend)
    // ...
}

Verification:

Rope swings realistically
Soft body maintains shape but deforms
No explosions or jitter

Phase 7: Polish and Optimization (Days 24-28)

Goal: Make it robust and fast

Tasks:

Add friction to collisions
Implement warm starting
Add broad-phase collision detection
Profile and optimize

Friction:

void apply_friction(Particle* a, Particle* b, Contact* c, float friction) {
    Vec2 rel_vel = vec2_sub(a->velocity, b->velocity);
    Vec2 tangent = vec2_sub(rel_vel, vec2_scale(c->normal, vec2_dot(rel_vel, c->normal)));
    float tangent_speed = vec2_length(tangent);

    if (tangent_speed < 0.001f) return;

    tangent = vec2_scale(tangent, 1.0f / tangent_speed);
    float jt = -vec2_dot(rel_vel, tangent) / (a->inv_mass + b->inv_mass);

    // Clamp to Coulomb friction
    jt = clamp(jt, -friction * c->normal_impulse, friction * c->normal_impulse);

    Vec2 friction_impulse = vec2_scale(tangent, jt);
    a->velocity = vec2_add(a->velocity, vec2_scale(friction_impulse, a->inv_mass));
    b->velocity = vec2_sub(b->velocity, vec2_scale(friction_impulse, b->inv_mass));
}

Verification:

Objects stack without sliding
System remains stable under stress
60 FPS with 100+ objects

Testing Strategy

Unit Tests

Vector operations:
- Dot product of perpendicular vectors = 0
- Normalized vector has length 1
- Reflection formula is correct
Collision detection:
- Overlapping circles detected
- Non-overlapping circles not detected
- Normal points from A to B
Integration:
- Constant velocity → position increases linearly
- Constant acceleration → velocity increases linearly
Constraints:
- Distance constraint maintains length
- Constraint solver converges

Physical Validation

Conservation of energy (elastic collision, no gravity)
Conservation of momentum (any collision)
Projectile motion (matches v₀t + ½gt²)
Pendulum period (matches 2π√(L/g))

Stress Tests

Many particles: 100+ without slowdown
Long chains: 50+ links stable
Stacking: 10+ objects stable
High speeds: No tunneling through walls

Common Pitfalls and Debugging

Explosion/Instability

Symptom: Particles fly off to infinity

Causes:

Time step too large
Forces computed incorrectly
Penetration not resolved

Fixes:

Reduce Δt (try 1/120 instead of 1/60)
Add position correction
Clamp velocities to maximum

Sinking Objects

Symptom: Objects slowly fall through floors

Causes:

Collision detected but not resolved
Bias/slop too large
Not enough solver iterations

Fixes:

Add position correction after velocity correction
Reduce bias (e.g., 0.01 instead of 0.1)
Increase iterations (10-20)

Jitter/Vibration

Symptom: Objects vibrate in place

Causes:

Over-correction
Conflicting constraints
No damping

Fixes:

Add damping: v *= 0.99
Reduce correction factor
Warm start solver

Constraint Stretching

Symptom: Springs stretch too much

Causes:

Not enough solver iterations
Stiffness too low
Integration error

Fixes:

More iterations (20+)
Use position-based dynamics
Smaller time step

Tunneling

Symptom: Fast objects pass through walls

Causes:

Discrete collision detection
Large time step
Small obstacles

Fixes:

Continuous collision detection (CCD)
Limit max velocity
Sub-step physics

Extensions and Challenges

Beginner Extensions

Add mouse interaction (drag particles)
Create preset scenes (bridge, rope, dominos)
Add visual effects (trails, colors)

Intermediate Extensions

Implement friction
Add rigid body rotation
Create hinge joints
Implement spatial hashing for broad phase

Advanced Extensions

Continuous collision detection
3D extension
Deformable bodies
Fluid simulation (SPH)

Real-World Connections

Game Engines

Every physics engine (Box2D, Bullet, PhysX) uses these concepts:

Constraint-based dynamics
Iterative solvers
Warm starting
Continuous collision detection

Robotics

Robot dynamics and control:

Forward/inverse kinematics
Constraint satisfaction
Force/torque computation

Computer Graphics

Physically-based animation:

Cloth simulation
Hair dynamics
Destruction effects

Engineering Simulation

Structural analysis:

Finite element methods
Stress/strain computation
Dynamic loading

Self-Assessment Checklist

Conceptual Understanding

I can explain why we need to solve Ax = b for physics
I understand the role of the Jacobian matrix
I can describe why iterative solvers work for constraints
I know what matrix conditioning means for stability
I understand the tradeoffs of different integration methods

Practical Skills

I implemented collision detection and response
I built a working constraint solver
My simulation runs stable at 60 FPS
I can create various constraint types
I can debug common physics bugs

Teaching Test

Can you explain to someone else:

How does impulse-based collision response work?
Why do constraints create systems of equations?
What makes a simulation stable vs. unstable?
How does Gauss-Seidel iteration work?

Resources

Primary

“Game Physics Engine Development” by Ian Millington
“Physics for Game Developers” by David Bourg
Box2D source code - Reference implementation

Reference

“Real-Time Collision Detection” by Christer Ericson
“Numerical Recipes in C” by Press et al. - Equation solvers
Erin Catto’s GDC presentations - Box2D author’s talks

Online

Box2D documentation - Algorithm explanations
Gaffer On Games (gafferongames.com) - Physics articles
Two-Bit Coding YouTube - Physics engine tutorials

When you complete this project, you’ll understand that physics simulation is applied linear algebra—solving equations, projecting vectors, and navigating numerical challenges. Every game, every simulation, every robot uses these techniques.