Animation & Physical Simulation

A keyframe is a snapshot of state — position, orientation, scale, or any other parameter — at a specific moment in time. The animator specifies a sparse set of these snapshots; the system must interpolate smoothly between them. This is the fundamental abstraction underlying all animation systems.

Why not linear interpolation? Linear interpolation (LERP) between keyframes produces motion that feels mechanical: it arrives at each pose at full speed and leaves immediately, creating visible velocity discontinuities. Real motion accelerates and decelerates. We need curves that are not just C0-continuous (no jumps) but C1-continuous (no velocity discontinuities) and ideally C2-continuous (no acceleration discontinuities, giving perceptually smooth motion).

Hermite cubics. A cubic Hermite spline is defined by two endpoints p0, p1 and two tangent vectors m0, m1. The curve is:

p(t) = (2t³-3t²+1)p0 + (t³-2t²+t)m0 + (-2t³+3t²)p1 + (t³-t²)m1

where t ∈ [0,1]. The four basis polynomials (H00, H10, H01, H11) are called Hermite basis functions. The key insight: the tangents give us direct control over velocity at each endpoint, making C1 continuity trivial to enforce by sharing tangents across adjacent segments.

Catmull-Rom splines automate tangent computation: the tangent at point pi is set to (pi+1 − pi−1) / 2. This gives a spline that passes through every control point with automatic C1 continuity — exactly what animators want, with no manual tangent editing.

Bezier curves (the de Casteljau construction) are the backbone of most vector graphics and many animation systems. A cubic Bezier is defined by four control points P0..P3. De Casteljau's algorithm evaluates it by repeated linear interpolation:

def de_casteljau(points, t):
    pts = list(points)
    while len(pts) > 1:
        pts = [lerp(pts[i], pts[i+1], t) for i in range(len(pts)-1)]
    return pts[0]

This is numerically stable and reveals the convex hull property: the curve lies entirely within the convex hull of its control points, enabling efficient bounding-box tests.

B-splines and NURBS give local control: moving one control point only affects nearby portions of the curve, unlike global Bezier. They also allow curves of arbitrary degree with guaranteed continuity. NURBS (Non-Uniform Rational B-Splines) extend this to rational (weighted) control points, allowing exact representation of conic sections — critical for CAD.

Practical considerations. In game engines, curves are stored as spline channels (one per DOF per bone) and evaluated thousands of times per frame. Pre-baking to fixed-timestep tables trades memory for speed. The choice of interpolation scheme is an artistic and technical tradeoff: Catmull-Rom for quick authoring, Bezier for fine control, B-splines when local editing is paramount.

import numpy as np

def catmull_rom(p0, p1, p2, p3, t):
    """Evaluate Catmull-Rom spline at parameter t in [0,1] for segment p1->p2."""
    t2 = t * t
    t3 = t2 * t
    return 0.5 * (
        (2*p1) +
        (-p0 + p2) * t +
        (2*p0 - 5*p1 + 4*p2 - p3) * t2 +
        (-p0 + 3*p1 - 3*p2 + p3) * t3
    )

# Example: interpolate 2D positions
points = np.array([[0,0],[1,2],[3,1],[4,3]], dtype=float)
t_vals = np.linspace(0, 1, 50)
for t in t_vals:
    p = catmull_rom(points[0], points[1], points[2], points[3], t)
    # p is smoothly interpolated between points[1] and points[2]

A character mesh has thousands of vertices. Replicating those vertices for every frame of animation is prohibitively expensive. Skeletal animation separates the problem: define a skeleton (a hierarchy of rigid bones connected by joints), animate the bones, then skin the mesh by deforming vertices based on nearby bone transforms.

The skeleton as a scene graph. Each bone has a local transform (translation + rotation, usually stored as a quaternion) relative to its parent. The global transform of a bone is the product of all local transforms from root to that bone. This is a standard tree traversal:

for each bone b in topological order:
    global[b] = global[parent[b]] * local[b]

