Inside Game Physics: Architecture, Numerical Solvers, and Pipelines

When a physics simulation in a competitive shooter drops from 144 FPS to 28 FPS during an explosive particle burst, a naive variable delta time loop instantly destabilizes. Rigid bodies gain phantom kinetic energy, ragdolls invert violently through level geometry, and constraint solvers diverge into numerical infinity. Game physics is fundamentally an exercise in real-time numerical stability under rigid CPU budgets, trading theoretical analytical perfection for bounded, predictable computational convergence.

A production-ready physics system does not simply compute textbook differential equations. It orchestrates a high-throughput pipeline: state extrapolation, broad-phase spatial culling, narrow-phase geometric intersection, manifold generation, iterative constraint relaxation, and network state reconciliation. Every stage must execute within a strict slice of a single frame, typically between 2.0 and 4.0 milliseconds on modern multi-threaded architectures.

This reference architecture dissects modern rigid body dynamics from the bare mathematics of state vectors to modern engine implementations in 2026. We will walk through the mechanical differences between explicit and implicit numerical integrators, analyze the collision detection pipeline, construct a sequential impulse contact solver, and evaluate performance characteristics across Jolt, PhysX 5, and Unreal Engine 5 Chaos.

Core Foundations of Game Physics and Simulation Loops

Interactive simulations represent rigid bodies as collections of dynamic properties: position vectors, linear momentum, orientation quaternions, and angular momentum. At runtime, the fundamental objective of game physics is to compute the evolution of these state variables over discrete intervals of time. In classical analytical mechanics, continuous differential equations govern acceleration and torque. In gaming physics, finite processor cycles mandate discrete numerical approximations, making the simulation loop the most critical architectural decision in an engine.

Binding physics updates directly to variable render frame times is a catastrophic architectural flaw. When delta time fluctuates, numerical integrators experience non-deterministic energy injection. A physics step taking 33 milliseconds injects significantly higher discrete momentum than two consecutive 16 millisecond steps, resulting in simulation jitter, unstable joint constraints, and game-breaking tunneling artifacts.

+-------------------------------------------------------------+
| Variable Render Loop |
| Accumulator += FrameTime (clamped to MaxAccumulatorTime) |
+-------------------------------------------------------------+
 |
 v
+-------------------------------------------------------------+
| Fixed Timestep Physics Sub-step |
| while (Accumulator >= FixedDeltaTime) |
+-------------------------------------------------------------+
| 1. Integrate Forces (Gravity, Thrusters) |
| 2. Broad-Phase Spatial Culling (BVH / Dynamic AABB Tree) |
| 3. Narrow-Phase Geometric Query (GJK / EPA / SAT) |
| 4. Solve Constraints & Contacts (Sequential Impulses) |
| 5. Integrate Velocities to Positions |
| 6. Accumulator -= FixedDeltaTime |
+-------------------------------------------------------------+
 |
 v
+-------------------------------------------------------------+
| Alpha Transform Interpolation |
| RenderState = Lerp(PreviousState, CurrentState, |
| Accumulator / FixedDeltaTime) |
+-------------------------------------------------------------+

To achieve determinism and absolute stability, engines employ a fixed timestep with an accumulator buffer. This pattern, popularized as the “Fix Your Timestep” model, buffers elapsed wall-clock time and consumes it in uniform sub-steps (such as 60 Hz or 120 Hz). Any leftover fractional time is stored in the accumulator and used to interpolate visual transforms between the previous and current physics state, decoupling visual refresh rates from mathematical integration.

#include <chrono>
#include <algorithm>

struct TransformState {
 float posX, posY, posZ;
 float rotX, rotY, rotZ, rotW;
};

class PhysicsSimulationLoop {
public:
 static constexpr double FixedDeltaTime = 1.0 / 60.0; // 60 Hz tick rate
 static constexpr double MaxFrameTime = 0.25; // Spiral-of-death clamp

