A physics engine is a numerical solver that advances a simulated world state \( \mathbf{x}(t) \) (positions, orientations, velocities, etc.) forward in time by applying physical laws plus approximations for events such as contact. In robotics, this includes:

• Rigid-body motion under forces and torques.
• Kinematic constraints (joints, closed chains).
• Collisions and contacts with friction.
• Actuation models (motors, controllers) at a chosen abstraction level.
• Sometimes soft bodies, fluids, and deformable contact (usually approximate).

flowchart TD
  S0["State at time k: q_k, v_k"] --> F["Compute forces/torques: gravity, motors, contact"]
  F --> CD["Collision detection"]
  CD --> C0["Build constraints (joints, contacts)"]
  C0 --> SOL["Solve constrained dynamics"]
  SOL --> INT["Integrate to get q_(k+1), v_(k+1)"]
  INT --> S1["State at time k+1"]
        

The engine chooses a model class (typically rigid-body with constraints) and a numerical time-stepping method. These choices determine realism, stability, and computational cost.

For a rigid body with position \( \mathbf{p}\in\mathbb{R}^3 \), orientation \( \mathbf{R}\in SO(3) \), linear velocity \( \mathbf{v} \), and angular velocity \( \boldsymbol{\omega} \), Newton–Euler equations are:

\[ m\dot{\mathbf{v}} = \mathbf{f}_{\text{ext}}, \qquad \mathbf{I}\dot{\boldsymbol{\omega}} + \boldsymbol{\omega}\times (\mathbf{I}\boldsymbol{\omega}) = \boldsymbol{\tau}_{\text{ext}}, \]

where \( m \) is mass, \( \mathbf{I} \) is inertia in the body frame, \( \mathbf{f}_{\text{ext}} \) and \( \boldsymbol{\tau}_{\text{ext}} \) include gravity, actuation, and contact impulses.

Engines typically stack many bodies into generalized coordinates \( \mathbf{q}\in\mathbb{R}^n \) and velocities \( \mathbf{v}=\dot{\mathbf{q}}\in\mathbb{R}^n \). The compact rigid-body system is:

\[ \mathbf{M}(\mathbf{q})\dot{\mathbf{v}} + \mathbf{C}(\mathbf{q},\mathbf{v})\mathbf{v} + \mathbf{g}(\mathbf{q}) = \boldsymbol{\tau} + \mathbf{J}(\mathbf{q})^\top \boldsymbol{\lambda}. \]

Here \( \mathbf{M} \) is the mass matrix, \( \mathbf{C}\mathbf{v} \) collects Coriolis/centrifugal terms, \( \mathbf{g} \) gravity forces, \( \boldsymbol{\tau} \) actuator torques/forces, and \( \mathbf{J}^\top\boldsymbol{\lambda} \) enforces constraints via Lagrange multipliers.

Joints are encoded as holonomic constraints \( \boldsymbol{\phi}(\mathbf{q})=\mathbf{0} \). Differentiating gives velocity and acceleration constraints:

\[ \mathbf{J}(\mathbf{q})\mathbf{v} = \mathbf{0}, \qquad \mathbf{J}(\mathbf{q})\dot{\mathbf{v}} + \dot{\mathbf{J}}(\mathbf{q},\mathbf{v})\mathbf{v}=\mathbf{0}, \]

where \( \mathbf{J}(\mathbf{q})=\partial\boldsymbol{\phi}/\partial\mathbf{q} \). Substituting into the rigid-body equations yields the saddle-point system:

\[ \begin{bmatrix} \mathbf{M} & -\mathbf{J}^\top \\ \mathbf{J} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \dot{\mathbf{v}} \\ \boldsymbol{\lambda} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\tau} - \mathbf{C}\mathbf{v} - \mathbf{g} \\ -\dot{\mathbf{J}}\mathbf{v} \end{bmatrix}. \]

Many engines solve this linearly each time step, or reduce it using Schur complements.

Interpretation: \( \boldsymbol{\lambda} \) are constraint forces; they do no virtual work along feasible motions. For any virtual displacement \( \delta\mathbf{q} \) with \( \mathbf{J}\delta\mathbf{q}=\mathbf{0} \),

\[ (\mathbf{J}^\top\boldsymbol{\lambda})^\top \delta\mathbf{q} = \boldsymbol{\lambda}^\top (\mathbf{J}\delta\mathbf{q}) = 0, \]