Bind pose and inverse bind matrices. At rigging time, the mesh is in a bind pose aligned with the skeleton. Each bone stores its inverse bind matrix M_b^{-1}: the transform that moves a vertex from world space into the bone's local space at bind time. At runtime, the skinning matrix for bone b is:

S_b = global[b] * M_b_inverse

This maps: world (bind) → bone local → world (current pose). Elegant.

Linear Blend Skinning (LBS). Each vertex v has a set of (bone, weight) pairs summing to 1. The skinned position is:

v' = sum_b w_b * S_b * v

LBS is fast (just a weighted average of matrix-transformed positions) and GPU-friendly (compute shader, 4 bones per vertex in hardware). Its weakness: the candy-wrapper artifact — when two bones at a joint rotate in opposite directions, LBS averages the matrices, collapsing volume in the middle.

Dual Quaternion Skinning (DQS) fixes this by blending dual quaternions (which encode both rotation and translation) instead of raw matrices. Dual quaternions can't represent scale, but for rigid bones they preserve volume perfectly, eliminating the candy-wrapper artifact at modest extra cost.

Forward Kinematics (FK) vs. Inverse Kinematics (IK). FK is the natural direction: given joint angles, compute the end-effector position. IK poses the inverse problem: given a desired end-effector position (e.g., the hand should touch a doorknob), find the joint angles. IK is solved iteratively. The classic technique is Jacobian-based IK: the Jacobian J maps joint-angle velocities to end-effector velocities. We solve ΔΘ = J^† Δx (using the pseudoinverse or transpose for speed) and iterate. FABRIK (Forward and Backward Reaching IK) is a simpler, faster alternative popular in games: alternate between pulling the chain forward (end-effector to target) and backward (base to its original position), iterating until convergence.

Blend trees and animation state machines. At runtime, characters blend multiple animations: a 'walk' clip, a 'run' clip, and an 'idle' clip might be blended based on a speed parameter. Blending is done at the quaternion level (SLERP between poses) rather than on final vertex positions. A state machine governs transitions between blend trees, triggered by gameplay events.

# Pseudocode: FABRIK IK solver
# joints: list of 3D positions [j0, j1, ..., jN], j0 = root, jN = end-effector
# bone_lengths: lengths[i] = dist(joints[i], joints[i+1])
# target: desired end-effector position

def fabrik(joints, bone_lengths, target, tolerance=1e-3, max_iter=10):
    root = joints[0].copy()
    for iteration in range(max_iter):
        # Forward pass: pull jN to target
        joints[-1] = target
        for i in range(len(joints)-2, -1, -1):
            direction = normalize(joints[i] - joints[i+1])
            joints[i] = joints[i+1] + direction * bone_lengths[i]
        # Backward pass: pull j0 back to root
        joints[0] = root
        for i in range(len(joints)-1):
            direction = normalize(joints[i+1] - joints[i])
            joints[i+1] = joints[i] + direction * bone_lengths[i]
        # Check convergence
        if dist(joints[-1], target) < tolerance:
            break
    return joints

Rotation is the hardest degree of freedom in animation. The naive approach — storing rotations as three Euler angles (pitch, yaw, roll) — suffers from gimbal lock: a configuration where two rotation axes align, collapsing a degree of freedom. Interpolating Euler angles also produces non-intuitive, curved paths in rotation space. Rotation matrices (3x3 orthogonal matrices with det=1) avoid gimbal lock but are redundant (9 numbers for 3 DOF) and interpolate poorly (linear interpolation of matrices does not stay on SO(3)).

Quaternions are four-component numbers q = w + xi + yj + zk, where i,j,k are imaginary units satisfying i²=j²=k²=ijk=-1. A unit quaternion (|q|=1) represents a rotation of angle θ around axis n as:

q = cos(θ/2) + sin(θ/2)(n_x i + n_y j + n_z k)

Rotating a vector v by q is computed as q * v * q^{-1} (quaternion product, treating v as a pure quaternion). Quaternion composition is just multiplication: q1 then q2 is q2*q1. This is more efficient than matrix multiplication and numerically stabler (re-normalizing a 4-vector is cheap).