 void Update(double realDeltaTime) {
 // Guard against long pauses (e.g. breakpoint or OS hitch)
 double frameTime = std:min(realDeltaTime, MaxFrameTime);
 accumulator += frameTime;

 while (accumulator >= FixedDeltaTime) {
 previousState = currentState;
 SimulateSubstep(currentState, FixedDeltaTime);
 accumulator -= FixedDeltaTime;
 }

 // Compute alpha for render interpolation: [0.0, 1.0]
 const double alpha = accumulator / FixedDeltaTime;
 renderState = InterpolateTransforms(previousState, currentState, static_cast<float>(alpha));
 }

private:
 double accumulator = 0.0;
 TransformState previousState{};
 TransformState currentState{};
 TransformState renderState{};

 void SimulateSubstep(TransformState& state, double dt) {
 // Step rigid body pipeline: Broadphase -> Narrowphase -> Solver -> Integration
 state.posY += static_cast<float>(-9.81 * dt * dt * 0.5);
 }

 TransformState InterpolateTransforms(const TransformState& prev, const TransformState& curr, float alpha) {
 TransformState out;
 out.posX = prev.posX + alpha * (curr.posX - prev.posX);
 out.posY = prev.posY + alpha * (curr.posY - prev.posY);
 out.posZ = prev.posZ + alpha * (curr.posZ - prev.posZ);
 // Quaternion slerp would be used here for rotation
 return out;
 }
};

Architecture Rule: Never mutate player or entity visual transforms directly inside the variable frame renderer. Render transforms must always be an interpolated projection between the previous completed physics state and the current unconsumed sub-step state.

Numerical Integration Schemes: From Euler to Verlet and RK4

At the center of physics in computer games lies numerical integration: the computational technique of calculating an entity position and velocity at time t + dt given its current state and acting forces. Because computer processors cannot evaluate infinitesimal differentials continuously, they rely on discrete approximations. The selection of an integration scheme directly impacts memory bandwidth, CPU overhead, and the conservation of mechanical energy.

Explicit methods use current state variables to compute future values. Implicit methods formulate an equation involving both current and future states, solving a system of equations at each step. While implicit integrators offer unconditional stability for stiff systems (such as high-strain cloth and complex hydraulics), their high compute cost makes them rare for general rigid bodies. Consequently, video game engines lean toward symplectic integrators.

Explicit Forward Euler computes the next position using the current velocity: x(t + dt) = x(t) + v(t) * dt. Because it projects forward tangentially along curves of motion, it continually introduces artificial kinetic energy into the system. An orbiting body integrated with Forward Euler will spiral outward into infinity. In contrast, Semi-Implicit Euler (also known as Euler-Cromer) evaluates velocity first using current accelerations, then immediately evaluates the new position using the updated velocity: v(t + dt) = v(t) + a(t) * dt followed by x(t + dt) = x(t) + v(t + dt) * dt.

struct RigidBodyState {
 float position[3] = {0.0f, 0.0f, 0.0f};
 float velocity[3] = {0.0f, 0.0f, 0.0f};
 float acceleration[3] = {0.0f, -9.81f, 0.0f};
 float inverseMass = 1.0f; // 0.0f indicates static object
};

// Production standard: Symplectic / Semi-Implicit Euler Integration
void IntegrateRigidBodySemiImplicit(RigidBodyState& body, float dt) {
 if (body.inverseMass <= 0.0f) return; // Skip infinite mass objects

 // 1. Advance linear velocity with current forces
 body.velocity[0] += body.acceleration[0] * dt;
 body.velocity[1] += body.acceleration[1] * dt;
 body.velocity[2] += body.acceleration[2] * dt;

 // 2. Advance position using the *newly computed* velocity
 body.position[0] += body.velocity[0] * dt;
 body.position[1] += body.velocity[1] * dt;
 body.position[2] += body.velocity[2] * dt;
}

Semi-Implicit Euler is symplectic: it preserves phase-space volume over long iterations. While it produces localized energy errors, these errors oscillate within a bounded envelope rather than diverging toward infinity, making it the practical choice for real-time physics engines.

The Collision Pipeline: Spatial Partitioning and Narrow-Phase Solvers

