↳ FRAME 05 · FIELD NOTES · MODULE 07

Numerical integration

A simulator does not solve the equations of motion; it steps them, and every step makes a small error. Which way those errors lean decides whether a swinging pendulum steadily gains energy, slowly loses it, or keeps it within a fixed band for ever. This module shows the three behaviours on the running example, and how to tell them apart before trusting a result.

11 formulas · 37/37 checks pass · 19 sources · in review, not yet fact-checked

Key ideas

  1. Two lines, in two orders

    Explicit Euler moves the position with the old velocity (F7.01); semi-implicit Euler moves it with the new one (F7.02). The code differs only in the order of two lines, and the behaviour differs completely. MuJoCo's integrator called Euler is the semi-implicit one, and it also treats joint damping implicitly.[1][2]

  2. Explicit Euler pumps energy in

    For a linear oscillator, each explicit Euler step multiplies the energy by 1+(ωh)21 + (\omega h)^2 (F7.09), whatever the step size: a smaller step slows the growth but never stops it. For the running example at 2 ms that factor is 1 + 2.37×10⁻⁴ per step, which compounds to ×10.7 over 10000 steps; the pendulum itself, slightly softer, gains ×10.5 from 5°.[3][4]

  3. Semi-implicit Euler keeps a nearby energy

    For a linear oscillator it conserves a slightly different quantity exactly (F7.06); for a single hinge such as the pendulum, closely enough that the true energy wobbles around a fixed value, by about ±ωh/2\pm\omega h/2, twice per swing, with no trend (F7.07). What remains is a frequency very slightly too high: 9.9 parts per million for the linearised example (F7.08).[2][5]

  4. RK4 is far more accurate, and slowly loses energy

    Classical Runge–Kutta evaluates the forces four times per step (F7.04) and is fourth-order accurate (F7.05). It conserves no energy-like quantity: for a linear oscillator each step removes a fraction of about (ωh)6/72(\omega h)^6/72 of the energy (F7.10), tiny at this step size, and a trend rather than a band.[2][6][7]

  5. Order says how fast the error shrinks, not which way it leans

    An order-pp method's error at a fixed end time falls by about 2p2^p when the step is halved (F7.05). The measured orders on the running example, from the angle error after 2 s of swing from 30° against a fine RK4 reference, are 1.01 for Euler (2 ms to 1 ms) and 4.02 for RK4 (4 ms to 2 ms). Halving the step is the cheapest test of whether a result depends on it.[7][8][9]

  6. Implicit variants are for stiff velocity-dependent forces

    implicit and implicitfast take each step with a linearisation of the velocity-dependent forces treated implicitly (F7.03); constraint forces are left out of it, and implicitfast also drops the Coriolis and centripetal terms. Without such forces both reduce to semi-implicit Euler, as the checks behind F7.03 confirm on the running example. MuJoCo recommends implicitfast for most models: similar cost to Euler, more stable. A fifth integrator, discrete, is meant for stiff springs, position servos or strong damping interacting with contacts, and is not covered here.[1]

  7. A single hinge is the easy case

    With one hinge the mass matrix is a constant, and semi-implicit Euler behaves as above. In a multi-link arm the mass matrix changes with the pose, MuJoCo's Euler is no longer symplectic, and energy can drift with no bug anywhere. A double pendulum released straight (elbow at 0°) from 60° gains 45 % of one link's mgLmgL in 20 s at 2 ms, 18 % at 1 ms and 7.9 % at 0.5 ms; under RK4, from the same release, it changes by −1×10⁻⁸ of mgLmgL. Released with the elbow bent at 90°, Euler's energy falls instead: it changes by −4.4 % of mgLmgL in 2 s at 2 ms and −2.4 % at 1 ms. Over 20 s the same release gives −8.2 % at 2 ms and −12 % at 1 ms: the runs at the two steps have followed different paths, and one halving no longer shrinks the error. MuJoCo's own documentation describes a pendulum mechanism that diverges under Euler or implicitfast and loses energy under implicit.[1][10]

  8. Even the stable methods have a step limit

    Semi-implicit Euler is stable only while ωh<2\omega h < 2 (F7.11); for the pendulum that allows steps up to 260 ms, far above 2 ms. A damping term applied explicitly has its own limit, hb/I<2h b/I < 2, and slightly less, 2−(ωh)2/22 - (\omega h)^2/2, once stiffness is included. That is the one a controller meets: a derivative gain applied as a torque by your own code must stay below about 2I/h2I/h, here 5.82 N·m·s/rad.[5]

