Multibody Dynamics · Module 9 of 10
Numerical Simulation of Multibody Systems
The equations rarely solve in closed form, so we march them forward in time. Choosing an integrator and a time step, and watching the errors, turns a model into a running simulation.
Readiness check
Tick only what you can do closed-notes before starting.
- Recall the explicit Euler update for a first-order system.
- State a state vector for a second-order system.
- Recall that step size controls accuracy and stability.
- Recall the natural frequency and period from Module 7.
- State why a stiff system needs small steps or an implicit method.
The core idea
Simulation writes the equations of motion as a first-order state system and marches it forward in small time steps. The integrator and the step size set the accuracy, stability, and cost of the run.
state: y = (q, q̇)explicit Euler: y_{n+1} = y_n + h f(y_n)step guide: h ≤ T / (tens of steps)The equations of motion are second order in the coordinates, so we rewrite them as a first-order system in the state y = (q, q̇): the derivative of position is velocity, and the derivative of velocity is the acceleration the dynamics solve returns. A time integrator advances this state. Explicit methods like forward Euler or Runge-Kutta evaluate the right side at known states and are simple and cheap, but they go unstable if the step is too large for the fastest motion in the system. Implicit methods solve a small equation each step and stay stable for stiff systems, at higher cost per step. For constrained models the integrator must also respect the algebraic constraints: it either solves the augmented DAE each step, with stabilization to prevent drift, or works in a reduced set of independent coordinates. The step size is chosen from the fastest timescale, the highest natural frequency or stiffest spring: resolve a period with tens of steps for accuracy, and stay under the stability limit set by the method.
The skills, taught in order
Five skills turn equations of motion into a stable run.
9.1 State-space form
Rewrite the second-order equations as a first-order system in (q, q̇). This is the form every general integrator expects.
9.2 Explicit integration
Use forward Euler or Runge-Kutta to step from known states. Simple and fast, but the step is capped by stability for stiff or high-frequency systems.
9.3 Implicit integration
Solve a small equation each step for stability on stiff systems (stiff springs, contact, fast actuators). It costs more per step but permits far larger steps.
9.4 DAE integration
For constrained models, integrate the augmented system with constraint stabilization or project onto the constraint manifold each step, so geometry stays satisfied.
9.5 Choosing the time step
Set h from the fastest timescale: resolve the shortest period with tens of steps and stay under the method’s stability limit. Then check that halving h barely changes the answer.
| Integrator | Cost per step | Best for |
|---|---|---|
| Explicit (Euler, RK) | low | smooth, non-stiff motion |
| Implicit | higher | stiff systems, contact |
| DAE solver | higher | constrained multibody models |
Match the integrator to the physics; a stiff contact model and a gently swinging linkage want different tools.
Engineering connection: a vehicle-dynamics simulation with stiff tire and bushing models uses an implicit or specialized solver, because an explicit step small enough to stay stable would make a real-time run impossible.
Worked example 1: one explicit Euler step
The linearized pendulum obeys θ̈ = −ω_n² θ with ω_n² = 19.62 (from L = 0.5 m). Starting from θ = 0.2 rad, θ̇ = 0, take one forward-Euler step of size h = 0.01 s.
- ProblemAdvance the pendulum state one Euler step as in Figure 1.
- Given / findω_n² = 19.62 s−2, state (θ, θ̇) = (0.2, 0), h = 0.01 s. Find the new state.
- AssumptionsSmall-angle linear model; forward (explicit) Euler.
- ModelState y = (θ, θ̇); f = (θ̇, −ω_n²θ). Update y_{n+1} = y_n + h f(y_n).
- Equationsθ̇_{new} = θ̇ + h(−ω_n² θ)θ_{new} = θ + h θ̇
- Solveθ̇_{new} = 0 + 0.01(−19.62 × 0.2) = 0.01(−3.924) = −0.0392 rad/s. θ_{new} = 0.2 + 0.01(0) = 0.2 rad.
- CheckThe angular velocity turns negative because the restoring acceleration points back toward zero; the angle barely moves in one small step, as expected.
- ConclusionAfter one step the state is (0.200, −0.0392). This is the raw explicit update; because forward Euler slowly adds energy to an undamped swing, a stable integrator such as semi-implicit Euler or Runge-Kutta is used to march the full motion.
Worked example 2: choosing a stable time step
The same pendulum has period T = 1.42 s. Choose a time step that resolves the motion with about fifty steps per period, and comment on stability.
- ProblemPick a time step for the pendulum in Figure 2.
- Given / findPeriod T = 1.42 s, target about 50 steps per period. Find h and check it is safe.
- AssumptionsThe fastest motion is this single oscillation; accuracy, not just stability, sets the step.
- ModelResolution rule: h = T / N with N steps per period. A conditionally stable explicit integrator (central-difference, semi-implicit Euler, or Runge-Kutta) stays stable for h of order 2/ω_n; plain forward Euler is unstable on an undamped oscillation at any step, so a stable scheme is used for the run.
- Equationsh = T / N = 1.42 / 50stable explicit scheme: h < 2/ω_n = 0.45 s
- Solveh = 1.42 / 50 = 0.0284 s, far below the 0.45 s stability ceiling.
- CheckAt 0.0284 s the step is set by accuracy, not stability, with comfortable margin, and simulating 5 s takes about 176 steps.
- ConclusionA step near 0.028 s resolves the swing accurately and safely. If a stiff spring were added, its higher frequency, not this one, would then dictate the step.
Misconceptions and diagnostics
| Mistake | Symptom | Diagnostic question | Correction |
|---|---|---|---|
| Sizing the step to the slow motion | Explicit run diverges on a stiff element | What is the fastest timescale here? | Set h from the highest frequency, not the visible motion. |
| Never checking step convergence | Trusting a result at one step size | Does halving h change the answer? | Refine the step until the solution stops moving. |
| Ignoring constraint drift in integration | Geometry slowly falls apart over a long run | Is Φ staying near zero? | Use a DAE solver with stabilization or projection. |
| No energy monitoring | Numerical damping or gain goes unnoticed | Is energy behaving as the physics says? | Track energy as a running validation signal. |
Practice ladder
For θ̈ = −19.62 θ, state (0.2, 0), take a step h = 0.02 s. Find the new angular velocity.
Show answer
θ̇_{new} = 0 + 0.02(−19.62 × 0.2) = −0.0785 rad/s.
A system’s fastest natural frequency is 100 rad/s. Estimate an explicit stability step ceiling.
Show answer
h < 2/ω = 2/100 = 0.02 s; use a fraction of this for accuracy.
A period is 0.5 s and you want 100 steps per period. Find the time step and the number of steps for a 2 s run.
Show answer
h = 0.5/100 = 0.005 s; a 2 s run is 2/0.005 = 400 steps.
For a model you know, state its fastest timescale, choose a time step, and describe how you would confirm the step is small enough.
What good work looks like
A good answer identifies the stiffest or highest-frequency element, sets h from it, and proposes a step-halving convergence and an energy or constraint-error check.
Working with AI, and proving it yourself
Ask AI to check your step choice, not to pick blindly
Portfolio task
Integrate one small model by hand for a few steps, then state the integrator, step size, and the convergence and energy checks you would run.
Retrieval and spaced review
Closed notes. Answer out loud, then reveal.
1. What is the state of a mechanical system?
The positions and velocities, y = (q, q̇).
2. Write the explicit Euler update.
y_{n+1} = y_n + h f(y_n).
3. When is an implicit method worth its cost?
For stiff systems where explicit steps must be tiny to stay stable.
4. What must a DAE integrator also do?
Keep the algebraic constraints satisfied, with stabilization or projection.
5. How is the time step chosen?
From the fastest timescale, with tens of steps per period and a stability margin.
Textbook mapping
Numerical treatment of the general formalism is developed by Wittenburg and computational-dynamics texts. Use these references to read further.
| Topic in this module | Where to read more |
|---|---|
| Equations of motion for numerical investigation | Wittenburg, Dynamics of Multibody Systems, ch. 5 |
| Integrators for multibody DAEs | Nikravesh, Computer-Aided Analysis of Mechanical Systems |
| ODE and stability basics | Chapra and Canale, Numerical Methods for Engineers |
Chapter references are to Wittenburg, Dynamics of Multibody Systems (Springer); integration methods draw on standard numerical-methods texts.