How to Use Lagrangian Mechanics to Build a Game Physics Engine

Building a physics engine from scratch is a formidable task. Most developers jump straight to Newtonian mechanics, dealing with forces, impulses, and explicit constraint solvers. While effective, this approach can become incredibly complex for systems with many interconnected parts or rigid constraints.

Enter Lagrangian Mechanics. Often overlooked in introductory game development, it offers an elegant, powerful, and sometimes simpler way to derive the equations of motion for complex systems. This post will demystify Lagrangian mechanics and show you how to leverage its power to build more robust and maintainable game physics.

Why Lagrangian Mechanics for Games?

Traditional Newtonian mechanics operates on forces and accelerations ($F=ma$). When you have objects connected by joints, rods, or moving on specific surfaces, you need to introduce constraint forces to enforce these rules. Calculating these constraint forces directly can be tricky and lead to numerical instability.

Lagrangian mechanics, however, offers a different perspective. Instead of forces, it focuses on energy. By defining the system’s kinetic and potential energies, you can derive the equations of motion without ever explicitly calculating constraint forces. This is particularly powerful for systems with “holonomic constraints” – constraints that can be expressed as an algebraic equation relating the coordinates.

The core advantages for game development are:

• Implicit Constraint Handling: Constraints are often handled naturally by choosing appropriate “generalized coordinates,” simplifying the problem significantly.
• Analytical Rigor: It guarantees energy conservation for conservative systems, leading to more stable simulations over long periods (unlike some force-based methods that can slowly “leak” energy).
• Simplified Derivation: For complex multi-body systems, deriving the equations of motion can be more straightforward than summing forces and torques.

The Core Idea: The Principle of Least Action

At the heart of Lagrangian mechanics is the Principle of Least Action. It states that the path a physical system takes between two points in time is the one that minimizes a quantity called “Action.” The Action ($S$) is the integral of the Lagrangian ($L$) over time:

$S = \int_{t_1}^{t_2} L(q, \dot{q}, t) , dt$

Where:

• $L$ is the Lagrangian.
• $q$ represents the generalized coordinates (positions, angles, etc.).
• $\dot{q}$ represents the generalized velocities (rates of change of generalized coordinates).
• $t$ is time.

Don’t worry, we won’t be directly minimizing this integral in our code. Instead, we’ll use a direct consequence of this principle: the Euler-Lagrange equation.

Defining the Lagrangian ($L$)

The Lagrangian ($L$) is simply defined as the difference between the system’s total Kinetic Energy ($T$) and its total Potential Energy ($V$):

$L = T - V$

• Kinetic Energy ($T$): The energy due to motion. For a particle of mass $m$ with velocity $v$, $T = \frac{1}{2}mv^2$. For a rotating body, it includes rotational kinetic energy.
• Potential Energy ($V$): The energy due to position or configuration in a force field. For gravity near the Earth’s surface, $V = mgh$ (mass * gravity * height). For a spring, $V = \frac{1}{2}kx^2$.

Let’s illustrate with a simple example: a particle falling under gravity. Let $y$ be the vertical position.

• Velocity: $\dot{y}$
• Kinetic Energy: $T = \frac{1}{2}m\dot{y}^2$
• Potential Energy: $V = mgy$ (assuming $y=0$ is the reference potential energy level)
• Lagrangian: $L = \frac{1}{2}m\dot{y}^2 - mgy$

The Euler-Lagrange Equation

The magic happens with the Euler-Lagrange equation. This equation provides a way to derive the equations of motion for each generalized coordinate without needing to sum forces. For each generalized coordinate $q_i$, the Euler-Lagrange equation is:

$\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} = 0$

Let’s apply this to our falling particle example where $q_1 = y$:

1. Calculate $\frac{\partial L}{\partial \dot{y}}$: $\frac{\partial}{\partial \dot{y}} \left( \frac{1}{2}m\dot{y}^2 - mgy \right) = m\dot{y}$

2. Calculate $\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{y}} \right)$: $\frac{d}{dt} (m\dot{y}) = m\ddot{y}$

3. Calculate $\frac{\partial L}{\partial y}$: $\frac{\partial}{\partial y} \left( \frac{1}{2}m\dot{y}^2 - mgy \right) = -mg$