The double-cover problem. Both q and -q represent the same rotation. When interpolating, you must check sign(q1·q2) and negate q2 if negative to avoid taking the 'long way around' (interpolating through 270° instead of 90°). This is called short-path SLERP.

SLERP (Spherical Linear Interpolation) interpolates between two quaternions along the great-circle arc on the unit 4-sphere S³:

SLERP(q0, q1, t) = q0 * (q0^{-1} * q1)^t

Equivalently:

SLERP(q0, q1, t) = sin((1-t)Ω)/sin(Ω) * q0 + sin(tΩ)/sin(Ω) * q1

where Ω = arccos(q0·q1). SLERP maintains constant angular velocity, which is exactly what animators want for smooth rotation.

SQUAD extends SLERP to cubic interpolation across a sequence of rotations, analogous to Catmull-Rom splines but on S³, by computing intermediate 'helper' quaternions using the De Casteljau construction on the sphere.

Converting between representations is a practical necessity: artists think in Euler angles (intuitive), rigs store quaternions (efficient), shaders use matrices (GPU-friendly). Every production engine implements conversions between all three, and the order of Euler angle application (XYZ, ZYX, etc.) must be tracked meticulously to avoid silent bugs.

import numpy as np

def slerp(q0, q1, t):
    """Spherical linear interpolation between unit quaternions."""
    q0 = np.array(q0, dtype=float)
    q1 = np.array(q1, dtype=float)
    dot = np.dot(q0, q1)
    # Short-path: negate q1 if dot < 0
    if dot < 0.0:
        q1 = -q1
        dot = -dot
    dot = min(dot, 1.0)  # clamp for numerical safety
    if dot > 0.9995:  # nearly identical: use LERP + normalize
        result = q0 + t * (q1 - q0)
        return result / np.linalg.norm(result)
    omega = np.arccos(dot)
    return (np.sin((1-t)*omega) * q0 + np.sin(t*omega) * q1) / np.sin(omega)

# Example: rotate from identity to 90 deg around Z
q_id  = np.array([1, 0, 0, 0], dtype=float)  # w,x,y,z
q_90z = np.array([np.cos(np.pi/4), 0, 0, np.sin(np.pi/4)], dtype=float)
q_45z = slerp(q_id, q_90z, 0.5)  # halfway => 45 deg around Z

A particle system is a collection of small, independent elements each obeying simple rules, whose aggregate produces complex visual phenomena: fire, smoke, sparks, rain, explosions, flocking birds. The power of the abstraction comes from emergence: rich behavior from simple per-particle rules.

The basic loop. Each particle has state: position x, velocity v, age, color, size, and optional user-defined fields. Each frame:

1. Emit new particles from an emitter (a point, sphere, mesh surface, etc.) with initial conditions drawn from probability distributions.
2. Update each particle: integrate equations of motion, update age, apply field forces.
3. Kill particles whose age exceeds their lifetime.
4. Render surviving particles (sprites, meshes, metaballs, volumetric density).

Forces and fields. Particles respond to gravity, wind, drag, vortex fields, attractor/repeller points, and user-defined force volumes. A common model:

a = (1/m) * (gravity + drag + F_user)
v += a * dt
x += v * dt

Soft particle rendering. Naive particle rendering shows harsh intersections where sprites clip into geometry. Soft particles fade opacity based on the difference between the particle's depth and the scene depth buffer, producing smooth integration with the environment — a technique now standard in every production engine.

GPU particles. Modern systems simulate millions of particles on the GPU using compute shaders. Each particle is an element of a structured buffer. The update kernel runs in parallel, and dead-particle management uses a concurrent freelist (GPU atomics). Rendering uses indirect draw calls — the CPU never reads back particle data.

Sub-emitters. Particles can spawn child emitters on death (sparks that emit smoke) or on collision (water droplets that splash), enabling compositional hierarchies of effects.