showing constraint forces do not change energy directly; they only restrict motion.

Collision detection finds candidate contact points. Contact modeling then enforces non-penetration. Let \( g(\mathbf{q}) \) be a signed gap function (positive when separated). Ideal rigid contact:

\[ g(\mathbf{q}) \ge 0,\qquad \lambda_n \ge 0,\qquad g(\mathbf{q})\lambda_n = 0. \]

These complementarity conditions mean either the bodies are apart (\( g>0 \Rightarrow \lambda_n=0 \)) or touching (\( g=0 \Rightarrow \lambda_n\ge 0 \)). Most engines approximate contact using:

• Penalty models (soft contact): \( \lambda_n = k\,\max(0,-g) \)
• Impulse / complementarity solvers (hard contact).

For Coulomb friction, tangential contact force \( \boldsymbol{\lambda}_t \) satisfies:

\[ \|\boldsymbol{\lambda}_t\| \le \mu \lambda_n, \quad \boldsymbol{\lambda}_t = -\mu \lambda_n \frac{\mathbf{v}_t}{\|\mathbf{v}_t\|} \ \text{when}\ \mathbf{v}_t \neq \mathbf{0}, \]

where \( \mathbf{v}_t \) is tangential relative velocity, and \( \mu \) is the friction coefficient. Numerically, this creates a nonlinear complementarity problem, often solved by iterative projection (PGS) or convex optimization methods.

Continuous dynamics are discretized with time step \( h \). A generic first-order system \( \dot{\mathbf{x}} = f(\mathbf{x},t) \) is integrated as:

\[ \mathbf{x}_{k+1} = \mathbf{x}_k + h\,\Phi(\mathbf{x}_k,\mathbf{x}_{k+1},t_k), \]

where \( \Phi \) defines the method:

• Explicit Euler: \( \Phi=f(\mathbf{x}_k,t_k) \)
• Implicit Euler: \( \Phi=f(\mathbf{x}_{k+1},t_{k+1}) \)
• Semi-implicit (symplectic) Euler: velocity updated explicitly, position updated using new velocity.

To connect with your Linear Control background, consider the scalar harmonic oscillator: \( \ddot{q} + \omega_0^2 q = 0 \). Define \( x_1=q \), \( x_2=\dot{q} \):

\[ \dot{\mathbf{x}} = \begin{bmatrix} 0 & 1 \\ -\omega_0^2 & 0 \end{bmatrix}\mathbf{x} \equiv \mathbf{A}\mathbf{x}. \]

Explicit Euler gives \( \mathbf{x}_{k+1}=(\mathbf{I}+h\mathbf{A})\mathbf{x}_k \). The eigenvalues of \( \mathbf{I}+h\mathbf{A} \) are \( 1 \pm j h\omega_0 \) with magnitude \( \sqrt{1+h^2\omega_0^2} > 1 \), so energy grows unphysically.

Claim (stability advantage of symplectic Euler): For the same oscillator, symplectic Euler preserves bounded energy (no exponential drift).

Proof sketch: Symplectic Euler updates: \( x_{2,k+1} = x_{2,k} - h\omega_0^2 x_{1,k} \), \( x_{1,k+1} = x_{1,k} + h x_{2,k+1} \). In matrix form:

\[ \mathbf{x}_{k+1} = \underbrace{\begin{bmatrix} 1-h^2\omega_0^2 & h \\ -h\omega_0^2 & 1 \end{bmatrix}}_{\mathbf{A}_s} \mathbf{x}_k. \]

The characteristic polynomial of \( \mathbf{A}_s \) is:

\[ \chi(\lambda)=\lambda^2 - 2\left(1-\tfrac{1}{2}h^2\omega_0^2\right)\lambda + 1. \]

Its roots satisfy \( \lambda_1\lambda_2=1 \). When \( h\omega_0 < 2 \), the discriminant is negative, so the roots are complex conjugates on the unit circle (\( |\lambda|=1 \)), giving bounded oscillations. Hence symplectic Euler avoids systematic energy explosion typical in explicit Euler.

This is why many real-time engines prefer symplectic or implicit methods: stability is crucial when contacts make dynamics stiff.

Physics engines are not full scientific simulators. Typical simplifications:

