How a simulator steps through time

Every lesson so far has used equations of motion: the rocket equation, the gravity turn's two equations for speed and flight-path angle, Kepler's orbits. The playground does not solve any of them. It does something much simpler and much more general. It takes the vehicle's state now — where it is, how fast it is moving, which way it points, how fast it is turning — works out every force and torque acting on it, and uses them to guess the state a fiftieth of a second later. Then it does it again, about 11,000 times for a hop and several million times for a mission that waits a day in orbit.

How good that guess is, how the errors add up, what they cost and how to tell whether they matter is the subject of numerical integration. It decides whether a simulated orbit closes on itself, whether a landing leg bounces or explodes, and whether two runs of the same flight give the same answer.

A flight as an initial-value problem

Write the state as one list of numbers, $\mathbf{y}$. For a vehicle in the playground it has thirteen: three of position, three of velocity, four for the attitude quaternion and three body rates. The laws of motion say how fast each of them is changing, given all of them and the time:

The function $\mathbf{f}$ is where all the physics lives: gravity with its $J_2$ term, the air and the wind, the engines, the fins, the legs. Given the state at the start, the equation fixes the whole flight. That is an initial-value problem, and a flight simulator is a machine for solving one.

Almost never can it be solved exactly. What a simulator does instead is choose a step $h$ and replace the continuous flight by a sequence of states ${\mathbf{y}}_0,{\mathbf{y}}_1,{\mathbf{y}}_2,\dots$ a step apart, each computed from the one before by a fixed recipe, an integrator. All the questions of this lesson are about the recipe and the step.

Euler, and why an orbit spirals out

The simplest recipe assumes the rates stay what they are now for the whole step:

This is explicit Euler. For an orbit it means: move along the present velocity for $h$ seconds, and change the velocity by the present acceleration.

Consider what one step does on a circular orbit, where the acceleration is at right angles to the velocity. The new position is along the tangent, and a tangent leads outwards: the radius grows from $r$ to $\sqrt{r^2+h^2{v}^2}$. The new velocity is the old one plus a small perpendicular piece, so the speed grows too, from $v$ to $\sqrt{v^2+h^2{a}^2}$. Both raise the orbital energy. Working it through, each step adds $h^2{a}^2$ to the energy per kilogram, which for a circular orbit with mean motion $n$ is a fraction $2(nh{)}^2$ of it, and over one orbit's $2\pi /(nh)$ steps

The error is proportional to the step: Euler is a first-order method. In a 200 km orbit, with $n$ = 1.18 × 10⁻³ per second, a 1 s step gains 1.5 % of the energy every orbit (1.44 % measured). A 10 s step gains 11.6 % in the first orbit, and the orbit visibly spirals out.

Semi-implicit Euler

Change one thing. Update the velocity first, and move with the new velocity:

This is semi-implicit (or symplectic) Euler. It costs exactly the same, one force evaluation per step, and it is still first order. But its errors have a different character. It preserves areas in position–velocity space, as the real motion does, and as a consequence it very nearly conserves a slightly altered energy for as long as it runs. The true energy then wobbles within a band instead of drifting: with a 10 s step in the 200 km orbit it swings by 1.4 × 10⁻⁴ of itself round the orbit and comes back to within 2 × 10⁻¹¹ at the end of it. The error goes into the phase instead. The vehicle keeps the right orbit but runs slightly ahead of where it should be: 4.1 km after one orbit, 12.3 km after three. For long simulations of planetary systems, where getting the shape of the orbits right for millions of revolutions matters more than the exact phase, symplectic methods are the standard tool.

Runge–Kutta

Euler uses one slope, the one at the start of the step. The classical fourth-order Runge–Kutta method samples four slopes across the step and averages them:

The weights are chosen so that the step agrees with the Taylor series of the true solution up to the $h^4$ term. The error of one step is then of order $h^5$, and after the $1/h$ steps of a fixed interval, of order $h^4$: RK4 is a fourth-order method. It costs four force evaluations per step, and buys an enormous amount with them. In the same 200 km orbit at 10 s it is 2 cm from the exact position after one orbit, and its energy decays by 4 × 10⁻¹¹ per orbit: RK4 is slightly dissipative, so its energy error is a slow, steady drift rather than a wobble.

