1. The Numerical Integration Challenge: Simulating Motion

1.1 Continuous Calculus vs Discrete Frames

Newton's second law of motion states that force equals mass times acceleration ($F = ma$). In physics simulations, we know the forces acting on a particle (like gravity or wind), allowing us to calculate its acceleration $a = F/m$. To find the particle's velocity and position, we must integrate acceleration over time. In continuous mathematics, velocity is the integral of acceleration, and position is the integral of velocity.

In a game engine, time advances in small discrete steps $\Delta t$ (e.g. $\Delta t \approx 16.67$ milliseconds for a 60 FPS target). We must use numerical approximations to update position and velocity values step-by-step. The choice of the integration algorithm determines the stability, accuracy, and performance of our physics engine.

1.2 Physical State Vectors vs Discrete Ticks

Continuous integration tracks infinite coordinates. Discrete ticks calculate positions using state vectors updated at interval boundaries.

Common Misconception — High frame rates guarantee physics engine stability: A common developer misconception is that a high frame rate alone guarantees physics stability. In reality, if your integration algorithm accumulates energy (like Explicit Euler), running at higher frame rates only delays the inevitable explosion. If the time-step $\Delta t$ increases temporarily during a lag spike, the error accumulates exponentially, causing objects to blast off out of the game world.

2. Euler Integration: The Traditional Approach

2.1 Explicit Euler Instability

The simplest integration method is **Explicit Euler**. At each step, we update position using the current velocity, and update velocity using the current acceleration:

Explicit Euler is popular because it requires only two lines of code. However, it is mathematically unstable (an first-order method). It assumes that velocity remains constant over the entire step $\Delta t$, which is false. In orbit simulations, Explicit Euler constantly overestimates position, adding kinetic energy to the system. Over time, a orbiting planet will spiral away from its sun, violating the conservation of energy.

Mermaid Diagram: The logical update sequence of the Explicit Euler integration model.

3. Verlet Integration: Velocity-Free Physics

3.1 Position-Based Physics

To resolve the instability of Euler, Loup Verlet developed **Verlet Integration** for molecular dynamics. Crucially, Verlet integration does not store velocity variables. Instead, it tracks the particle's **Current Position ($x_t$)** and **Previous Position ($x_{t-\Delta t}$)**.

The velocity is calculated implicitly by subtracting the previous position from the current position. Because velocity is derived directly from actual displacement, the system is self-correcting: it is physically impossible for the velocity to drift away from the coordinates, ensuring excellent energy conservation.

4. The Mathematical Derivation of Verlet Integration

4.1 Taylor Series Expansion

To understand why Verlet is stable, we derive it using Taylor series expansions for the position function $x(t)$ in both the forward and backward directions:

Where $x'(t)$ is velocity $v(t)$, and $x''(t)$ is acceleration $a(t)$. Adding these two equations cancels out the odd-order derivative terms (such as velocity and jerk):

Rearranging the terms to solve for the next position $x(t+\Delta t)$ yields the standard Verlet update equation:

Because the third-order derivative terms cancel out, Verlet integration is a second-order accurate method. The local truncation error is $O(\Delta t^4)$, making it far more accurate than Euler integration with zero additional computational cost.

5. Constraint Solving in Verlet Systems

5.1 Distance Constraint Relaxation

In physics engines, constraint solving coordinates relationships between particles (e.g. keeping two rope segments linked by a fixed distance $L$). In Euler systems, resolving constraints requires calculating complex impulse forces, which often leads to solver failures.

