M2-T5 — Simulation You Can Trust

Learning objectives

After this lesson, you should be able to:

• Distinguish verification (solving the equations right) from validation (having the right equations).
• Measure the order of accuracy of an integrator from a convergence study.
• Compute the stability limit of the step for an explicit integrator and choose a step for a given fast dynamics.
• Test a simulator with conserved quantities, known solutions and quaternion normalization.

The intuition

A number from a simulator can be wrong for three reasons, and each needs a different test.

1. The model is wrong: the equations leave out something that matters, or the parameters are off. Finding this is validation, comparing the model with the real vehicle. The last four lessons built the model and said what was left out; flight tests and system identification validate it.
2. The equations are solved badly: the integrator is too crude, the step too large. This is a numerical error, and finding it is verification. It is the subject of this lesson.
3. The code does not implement the equations: a wrong sign, a degrees-radians mistake, a frame mix-up. Verification catches this too, because the checks below test the whole chain from equation to output.

The remarkable thing about verification is that it needs no data from the real vehicle. The physics supplies the tests: some quantities must be conserved, some solutions are known exactly, and a correct integrator has an order of accuracy that can be measured.

Five tests, cheapest first

• Convergence. Run the same problem with the step halved, again and again. The error should fall by a fixed factor, $2^p$ for an integrator of order $p$. If it does not, the integrator or the code is wrong, and the observed order tells which.
• Stability of the step. An explicit integrator blows up if the step is larger than a limit set by the fastest dynamics in the model. For Hopper that is the gimbal servo at 40 rad/s.
• Conserved quantities. With no thrust and no torque, the energy and the angular momentum of a rigid body must not change. Any drift is the integrator's.
• Known solutions. A thrust-free fall is a parabola, a harmonic oscillator is a cosine. Compare with the formula.
• Invariants of the representation. A quaternion must have unit length, a rotation matrix must be orthogonal. They drift if the update does not respect them.

Accuracy and stability are different limits

A step can be too large in two ways. Accuracy says how close the answer is to the truth, and improves smoothly as the step shrinks, at a rate set by the order. Stability says whether the computation stays bounded at all, and is a cliff: below a limit it is fine, above it the numbers grow without bound within a few steps. The limit depends on the fastest decaying mode in the model. Hopper's gimbal servo has a time constant of 25 ms, a rate of 40 rad/s, so it sets the limit even if nothing else in the model moves quickly.

Test a simulator the way you test a controller: against things whose answer you know. If you cannot name the test that would catch a particular bug, that bug is waiting.

Definitions

Verification and validation

Verification: does the code solve the stated equations correctly? Validation: do the equations describe the real vehicle? Only verification can be done without hardware.

Step $h$ and integrator

The fixed time increment, in seconds, and the rule that advances the state by one step. Explicit Euler: $\mathbf{x}_{k+1} = \mathbf{x}_k + h\,\mathbf{f}(\mathbf{x}_k)$. RK4 combines four evaluations of $\mathbf{f}$ per step.

Order of accuracy $p$

The global error shrinks as $h^p$: Euler has $p = 1$, RK4 has $p = 4$.

Stability limit

The largest value of $|\lambda|\,h$ (the fastest decay rate times the step) for which an explicit method does not grow a decaying solution.

Invariant

A quantity that the exact solution leaves unchanged under the stated conditions: energy, angular momentum, the norm of a quaternion.

A convergence study fits a straight line to the logarithm of the error against the logarithm of the step. The slope is the observed order:

For the test equation $\dot{x} = \lambda x$ with $\lambda < 0$, one step of an explicit method multiplies the state by a factor $R(z)$, with $z = \lambda h$. The method is stable when $|R(z)| \le 1$:

A torque-free rigid body must keep its kinetic energy and the magnitude of its angular momentum:

Derivation: where the stability limits and the energy drift of Euler come from

The limits that the worked example finds by experiment each have a short derivation.

The step limit of explicit Euler and RK4

1. $$x_{k+1} = x_k + h\lambda x_k = (1 + \lambda h)\,x_k = R(z)\,x_k$$

   Apply explicit Euler to $\dot{x} = \lambda x$. Each step multiplies the state by the same factor.

2. $$|1 + z| \le 1 \;\Longleftrightarrow\; z \ge -2 \;\Longleftrightarrow\; h \le \frac{2}{|\lambda|}$$

   The state stays bounded only if the factor has magnitude at most 1. For a real negative $z$ that means $z \ge -2$. For the gimbal servo, $|\lambda| = 40$ gives $h \le 0.05$ s.