Formulas

The same entries as the formula sheet. Open an entry's details for its symbols and the checks behind it.

  1. Explicit (forward) Euler

    Not offered by MuJoCo; included because it is the textbook baseline.

    vt+h=vt+h at,qt+h=qt+h vtv_{t+h} = v_t + h\,a_t,\qquad q_{t+h} = q_t + h\,v_t
    F7.01✓ 2/2 checked[1]
    Symbols and checks
    SymbolMeaningUnit
    hhtimesteps
    aaaccelerationrad/s²
    • ✓ the page's two update loops compile (-Wall -Wextra -Wconversion -Werror)expected 0 · got 0 · absolute error 0
    • ✓ C++ explicit loop: energy ratio after 20 s, vs the Python runexpected 10.5334 · got 10.5334 · relative error 0

  2. Semi-implicit Euler, MuJoCo's default

    Position is updated with the new velocity. Symplectic when the mass matrix is constant (a single hinge, as in the running example); in multi-link models the mass matrix depends on the pose, MuJoCo's Euler is not symplectic, and the energy can drift either way. MuJoCo also treats joint damping implicitly unless the eulerdamp flag is disabled.

    vt+h=vt+h at,qt+h=qt+h vt+hv_{t+h} = v_t + h\,a_t,\qquad q_{t+h} = q_t + h\,v_{t+h}
    F7.02✓ 8/8 checked[1][2][10][11]
    Symbols and checks
    • ✓ modified energy drift over 1e5 steps (relative)expected 0 · got 4.20402×10⁻¹⁴ · absolute error 4.2×10⁻¹⁴
    • ✓ straight from 60 deg, Euler at 2 ms: the 2 s window means rise steadily (a drift, not a band)observed as stated9 of 9 successive windows higher; first +0.014, last +0.409 m g L
    • ✓ straight from 60 deg: halving the step at least halves the gain in 20 s, twice overobserved as stated+0.450, +0.177, +0.079 m g L at 2, 1, 0.5 ms
    • ✓ elbow bent 90 deg, 2 s: Euler loses energy, and each halved step about halves the lossobserved as stated-0.0435, -0.0238, -0.0125 m g L at 2, 1, 0.5 ms
    • ✓ straight from 60 deg under RK4 at 2 ms: change below 1e-6 m g Lobserved as statedchange -1.2e-08 m g L
    • ✓ semi-implicit loop coded by hand vs MuJoCo Euler after 20 s (rad)expected 0 · got 4.44089×10⁻¹⁶ · absolute error 4.4×10⁻¹⁶
    • ✓ the page's two update loops compile (-Wall -Wextra -Wconversion -Werror)expected 0 · got 0 · absolute error 0
    • ✓ C++ semi-implicit loop vs MuJoCo Euler after 20 s (rad)expected 0 · got 4.44089×10⁻¹⁶ · absolute error 4.4×10⁻¹⁶

  3. Implicit-in-velocity update (MuJoCo implicit, implicitfast)

    With no velocity-dependent forces (D = 0) this reduces to semi-implicit Euler. implicitfast drops the Coriolis and centripetal terms from D. In both, D leaves out the constraint forces J^T f, as MuJoCo's documentation states.

    (M+hD) vt+h=(M+hD) vt+h Ma(vt),D≡−∂∂v(τ−c)(M + hD)\,v_{t+h} = (M+hD)\,v_t + h\,M a(v_t),\qquad D \equiv -\dfrac{\partial}{\partial v}\big(\tau - c\big)
    F7.03✓ 3/3 checked[1]
    Symbols and checks
    SymbolMeaningUnit
    DDvelocity derivative of the forces (damping-like)N·m·s/rad
    • ✓ undamped: max |q_Euler - q_implicitfast| (rad)expected 0 · got 0 · absolute error 0
    • ✓ joint damping only: max difference (rad)expected 0 · got 0 · absolute error 0Euler already treats joint damping implicitly
    • ✓ undamped: max |q_Euler - q_implicit| (rad)expected 0 · got 5.55112×10⁻¹⁶ · absolute error 5.6×10⁻¹⁶

  4. Classical fourth-order Runge–Kutta (RK4)

    Four force evaluations per step, so about four times the cost per step of Euler (MuJoCo docs, timestep).

    k1=f(yt), k2=f(yt+h2k1), k3=f(yt+h2k2), k4=f(yt+hk3),yt+h=yt+h6 (k1+2k2+2k3+k4)k_1 = f(y_t),\ k_2 = f\big(y_t+\tfrac h2 k_1\big),\ k_3 = f\big(y_t+\tfrac h2 k_2\big),\ k_4 = f(y_t+h k_3),\quad y_{t+h} = y_t + \tfrac h6\,(k_1 + 2k_2 + 2k_3 + k_4)
    F7.04✓ 1/1 checked[8][11]
    Symbols and checks
    • ✓ same cost (RK4 at 2 ms, Euler at 0.5 ms): RK4's energy error is smallerobserved as statedRK4 1.8e-09, Euler 1.9e-03

  5. Order of accuracy and what halving h buys

    Asymptotic: h small enough, and the solution smooth.

    e(h)≈C hp ⇒ e(h)e(h/2)≈2p,pmeasured=log⁡2 ⁣e(h)e(h/2)e(h) \approx C\,h^{p}\ \Rightarrow\ \dfrac{e(h)}{e(h/2)} \approx 2^{p},\qquad p_{\text{measured}} = \log_2\!\dfrac{e(h)}{e(h/2)}
    F7.05✓ 9/9 checked[7][8]
    Symbols and checks
    SymbolMeaningUnit
    eeglobal error at a fixed end timeunits of the quantity
    pporder (1 for Euler and 4 for RK4)dimensionless
    • ✓ measured order, MuJoCo Eulerexpected 1 · got 1.01198 · relative error 0.012
    • ✓ measured order, MuJoCo RK4expected 4 · got 4.02344 · relative error 0.0059
    • ✓ Euler order from differences of three runs, no referenceexpected 1 · got 1.03436 · relative error 0.034
    • ✓ straight from 60 deg: halving the step at least halves the gain in 20 s, twice overobserved as stated+0.450, +0.177, +0.079 m g L at 2, 1, 0.5 ms
    • ✓ elbow bent 90 deg, 2 s: Euler loses energy, and each halved step about halves the lossobserved as stated-0.0435, -0.0238, -0.0125 m g L at 2, 1, 0.5 ms
    • ✓ the same release over 20 s: halving the step from 2 to 1 ms does not shrink the lossobserved as stated-0.0821 then -0.1219 m g L: the runs have parted ways
    • ✓ explicit Euler: halving h halves the log of the growthexpected 2 · got 1.99473 · relative error 0.0026
    • ✓ MuJoCo Euler: halving h halves the excursionexpected 2 · got 2.00774 · relative error 0.0039
    • ✓ RK4: halving h cuts the loss by 2^5expected 32 · got 32.0025 · relative error 8×10⁻⁵per step (wh)^6, and twice as many steps

  6. What semi-implicit Euler conserves (linear oscillator)

    Exact for q̈ = −ω²q with the update in F7.02. True energy oscillates around it; at turning points (v = 0) the two agree.

    H~=v2+ω2q2−h ω2q v=const\tilde H = v^2 + \omega^2 q^2 - h\,\omega^2 q\,v = \text{const}
    F7.06✓ 1/1 checked[2][5]
    Symbols and checks
    SymbolMeaningUnit
    ωωoscillator angular frequencyrad/s
    qqpositionrad
    vvvelocityrad/s
    • ✓ modified energy drift over 1e5 steps (relative)expected 0 · got 4.20402×10⁻¹⁴ · absolute error 4.2×10⁻¹⁴

  7. Energy wobble of semi-implicit Euler

    Linear oscillator, ωh ≪ 1; a nonlinear pendulum at large amplitude wobbles somewhat less.

    ΔEE≈±ωh2at twice the swing frequency, no trend\dfrac{\Delta E}{E} \approx \pm\dfrac{\omega h}{2}\quad\text{at twice the swing frequency, no trend}
    F7.07✓ 5/5 checked[2][5]
    Symbols and checks
    • ✓ Euler energy wobble amplitude at 5 degexpected 0.00769125 · got 0.0076876 · relative error 4.7×10⁻⁴
    • ✓ no trend: mean of last 10 % minus first 10 %expected 0 · got 5.21858×10⁻⁵ · absolute error 5.2×10⁻⁵
    • ✓ MuJoCo Euler, 5 deg, 2 ms: largest energy excursion vs wh/2expected 0.00769125 · got 0.0077471 · relative error 0.0073
    • ✓ MuJoCo Euler: halving h halves the excursionexpected 2 · got 2.00774 · relative error 0.0039
    • ✓ from 30 deg, MuJoCo Euler's largest excursion vs wh/2expected 0.00769125 · got 0.00761843 · relative error 0.0095

  8. Frequency error of semi-implicit Euler

    Linear oscillator.

    cos⁡(ω~h)=1−(ωh)22 ⇒ ω~≈ω(1+(ωh)224)\cos(\tilde\omega h) = 1 - \dfrac{(\omega h)^2}{2}\ \Rightarrow\ \tilde\omega \approx \omega\Big(1 + \dfrac{(\omega h)^2}{24}\Big)
    F7.08✓ 2/2 checked[5]
    Symbols and checks
    • ✓ numerical frequency, exact vs seriesexpected 7.69132 · got 7.69132 · relative error 2.6×10⁻¹⁰
    • ✓ measured frequency of the recursionexpected 7.69132 · got 7.69132 · relative error 2.9×10⁻¹¹

  9. Explicit Euler pumps energy in

    Linear oscillator: energy grows every step, for any h.

    En+1En=1+(ωh)2\dfrac{E_{n+1}}{E_n} = 1 + (\omega h)^2
    F7.09✓ 6/6 checked[3]
    Symbols and checks
    • ✓ explicit Euler energy ratio per stepexpected 1.00024 · got 1.00024 · relative error 0
    • ✓ explicit Euler, 5 deg, 20 s at 2 ms: energy ratio vs (1 + (wh)^2)^nexpected 10.654 · got 10.5334 · relative error 0.011small-swing prediction; the pendulum is slightly softer
    • ✓ explicit Euler: halving h halves the log of the growthexpected 2 · got 1.99473 · relative error 0.0026
    • ✓ the same loop with sin q replaced by q: energy ratio vs (1 + (wh)^2)^nexpected 10.654 · got 10.654 · relative error 5.5×10⁻¹³so the shortfall at 5 deg is the pendulum's softness
    • ✓ from 30 deg, explicit Euler gains less than from 5 degobserved as statedx7.52 against x10.5
    • ✓ C++ explicit loop: energy ratio after 20 s, vs the Python runexpected 10.5334 · got 10.5334 · relative error 0

  10. RK4 numerical damping of an oscillator

    Linear oscillator. Equivalent damping ratio ≈ (ωh)⁵/144, which is negligible at small ωh, but builds up over very long runs.

    ∣R(iωh)∣2=1−(ωh)672+(ωh)8576|R(i\omega h)|^2 = 1 - \dfrac{(\omega h)^6}{72} + \dfrac{(\omega h)^8}{576}
    F7.10✓ 3/3 checked[6][7]
    Symbols and checks
    • ✓ |R(i w h)|^2 exact vs formulaexpected 1 · got 1 · relative error 1.1×10⁻¹⁶
    • ✓ RK4, 5 deg, 20 s at 2 ms: energy lost vs n (wh)^6 / 72expected 1.84004×10⁻⁹ · got 1.83447×10⁻⁹ · relative error 0.003
    • ✓ RK4: halving h cuts the loss by 2^5expected 32 · got 32.0025 · relative error 8×10⁻⁵per step (wh)^6, and twice as many steps

  11. Stability limits

    Linear system. The second is why a D-gain applied as a torque by your own code has a ceiling, K_d < 2I/h, lower once stiffness is included.

    semi-implicit Euler: ωh<2,explicit damping alone: h bI<2\text{semi-implicit Euler: } \omega h < 2,\qquad \text{explicit damping alone: } h\,\dfrac{b}{I} < 2
    F7.11✓ 5/5 checked[5]
    Symbols and checks
    • ✓ symplectic Euler: spectral radius at wh=1.99expected 1 · got 1 · relative error 3.3×10⁻¹⁶
    • ✓ symplectic Euler: unstable at wh=2.01 (rho>1)observed as stated
    • ✓ explicit damper alone (wh -> 0): eigenvalue 1 - hb/I, stable iff hb/I < 2observed as stated
    • ✓ damper plus stiffness at wh = 1: limit hb/I = 2 - (wh)^2/2 = 1.5observed as stated
    • ✓ the pendulum at 2 ms: unstable between 2 - (wh)^2/2 and 2observed as statedso stiffness lowers the damping limit, slightly at this step