A naive collision detection system checks every object against every other object in the scene. For N entities, this scales at O(N^2) complexity. Testing 10,000 objects would require 50 million mathematical checks per frame, overwhelming modern CPUs. Modern video game physics circumvents this bottleneck by decomposing collision detection into a multi-tiered pipeline: spatial acceleration broad-phase culling followed by exact narrow-phase geometry tests.

Scene Primitives (Meshes, Convex Hulls, Boxes, Capsules)
 |
 v
+-------------------------------------------------------------+
| Broad-Phase Spatial Partitioning |
| Dynamic BVH / Bounding Volume Tree (AABB) |
| Prunes Non-Intersecting Distant Pairs |
+-------------------------------------------------------------+
 |
 Candidate Pairs (AABB Overlaps)
 v
+-------------------------------------------------------------+
| Mid-Phase Primitive Decomposition |
| Compound Shapes / Triangle BVH Mesh Culling |
+-------------------------------------------------------------+
 |
 Convex Primitive Pairs
 v
+-------------------------------------------------------------+
| Narrow-Phase Geometry Tests |
| 1. Separating Axis Theorem (SAT) for Polytopes |
| 2. Gilbert-Johnson-Keerthi (GJK) for Convex Distance |
| 3. Expanding Polytope Algorithm (EPA) for Penetration |
+-------------------------------------------------------------+
 |
 v
+-------------------------------------------------------------+
| Contact Manifold Generation |
| Computes Penetration Depth, Normal, and Contacts |
+-------------------------------------------------------------+

The execution of this pipeline follows four rigorous stages:

1. Broad-Phase Bounding Volume Hierarchy (BVH): Objects are wrapped in loose Axis-Aligned Bounding Boxes (AABBs). A dynamic bounding volume tree organizes these boxes into a self-balancing binary tree structure. As objects move, their leaf nodes update, and tree traversal yields a small subset of candidate pairs whose bounding volumes overlap, executing in O(N log N) or near O(N) time.
2. Mid-Phase Mesh Traversal: For complex triangle meshes, the broad-phase candidate undergoes mid-phase evaluation. The candidate AABB queries an internal BVH representing the sub-triangles of the static terrain mesh, culling unimpacted geometry before expensive geometric tests run.
3. Narrow-Phase Convex Solving (GJK Algorithm): The Gilbert-Johnson-Keerthi (GJK) algorithm determines whether two arbitrary convex shapes intersect. GJK computes the Minkowski difference of the two sets: A (-) B = { a - b | a in A, b in B }. If the geometric shapes intersect, their Minkowski difference contains the coordinate origin. Using a directed support mapping function, GJK iteratively constructs a simplex (a tetrahedron in 3D) inside the Minkowski difference to enclose the origin.
4. Penetration Depth Estimation (EPA): If GJK detects an intersection, it cannot natively compute how deep the objects have penetrated or which surface normal will separate them. Control transfers to the Expanding Polytope Algorithm (EPA). EPA takes the simplex constructed by GJK, expands its faces outward toward the boundary of the Minkowski difference using new support points, and identifies the absolute closest face to the origin. This yields the exact penetration depth and contact normal.

Optimization Insight: Real-time engines bypass full GJK/EPA sweeps whenever analytical solutions exist. Sphere-to-sphere, sphere-to-capsule, and box-to-box intersections rely on direct mathematical formulas, routing to GJK and EPA strictly when irregular convex hulls collide.

Constraint Solving, Contact Manifolds, and Tunneling Artifacts

Detecting collision points is merely the first half of the simulation challenge. The physics engine must subsequently resolve overlaps and transfer momentum without introducing visual jitter, unstable oscillations, or unnatural energy generation. Modern production systems format physical contacts, hinges, ragdoll joints, and friction as constraint systems solved via Sequential Impulses.

A mechanical constraint is formulated as a velocity kinematic equation: J * v + b = 0, where J represents the Jacobian matrix (the partial derivatives of the constraint with respect to entity degrees of freedom), v is the generalized system velocity vector, and b represents the constraint bias (often incorporating Baumgarte Stabilization to correct position penetration). Rather than solving enormous, computationally prohibitive global matrix equations via Linear Complementarity Problem (LCP) direct solvers, real-time engines use Projected Gauss-Seidel (PGS) iterations, processing constraints sequentially until impulses stabilize.