• Ideal rigidity: bodies do not deform unless a soft-body module is used.
• Simple friction: Coulomb or regularized variants; no micro-slip modeling.
• No detailed fluid/air interaction unless specialized plugins exist.
• Motor models abstracted: torque sources, PD servos, or ideal velocity drives.
• Discrete contacts: time-step contacts instead of continuous collision evolution.

flowchart LR
  R["Rigid bodies"] -->|modeled well| DYN["Newton-Euler / M(q)vdot"]
  C["Contacts"] -->|approx.| LCP["Complementarity / penalty"]
  F["Friction"] -->|approx.| MU["Coulomb cone"]
  S["Soft bodies, fluids"] -->|often ignored or coarse| APX["Extra modules"]
        

For digital twins, these approximations must be aligned with the real robot: if flexibility or fluid drag matters for your task, a basic rigid engine may be insufficient.

We simulate a 1D point-mass "robot foot" impacting the ground using a soft penalty model. This mirrors how many engines stabilize contacts in real time.

import numpy as np

# Parameters
m = 2.0                 # mass (kg)
g = 9.81                # gravity (m/s^2)
k = 5000.0              # contact stiffness (N/m)
d = 50.0                # contact damping (N*s/m)
h = 1e-3                # time step (s)
T = 1.0                 # total time

# State: position q (m), velocity v (m/s)
q, v = 0.3, -1.0        # start above ground, moving down

def contact_force(q, v):
    # ground at q=0; gap g(q)=q
    if q >= 0:
        return 0.0
    # penalty: lambda_n = k(-q) - d v (only when penetrating)
    return k * (-q) - d * v

qs, vs, ts = [], [], []
t = 0.0
while t <= T:
    f_ext = -m*g + contact_force(q, v)
    a = f_ext / m

    # symplectic Euler
    v = v + h * a
    q = q + h * v

    qs.append(q); vs.append(v); ts.append(t)
    t += h

print("Final position, velocity:", q, v)
      

Even this toy model highlights why contact stiffness and step size matter: large \( k \) with large \( h \) causes instability.

// Minimal 3D rigid-body translational integrator (no rotation)
// Requires Eigen: https://eigen.tuxfamily.org/
#include <iostream>
#include <Eigen/Dense>

struct Body {
  double m;
  Eigen::Vector3d p;  // position
  Eigen::Vector3d v;  // velocity
};

int main() {
  Body b;
  b.m = 1.5;
  b.p = Eigen::Vector3d(0,0,1);
  b.v = Eigen::Vector3d(0,0,0);

  Eigen::Vector3d g(0,0,-9.81);
  double h = 0.001;

  for (int k=0; k<10000; ++k) {
    Eigen::Vector3d f_ext = b.m * g; // gravity only
    Eigen::Vector3d a = f_ext / b.m;

    // semi-implicit Euler
    b.v += h * a;
    b.p += h * b.v;
  }

  std::cout << "p=" << b.p.transpose() << " v=" << b.v.transpose() << std::endl;
  return 0;
}
      

Real engines add rotational dynamics, constraints, and collision impulses, but the integration pattern is the same.

A basic impulse update for two disks colliding along normal \( \mathbf{n} \).

class Vec2 {
    public double x, y;
    Vec2(double x, double y){ this.x=x; this.y=y; }
    Vec2 add(Vec2 o){ return new Vec2(x+o.x, y+o.y); }
    Vec2 sub(Vec2 o){ return new Vec2(x-o.x, y-o.y); }
    Vec2 mul(double s){ return new Vec2(s*x, s*y); }
    double dot(Vec2 o){ return x*o.x + y*o.y; }
    double norm(){ return Math.sqrt(x*x + y*y); }
    Vec2 unit(){ double n=norm(); return new Vec2(x/n, y/n); }
}

public class DiskCollision {
    static void collide(
        double m1, Vec2 v1,
        double m2, Vec2 v2,
        Vec2 n, double restitution
    ){
        // relative normal speed
        double vn = v1.sub(v2).dot(n);
        if (vn >= 0) return; // separating

        double j = -(1.0 + restitution) * vn / (1.0/m1 + 1.0/m2);
        Vec2 impulse = n.mul(j);

        // apply impulses
        v1.x += impulse.x / m1;  v1.y += impulse.y / m1;
        v2.x -= impulse.x / m2;  v2.y -= impulse.y / m2;
    }
}
      

Engines generalize this to many contacts and masses via a global solver.