Stability. The naive Euler integrator accumulates error and can explode for stiff forces (e.g., strong springs). For particles under moderate forces, Euler is usually fine. For cloth and springs, switch to Verlet or implicit integration (see the Physics Simulation concept below).

Flocking (Boids). A deceptively simple particle system with three rules per agent:

1. Separation: steer away from nearby agents.
2. Alignment: steer toward average heading of neighbors.
3. Cohesion: steer toward average position of neighbors. This produces lifelike schooling, flocking, and swarming from purely local interactions — a canonical example of emergent behavior and relevant to multi-agent simulation.

# Minimal particle system in Python (NumPy-vectorized)
import numpy as np

N = 10_000
pos = np.random.randn(N, 3).astype(np.float32)  # initial positions
vel = np.random.randn(N, 3).astype(np.float32) * 0.1
age = np.zeros(N, dtype=np.float32)
lifetime = np.random.uniform(1.0, 5.0, N).astype(np.float32)

gravity = np.array([0.0, -9.8, 0.0], dtype=np.float32)
drag = 0.98  # velocity damping per frame
dt = 1.0 / 60.0

def update(pos, vel, age, lifetime, dt):
    vel += gravity * dt
    vel *= drag
    pos += vel * dt
    age += dt
    alive = age < lifetime
    # Respawn dead particles at origin with new velocities
    dead = ~alive
    pos[dead] = 0.0
    vel[dead] = np.random.randn(dead.sum(), 3).astype(np.float32) * 2.0
    age[dead] = 0.0
    return pos, vel, age

# Simulation loop
for frame in range(300):
    pos, vel, age = update(pos, vel, age, lifetime, dt)

The core problem of physics simulation is solving ordinary differential equations (ODEs) of the form dx/dt = f(x, t) forward in time with finite discrete steps. The choice of integrator determines accuracy, stability, and energy behavior — the tradeoffs are subtle and consequential.

Explicit (Forward) Euler is the simplest:

v(t+dt) = v(t) + a(t) * dt
x(t+dt) = x(t) + v(t) * dt

Error is O(dt²) per step, O(dt) globally. Crucially, Euler is not energy-conserving: for an oscillator, energy grows unboundedly over time, causing the simulation to explode. Never use pure Euler for anything involving springs or sustained oscillation.

Symplectic Euler (Semi-implicit Euler) uses the new velocity to update position:

v(t+dt) = v(t) + a(t) * dt
x(t+dt) = x(t) + v(t+dt) * dt  # <-- uses updated v

This tiny change makes it a symplectic integrator: it exactly conserves a slightly modified Hamiltonian, so energy remains bounded. It is first-order accurate but energy-stable, which is why almost every game engine uses it for rigid bodies and particles.

Verlet integration maintains position at two consecutive timesteps and derives velocity implicitly:

x(t+dt) = 2*x(t) - x(t-dt) + a(t) * dt²
v(t) = (x(t+dt) - x(t-dt)) / (2*dt)  # central difference

Verlet is second-order accurate and time-reversible. Velocity Verlet is equivalent but stores velocity explicitly:

x(t+dt) = x(t) + v(t)*dt + 0.5*a(t)*dt²
a(t+dt) = f(x(t+dt)) / m
v(t+dt) = v(t) + 0.5*(a(t) + a(t+dt))*dt

Verlet is the standard for cloth simulation, molecular dynamics, and position-based dynamics.

Runge-Kutta 4 (RK4) achieves fourth-order accuracy by evaluating the derivative at four sub-steps:

k1 = f(t, x)
k2 = f(t+dt/2, x + dt/2 * k1)
k3 = f(t+dt/2, x + dt/2 * k2)
k4 = f(t+dt, x + dt * k3)
x(t+dt) = x + dt/6 * (k1 + 2*k2 + 2*k3 + k4)

RK4 is excellent for offline simulation and cutscene animation where accuracy matters more than speed. For real-time games, symplectic Euler or Verlet usually suffice.