// Production implementation of a 1D Contact Impulse Solver with Baumgarte stabilization
struct ContactConstraint {
 float normal[3]; // Contact normal pointing from BodyA to BodyB
 float contactPoint[3]; // World space contact position
 float penetration; // Penetration depth (scalar)
 float impulseSum; // Accumulated impulse clamped for friction/normal
 float effectiveMass; // Inverse of (J * M^-1 * J^T)
 float bias; // Baumgarte stabilization bias term
};

void PrepareContactConstraint(ContactConstraint& c, float dt, float restitution) {
 const float baumgarteBeta = 0.2f; // Stability factor: absorbs 20% penetration per frame
 const float slop = 0.005f; // Penetration threshold tolerance (5mm)

 // Baumgarte stabilization term prevents objects from sinking into floors
 float penetrationError = std:max(0.0f, c.penetration - slop);
 c.bias = -(baumgarteBeta / dt) * penetrationError;
 c.impulseSum = 0.0f; // Reset or warm-start with persistent manifold cache
}

void ResolveContactVelocity(ContactConstraint& c, float* velA, float* velB, float invMassA, float invMassB) {
 // Relative velocity along the normal: (velB - velA). normal
 float relVel = (velB[0] - velA[0]) * c.normal[0] +
 (velB[1] - velA[1]) * c.normal[1] +
 (velB[2] - velA[2]) * c.normal[2];

 // Calculate raw Lagrange multiplier impulse delta
 float deltaImpulse = -(relVel + c.bias) * c.effectiveMass;

 // Clamp accumulated impulse (contacts cannot pull objects together, normal force >= 0)
 float oldImpulse = c.impulseSum;
 c.impulseSum = std:max(oldImpulse + deltaImpulse, 0.0f);
 float appliedImpulse = c.impulseSum - oldImpulse;

 // Apply impulse to velocities
 velA[0] -= appliedImpulse * invMassA * c.normal[0];
 velA[1] -= appliedImpulse * invMassA * c.normal[1];
 velA[2] -= appliedImpulse * invMassA * c.normal[2];

 velB[0] += appliedImpulse * invMassB * c.normal[0];
 velB[1] += appliedImpulse * invMassB * c.normal[1];
 velB[2] += appliedImpulse * invMassB * c.normal[2];
}

High velocities introduce an artifact known as tunneling: fast-moving entities (such as projectiles or high-speed vehicles) pass completely through solid objects between frames because their displacement per sub-step exceeds the width of the target collider. Resolving this issue requires distinct strategies based on performance cost and accuracy requirements.

• Raycast Sweeping: Project a mathematical line segment from the old position to the new position. Effective for bullets, but fails to handle tumbling angular volumes like flying debris.
• Convex Hull Sweeping: Extrude the convex hull of the moving body along its velocity vector to form a swept volume, then test for intersection against the static scene.
• Conservative Advancement: Continuously step time forward by calculating upper-bound distances to the nearest obstacle, advancing the object safely until time expires or impact occurs.
• Speculative Contacts: Expand the bounding box along velocity vectors and generate negative-distance contact constraints during broad-phase, causing the velocity solver to naturally brake the body before penetration happens.

Evaluating Modern Video Game Physics Engines in 2026

Modern studios rarely build complete physics pipelines from scratch. Selecting an off-the-shelf video game physics engine dictates memory layouts, SIMD acceleration strategies, and execution performance. By 2026, the landscape has coalesced around three distinct industry runtimes: NVIDIA PhysX 5, Unreal Engine 5 Chaos Physics, and the open-source Jolt Physics engine.

Jolt Physics, created by Jorrit Rouwe, has gained massive traction across the gaming industry due to its clean memory efficiency and fully multi-threaded, data-oriented design. Unlike legacy engines weighed down by decades of abstractions, Jolt organizes rigid bodies in memory-dense structures of arrays (SoA), executing collision queries and constraint updates across SIMD-vectorized execution pipelines with zero thread contention.

