Regardless of platform, a robot simulation advances a continuous-time state model \( \mathbf{x}(t) \in \mathbb{R}^n \) driven by controls \( \mathbf{u}(t) \in \mathbb{R}^m \):

\[ \dot{\mathbf{x}}(t) = \mathbf{f}(\mathbf{x}(t), \mathbf{u}(t), t), \quad \mathbf{x}(0)=\mathbf{x}_0. \]

For rigid-body robots (Lesson 2), a typical second-order form uses generalized coordinates \( \mathbf{q}\in\mathbb{R}^d \) and velocities \( \dot{\mathbf{q}}\in\mathbb{R}^d \):

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

where \( \mathbf{M} \) is the mass matrix, \( \mathbf{C} \) collects Coriolis/centrifugal terms, \( \mathbf{g} \) gravity, \( \boldsymbol{\tau} \) actuator torques, and constraints/contacts inject forces via Lagrange multipliers \( \boldsymbol{\lambda} \) and constraint Jacobian \( \mathbf{J} \).

Simulators differ mainly in:

• Constraint/contact model: penalty vs. complementarity / impulse.
• Integrator: explicit, semi-implicit (symplectic), or implicit schemes.
• System representation: maximal coordinates (world-frame bodies) vs. generalized coordinates.
• Robotics tooling: ROS, synthetic sensors, GPU acceleration, scripting APIs.

Discretization produces a map \( \mathbf{x}_{k+1} = \Phi_h(\mathbf{x}_k,\mathbf{u}_k) \) with time step \( h \). Two common families:

Explicit Euler (simple, cheaper, less stable):

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

Semi-implicit (symplectic) Euler (better energy behavior):

\[ \dot{\mathbf{q}}_{k+1} = \dot{\mathbf{q}}_{k} + h\,\ddot{\mathbf{q}}(\mathbf{q}_k,\dot{\mathbf{q}}_k), \quad \mathbf{q}_{k+1} = \mathbf{q}_k + h\,\dot{\mathbf{q}}_{k+1}. \]

Stability fact (linear systems). For a stable continuous LTI system \( \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} \) (students know linear control), explicit Euler yields \( \mathbf{x}_{k+1} = (\mathbf{I}+h\mathbf{A})\mathbf{x}_k \). Discrete stability requires all eigenvalues of \( \mathbf{I}+h\mathbf{A} \) lie inside the unit disk: \( |1+h\lambda_i(\mathbf{A})| < 1 \).

Proof sketch.

Diagonalize (or Jordan-form) \( \mathbf{A} \) so each mode evolves as \( x_{k+1} = (1+h\lambda)x_k \). Convergence to zero occurs iff \( |1+h\lambda| < 1 \). For real negative \( \lambda=-\alpha,\alpha>0 \), this gives \( -1 < 1-h\alpha < 1 \Rightarrow 0 < h < 2/\alpha \). Hence stiff dynamics (large \( \alpha \)) demand small \( h \).

Engines emphasizing high-speed contact/legged locomotion typically adopt semi-implicit or implicit solvers to allow larger stable steps. This axis will help interpret differences below.

3.1 Gazebo (Classic / Ignition / Gazebo Sim)

Gazebo is a robotics-first simulator designed around modular physics engines. It uses maximal coordinate rigid-body dynamics and supports multiple backends (ODE, Bullet, DART, Simbody), so the actual dynamics map \( \Phi_h \) depends on the chosen backend.

Mathematically, Gazebo solves the constrained dynamics of Section 1 via either:

• Penalty contacts: approximate \( \mathbf{J}\ddot{\mathbf{q}}+\dot{\mathbf{J}}\dot{\mathbf{q}} \ge 0 \) with springs/dampers, giving smoother but potentially “soft” contact.
• Impulse/constraint methods: solve for \( \boldsymbol{\lambda} \) explicitly per step.

Strengths: tight ROS/ROS2 integration, large model zoo, sensor plugins. Typical use: system-level robotics testing and digital-twin-style pipelines.

3.2 Isaac Sim

Isaac Sim is a high-fidelity simulator built on NVIDIA PhysX, emphasizing GPU-accelerated rigid-body/contact and photorealistic rendering. PhysX uses a projected Gauss–Seidel style constraint solver:

\[ \text{find } \boldsymbol{\lambda} \text{ such that } \mathbf{A}\boldsymbol{\lambda} + \mathbf{b} \ge 0,\; \boldsymbol{\lambda}\ge 0,\; \boldsymbol{\lambda}^T(\mathbf{A}\boldsymbol{\lambda}+\mathbf{b})=0, \]

which is a discrete complementarity form of contact (Lesson 2). This supports stable stacking and frictional contact at larger time steps.

Strengths: scalability, synthetic sensor data, domain randomization for sim-to-real experiments (conceptual only here), and ROS bridges.