4. Substitute into Euler-Lagrange equation: $m\ddot{y} - (-mg) = 0$ $m\ddot{y} + mg = 0$ $\ddot{y} = -g$

This is exactly what we expect from Newtonian mechanics: the acceleration due to gravity. Simple, right? The real power emerges when systems become more complex.

Generalized Coordinates

Choosing the right generalized coordinates ($q_i$) is crucial. They are the minimum set of independent variables required to describe the system’s configuration. They implicitly handle constraints.

• Example 1: Particle on a wire: If a bead is constrained to move along a circular wire, instead of using Cartesian $(x,y)$ coordinates and then enforcing $x^2+y^2=R^2$ with constraint forces, you can simply use a single angle $\theta$ as your generalized coordinate. The constraint is gone!
• Example 2: Simple Pendulum: Instead of $(x,y)$ coordinates for the bob and then a constraint $x^2+y^2=L^2$ (where $L$ is rod length), use the angle $\theta$ from the vertical.

Applying to Game Physics: Simple Pendulum Example

Let’s build a simple physics engine for a pendulum. We’ll use Python with sympy for symbolic differentiation and scipy for numerical integration.

System Setup

• A point mass $m$ at the end of a rigid, massless rod of length $L$.
• The rod is pivoted at the origin $(0,0)$.
• Gravity acts downwards ($g$).
• Generalized coordinate: $\theta$, the angle the rod makes with the positive y-axis (downwards).

1. Define Kinetic Energy ($T$)

The position of the mass is $(x, y) = (L \sin\theta, -L \cos\theta)$. The velocity components are: $\dot{x} = L \cos\theta \dot{\theta}$ $\dot{y} = L \sin\theta \dot{\theta}$

Kinetic Energy: $T = \frac{1}{2}m(\dot{x}^2 + \dot{y}^2)$ $T = \frac{1}{2}m((L \cos\theta \dot{\theta})^2 + (L \sin\theta \dot{\theta})^2)$ $T = \frac{1}{2}m(L^2 \cos^2\theta \dot{\theta}^2 + L^2 \sin^2\theta \dot{\theta}^2)$ $T = \frac{1}{2}mL^2\dot{\theta}^2(\cos^2\theta + \sin^2\theta)$ $T = \frac{1}{2}mL^2\dot{\theta}^2$

2. Define Potential Energy ($V$)

Potential Energy: $V = mgy$. If we take $y=0$ at the pivot point, then $y = -L \cos\theta$. $V = mg(-L \cos\theta) = -mgL \cos\theta$

3. Define the Lagrangian ($L$)

$L = T - V = \frac{1}{2}mL^2\dot{\theta}^2 - (-mgL \cos\theta)$ $L = \frac{1}{2}mL^2\dot{\theta}^2 + mgL \cos\theta$

4. Apply Euler-Lagrange Equation

For $q_1 = \theta$: $\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}} \right) - \frac{\partial L}{\partial \theta} = 0$

1. Calculate $\frac{\partial L}{\partial \dot{\theta}}$: $\frac{\partial}{\partial \dot{\theta}} \left( \frac{1}{2}mL^2\dot{\theta}^2 + mgL \cos\theta \right) = mL^2\dot{\theta}$

2. Calculate $\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\theta}} \right)$: $\frac{d}{dt} (mL^2\dot{\theta}) = mL^2\ddot{\theta}$

3. Calculate $\frac{\partial L}{\partial \theta}$: $\frac{\partial}{\partial \theta} \left( \frac{1}{2}mL^2\dot{\theta}^2 + mgL \cos\theta \right) = -mgL \sin\theta$

4. Substitute into Euler-Lagrange equation: $mL^2\ddot{\theta} - (-mgL \sin\theta) = 0$ $mL^2\ddot{\theta} + mgL \sin\theta = 0$

This is the equation of motion for the simple pendulum. We need to solve for $\ddot{\theta}$: $\ddot{\theta} = -\frac{mgL \sin\theta}{mL^2} = -\frac{g}{L} \sin\theta$

This is a second-order ordinary differential equation (ODE). For numerical simulation, we typically convert second-order ODEs into a system of first-order ODEs. Let $q_0 = \theta$ and $q_1 = \dot{\theta}$. Then: $\dot{q}_0 = q_1$ $\dot{q}_1 = -\frac{g}{L} \sin(q_0)$