Implicit integration evaluates forces at the next timestep, requiring a solve at each step. For stiff systems (cloth, soft bodies) implicit methods allow large timesteps without explosion at the cost of solving a linear system. Position-Based Dynamics (PBD) takes a different approach: instead of integrating forces, it directly projects positions to satisfy constraints, which is unconditionally stable and highly amenable to parallelism — the approach used in Nvidia PhysX cloth and Unreal's Chaos.

Stability and stiffness. A spring with stiffness k and mass m has natural frequency ω = √(k/m). The explicit Euler stability condition is dt < 2/ω. For a stiff spring (large k), this forces tiny timesteps. The alternative is implicit integration or constraint-based methods.

# RK4 integrator for a harmonic oscillator: x'' = -omega^2 * x
import numpy as np

omega = 2.0  # rad/s
dt = 0.01
t_end = 10.0

def deriv(state, t):
    # state = [x, v]
    x, v = state
    return np.array([v, -omega**2 * x])

def rk4_step(state, t, dt):
    k1 = deriv(state, t)
    k2 = deriv(state + 0.5*dt*k1, t + 0.5*dt)
    k3 = deriv(state + 0.5*dt*k2, t + 0.5*dt)
    k4 = deriv(state + dt*k3, t + dt)
    return state + (dt/6.0) * (k1 + 2*k2 + 2*k3 + k4)

state = np.array([1.0, 0.0])  # x=1, v=0
t = 0.0
while t < t_end:
    state = rk4_step(state, t, dt)
    t += dt
# state[0] ≈ cos(omega * t_end) -- check energy conservation!

A rigid body is an idealized solid that does not deform. Its motion is governed by two coupled ODEs: one for the center of mass (linear momentum) and one for rotation (angular momentum).

State representation. The full state of a rigid body is: position x (center of mass), linear velocity v, orientation q (unit quaternion), angular velocity ω (in world space), linear momentum p = mv, and angular momentum L = Iω. Storing both momentum and velocity is redundant but convenient.

Equations of motion:

dx/dt = v
dv/dt = F/m
dq/dt = 0.5 * [0, ω] * q   (quaternion derivative)
dω/dt = I^{-1} * (τ - ω × (I*ω))

where F is net force, τ is net torque, and I is the inertia tensor — a 3x3 matrix encoding how mass is distributed. For a solid sphere of mass m and radius r, I = (2/5)mr²·Identity. The ω × (Iω) term is the gyroscopic term, responsible for the precession of a spinning top.

The inertia tensor in world space. I is typically computed in body (local) frame where it is diagonal. In world frame: I_world = R * I_body * R^T, where R is the rotation matrix derived from q. At each step, we recompute I_world^{-1} from I_body^{-1} and the current rotation.

Impulse-based response. When two rigid bodies collide, we apply instantaneous impulses rather than forces (which would require infinitesimal timesteps). The impulse magnitude j is derived from the constraint that post-collision relative velocity at the contact point is -e * pre-collision relative velocity (e = coefficient of restitution). The formula:

j = -(1+e) * v_rel · n
    / (1/m1 + 1/m2 + n·(I1^{-1}(r1×n))×r1 + n·(I2^{-1}(r2×n))×r2)

where n is the contact normal, r1, r2 are vectors from the bodies' centers of mass to the contact point, and v_rel is the relative velocity at the contact point.

Constraints and joints. Articulated structures (ragdolls, vehicle suspensions) require joints that constrain relative motion between bodies. Constraint solving is formulated as a linear complementarity problem (LCP) or handled by sequential impulses (Projected Gauss-Seidel), which iteratively applies impulses to reduce constraint error. Unreal Engine's Chaos and PhysX both use sequential impulse solvers.

Sleeping. Bodies below a velocity/kinetic-energy threshold are put to sleep: they are excluded from integration and collision detection until a collision wakes them. This is a critical optimization — a scene with thousands of sleeping objects remains efficient.

# Rigid body linear + angular update (symplectic Euler, no collision)
import numpy as np
from scipy.spatial.transform import Rotation

