mathematics//numerical methods//numerical integration//Runge-Kutta method

A Runge-Kutta method is a scheme of numerical integration that evaluates the dynamics at several points inside each step and combines those slopes into one weighted average, and its classic fourth-order version, RK4, is the default when a simulation must be accurate with a reasonable step: a drone simulator, an offline trajectory study, the plant model in a test bench. With \(\dot x=f(x)\):


A Runge-Kutta method is a scheme of numerical integration that evaluates the dynamics at several points inside each step and combines those slopes into one weighted average, and its classic fourth-order version, RK4, is the default when a simulation must be accurate with a reasonable step: a drone simulator, an offline trajectory study, the plant model in a test bench. With x˙=f(x)\dot x=f(x)x˙=f(x):

k1=f(xk),k2=f ⁣(xk+Δt2k1),k3=f ⁣(xk+Δt2k2),k4=f(xk+Δt k3),xk+1=xk+Δt6 (k1+2k2+2k3+k4).\begin{aligned} k_1&=f(x_k), & k_2&=f\!\left(x_k+\tfrac{\Delta t}{2}k_1\right),\\ k_3&=f\!\left(x_k+\tfrac{\Delta t}{2}k_2\right), & k_4&=f(x_k+\Delta t\,k_3),\\ x_{k+1}&=x_k+\tfrac{\Delta t}{6}\,(k_1+2k_2+2k_3+k_4). \end{aligned}k1​k3​xk+1​​=f(xk​),=f(xk​+2Δt​k2​),=xk​+6Δt​(k1​+2k2​+2k3​+k4​).​k2​k4​​=f(xk​+2Δt​k1​),=f(xk​+Δtk3​),​

k1k_1k1​ is the slope at the start, k2k_2k2​ and k3k_3k3​ two estimates at the midpoint, k4k_4k4​ the slope at the end; the middle ones weigh double.

The payoff is the order. RK4 matches the Taylor polynomial of the true solution up to the fourth power of the step, so its accumulated error is proportional to Δt4\Delta t^4Δt4: halving the step divides the error by sixteen, where the Euler method only halves it. Each step costs four evaluations of fff instead of one, and in exchange it allows steps many times larger for the same error, which is usually the better deal.

It is more stable than Euler, and still bounded. On a decaying mode x˙=λx\dot x=\lambda xx˙=λx RK4 stays stable up to about Δt<2.79/∣λ∣\Delta t<2.79/|\lambda|Δt<2.79/∣λ∣, against 2/∣λ∣2/|\lambda|2/∣λ∣ for explicit Euler: with λ=−50\lambda=-50λ=−50 per second, 56 ms instead of 40 ms. More margin, but a ceiling all the same, so a system with a very fast mode still forces a small step (numerical stability, stiffness). On an undamped oscillator with large steps it slowly loses energy instead of gaining it, the safer direction.

Embedded variants estimate their own error. RK45 (Dormand-Prince, the default of SciPy's solve_ivp and MATLAB's ode45) computes a fourth and a fifth order solution from the same evaluations, uses the difference as an error estimate and adapts the step to meet a tolerance. That makes it efficient offline and unsuitable where every step has a deadline.

The cost is modest on today's hardware. One RK4 step of a 12-state rigid-body model is on the order of a thousand floating-point operations, microseconds on a microcontroller with an FPU, which is why a planar drone simulated in Python with RK4 is a good first project for learning control (planar drone).