3. $$|R_{\text{RK4}}(-z)| \le 1 \;\Longleftrightarrow\; z \le 2.785$$

   Do the same with the RK4 polynomial. There is no neat closed form: the limit is the root of the polynomial equation, found numerically (the program below finds it by bisection). For the servo $h \le 0.0696$ s.

Staying below the stability limit is not enough for accuracy. At $z = -0.8$ (a 20 ms step for the servo), Euler's factor is $1 - 0.8 = 0.20$ while the exact decay over the step is $e^{-0.8} = 0.45$. The simulation is stable and decays more than twice as fast as the servo does.

Why explicit Euler gains energy

1. $$\begin{pmatrix} x \\ v/\omega \end{pmatrix}_{k+1} = \begin{pmatrix} 1 & h\omega \\ -h\omega & 1 \end{pmatrix}\begin{pmatrix} x \\ v/\omega \end{pmatrix}_{k}$$

   Apply Euler to the oscillator $\ddot{x} = -\omega^2 x$ in the scaled variables $(x, v/\omega)$. The exact solution rotates this vector without changing its length.

2. $$\left\lVert \mathbf{s}_{k+1} \right\rVert^2 = \bigl(1 + (h\omega)^2\bigr)\,\lVert\mathbf{s}_k\rVert^2$$

   The matrix is a rotation scaled by $\sqrt{1 + (h\omega)^2}$. The squared length, which is proportional to the energy, is multiplied by $1 + (h\omega)^2$ on every step.

3. $$\frac{E(T)}{E(0)} = \bigl(1 + (h\omega)^2\bigr)^{T/h} \approx e^{\omega^2 h\,T}$$

   After $T/h$ steps. For small $h$ the exponent is proportional to $h$: the energy drift of Euler is first order in the step, so it is never zero and grows with the length of the run.

Euler does not merely add a small error: it adds one with a sign. Every step pumps in energy, so a long simulation of an oscillating or rotating system is wrong in a systematic direction.

A worked example

Run the five tests on the tools of this course.

1. Convergence. Integrate $\dot{x} = -x$ to $t = 2$ with steps 0.2, 0.1, 0.05 and 0.025 and compare with $e^{-2}$. Euler errors fall from 2.8×10⁻² to 3.4×10⁻³, an observed order of 1.01. RK4 errors fall from 4.3×10⁻⁶ to 9.0×10⁻¹⁰, an order of 4.07. The library implements what it claims.
2. Stability. Integrate $\dot{x} = -40x$ (the servo). At $h = 0.04$ both methods decay. At 0.06 Euler blows up and RK4 decays. At 0.07 both blow up. The limits of $|\lambda| h$, found by bisection, are 2.0000 for Euler and 2.7853 for RK4.
3. Conserved quantities. Integrate a torque-free body with inertias (1, 2, 3) kg m² and rates (0.1, 2.0, 0.1) rad/s for 20 s. Explicit Euler with a step of 0.05 s changes the energy by 46% and the angular momentum by 21%. RK4 with the same step changes them by 7×10⁻⁸ and 4×10⁻⁸, and with a step of 1 ms by 10⁻¹⁴.
4. Quaternions. Rotate for 100 s at 1°/s about all three axes with a step of 10 ms. The exponential update keeps the norm at 1 to twelve places. The naive first-order update, $\mathbf{q} \leftarrow \mathbf{q} + \tfrac12 h\,\mathbf{q}\otimes(0, \boldsymbol{\omega})$, ends at 1.000115, and the error keeps growing with time.
5. Known solution. Hopper in a thrust-free fall from 100 m with a 3 m/s sideways speed. After 2 s the height is 80.3800 m, the exact value, and the energy changes by 10⁻¹⁴ (relative). Then the oscillator: explicit Euler multiplies the energy by 1.4889 over two seconds at $\omega = 2$ and $h = 0.05$, exactly the predicted factor $(1 + (h\omega)^2)^{T/h}$.
6. Choose the step. The fastest dynamics in Hopper is the servo at 40 rad/s. RK4 is stable up to 70 ms, but its error per step relative to the exact decay is 0.5% at 20 ms, 1×10⁻⁴ at 10 ms, 3×10⁻⁶ at 5 ms and 9×10⁻¹⁰ at 1 ms. A step of 1 ms leaves a wide margin, costs nothing noticeable, and is what the rest of the course uses.

import math
from gnckit import (integrate, euler_step, rk4_step, convergence_order, rigid_body_rates, rigid_body_invariants,
                    quat_integrate, quat_integrate_naive, quat_norm, HOPPER, hopper_deriv, hopper_energy)