Code Example: Symbolic Derivation with SymPy

Let’s use sympy to perform these symbolic derivations automatically.

import sympy as sp

# Define symbolic variables
t = sp.symbols('t')
m, L, g = sp.symbols('m L g', real=True, positive=True)

# Define generalized coordinate as a function of time
theta = sp.Function('theta')(t)

# Define generalized velocities
theta_dot = sp.diff(theta, t)

# Define Kinetic Energy (T)
T = sp.Rational(1, 2) * m * L**2 * theta_dot**2

# Define Potential Energy (V)
V = -m * g * L * sp.cos(theta)

# Define Lagrangian (L)
L_expr = T - V

# Euler-Lagrange equation: d/dt (dL/d(theta_dot)) - dL/d(theta) = 0
dL_d_theta_dot = sp.diff(L_expr, theta_dot)
dt_dL_d_theta_dot = sp.diff(dL_d_theta_dot, t)
dL_d_theta = sp.diff(L_expr, theta)

# Equation of motion
eom = dt_dL_d_theta_dot - dL_d_theta

# Solve for the second derivative (acceleration)
theta_ddot = sp.diff(theta_dot, t)
solved_eom = sp.solve(eom, theta_ddot)

print("Lagrangian (L):")
print(L_expr)
print("\nEquation of Motion (before solving for theta_ddot):")
print(eom)
print("\nSolved for theta_ddot (angular acceleration):")
# solved_eom is a list, take the first element
print(solved_eom[0])

Lagrangian (L):
L*g*m*cos(theta(t)) + L**2*m*Derivative(theta(t), t)**2/2

Equation of Motion (before solving for theta_ddot):
L**2*m*Derivative(theta(t), (t, 2)) + L*g*m*sin(theta(t))

Solved for theta_ddot (angular acceleration):
-g*sin(theta(t))/L

The sympy output L**2*m*Derivative(theta(t), (t, 2)) + L*g*m*sin(theta(t)) corresponds to $mL^2\ddot{\theta} + mgL \sin\theta = 0$, and sp.solve correctly gives us -g*sin(theta(t))/L which is $\ddot{\theta} = -\frac{g}{L} \sin\theta$. This is fantastic! SymPy handles the tedious calculus.

Code Example: Numerical Simulation with SciPy

Now that we have the equation of motion ($\ddot{\theta} = -\frac{g}{L} \sin\theta$), we can numerically integrate it over time. We’ll set up the system of first-order ODEs and use scipy.integrate.odeint.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# --- Physics Constants ---
g = 9.81  # m/s^2 (acceleration due to gravity)
L = 1.0   # m (length of the pendulum rod)
m = 1.0   # kg (mass of the pendulum bob - not used in EOM, but good to define)

# --- Initial Conditions ---
# [theta, theta_dot]
initial_state = [np.radians(30), 0.0] # 30 degrees, no initial angular velocity

# --- Time Steps ---
t_start = 0.0
t_end = 10.0
num_points = 500
t_span = np.linspace(t_start, t_end, num_points)

# --- Define the system of ODEs ---
# state = [theta, theta_dot]
# dstate/dt = [theta_dot, theta_ddot]
def pendulum_eom(state, t, L_val, g_val):
    theta, theta_dot = state
    theta_ddot = - (g_val / L_val) * np.sin(theta)
    return [theta_dot, theta_ddot]

# --- Solve the ODE ---
# The odeint function takes:
# 1. The function defining the derivatives (pendulum_eom)
# 2. The initial state
# 3. The time points to evaluate
# 4. Optional arguments passed to the derivative function
solution = odeint(pendulum_eom, initial_state, t_span, args=(L, g))

# Extract results
theta_rad = solution[:, 0]
theta_deg = np.degrees(theta_rad)
theta_dot_rad_per_sec = solution[:, 1]