RK4 is what the playground integrates every vehicle with, at a step of 20 ms. At that step a 200 km orbit's error per revolution is far below the rounding of the numbers themselves, which is why the audit of the simulator found a coasting body's energy and angular momentum conserved to about 10⁻¹³.

Start with explicit Euler at 10 s for three orbits: the orbit spirals out to 8,600 km and the energy error climbs past 20 %. Switch to semi-implicit Euler at the same step and cost: the orbit stays on the dashed line, the energy error wobbles below 10⁻³ and returns, and the readout shows the price — 12 km of along-track error. RK4 at 10 s: 7 cm after three orbits, and an energy error growing by the same small amount each orbit. Then switch to the elliptical transfer orbit, from 200 km up to geostationary height, and watch how badly a fixed step fits it: near the low point the vehicle moves at 10 km/s and turns fast, near the high point it barely moves, and the same step size serves neither.

Error order: halve the step

A method of order $p$ has a global error that behaves, for small steps, like

so halving the step divides the error by $2^p$: by 2 for Euler, by 16 for RK4. That is the most useful fact in numerical integration, because it can be checked. Run the same problem at $h$ and at $h/2$; if the ratio of the errors — or of the differences between successive runs, when there is no exact answer — is close to $2^p$, the step is small enough that the error is behaving as the theory says, and its size can be estimated from the difference. If the ratio is something else, something other than the integrator's truncation is dominating the error.

After one and a quarter orbits, halving the step from 10 s to 5 s divides Euler's error by 1.9 (it approaches 2 as the step shrinks), semi-implicit Euler's by 2.0 and RK4's by 16.8. Now measure after exactly one orbit instead. Semi-implicit Euler suddenly divides its error by 4.0 — it seems second order. It is not: its first-order error is the wobble, and the wobble comes back to nearly nothing at the end of every revolution, leaving only a smaller second-order error behind. A convergence test made at a special moment can be fooled by a special cancellation. Make it where nothing special happens.

The fifth-order formula of the Dormand–Prince pair (next section), flown at a fixed step on the transfer orbit, divides its error by 35 to 38 per halving, close to $2^5=32$. On the circular orbit it does not behave so neatly: its errors there are small enough to cancel partly and to meet the rounding floor, around 10⁻⁷ m, within a few halvings.

Plot the errors against the work — force evaluations per orbit — and the comparison that matters appears. For anything better than kilometres, the fourth-order method wins by orders of magnitude, even at four evaluations a step.

Adaptive steps: an embedded pair

A fixed step has to be small enough for the hardest part of the trajectory, and then wastes effort everywhere else. The alternative is to let the integrator estimate its own error and choose each step to keep that error at a set tolerance.

The standard way is an embedded pair: two Runge–Kutta formulas of different order that share their slope evaluations, so the second costs almost nothing. Dormand and Prince's 5(4) pair of 1980 takes seven slopes per step and combines them twice, once into a fifth-order solution and once into a fourth-order one. Their difference estimates the error of the lower-order one. If it is within the tolerance the step is accepted — and the fifth-order solution kept, since it is the better one — and if not, the step is retried shorter. Either way the next step is scaled by

where the exponent is one over the order of the step's error, and the 0.9 and the limits keep one bad estimate from collapsing or exploding the step. The seventh slope is evaluated at the accepted solution, so it is the first slope of the next step (first same as last), and an accepted step costs six new force evaluations.

On the transfer orbit at a tolerance of 10⁻⁹, the steps run from 42 s near the low point to 584 s near the high point. One orbit takes 158 steps and 955 force evaluations and ends 3 m from the exact position. RK4 at a fixed 10 s takes 15,148 evaluations for 0.56 m; stretched to 60 s, where it still costs more than twice as much as the adaptive pair, it ends 892 m off. The more the pace of the problem varies, the more an adaptive step is worth.

In the playground