In Verlet systems, we use **Distance Constraint Relaxation** (or Jakobsen's method). Instead of computing forces, we look at the actual distance $D$ between two particles. If $D \neq L$, we calculate the displacement error and shift the particle coordinates directly to restore the target distance $L$. Because velocity is derived from position shifts, this coordinate change automatically adjusts the implicit velocities, resolving constraints stably.

6. Managing Time-Steps: The Dangers of Variable delta t

6.1 Time-Step Glitch Explosions

The standard Verlet update equation assumes that the time-step $\Delta t$ is constant. If $\Delta t$ changes dynamically between frames (due to garbage collection spikes or render lags), the implicit velocity calculation breaks down, causing objects to accelerate or decelerate incorrectly.

To prevent this, physics engines use **Time-Step Correction**. We multiply the displacement term by the ratio of the current time-step to the previous time-step, or use a fixed physics accumulator (such as `FixedUpdate` in Unity) to execute physics updates in uniform ticks.

6.2 Time-Step Correction Formulas and fixed Accumulator Loops

Let $\Delta t$ be the current frame time-step, and let $\Delta t_{\text{prev}}$ be the time-step of the previous frame. If they are unequal, the standard Verlet position update will multiply the previous velocity term incorrectly. The **Time-Step Corrected Verlet** update adjusts the displacement factor by their ratio:

While this correction scale factor reduces glitches, high ratios (e.g. during a major lag spike where $\Delta t / \Delta t_{\text{prev}} > 5$) will still cause numerical integration errors. The industry standard solution is to decouple the rendering frame rate from the physics engine update rate using a **Fixed Time-Step Accumulator**.

Instead of integrating with variable frame times, the engine queues elapsed time into an accumulator buffer. The physics solver then runs updates using a fixed time-step $\Delta t_{\text{physics}}$ (typically 10ms or 16.6ms) as long as the accumulator contains enough time credit:

This fixed accumulator loop guarantees that the physics engine executes identical numerical updates on all machines, regardless of frame rate fluctuations. Below is a comparison table of time integration parameters:

Pitfall — The Spiral of Death: If the physics integration phase takes longer to compute than $\Delta t_{\text{physics}}$, the accumulator queue will grow on every frame. This forces the while loop to run more times in the next frame, further increasing render latency, locking the game. Prevent this by clamping the maximum accumulated time to a safety limit (e.g. 100ms).

7. Advanced: Verlet integration with angular forces and rigid bodies

7.1 Rigid Body Constraints

While Verlet is inherently particle-based, we can model rigid bodies by grouping particles and linking them with rigid distance constraints (sticks). For example, a rigid square can be modeled using 4 particles at the corners and 6 sticks (including diagonals) to prevent shearing.

This eliminates the need for expensive angular velocity vectors and rotation matrices, making Verlet engines ideal for ragdoll systems and destructible structures.

7.2 Distance Constraint Relaxation Vector Mathematics

To solve a distance constraint between two particles at positions $\vec{p}_1 = (x_1, y_1)$ and $\vec{p}_2 = (x_2, y_2)$ linked by a target rest length $L$, we first calculate the difference vector $\vec{\Delta}$ and its Euclidean length $D$:

If $D \neq L$, the constraint is violated. The scalar error difference is $\text{diff} = D - L$. To resolve this error, we shift both particles along the line of action $\vec{\Delta} / D$. If both particles have equal mass, we shift each by half the error:

The updated coordinates are obtained by shifting particle 1 forward and particle 2 backward:

If the particles have different masses $m_1$ and $m_2$, we distribute the displacement shifts inversely proportional to their masses. Let $w_1 = 1/m_1$ and $w_2 = 1/m_2$ be the inverse masses. The displacement scale is normalized by the total inverse mass $w_{\text{total}} = w_1 + w_2$:

In a complex network of constraints (like a cloth grid or a ragdoll skeleton), resolving one constraint will pull other particles, violating neighboring sticks. We solve this by running multiple relaxation sweeps (iterations) per frame. Below is a comparison table of Verlet structural components:

Pitfall — Solver locking via circular constraints: If a loop of rigid distance sticks is over-constrained (e.g. attempting to force diagonals of a triangle to be longer than the sum of sides), the relaxation loop will oscillate infinitely, causing the structure to jitter or explode. Always ensure constraints are physically possible before execution.

8. Advanced: Comparison of Numerical Integration Methods

8.1 Integration Strategy Matrix

The table below compares the mathematical accuracy, stability, and applications of different integration algorithms:

8.2 Symplectic Geometry and Solver Trade-offs

In physics engine design, choosing an integration algorithm requires balancing accuracy, stability, and execution speed. Let $N$ be the number of particles in the simulation. The computational overhead and performance characteristics differ significantly:

• **Explicit Euler**: Requires 1 force evaluation per step, but drifts rapidly, requiring very small $\Delta t$ to avoid explosions.
• **Semi-Implicit Euler**: Symplectic integrator of 1st order. Requires 1 force evaluation. It is highly stable for rigid bodies and character controllers, and is the default choice for general 2D/3D physics engines (like Box2D or PhysX).
• **Verlet Integration**: Symplectic integrator of 2nd order. Requires 1 force evaluation. Position-based shifts make constraint solving extremely stable, making it the standard choice for cloth and soft body dynamics.
• **Runge-Kutta 4 (RK4)**: 4th-order accurate, but requires 4 force evaluations per step. It is computationally expensive, making it unsuitable for hundreds of particles in real-time games, but ideal for high-precision orbit trajectory tracking.

Pitfall — Ignoring collision boundaries during relaxation: When resolving constraint overlaps, shifting a particle's coordinates directly might push it through a thin static collision boundary (like a wall). Always verify boundaries *after* resolving internal constraints, clamping coordinates to valid spaces to prevent clipping glitches.

9. Complete Verlet Physics Engine in JavaScript

9.1 The Class Code

The following JavaScript code implements a basic Verlet particle and distance stick solver, executing coordinate relaxation sweeps:

10. Interactive: Verlet Rope and Particle Simulator

Click "Apply Wind Force" to push the rope particles. Click "Resolve Constraints" to relax the distance offsets back to their target length:

Energy Drift Comparison: Explicit Euler vs Verlet

The chart below compares the total physical energy of a simple pendulum over time using Explicit Euler vs Verlet integration, illustrating how Euler gains energy and explodes:

12. Frequently Asked Questions