Worked example: twenty seconds of free swing, three ways

The running example released from 5° at rest, with no damping, for 20 s: 10000 steps at 2 ms. Energy is computed by hand from the position and velocity (F6.10), at every step for Euler's band and at the end for the other two, not read from d->energy, which after mj_step can describe an earlier state (F5.06). Every number below comes from tools/robotics/verify.py.

  1. Predict first

    Before reading on, write down for each method whether the energy will rise, fall or stay put after 20 s, roughly by how much, and what halving the step will do. Module 00 explains why the guess is worth making even if it is wrong.[12]

  2. Explicit Euler, by hand

    MuJoCo does not offer it, so it is coded directly on the pendulum's equation of motion, θ¨=−(k/I)sin⁡θ\ddot\theta = -(k/I)\sin\theta with k=mgdk = mgd. The position moves with the velocity from before the step.[1]

    const double a = -(k / I) * std::sin(q);
q += h * v;  // the old velocity
v += h * a;

  3. Swap two lines

    Update the velocity first and move the position with the new one. For this model the result differs from MuJoCo's Euler by 4×10⁻¹⁶ rad after 20 s, which is rounding: MuJoCo's Euler is the semi-implicit kind.[4]

    const double a = -(k / I) * std::sin(q);
v += h * a;
q += h * v;  // the new velocity

  4. Run all three, then halve the step

    Explicit Euler keeps gaining energy, and halving the step halves the logarithm of the gain. Semi-implicit Euler never strays more than a fixed band, which halves with the step. RK4 loses a tiny, steady amount, and halving the step cuts it by 32, close to 252^5: twice as many steps, each losing 262^6 times less.[4]

    Energy over 20 s of free swing from 5°: ratio at the end, largest excursion, or change at the end
    MethodMeasureStep 2 msStep 1 msBehaviour
    Explicit Euler (by hand)ratio at 20 s×10.5×3.26grows every step
    MuJoCo Eulerlargest excursion±0.77 %±0.39 %wobbles, no trend
    MuJoCo RK4change at 20 s−1.8×10⁻⁷ %−5.7×10⁻⁹ %loses a little

  5. Check the numbers against the formulas

    The small-swing formulas predict all three at 2 ms. Explicit Euler comes in slightly under its prediction because a real pendulum is a little softer than its linear approximation: the same loop with sin⁡θ\sin\theta replaced by θ\theta matches the formula to nine figures.[4]

    Small-swing prediction versus the run, 5°, 2 ms
    MethodFormulaPredictedMeasured
    Explicit EulerF7.09: (1+(ωh)2)n(1+(\omega h)^2)^n×10.7×10.5
    MuJoCo EulerF7.07: ±ωh/2\pm\omega h/2±0.77 %±0.77 %
    MuJoCo RK4F7.10: n(ωh)6/72n(\omega h)^6/72−1.8×10⁻⁷ %−1.8×10⁻⁷ %

  6. Compare at equal cost, and at a bigger swing

    RK4 does four force evaluations per step, so the fair comparison is with Euler at a quarter of the step: at 0.5 ms Euler's band is still ±0.19 %, many orders of magnitude above RK4's loss at 2 ms. From 30° instead of 5°, Euler's band is about the same, ±0.76 %, while explicit Euler gains less, ×7.52, because the larger swing is softer still.[1][4]