# 1. order of accuracy: x' = -x to t = 2
dts = [0.2, 0.1, 0.05, 0.025]
for step in (euler_step, rk4_step):
    errs = [abs(integrate(lambda t, x: [-x[0]], [1.0], 2.0, dt, step)[1][-1][0] - math.exp(-2.0)) for dt in dts]
    print(f'{step.__name__:10s} errors {["%.1e" % e for e in errs]}, observed order {convergence_order(errs, dts):.2f}')

# 2. stability limit for x' = -40 x (the gimbal servo)
for dt in (0.04, 0.06, 0.07):
    out = []
    for step in (euler_step, rk4_step):
        x = integrate(lambda t, x: [-40.0 * x[0]], [1.0], 3.0, dt, step)[1][-1][0]
        out.append(f'{step.__name__} {"blows up" if abs(x) > 1 else "decays"}')
    print(f'dt = {dt:.2f} s (40 dt = {40 * dt:.1f}): ' + ', '.join(out))

def limit(R):                                         # largest lam*dt for which the one-step factor R(-lam*dt) has magnitude at most 1
    lo, hi = 0.0, 10.0
    for _ in range(60):
        mid = (lo + hi) / 2
        lo, hi = (mid, hi) if abs(R(-mid)) <= 1 else (lo, mid)
    return lo
print(f'stability limit of lam*dt: euler {limit(lambda z: 1 + z):.4f}, rk4 {limit(lambda z: 1 + z + z * z / 2 + z ** 3 / 6 + z ** 4 / 24):.4f}')

# 3. torque-free rigid body: energy and momentum must not change
J = (1.0, 2.0, 3.0)
w0 = [0.1, 2.0, 0.1]
e0, h0 = rigid_body_invariants(J, w0)
for label, step, dt in (('euler 0.05 s', euler_step, 0.05), ('rk4 0.05 s', rk4_step, 0.05), ('rk4 0.001 s', rk4_step, 0.001)):
    w = integrate(lambda t, x: rigid_body_rates(J, x), w0, 20.0, dt, step)[1][-1]
    e1, h1 = rigid_body_invariants(J, w)
    print(f'{label:13s} energy drift {abs(e1 - e0) / e0:9.2e}, momentum drift {abs(h1 - h0) / h0:9.2e}')

# 4. quaternion norm after 100 s at 1 deg/s about each axis, step 0.01 s
w = [0.0175, 0.0175, 0.0175]
qa = qb = [1.0, 0.0, 0.0, 0.0]
for _ in range(10000):
    qa = quat_integrate(qa, w, 0.01)
    qb = quat_integrate_naive(qb, w, 0.01)
print(f'quaternion norm: exact update {quat_norm(qa):.12f}, naive update {quat_norm(qb):.6f}')

# 5. Hopper in free fall: energy
s = [0.0, 100.0, 0.0, 3.0, 0.0, 0.0, 100.0]
E0 = hopper_energy(s)
xs = integrate(lambda t, x: hopper_deriv(x, 0.0, 0.0), s, 2.0, 0.001)[1]
print(f'thrust-free fall: energy changes by {abs(hopper_energy(xs[-1]) - E0) / E0:.1e}, height {xs[-1][1]:.4f} m (exact {100 - 0.5 * 9.81 * 4:.4f} m)')

# 6. explicit Euler on an oscillator: the energy is multiplied by (1 + (w dt)^2) at every step
w, dt, T = 2.0, 0.05, 2.0
xs = integrate(lambda t, x: [x[1], -w * w * x[0]], [1.0, 0.0], T, dt, euler_step)[1]
E = lambda s: 0.5 * s[1] ** 2 + 0.5 * w * w * s[0] ** 2
print(f'oscillator, explicit Euler: energy x{E(xs[-1]) / E(xs[0]):.4f}, predicted (1 + (w dt)^2)^(T/dt) = x{(1 + (w * dt) ** 2) ** (T / dt):.4f}')

euler_step errors ['2.8e-02', '1.4e-02', '6.8e-03', '3.4e-03'], observed order 1.01
rk4_step   errors ['4.3e-06', '2.5e-07', '1.5e-08', '9.0e-10'], observed order 4.07
dt = 0.04 s (40 dt = 1.6): euler_step decays, rk4_step decays
dt = 0.06 s (40 dt = 2.4): euler_step blows up, rk4_step decays
dt = 0.07 s (40 dt = 2.8): euler_step blows up, rk4_step blows up
stability limit of lam*dt: euler 2.0000, rk4 2.7853
euler 0.05 s  energy drift  4.63e-01, momentum drift  2.09e-01
rk4 0.05 s    energy drift  7.28e-08, momentum drift  3.63e-08
rk4 0.001 s   energy drift  1.33e-14, momentum drift  5.76e-15
quaternion norm: exact update 1.000000000000, naive update 1.000115
thrust-free fall: energy changes by 9.7e-15, height 80.3800 m (exact 80.3800 m)
oscillator, explicit Euler: energy x1.4889, predicted (1 + (w dt)^2)^(T/dt) = x1.4889

