Numerical integration

With the accelerations right, longer runs exposed the next problem: semi-implicit Euler drifts. A frictionless double pendulum should conserve energy, yet over 10 seconds at a 2 ms step Euler's total energy wanders by up to 0.56 J. Euler updates the velocity from the current acceleration, then moves the position with that new velocity, which is simple but only first-order accurate. Velocity Verlet splits the velocity update into two half-steps around the position update and recomputes the accelerations in between, halving the drift to 0.28 J. Its weakness is that velocity-dependent forces like Coriolis are evaluated at an estimated velocity. Implicit midpoint evaluates the dynamics at the average of the old and new states, making each step time-symmetric. Because the new state isn't known in advance, it iterates four times to converge, and it stays within 0.0006 J. RK4 instead samples the dynamics four times across the step and blends the results, reaching fourth-order accuracy and staying within 0.0001 J.

V2's dynamics were wrong in ways only a reference could show: brought back to life next to MuJoCo, a one-link arm swung exactly 4× too fast (it used the free-body inverse inertia, dropping the parallel-axis term), and a two-link arm moved in the wrong shape. V3 runs the full Articulated Body Algorithm, with Coriolis and centrifugal bias terms, parent-to-body transforms and gravity as a fictitious base acceleration, and steps it forward with semi-implicit Euler, the same default MuJoCo uses. Here all three drive the same arm with the same shoulder torque and no gravity, since V2 never had any. V3 stays within 3×10⁻⁵ rad of MuJoCo for all 4 seconds, while V2 has swung the shoulder 4.5 times as far and folded the elbow forward instead of letting it trail behind.

MuJoCo is only a reference, so I also checked V3 against an exact answer. A double pendulum is chaotic and has no closed-form trajectory: two simulations that differ by a single rounding error eventually go completely different ways. Its equations of motion do have a closed form, though. Lagrangian mechanics reduces them to a 2×2 linear system that gives the exact joint accelerations at any given position and velocity. I derived that system independently of the engine's code, then compared it with V3's accelerations at every state of a 10-second run. Because each state is checked on its own, errors can't build up and chaos can't amplify them. The two agree to within 3×10⁻⁶, which is float32 round-off, and V3 matches MuJoCo to the same 3×10⁻⁶.