Order of accuracy tells you how fast the error shrinks; it does not tell you which way the energy leans. Check both: halve the step, and watch the energy of an undriven run.

Best practice

Questions

Is MuJoCo's Euler the textbook Euler method?
No. It is semi-implicit: the position moves with the new velocity (F7.02), and joint damping is treated implicitly. The check behind this page reproduces it with a hand-written semi-implicit loop, to within rounding over 20 s.[1][4]
RK4 is fourth order. Why does its energy drift at all?
Order measures how fast the error shrinks with the step, not whether any quantity is conserved. RK4 is neither symplectic nor symmetric, so nothing stops its energy error growing with the length of the run; for a linear oscillator it removes a small, steady fraction per step (F7.10). Semi-implicit Euler is less accurate per step but, on a single hinge, nearly conserves a nearby energy (exactly, for a linear oscillator, F7.06).[2][6]
Why does the energy wobble twice per swing?
The quantity semi-implicit Euler conserves differs from the true energy by a term proportional to h ω2qvh\,\omega^2 q v (F7.06). Position times velocity goes through two cycles for every swing, so the difference does too. At the turning points and at the bottom of the swing, where vv or qq is zero, the two agree.[5]
Halving the step should cut the error by 2p+12^{p+1}, not 2p2^p. Which is it?
Both, for different errors. The error made in a single step falls by 2p+12^{p+1}; the error at a fixed end time, after twice as many steps, falls by 2p2^p (F7.05). The description of Moler's video states the first; this page's measured orders use the second. RK4's energy figure in the worked example falls by 252^5 for a different reason: its energy error per step is of order h6h^6 (F7.10).[7][16]
Then why not always use RK4?
It costs four force evaluations per step, its energy still drifts in very long runs, and when large velocity-dependent forces are present a single-step method that treats them implicitly can be significantly more stable, as MuJoCo's documentation notes.[1][6]
Does a smaller step make explicit Euler acceptable?
It makes it slower to go wrong, but never stops it: the energy grows by a factor 1+(ωh)21 + (\omega h)^2 every step for any positive hh (F7.09), so a long enough run always gains without bound. Halving the step from 2 ms to 1 ms took the 20-second gain from ×10.5 to ×3.26.[3][4]