In MATLAB, rigid multibody simulation is typically done with Simscape Multibody. Conceptually, it solves:

\[ \mathbf{M}(\mathbf{q})\dot{\mathbf{v}}=\boldsymbol{\tau}_{\text{applied}}+\mathbf{J}^\top\boldsymbol{\lambda}, \quad \boldsymbol{\phi}(\mathbf{q})=\mathbf{0}, \]

using variable-step implicit integrators for stiffness. A simple discrete-time script version of symplectic Euler is:

m = 2.0; g = 9.81; h = 1e-3; T = 1.0;
q = 0.3; v = -1.0;
k = 5000; d = 50;

contact_force = @(q,v) (q >= 0) * 0 + (q < 0) * (k*(-q) - d*v);

ts = 0:h:T;
qs = zeros(size(ts)); vs = zeros(size(ts));
for i=1:length(ts)
    f_ext = -m*g + contact_force(q,v);
    a = f_ext/m;

    v = v + h*a;  % symplectic Euler
    q = q + h*v;

    qs(i)=q; vs(i)=v;
end
disp([q v])
      

Problem 1 (Constraint Forces): Consider a particle of mass \( m \) constrained to a circle of radius \( R \) in the plane. Let \( \boldsymbol{\phi}(\mathbf{q}) = x^2+y^2-R^2=0 \). Derive the Lagrange multiplier equation and show that the constraint force acts radially.

Solution: We have \( \mathbf{q}=[x\ y]^\top \) and \( \mathbf{J}=\nabla\phi = [2x\ 2y] \). The constrained dynamics:

\[ m\ddot{\mathbf{q}} = \mathbf{f}_{\text{ext}} + \mathbf{J}^\top\lambda. \]

Since \( \mathbf{J}^\top\lambda = \lambda[2x\ 2y]^\top = 2\lambda\,\mathbf{q} \), the constraint force is proportional to \( \mathbf{q} \), i.e., purely radial. Its role is to cancel the radial component of \( \mathbf{f}_{\text{ext}} \) and keep \( x^2+y^2=R^2 \).

Problem 2 (Explicit Euler Instability): For the harmonic oscillator \( \ddot{q}+\omega_0^2 q=0 \), show that explicit Euler causes energy growth for any \( h>0 \).

Solution: As in Section 5, \( \mathbf{x}_{k+1}=(\mathbf{I}+h\mathbf{A})\mathbf{x}_k \), with eigenvalues \( 1\pm jh\omega_0 \). Their magnitude is:

\[ |1\pm jh\omega_0|=\sqrt{1+h^2\omega_0^2} > 1. \]

Thus \( \|\mathbf{x}_k\| \) grows approximately like \( (\sqrt{1+h^2\omega_0^2})^k \) implying energy drift upward.

Problem 3 (Complementarity Contact): A mass is above a rigid ground with gap \( g=q \). Write the complementarity conditions and interpret each case physically.

Solution: The conditions are:

\[ g(q)\ge 0,\quad \lambda_n\ge 0,\quad g(q)\lambda_n=0. \]

If \( g>0 \) (separated), then necessarily \( \lambda_n=0 \) (no normal force). If \( g=0 \) (touching), then \( \lambda_n\ge 0 \), preventing penetration. This matches ideal rigid contact.

Problem 4 (Penalty Contact Tuning): In the penalty model \( \lambda_n = k(-g) - d v \) for \( g<0 \), give a sufficient condition on \( h \) to avoid numerical oscillation in a 1D mass–spring contact.

Solution: When penetrating, the dynamics approximate a damped spring: \( m\ddot{q}+d\dot{q}+kq=0 \) (with \( q<0 \)). Its natural frequency is \( \omega_c=\sqrt{k/m} \). A conservative stability guideline for explicit/semi-implicit steps is:

\[ h\omega_c < 2 \quad \Rightarrow \quad h < 2\sqrt{\frac{m}{k}}. \]

This ensures the discrete eigenvalues remain inside or on the unit circle, preventing spurious high-frequency blow-up.

Physics engines approximate robot-world interaction by combining rigid-body dynamics, constraint forces, collision/contact models, friction laws, and stable time-stepping. Different engines mainly differ in contact/constraint solvers and integrators, which explains why simulated behaviors can vary. Understanding these models lets you choose, tune, and trust simulations when building digital twins.