The playground's nbody.ts contains this integrator, with the controller above and an error measured relative to the state's own size, $\epsilon =\max (|\delta \mathbf{r}|/|\mathbf{r}|,|\delta \mathbf{v}|/|\mathbf{v}|)$. Its role is narrower than it might seem. The cruise rails that carry a vehicle between the planets — a program asks for them with fc.cruise() — do not use it. They choose each step from the geometry instead: at most an hour, at most as long as the trajectory takes to turn through 0.5°, and never more than half the time it would take to reach a 2,000 km clearance shell round any body or a sphere of influence. Each such step is then taken by classical RK4 through the gravity of the Sun, the Earth, the Moon and, on a Mars mission, Mars. The reason given in the code is cost and generality: an adaptive solver would re-evaluate the planets' positions at every stage of every attempt, and it is written for a frame the Moon- and Mars-centred flights do not use. The Dormand–Prince integrator is the reference the tests measure the cruise against.

The comparison is instructive. A vehicle sent out of a 300 km orbit at 1.4 times circular speed flies its first 11 minutes, until it is 2,000 km up, in 34,184 steps of 20 ms, then 673 cruise steps with a median of 375 s. After 72 hours and 310,000 km it is 0.39 km from the Dormand–Prince solution at a tolerance of 10⁻¹², which covers the same 72 hours in 455 steps. A step rule taken from the physics is adaptive too; it simply does not measure its own error.

Why the flight itself runs at a fixed 50 Hz

For the flight itself — everything with engines, air, fins or legs — the playground does the opposite of adapting. Every vehicle advances in steps of exactly 20 ms, and the flight program's update(fc) is called once before each of them, fifty times per simulated second. There are four reasons, and only one of them is about accuracy.

The flight computer is a sampled system. Real flight software runs in fixed cycles, reads its sensors once a cycle and holds its commands until the next. The playground's autopilot and actuators work the same way. That makes the control loop part of the model: change the step and you change the controller, not only the integration. The simulator shows this plainly. A burn in vacuum with the engines alone, integrated with the thrust and mass evaluated at every Runge–Kutta stage, gives the same state to within a micrometre at every step from 40 ms to 1.25 ms: the integration is exact to rounding. Let the autopilot turn the vehicle 2° during the same burn and the result moves by 0.43 m between 20 ms and the finest step, and converges unevenly, at roughly first order. That is not an error to be integrated away. It is what a 50 Hz flight computer does.

Events happen on ticks. Ignition, shutdown, staging and touchdown happen at step boundaries. An adaptive integrator meeting an engine cut or a foot striking the ground would see its error estimate explode and cut its step again and again to resolve the corner; a fixed step simply takes it.

Determinism. The same vehicle, mission and program give the same flight, bit for bit, however the flight is watched. The app does not step the simulation per frame. Each animation frame it adds the real time elapsed, multiplied by the time warp, to an account of simulated time owed, then runs whole 20 ms steps in chunks of about a millisecond of computing until the account is paid or 11 ms of the frame have gone. The step is never lengthened: a warp the computer cannot keep up with makes the simulation run slower, not coarser. The picture is interpolated between the last two states. So a hop flown at 60 or 144 frames a second, at ×1, ×10 or ×100, or as fast as the machine allows, performs the same 11,438 steps and ends in the same bits — checked here by hashing the final state of each.

A budget. A fixed step has a fixed cost, which is what a real-time simulation needs. One step of a vehicle takes a fraction of a millisecond; during an ascent, with a booster flying home and running its predictors, about half a millisecond on the machine these figures were measured on. The 11 ms of a 60 Hz frame then buy about twenty steps, which is why the warp that the simulation can actually keep up with is lower during a busy ascent than on a quiet coast.

What is held for a whole step is short to list: the commands, as on any flight computer; the wind, sampled once per step at the vehicle's position; and, on the classic equatorial range, which reproduces the original model, the thrust and the mass as well. That last makes powered flight first order there: 9.4 m of error after 70 s of burn at 20 ms, against 0.27 m with the autopilot flying on a real site.