Watch

Read

Sources

Numbered in order of first use on this page. read: checked against the claim. listed: exists and covers the topic, wording not checked. derived: worked out here and checked numerically.

  1. [1]

    MuJoCo Computation chapter, Numerical integration (3.17.0) (opens in a new tab)

    Google DeepMind

    explicit, semi-implicit and implicit-in-velocity Euler; D matrix; implicitfast drops Coriolis and centripetal terms; discrete integrator; §Friction loss (dry friction modelled as a constraint whose effect is to drive the joint velocity to zero); §Reproducibility; §Numerical integration, RK4 and the tip “Choosing timestep and integrator”

    The exact update equations of every MuJoCo integrator, and how friction loss is modelled. The docs recommend implicitfast for most models (similar cost to Euler, more stable), keep Euler for compatibility, and call RK4 best for energy-conserving or nearly conserving systems; they report RK4 qualitatively better than single-step methods on an energy-conserving example even at a quarter of the timestep, and single-step methods with implicit velocity-dependent forces significantly more stable than RK4 when those forces are large.

    read
  2. [2]

    Geometric numerical integration illustrated by the Störmer–Verlet method (opens in a new tab)

    Ernst Hairer, Christian Lubich, Gerhard Wanner · Acta Numerica 2003, Cambridge University Press, pp. 399–450. doi:10.1017/CBO9780511550157.006

    §1 (symplectic Euler defined); §3 (Kepler example — no drift for symplectic, drift for explicit Runge–Kutta); §4 (modified Hamiltonian, backward error analysis); §5.1 (energy approximately preserved over very long times); §5.4 (oscillations in H of size proportional to h)

    Why symplectic Euler’s energy error stays bounded, and why explicit Runge–Kutta methods drift (the paper’s general statement; it does not single out RK4).

    read
  3. [3]

    Energy growth of explicit Euler on q̈ = −ω²q

    F7.09.

    derived
  4. [4]

    Module 07 worked example (tools/robotics/verify.py)

    Twenty seconds of free swing of the running example from 5°: explicit Euler coded by hand (and on the linearised pendulum), MuJoCo’s Euler and RK4, at 2, 1 and 0.5 ms; Euler’s band from the energy at every step, the other figures from the state at the end; compared with F7.07, F7.09 and F7.10. From 30°, explicit Euler and MuJoCo’s Euler at 2 ms. The page’s two C++ update loops are compiled and run as printed, and the semi-implicit one matches MuJoCo’s Euler. Checks: integration_example.

    derived
  5. [5]

    Invariant, energy wobble, frequency shift and stability of semi-implicit Euler on q̈ = −ω²q

    F7.06–F7.08, F7.11. Consistent with Hairer, Lubich & Wanner (2003) §4–5, which give the general result (bounded energy error of size O(h)).

    derived
  6. [6]

    RK4 amplification factor on q̈ = −ω²q from its stability function R(z) = 1 + z + z²/2 + z³/6 + z⁴/24

    F7.10.

    derived
  7. [7]

    Solving Ordinary Differential Equations I — Nonstiff Problems, 2nd revised edition (opens in a new tab)

    Ernst Hairer, Syvert P. Nørsett, Gerhard Wanner · Springer, 1993. ISBN 978-3-540-56670-0

    Order of accuracy, the classical RK4 method and its stability function.

    listed
  8. [8]

    Numerical Methods for Ordinary Differential Equations, 3rd edition (opens in a new tab)

    J. C. Butcher · Wiley, 2016. ISBN 978-1-119-12150-3

    Runge–Kutta methods, order, stability functions.

    listed
  9. [9]

    Measured order of MuJoCo’s Euler and RK4 (tools/robotics/verify.py)

    The running example released from 30°: angle after 2 s at steps of 4, 2 and 1 ms, against an RK4 reference at 10 µs (Euler from 2 to 1 ms, RK4 from 4 to 2 ms), and for Euler also from the differences of the three runs with no reference. Checks: convergence_orders.

    derived
  10. [11]

    MuJoCo XML Reference (3.17.0) (opens in a new tab)

    Google DeepMind

    option-timestep (default 0.002 s; RK4 costs about four times per step); option-gravity (default 0 0 −9.81); option-integrator; option-flag-energy; option-flag-eulerdamp (implicit joint damping in Euler, on by default); body-joint-armature (reflected inertia, gear ratio squared; positive armature improves stability); actuator-general-armature (rotor inertia, scaled by gear squared); body-joint-damping (Euler treats joint damping implicitly); body-joint-frictionloss; compiler-inertiafromgeom; body-geom-density (default 1000 kg/m³, used only when mass is not given)

    MuJoCo defaults, energy definitions, armature, damping and friction-loss behaviour. Read from the 3.17.0 source of the docs (doc/XMLreference.rst).

    read
  11. [12]

    The pretesting effect: Do unsuccessful retrieval attempts enhance learning? (opens in a new tab)

    Lindsey E. Richland, Nate Kornell, Liche Sean Kao · Journal of Experimental Psychology: Applied 15(3), 243–257 (2009). doi:10.1037/a0016496

    Taking a test before studying, even one mostly failed (5–22 % correct), improved later recall of the pretested information more than spending the same time studying; untested information did not differ from the extended-study condition. The material was a factual text on vision.

    read
  12. [13]

    MuJoCo Modeling chapter, solver parameters and friction (3.17.0) (opens in a new tab)

    Google DeepMind

    Solver parameters: solref default "0.02 1" (time constant, damping ratio 1: critically damped); timeconst at least twice the timestep (enforced unless refsafe is off); Restitution (zero damping in the direct format); friction loss is first-order, purely damped, no exact stick (slow slippage)

    Constraints are soft mass-spring-dampers, critically damped by default; restitution needs the direct format with zero damping, and even then energy only fluctuates around its initial value. Why simulated friction loss creeps instead of sticking.

    read
  13. [14]

    mj_RungeKutta and mj_forwardSkip, src/engine/engine_forward.c (MuJoCo 3.17.0) (opens in a new tab)

    Google DeepMind, GitHub

    RK4 evaluates forward dynamics at each intermediate stage with skipsensor set and a comment saying energy is not recomputed, but mj_forwardSkip recomputes energy whenever its own flags are clear, so energy is recomputed at every stage; the last evaluation is at the final stage. The Computation chapter says the optional computations run only once.

    read
  14. [15]

    Compiled and run checks for module 05 (tools/robotics/verify.py)

    Module 05’s checks in tools/robotics/verify.py: MuJoCo 3.17.0 behaviour (the step loop and its split under every integrator, energy timing including RK4’s stages and energy sensors, joint sizes, angle units, default density, joint and actuator armature, determinism); the worked example’s two XML models compiled and read back field by field and with mj_printFormattedModel; its GoogleTest, extracted from the page, compiled against the C library from the same release and run against both models and a warning case; mj_loadXML probed with a clean model, one with a compiler warning and one that fails; and five of the sheet’s six C++ entries (F5.05, F5.06, F5.10–F5.12) compiled with -Wconversion -Wsign-conversion -Werror.

    derived
  15. [16]

    Order, Naming Conventions (opens in a new tab)

    Cleve Moler · MIT OpenCourseWare RES.18-009, Learn Differential Equations: Up Close with Gilbert Strang and Cleve Moler (2015)

    OCW page description (video ID 6O9D6am_RK4)

    “A method is said to have order p if cutting the step size in half reduces the error in one step by a factor of two to the power p+1.” Read from the OCW page’s description; the video itself was not watched for this page.

    read
  16. [17]

    Euler, ODE1 (opens in a new tab)

    Cleve Moler · MIT OpenCourseWare RES.18-009, Learn Differential Equations: Up Close with Gilbert Strang and Cleve Moler (2015)

    Introduces numerical methods for ODEs through Euler’s method (MATLAB’s ODE1), with exponential growth and compound interest as examples. The first of the 12 videos in the OCW gallery “Solving ODEs in MATLAB”. Video ID and description from the OCW page; not watched for this page.

    listed
  17. [18]

    Classical Runge-Kutta, ODE4 (opens in a new tab)

    Cleve Moler · MIT OpenCourseWare RES.18-009, Learn Differential Equations: Up Close with Gilbert Strang and Cleve Moler (2015)

    The classical Runge–Kutta method, described as the most widely used numerical method for ODEs; its drawback is the lack of an error estimate. Video ID and description from the OCW page; not watched for this page.

    listed
  18. [19]

    Differential equations, studying the unsolvable (DE1) (opens in a new tab)

    Grant Sanderson · 3Blue1Brown, 31 March 2019

    An animated introduction to differential equations built around the pendulum equation.

    listed