3.3 MuJoCo

MuJoCo uses generalized coordinates with an efficient articulated-body algorithm. Contacts are modeled with smooth soft constraints rather than hard impacts, often described by a penalty potential \( V_c(\mathbf{q}) \) producing forces \( \mathbf{f}_c = -\nabla_{\mathbf{q}}V_c \).

\[ \mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} = \boldsymbol{\tau} - \nabla_{\mathbf{q}} V_c(\mathbf{q}) + \cdots \]

This yields differentiable dynamics and stable control/optimization experiments, especially for legged systems, but contacts may feel more compliant.

Strengths: fast, numerically robust, excellent for control prototyping and research on contact-rich locomotion.

3.4 PyBullet

PyBullet is a Python API to the Bullet engine. Bullet uses a sequential impulse constraint method per time step:

\[ \Delta\dot{\mathbf{q}} = \mathbf{M}^{-1}\mathbf{J}^T\boldsymbol{\Lambda}, \quad \boldsymbol{\Lambda} = \arg\min_{\boldsymbol{\lambda}\ge 0} \| \mathbf{J}\dot{\mathbf{q}} + \mathbf{b}\|^2, \]

where \( \boldsymbol{\Lambda} \) is the impulse vector approximated iteratively. PyBullet prioritizes accessibility and quick scripting over photorealism.

Strengths: simple Python experiments, URDF support, good for education and rapid tests.

3.5 Quick Comparison

flowchart TD
  S["Start: define simulation goal"] --> G1["Need ROS-based system test?"]
  G1 -->|yes| ZB["Prefer Gazebo"]
  G1 -->|no| G2["Need photorealism / GPU scale?"]
  G2 -->|yes| IS["Prefer Isaac Sim"]
  G2 -->|no| G3["Need fast contact-rich control research?"]
  G3 -->|yes| MJ["Prefer MuJoCo"]
  G3 -->|no| PB["Prefer PyBullet for rapid scripting"]
        

This choice is not absolute. For example, a digital twin workflow might prototype algorithms in PyBullet or MuJoCo, then validate with Gazebo or Isaac depending on deployment needs.

The following snippets all implement the same idea: load a simple rigid body model, step the simulator, and read state. These are intentionally minimal so you can map the common loop: initialize → load model → for k: apply u_k → step → read x_k.

5.1 Python — PyBullet

import pybullet as p
import pybullet_data
import time

cid = p.connect(p.GUI)  # or p.DIRECT
p.setAdditionalSearchPath(pybullet_data.getDataPath())
p.setGravity(0, 0, -9.81)

plane = p.loadURDF("plane.urdf")
box = p.loadURDF("r2d2.urdf", [0, 0, 0.2])

dt = 1.0/240.0
p.setTimeStep(dt)

for k in range(1000):
    # Example control: zero torques (free fall)
    p.stepSimulation()
    pos, orn = p.getBasePositionAndOrientation(box)
    if k % 120 == 0:
        print("t=", k*dt, "pos=", pos)
    time.sleep(dt)
p.disconnect()

5.2 C++ — Gazebo (plugin skeleton)

// Minimal Gazebo model plugin (classic-style API)
#include <gazebo/gazebo.hh>
#include <gazebo/physics/physics.hh>

namespace gazebo {
class HelloSim : public ModelPlugin {
  physics::ModelPtr model;
  event::ConnectionPtr updateConn;

public:
  void Load(physics::ModelPtr _model, sdf::ElementPtr) {
    model = _model;
    updateConn = event::Events::ConnectWorldUpdateBegin(
      std::bind(&HelloSim::OnUpdate, this));
  }

  void OnUpdate() {
    // Read pose each sim step
    auto pose = model->WorldPose();
    std::cout << "z=" << pose.Pos().Z() << std::endl;
  }
};

GZ_REGISTER_MODEL_PLUGIN(HelloSim)
}

5.3 Java — Bullet (via jBullet / libGDX Bullet bindings)

// Pseudocode for Java Bullet bindings: set up world and step
DynamicsWorld world = new DiscreteDynamicsWorld(...);
world.setGravity(new Vector3(0f, 0f, -9.81f));

// Create a rigid body (box)
RigidBody box = RigidBodyFactory.createBox(mass, size, initialTransform);
world.addRigidBody(box);

float dt = 1f/240f;
for(int k=0; k<1000; k++){
    world.stepSimulation(dt, 10);
    Transform tr = box.getWorldTransform();
    System.out.println("t=" + k*dt + " z=" + tr.origin.z);
}

5.4 Matlab/Simulink — stepping a rigid-body plant (conceptual)

% Simple discrete simulation loop consistent with engine stepping
dt = 1/240;
T  = 4;                 % seconds
N  = round(T/dt);