# --- Plotting (Conceptual Output for Blog) ---
print(f"Simulation of Simple Pendulum (L={L}m, g={g}m/s^2)")
print(f"Initial Angle: {initial_state[0]:.2f} rad ({np.degrees(initial_state[0]):.2f} deg)")
print(f"Initial Angular Velocity: {initial_state[1]:.2f} rad/s")
print(f"Simulation Duration: {t_end} seconds")
print(f"\nFirst 5 data points (Time, Angle [deg], Angular Velocity [rad/s]):")
for i in range(5):
    print(f"t={t_span[i]:.2f}, theta={theta_deg[i]:.2f}, theta_dot={theta_dot_rad_per_sec[i]:.2f}")
print("...")
print(f"\nLast 5 data points (Time, Angle [deg], Angular Velocity [rad/s]):")
for i in range(num_points - 5, num_points):
    print(f"t={t_span[i]:.2f}, theta={theta_deg[i]:.2f}, theta_dot={theta_dot_rad_per_sec[i]:.2f}")

# Note: In a real environment, you'd show a plot here.
# plt.figure(figsize=(10, 6))
# plt.plot(t_span, theta_deg, label='Angle (degrees)')
# plt.plot(t_span, theta_dot_rad_per_sec, label='Angular Velocity (rad/s)')
# plt.title('Simple Pendulum Simulation')
# plt.xlabel('Time (s)')
# plt.ylabel('Value')
# plt.grid(True)
# plt.legend()
# plt.show()

Simulation of Simple Pendulum (L=1.0m, g=9.81m/s^2)
Initial Angle: 0.52 rad (30.00 deg)
Initial Angular Velocity: 0.00 rad/s
Simulation Duration: 10.0 seconds

First 5 data points (Time, Angle [deg], Angular Velocity [rad/s]):
t=0.00, theta=30.00, theta_dot=0.00
t=0.02, theta=29.99, theta_dot=-0.04
t=0.04, theta=29.96, theta_dot=-0.07
t=0.06, theta=29.90, theta_dot=-0.11
t=0.08, theta=29.82, theta_dot=-0.14
...

Last 5 data points (Time, Angle [deg], Angular Velocity [rad/s]):
t=9.92, theta=27.95, theta_dot=0.82
t=9.94, theta=28.32, theta_dot=0.74
t=9.96, theta=28.66, theta_dot=0.66
t=9.98, theta=28.97, theta_dot=0.58
t=10.00, theta=29.25, theta_dot=0.50

The output shows the angle oscillating, as expected for a pendulum. scipy.integrate.odeint is a robust solver for this type of problem. For a real-time game engine, you’d likely implement a simpler fixed-step integrator like Runge-Kutta 4 (RK4) or Verlet.

More Complex Systems: Double Pendulum

The real power of Lagrangian mechanics shines with systems like the double pendulum. This system is notoriously chaotic and difficult to derive using Newtonian methods because of the complex interaction forces. With Lagrangian mechanics, you define two generalized coordinates (e.g., two angles, $\theta_1, \theta_2$) and apply the Euler-Lagrange equation for each.

The process remains the same:

1. Define the Cartesian coordinates of each mass in terms of $\theta_1, \theta_2$.
2. Calculate the kinetic energy for each mass, then sum them for total $T$.
3. Calculate the potential energy for each mass, then sum them for total $V$.
4. Form the Lagrangian $L = T - V$.
5. Apply the Euler-Lagrange equation for $\theta_1$ and then for $\theta_2$.
6. You’ll get two coupled second-order ODEs. Solve them simultaneously using sympy and then numerically.

The derivation becomes algebraically messy by hand, but SymPy handles it like a champ. This is a critical workflow for using Lagrangian mechanics in game development: derive symbolically, implement numerically.

Constraints and Lagrange Multipliers (Briefly)

What if your constraints aren’t naturally handled by generalized coordinates? For example, a non-holonomic constraint (like a wheel rolling without slipping) or an inequality constraint (an object resting on a surface). In these cases, you might introduce Lagrange Multipliers ($\lambda_k$) into the Euler-Lagrange equations.

The modified Euler-Lagrange equation becomes:

$\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{q}_i} \right) - \frac{\partial L}{\partial q_i} + \sum_k \lambda_k \frac{\partial f_k}{\partial q_i} = 0$

Where $f_k(q_1, …, q_n, t) = 0$ are the constraint equations. The $\lambda_k$ values represent the forces required to enforce these constraints.