NVIDIA PhysX 5 remains a formidable runtime, offering robust support for finite element method (FEM) soft body physics, liquids, and complex vehicle dynamics. However, its tight integration with hardware acceleration libraries and larger binary footprint can complicate custom lightweight engine architectures.

Chaos Physics is built deeply into Unreal Engine 5. While early iterations experienced performance regressions compared to UE4 legacy PhysX integration, its modern 2026 iterations provide industry-leading procedural destruction networks (Chaos Destruction) and tight binding to Unreal Niagara visual VFX pipeline.

Selection Metric: For greenfield custom C++ game engines or large open worlds requiring scalable determinism, Jolt Physics is the architectural gold standard. For projects embedded in the Unreal Engine ecosystem, Chaos Physics provides unparalleled authoring workflows despite higher base CPU overhead.

Network Synchronization and Determinism in Multiplayer Physics

Synchronizing physical entities across a distributed, non-zero latency network is among the hardest problems in software engineering. When two players view the same falling crate, sending continuous 60 Hz position packets saturates bandwidth. Instead, network programmers rely on mathematical determinism, client-side prediction, and authoritative server rollbacks.

Server Timeline: (Tick 104 Auth State) ----------------------------> Output Packet
 |
 Network Latency (RTT / 2) |
 v
Client Timeline: [Server Packet In]
(Tick 108 Local State) <-- Prediction Buffer Ahead |
 v
 +------------------------------------+
 | Tick 104 Discrepancy Detected? |
 +------------------------------------+
 / \
 [NO] [YES]
 / \
 Keep Simulating 1. Rewind to Tick 104
 2. Snap to Server State
 3. Resimulate Inputs
 (Ticks 105 -> 108)
 4. Smooth Error Delta

Multiplayer physics synchronization follows strict architectural patterns:

1. Client Input Streaming: The client tags user inputs with sequential tick IDs, applies them locally to the client physics simulation immediately, and transmits the input packet to the server.
2. Input Ring Buffering: The local machine maintains a circular ring buffer containing input history, internal velocity states, and collision bounding boxes for the past 64 to 128 ticks.
3. Authoritative Server Simulation: The server executes the identical sub-step simulation loop using authenticated client inputs, broadcasting lightweight compressed snapshot updates containing physical transforms and angular momentum.
4. State Verification and Rollback: When the client receives the server snapshot for Tick K, it compares the authoritative server state with the saved entry in its circular ring buffer. If positional discrepancies exceed a narrow tolerance threshold, the client executes a rollback: it replaces its local state at Tick K with the server state, then reapplies all stored player inputs from Tick K+1 up to current Tick K+N within a single frame.
5. Visual Error Smoothing: Snapping the physical mesh instantaneously to the corrected transform causes noticeable visual pop. Production code isolates the visual render transform from the collision transform, applying an exponential decay offset to smooth positional corrections over 100 milliseconds.

Executing synchronized physics requires addressing hardware-level non-determinism. Standard IEEE 754 floating-point operations produce slightly divergent values across Intel, AMD, and ARM architectures due to differing microcode, compiler vectorization optimizations (such as fused multiply-add instructions), and variable rounding modes. High-precision lockstep titles eliminate this drift by compiling with strict floating-point flags, using fixed-point math wrappers, or relying on periodic authoritative full-state reconciliation.

Modern game physics operates at the critical intersection of applied mathematics, high-performance computing, and visual art. From designing frame-rate-independent simulation loops using Semi-Implicit Euler integration to managing multi-tiered collision hierarchies with dynamic BVH trees and GJK solvers, maintaining computational stability requires absolute technical precision.

As game worlds expand in scale and physical complexity, modern runtimes prioritize data-oriented design, SIMD execution efficiency, and robust network prediction architectures over monolithic single-threaded calculations. By understanding the mechanical pipeline of numerical integration, constraint solving, and synchronization, engineering teams can build dynamic, reactive interactive worlds that run with uncompromising frame-time stability.