% Example: 1D point mass with gravity (as a toy proxy)
m = 1; g = 9.81;
x = zeros(N,1); v = zeros(N,1);
x(1) = 0.2; v(1) = 0;

for k = 1:N-1
    a = -g;                 % no control
    v(k+1) = v(k) + dt*a;   % semi-implicit Euler
    x(k+1) = x(k) + dt*v(k+1);
end

plot((0:N-1)*dt, x), xlabel('t'), ylabel('height')

Note how the Matlab loop mirrors simulator structure. Real robotics models are loaded via URDF/SDF or engine-specific model formats, but the numerical logic is the same.

Digital twins (Lesson 5) refine \( \Phi_h \) and sensor models to match the real robot. At the loop level, every simulator provides:

• State access \( \mathbf{x}_k \)
• Control injection \( \mathbf{u}_k \)
• Sensor outputs \( \mathbf{y}_k = \mathbf{h}(\mathbf{x}_k) \)
• Parameter hooks to tune \( \mathbf{M},\mathbf{C},\mathbf{g}, \) and contact terms

flowchart TD
  R["Real robot"] -->|measure y_k| ID["Parameter identification"]
  ID --> SIM["Simulator (Gazebo/Isaac/MuJoCo/PyBullet)"]
  SIM -->|predict x_k, y_k| VAL["Validation vs. real data"]
  VAL -->|update params| SIM
  SIM --> CTRL["Controller / software under test"]
  CTRL -->|u_k| SIM
        

Problem 1 (Euler stability bound): Consider the scalar stable system \( \dot{x} = -\alpha x \) with \( \alpha>0 \). Using explicit Euler, find the largest step size \( h_{\max} \) that preserves stability.

Solution:

Explicit Euler gives \( x_{k+1} = (1-h\alpha)x_k \). Stability requires \( |1-h\alpha| < 1 \). Hence \( -1 < 1-h\alpha < 1 \Rightarrow 0 < h < 2/\alpha \). Therefore \( h_{\max} = 2/\alpha \).

Problem 2 (Energy drift comparison): A unit-mass spring system obeys \( \ddot{x} + \omega^2 x = 0 \). Show that explicit Euler causes energy to increase for any nonzero step, whereas semi-implicit Euler does not.

Solution:

Let state \( \mathbf{x}=[x\;\;v]^T \), with \( \dot{x}=v,\;\dot{v}=-\omega^2 x \). Explicit Euler: \( v_{k+1}=v_k-h\omega^2 x_k,\; x_{k+1}=x_k+h v_k \). Mechanical energy \( E_k=\tfrac12(v_k^2+\omega^2 x_k^2) \). Substitute and expand:

\[ E_{k+1} = E_k + \tfrac12 h^2\omega^2(v_k^2+\omega^2 x_k^2) = E_k\left(1+h^2\omega^2\right) > E_k \quad \text{for } h\neq 0. \]

So energy strictly grows. Semi-implicit Euler: \( v_{k+1}=v_k-h\omega^2 x_k,\; x_{k+1}=x_k+h v_{k+1} \). A similar expansion yields bounded oscillatory energy with \( E_{k+1}-E_k = \mathcal{O}(h^3) \), so no systematic growth. This explains why many engines prefer semi-implicit schemes.

Problem 3 (Normal impulse with restitution): Two bodies collide along a contact normal with pre-impact relative normal velocity \( v_n^- < 0 \). With restitution coefficient \( e\in[0,1] \), the post-impact normal velocity satisfies \( v_n^+ = -e v_n^- \). Derive the impulse magnitude \( \Lambda \) assuming effective mass \( m_\text{eff} \).

Solution:

Impulse changes normal velocity by \( v_n^+ = v_n^- + \Lambda/m_\text{eff} \). Enforce restitution: \( -e v_n^- = v_n^- + \Lambda/m_\text{eff} \). Hence \( \Lambda = -(1+e)m_\text{eff} v_n^- \). Sequential-impulse engines (Bullet/PyBullet) approximate this iteratively each step.

Problem 4 (Simulator choice reasoning): You need to test a ROS2 navigation stack with noisy LiDAR and camera sensors, but not photorealistic rendering. Which simulator is most appropriate and why?

Solution:

Gazebo is most appropriate. It is ROS-native, supports standard robot model formats (URDF/SDF), and provides configurable sensor plugins. Since photorealism is not required, the extra GPU rendering strengths of Isaac Sim are unnecessary; Gazebo’s tooling aligns directly with the software-in-the-loop goal.

We compared Gazebo, Isaac Sim, MuJoCo, and PyBullet by expressing all of them as discrete solvers of the same constrained robot dynamics. Differences arise from coordinate representations, contact formulations, and integrator choices, which determine stability and realism for a given step size. You now have a principled basis to select a simulator and to interpret simulation errors as consequences of \( \Phi_h \), not “mystery behavior.”