class RigidBody:
    def __init__(self, mass, inertia_body):
        self.mass = mass
        self.inv_mass = 1.0 / mass
        self.inertia_body_inv = np.linalg.inv(inertia_body)
        self.pos = np.zeros(3)
        self.vel = np.zeros(3)
        self.quat = np.array([1,0,0,0], dtype=float)  # w,x,y,z
        self.omega = np.zeros(3)  # world-space angular velocity

    def apply_force_torque(self, F, tau, dt):
        # Linear: symplectic Euler
        self.vel += (F * self.inv_mass) * dt
        self.pos += self.vel * dt
        # Angular: compute world-space inertia inverse
        R = Rotation.from_quat([*self.quat[1:], self.quat[0]]).as_matrix()
        I_inv_world = R @ self.inertia_body_inv @ R.T
        # Integrate angular momentum
        self.omega += I_inv_world @ tau * dt
        # Update quaternion: dq/dt = 0.5 * [0,omega] * q
        omega_quat = np.array([0, *self.omega])
        dq = 0.5 * quat_mul(omega_quat, self.quat)
        self.quat += dq * dt
        self.quat /= np.linalg.norm(self.quat)  # renormalize

Collision detection asks: are these objects overlapping, and if so, where? It is typically split into two phases: broad phase (quickly eliminate pairs that definitely do not collide) and narrow phase (precisely compute contact manifold for candidate pairs).

Broad phase: spatial data structures. The naive O(n²) pairwise check is untenable for n > a few hundred. Common solutions:

• Bounding Volume Hierarchies (BVH): Build a tree where each node's bounding box contains all descendant geometry. Test two BVHs by recursively descending whenever bounding boxes overlap. Construction (top-down with median split, or SAH — Surface Area Heuristic) is O(n log n); query is O(log² n) for non-overlapping objects.
• Sweep and Prune (Sort and Sweep): Sort AABBs along each axis and find overlapping intervals. Update via insertion sort each frame (near-O(n) for temporally coherent scenes).
• Spatial Hashing: Partition space into a fixed grid; hash each cell. O(1) average lookup, excellent for uniform-density particle systems.

Narrow phase: GJK algorithm. The Gilbert-Johnson-Keerthi (GJK) algorithm determines whether two convex shapes overlap using only a support function (the point on the shape furthest in a given direction). GJK iteratively builds a simplex (point, line, triangle, tetrahedron in 3D) that approximates the Minkowski difference A⊖B. If the origin lies inside A⊖B, the shapes overlap. GJK is extremely fast in practice and works for any convex shape (spheres, capsules, convex hulls) via a uniform interface.

EPA (Expanding Polytope Algorithm) extends GJK to compute the penetration depth and contact normal when objects do overlap, by expanding the simplex toward the origin boundary.

SAT (Separating Axis Theorem). For convex polyhedra, if no separating axis exists among the face normals and edge cross-products of both shapes, they intersect. Efficient for low-polygon convex hulls; outperforms GJK for simple cases.

Continuous Collision Detection (CCD). Discrete collision detection can miss fast-moving objects that tunnel through thin geometry between frames. CCD sweeps the object's path over the timestep and finds the first time of impact (TOI) analytically or via bisection, preventing tunneling. Essential for projectiles and thin obstacles.

Deformable and triangle-mesh collision. Non-convex meshes require decomposition into convex parts (V-HACD) or specialized algorithms (BVH traversal for triangle-mesh pairs). This is computationally expensive and is a major reason game collision meshes are simplified convex proxies of the visual mesh.

# Broad phase: Axis-Aligned Bounding Box overlap test and simple spatial hash
import numpy as np
from collections import defaultdict

def aabb_overlap(min_a, max_a, min_b, max_b):
    """True if AABB [min_a, max_a] overlaps [min_b, max_b] in 3D."""
    return (np.all(min_a <= max_b) and np.all(min_b <= max_a))