Long coasts get the same step with a cheaper sum. A vehicle whose program is asleep above 151 km is advanced on rails: still 20 ms, but with only gravity and, on real sites, point-mass drag to evaluate. Only the cruise between planets changes the step, and only when nothing can act on the vehicle.

Stiffness: when stability, not accuracy, sets the step

Some forces act much faster than anything the flight cares about. A landing leg's spring, the friction that holds a foot where it stands, a crush core yielding in milliseconds. Such a system is stiff: it contains a mode so fast that the step must be tiny to follow it, even when the motion that matters is slow.

To see why, take one mode, a spring and a damper, $\ddot{x}+2\zeta \omega \dot{x}+{\omega }^{2}x=0$. Its solutions behave like $e^{\lambda t}$ with $\lambda =\omega \left(-\zeta \pm \sqrt{{\zeta }^2-1}\right)$, a pair of complex numbers when it oscillates and two real ones when it is overdamped. One step of a method multiplies the mode by a fixed factor $R(\lambda h)$:

The step is stable only if $|R(\lambda h)|\le 1$: the product of the mode's eigenvalue and the step must lie inside the method's stability region. For explicit Euler that is a disc of radius 1 centred on −1. For RK4 it is larger, reaching −2.785 along the real axis and ±2.83 along the imaginary one. Outside it, the numerical solution grows by $|R|$ every step, however strongly the real one decays. And for an undamped oscillation explicit Euler is unstable at any step at all.

With a damping ratio of 0.7, RK4 at 20 ms is stable only for modes below 21.5 Hz, explicit Euler below 11.2 Hz. A leg's bounce, at 1.5 Hz in the playground, is no trouble: $\omega h$ = 0.19. The friction of the feet is another matter. In the realistic model a foot that has stuck to the ground is held by a tangential spring of 3 Hz, damped at 0.7, on a quarter of the mass the legs are built for. The four feet sit 8.2 m from the axis of a vehicle whose roll inertia is small, and on the landed Hopper they would hold its roll with a mode of about 22.5 Hz and a damping ratio of 5.3. That damping is so heavy that the mode's fast part decays at about 1,480 per second: $\lambda h\approx -30$ at 20 ms, and −7.4 even in 5 ms sub-steps. Both are far outside RK4's region.

The figure opens on that case: RK4 in four 5 ms sub-steps, the fast eigenvalue off the plot to the left, and the numerical response gone after one step while the exact one settles smoothly. Switch to implicit Euler, which evaluates the force at the end of the step and so has to solve an equation to take it: it is stable at any step, because its region is everything outside a disc on the right. Bring the damping ratio down to 0.7 and the frequency to 15 Hz and RK4 is comfortable again.

The playground stays explicit and does two things instead. Within 2 m of the ground it splits each 20 ms step into four sub-steps of 5 ms, which moves every eigenvalue four times closer to the origin. And it limits the stiffness of the contact to what that sub-step can take: the feet's friction is capped so that its damping rate times the sub-step is at most 1, and the stuck-foot spring so that $\omega h$ is at most 0.5. On the landed Hopper that softens the stuck feet's roll spring from about 22.5 Hz to 15.9 Hz. It is a physical simplification made for a numerical reason — the kind every real-time simulator makes somewhere — and it is written down in contact.ts next to the constants it changes.

How to know a result is converged

A number from a simulation has three kinds of error: the model's (the physics left out or simplified), the integration's (the step), and the arithmetic's (rounding). The last two can be measured.

Sweep the step. Run the same flight at 40, 20, 10, 5 and 2.5 ms and look at the quantity you care about. The playground allows it: a headless new Simulation({ dt }) takes any step up to 0.1 s, and everything that depends on the step — the autopilot's filters, the predictors' timing, the program's fc.dt — follows it. The Hopper's hop from Cape Canaveral on a calm day lands 0.25, 0.20, 0.17, 0.16 and 0.16 m from the centre of the pad: converged to a few centimetres. Its touchdown speed does not converge: 1.21, 1.33, 1.21, 1.18 and 1.70 m/s. It is set by the tick on which the program cuts the engine, a discrete event, and a step half the size moves the cut to a different tick. The upper stage of the full mission, cut off when a threshold on the orbit is crossed while accelerating at 5 g, ends with an apoapsis of 244.0, 241.2, 237.2 and 235.9 km at 40, 20, 10 and 5 ms: every 20 ms step at that moment adds 1 m/s, which is 3.4 km of apoapsis, and the cut falls wherever the threshold does within a step. A quantity that depends on when a discrete event happens converges, if at all, at first order.