In robotics

These tests are standard practice, and the stakes grow with the programme.

• Monte Carlo campaigns. Module 12 flies thousands of dispersed simulated flights to show that a requirement is met. A simulator that has not been verified turns that statistical evidence into a statistical fiction, and a biased integrator biases every run the same way, which is far worse than noise because averaging does not remove it.
• Fixed-step, real-time simulation. Software-in-the-loop and hardware-in-the-loop setups (see UAV M6-T14: Software-in-the-Loop simulation) run at a fixed step tied to the controller's clock. The step cannot be chosen for accuracy alone: it is the controller period, and the model must be integrated with sub-steps if the controller's period is too long for the vehicle's fast dynamics.
• Orbit propagation. Over thousands of orbits a small energy drift becomes a wrong orbit. Propagators therefore use high-order variable-step integrators or symplectic integrators, which conserve energy by construction. The same discipline applies to any long, oscillatory simulation.
• Estimators. An inertial navigation filter integrates gyroscope and accelerometer data at a few hundred hertz. The attitude update there is the exponential update of this lesson, and module 6 adds corrections for coning and sculling, the errors of integrating rotation in steps.
• Credibility of models. Space programmes judge a simulation by documented verification and validation evidence, as set out in NASA's standard for models and simulations (NASA-STD-7009, listed in the resources). The five tests above are the verification half of that evidence.

The verification and validation process as a whole is the subject of UAV M8-T3: Verification and validation. The labs of this course use the same discipline in miniature: lab M2-L2 hands you a simulator with three defects and asks you to find them with these tests.

Common misconceptions

Key takeaways

1. Verification (solving the equations right) needs no flight data; validation (the right equations) does. This lesson is verification.
2. The observed order of accuracy, the slope of log error against log step, tells whether an integrator does what it should: 1 for Euler, 4 for RK4.
3. An explicit method is stable only if $|\lambda| h$ is below a limit: 2 for Euler, 2.785 for RK4. Hopper's servo at 40 rad/s requires $h < 0.05$ s for Euler and 0.07 s for RK4, and accuracy needs far less.
4. Explicit Euler multiplies the energy of an oscillator by $1 + (h\omega)^2$ at every step: a systematic gain.
5. Conserved quantities (energy, angular momentum), exact solutions and unit norms are free tests of a simulator.
6. A simulator is only as good as the tests it has passed; keep them and run them every time it changes.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• NASA, NASA-STD-7009: Standard for Models and Simulations . What must be shown before a simulation supports a decision. Free to read.
• R. Tedrake, Underactuated Robotics , the appendix on simulation and numerical integration, with notebooks. Free to read; no licence statement was found.
• B. L. Stevens, F. L. Lewis and E. N. Johnson, Aircraft Control and Simulation, 3rd edition, Wiley, 2016, the sections on numerical integration of the equations of motion. Paid book.
• E. Hairer, S. P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd edition, Springer, 1993, the chapters on Runge-Kutta methods and on stability. The standard reference. Paid book.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

• Before
  M2-T4: Trim and Linearization

  Find an equilibrium of a nonlinear model, linearize about it to a state-space model, and state the range over which the linear model is trusted.

• Lab
  M2-L2: Find the Simulator Bug

  A simulator gains energy with nothing driving it. Use energy, momentum and step-size checks to locate the three defects and correct them.

Other courses on this site treat neighbouring topics from their own angle:

• UAV M1-T4: UAV coordinate frames
• UAV M1-T5: Six-degree-of-freedom UAV dynamics
• UAV M5-T1: UAV equations of motion
• UAV M5-T2: Quadrotor dynamic modelling
• Linear Algebra M7-T3: Three-Dimensional Rotation Matrices
• Linear Algebra M7-T8: Euler Angles and Quaternions
• Linear Algebra M8-T4: Jacobian Matrices
• Navigation M1-T3: Robot Motion Models for Navigation
• Robot Learning M7-T2: System Identification

This course leaves the following to the pages that already teach them:

• Rotation matrices, Euler angles and quaternions from first principles
• The quadrotor model, motor mixing and thrust generation
• Unicycle and car-like kinematic models