Note: While powerful, implementing Lagrange multipliers in a real-time game engine can be computationally expensive and complex due to the need to solve a coupled system of differential-algebraic equations (DAEs). Most commercial game engines use iterative force-based constraint solvers (like sequential impulse or iterative projection) for rigid body dynamics, which are often more robust and performant for collision-rich environments. Lagrangian mechanics shines when constraints are easily expressed via generalized coordinates, avoiding the need for multipliers altogether.

Advantages of Lagrangian Mechanics for Games

• Natural Constraint Handling: As shown, choosing good generalized coordinates simplifies constraint management dramatically.
• Energy Conservation: For conservative forces, the derived equations inherently conserve energy, leading to more stable long-term simulations.
• Modularity: Once the Lagrangian is defined, the system’s equations of motion can be derived mechanically (even automatically with sympy), which helps manage complexity.
• Reduced DOF: By using generalized coordinates, you work with the minimum degrees of freedom, reducing the size of the system you need to solve.

Disadvantages / Challenges

• Derivation Complexity: While sympy helps, defining $T$ and $V$ for very complex multi-body systems (e.g., a human skeleton) can still be non-trivial.
• Numerical Stability: Integrating the derived ODEs still requires careful selection of numerical integrators and step sizes to avoid instability, especially for stiff systems or those with high frequencies.
• Collision Detection/Response: Lagrangian mechanics focuses on the dynamics (how things move). Collision detection (when objects touch) and response (how they bounce/slide) are entirely separate problems that need to be handled by other means (e.g., GJK/EPA for detection, impulse-based methods for response).
• Real-time Performance: The solutions derived are often non-linear ODEs, requiring iterative numerical solvers that can be computationally intensive for many degrees of freedom at high update rates. It’s often suitable for specific, high-fidelity parts of a game (e.g., character ragdolls, specific mechanical puzzles) rather than an entire rigid body engine for many interacting objects.

Practical Considerations for a Game Engine

1. Symbolic vs. Numerical:
   • Symbolic (Offline): Use sympy (or Mathematica/Maple) to derive the equations of motion for your specific mechanical systems (pendulums, wheeled vehicles, robot arms, etc.).
   • Numerical (Runtime): Take the derived equations (e.g., $\ddot{\theta} = -\frac{g}{L} \sin\theta$) and implement them in your game engine’s C++/C# code.
2. Integrators:
   • For game physics, simple integrators like Euler (explicit) are easy but prone to instability and energy loss.
   • Runge-Kutta 4 (RK4) is a good general-purpose choice, offering good accuracy and stability for a fixed step size.
   • Verlet integration is excellent for position-based dynamics and naturally conserves energy for conservative forces.
   • Implicit solvers (e.g., implicit Euler) offer better stability for stiff systems but are more complex to implement as they require solving a system of equations at each step.
3. State Management: Your simulation state will typically be a vector of generalized coordinates and their velocities (e.g., [theta, theta_dot, phi, phi_dot, ...]).
4. Collision Handling: This is the biggest gap. Lagrangian mechanics doesn’t inherently handle collisions. You’ll need:
   • Collision detection library (e.g., Bullet or Havok for complex shapes, or simpler AABB/sphere tests).
   • Collision response: Once a collision is detected, you can apply impulses (Newtonian style) or modify the generalized velocities to reflect the bounce/slide. This often requires temporarily switching away from the Lagrangian formulation or incorporating non-conservative forces into the Lagrangian framework.

Conclusion

Lagrangian mechanics offers an elegant and powerful analytical framework for understanding and simulating mechanical systems. For game developers, it’s particularly valuable for systems where constraints are easily absorbed into generalized coordinates, simplifying the derivation of complex equations of motion.

While not a complete replacement for traditional force-based rigid body engines (especially for broad-phase collision detection and contact resolution), Lagrangian mechanics excels in niche areas:

• High-fidelity simulations of specific mechanical puzzles or contraptions.
• Character ragdolls or inverse kinematics systems where articulated joints are essentially constraints.
• Educational tools and physics demonstrators where accuracy and conceptual clarity are paramount.

By combining the symbolic power of sympy with efficient numerical integrators, you can build surprisingly robust and stable physics simulations for your games, allowing you to focus on the fun parts of game development instead of wrestling with tricky constraint forces. Give it a try, and you might find it simplifies some of your toughest physics challenges.