Expect the sweep to find things. The playground's own audit swept a hop in wind at 20, 10 and 5 ms and found it landing 2.5, 3.6 and 34.7 m off the pad: the step was hard-coded in the autopilot's rate estimate and several filters, so a smaller step changed the controller wrongly. With every step-dependent constant following the simulation's step, the same comparison now agrees to under a metre, and a test holds it there.

Check what must be conserved. A coasting body's energy and angular momentum are known exactly; a burn in vacuum must give the rocket equation's $\Delta v$. The playground's audit checks both, and an unpowered descent through the air converges at the fourth order it should: the error falls by 13 and then 18 times per halving.

Compare with an independent method. The cruise rails against Dormand–Prince; a hand-written reference integrator against the vehicle model. Agreement between two different methods is much stronger evidence than a small change between two runs of one.

Floating point, and the same answer twice

A double-precision number carries 53 bits of mantissa, about sixteen decimal digits. At the Earth's radius, 6.4 × 10⁶ m, adjacent doubles are 0.93 nm apart: that is the finest position the simulator can represent. And arithmetic on them is not associative:

So the same sums done in a different order give different bits. Even the clock is not exact: 0.02 has no finite binary expansion, and adding it to itself through a 30-hour mission, 5.4 million times, drifts by 4 µs.

Bit-for-bit determinism therefore needs more than a fixed step. It needs every operation done in the same order every time. The playground takes care of that: random draws come from seeded streams, one per quantity; vehicles are stepped in the order they were created; the other bodies' positions are updated on a fixed schedule of mission time; and reading telemetry for the screen never touches the numbers the program sees. The language adds one caveat of its own: JavaScript does not specify the last bit of functions such as Math.sin or Math.exp, so identical results are guaranteed on one engine, not across all of them.

How much do the bits matter? Ten seconds into the Hopper's hop, change one coordinate of its position by the smallest amount a double can change it there, 0.12 nm, and fly on. The two flights never again agree in every bit; they drift apart by a few hundredths of a micrometre and no more, and land at the same point to a micrometre, because the guidance steers both onto the same pad. Determinism is not accuracy. The step sweep says that landing point is known to a few centimetres. What determinism buys is that a flight can be shared as a link, replayed, and argued about, and that a difference between two runs always means something changed.

Try it

Choose the Hopper preset, launch from Cape Canaveral with the weather set to Custom — wind, gusts and turbulence at zero — and Hardware dispersions off under Failures & sensors, and fly the Hop example at normal speed. Note the flight report's landing: 0.2 m from the centre at 1.3 m/s, with 10,144 kg of propellant left. Fly it again with the time warp at its highest setting. The report is identical to the last digit: the warp changed how many 20 ms steps each frame computed, and nothing else.

Add fc.log(fc.t, fc.dt) at the top of update, fly a second or two, and read the console: the program is called every 0.02 s of mission time, whatever the frame rate. Then fly the Full mission example at the highest warp and watch the achieved warp in the top bar: modest while the booster is flying home and running its predictors, far higher once the upper stage is asleep in orbit. The step is the same 20 ms throughout; a step on rails simply costs much less to compute.

What carries forward

Every number in these courses — the losses of an ascent, the orbits, the landings — came out of this machinery: a fixed 20 ms step of RK4 for anything under its own power, cheaper sums for a coast, and a looser step for the space between planets. Knowing where it is exact, where it is only first order, and where it was softened to stay stable is what lets a number from the simulator be trusted, and says how far. For what the model itself leaves out, see fidelity and its limits; for the equations it integrates, the physics model; and for the orbits it checks itself against, orbits as conic sections and Lambert's problem.