class SpatialHash:
    def __init__(self, cell_size):
        self.cell_size = cell_size
        self.table = defaultdict(list)

    def _hash(self, x, y, z):
        return (int(x)//self.cell_size, int(y)//self.cell_size, int(z)//self.cell_size)

    def insert(self, obj_id, pos):
        key = self._hash(*pos)
        self.table[key].append(obj_id)

    def query_neighbors(self, pos):
        cx, cy, cz = self._hash(*pos)
        neighbors = []
        for dx in (-1, 0, 1):
            for dy in (-1, 0, 1):
                for dz in (-1, 0, 1):
                    neighbors += self.table[(cx+dx, cy+dy, cz+dz)]
        return neighbors

Procedural generation replaces hand-authored content with algorithms that produce infinite, varied, plausible results from compact parameters. This is not just a convenience — it reflects deep mathematical structure in natural phenomena.

Perlin noise (and its improved variant, Simplex noise) is the workhorse of procedural texture and terrain. Perlin noise is a gradient noise: at each integer lattice point, assign a random gradient vector; interpolate smoothly between the dot products of gradient and offset vectors. The result is a smooth, band-limited noise with a characteristic frequency. Key properties:

• Coherence: nearby points have similar values (unlike white noise).
• Zero DC: the average is zero.
• Isotropy (approximately): looks similar in all directions.

Fractional Brownian Motion (fBm) sums octaves of noise at increasing frequencies and decreasing amplitudes:

fbm(x) = sum_i (1/2^i) * noise(2^i * x)

This produces the self-similar, fractal-like roughness of mountains, clouds, and coastlines — matching Mandelbrot's fractal geometry of nature.

Domain warping feeds the output of one fBm back as the input coordinates of another, producing the swirling, folded patterns of marble and turbulent smoke. A technique pioneered by Inigo Quilez and widely used in shader art.

L-Systems (Lindenmayer Systems) model plant growth using formal grammars. Start with an axiom (a string of symbols), apply rewriting rules repeatedly, then interpret the resulting string as turtle-graphics drawing commands:

Axiom: F
Rule: F -> F[+F]F[-F]F

After a few iterations, this draws a branching plant. L-systems capture the recursive, self-similar structure of botanical morphology. They underlie most procedural vegetation in game engines.

Fractals and the Mandelbrot set. The Mandelbrot set is defined by the iteration z_{n+1} = z_n² + c, marking which complex c values keep z bounded. Its fractal boundary has Hausdorff dimension ~2 and infinite complexity — an icon of emergent complexity from simple rules. Julia sets (fixed c, varying initial z) produce the organic, decorative forms used in procedural art.

Voronoi diagrams partition space into cells, each containing all points closest to a given seed. Voronoi noise is used for cracked earth, cell patterns, and stylized rendering. The dual structure (Delaunay triangulation) is used for mesh generation and spatial subdivision.

Wave Function Collapse (WFC) is a modern constraint-based procedural algorithm: given examples of locally valid tile patterns, it generates globally consistent tilings by propagating constraints from a partially-collapsed state, analogous to quantum measurement. Used in Caves of Qud and many indie games for dungeon/world generation.

Practical pipeline. In a game like Minecraft or No Man's Sky, the world is generated by composing: terrain height from fBm, biome selection from low-frequency noise, feature placement via Poisson disk sampling, structure generation via template instantiation, and surface decoration via scatter placement on the terrain mesh — all without storing the world in advance.

# Simplex-style value noise (2D) and fBm terrain height
import numpy as np

def smoothstep(t):
    return t * t * (3 - 2 * t)

def value_noise_2d(x, y, seed=0):
    """Simple value noise via hashed lattice + bilinear interpolation."""
    rng = np.random.default_rng(seed)
    def rand_val(ix, iy):
        rng2 = np.random.default_rng((int(ix) * 1664525 + int(iy) * 1013904223) & 0xFFFFFFFF)
        return rng2.random()
    ix, iy = int(np.floor(x)), int(np.floor(y))
    fx, fy = x - ix, y - iy
    ux, uy = smoothstep(fx), smoothstep(fy)
    v00 = rand_val(ix,   iy  )
    v10 = rand_val(ix+1, iy  )
    v01 = rand_val(ix,   iy+1)
    v11 = rand_val(ix+1, iy+1)
    return (v00*(1-ux)*(1-uy) + v10*ux*(1-uy) +
            v01*(1-ux)*uy     + v11*ux*uy)

def fbm(x, y, octaves=6, lacunarity=2.0, gain=0.5):
    """Fractional Brownian Motion: sum of noise octaves."""
    value, amplitude, frequency = 0.0, 0.5, 1.0
    for _ in range(octaves):
        value += amplitude * value_noise_2d(x * frequency, y * frequency)
        amplitude *= gain
        frequency *= lacunarity
    return value

# Generate 256x256 terrain heightmap
size = 256
heightmap = np.array([[fbm(x/64, y/64) for x in range(size)] for y in range(size)])

Cloth simulation extends particle systems with structural constraints between particles, enforcing that connected particles maintain (approximately) fixed distances. This moves us from independent particle dynamics to a coupled constraint system.

Mass-spring model. Represent the cloth as a grid of mass particles connected by springs: structural springs (horizontal/vertical neighbors), shear springs (diagonal neighbors), and bend springs (two-hop neighbors that resist out-of-plane bending). Each spring contributes a Hooke's law force: F = -k(|x1-x2| - rest_length) * normalize(x1-x2).

The problem: stiff springs (high k) demand tiny timesteps with explicit Euler. Solutions:

1. Implicit integration (Baraff-Witkin 1998): form a linear system at each step, solve it. Unconditionally stable but expensive to factor.
2. Position-Based Dynamics (PBD): After integrating positions with no forces, directly project each constraint to its satisfied state, iterating for a fixed number of Gauss-Seidel iterations. Stable, parallelizable, and tunable via iteration count. Used in Unreal's Chaos cloth.
3. XPBD (Extended PBD): Adds compliance (inverse stiffness) to PBD, making constraint stiffness timestep-independent and physically meaningful.

Self-collision. A falling cloth can fold onto itself. Detecting self-collision requires BVH-BVH queries over the cloth mesh at each frame — the dominant cost in high-quality cloth simulation.

Finite Element Method (FEM) for soft bodies models material deformation using continuum mechanics: the strain tensor measures local deformation, the stress tensor gives internal forces, and the constitutive model (e.g., Neo-Hookean, Mooney-Rivlin) relates the two. FEM is the gold standard for medical simulation and film VFX but is too expensive for real-time games without significant approximation.

Practical notes. Most game cloth (capes, hair, foliage) uses PBD/XPBD with 2-10 constraint iterations, simulated at 60-120Hz, with wind forces sampled from a turbulence field. Self-collision is often skipped or approximated with sphere-mesh collision against the character body.

# Position-Based Dynamics: cloth simulation (core constraint projection)
import numpy as np

N = 10  # grid size
rest = 1.0 / (N - 1)  # rest length between adjacent particles

# Initialize grid positions (flattened N*N array)
pos = np.array([[x/(N-1), 0.0, y/(N-1)] for y in range(N) for x in range(N)], dtype=float)
old_pos = pos.copy()
vel = np.zeros_like(pos)

# Build constraint list: (i, j, rest_length)
constraints = []
for y in range(N):
    for x in range(N):
        idx = y*N + x
        if x+1 < N: constraints.append((idx, y*N+(x+1), rest))
        if y+1 < N: constraints.append((idx, (y+1)*N+x, rest))

def pbd_step(pos, old_pos, dt, iterations=5):
    gravity = np.array([0, -9.8, 0])
    # Verlet integrate
    new_pos = 2*pos - old_pos + gravity * dt**2
    old_pos = pos.copy()
    # Satisfy constraints
    for _ in range(iterations):
        for (i, j, d0) in constraints:
            delta = new_pos[j] - new_pos[i]
            dist = np.linalg.norm(delta)
            if dist < 1e-8: continue
            correction = 0.5 * (dist - d0) * delta / dist
            new_pos[i] += correction
            new_pos[j] -= correction
    # Pin top row
    for x in range(N): new_pos[x] = pos[x]
    return new_pos, old_pos