Simulating Chaotic Dynamical Systems & N-Body Orbits with Runge-Kutta Methods

By Aman Sahani•

2026-09-08

•

5 min read

⚡ Key Takeaway (TL;DR)

Phase-space analysis, Lyapunov exponents, and numerical integration techniques for non-linear physical equations on ARM64.

Simulating Chaotic Dynamical Systems & N-Body Orbits with Runge-Kutta Methods

Nonlinear physical systems, from the Lorenz atmospheric attractor to gravitational $N$ -body celestial mechanics, exhibit sensitive dependence on initial conditions.

The Lorenz Attractor Equations

The continuous 3D phase-space system is governed by:

Mathematical Model / Equation

\begin{aligned} \frac{dx}{dt} &= \sigma (y - x) \\ \frac{dy}{dt} &= x (\rho - z) - y \\ \frac{dz}{dt} &= xy - \beta z \end{aligned}

With canonical Rayleigh parameters $\sigma = 10$ , $\rho = 28$ , $\beta = \frac{8}{3}$ .

Fourth-Order Runge-Kutta (RK4) Integrator

To integrate $\frac{d\mathbf{y}}{dt} = f(t, \mathbf{y})$ with step size $h$ :

Mathematical Model / Equation

\mathbf{y}_{n+1} = \mathbf{y}_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)

Where:

Mathematical Model / Equation

k_1 = f(t_n, \mathbf{y}_n), \quad k_2 = f\left(t_n + \frac{h}{2}, \mathbf{y}_n + \frac{h}{2}k_1\right)

Mathematical Model / Equation

k_3 = f\left(t_n + \frac{h}{2}, \mathbf{y}_n + \frac{h}{2}k_2\right), \quad k_4 = f(t_n + h, \mathbf{y}_n + hk_3)

cpp

struct State { double x, y, z; };

State rk4_step(const State& s, double dt, double sigma, double rho, double beta) {
    auto deriv = [&](const State& st) -> State {
        return { sigma * (st.y - st.x), st.x * (rho - st.z) - st.y, st.x * st.y - beta * st.z };
    };
    State k1 = deriv(s);
    State s2 = { s.x + 0.5*dt*k1.x, s.y + 0.5*dt*k1.y, s.z + 0.5*dt*k1.z };
    State k2 = deriv(s2);
    State s3 = { s.x + 0.5*dt*k2.x, s.y + 0.5*dt*k2.y, s.z + 0.5*dt*k2.z };
    State k3 = deriv(s3);
    State s4 = { s.x + dt*k3.x, s.y + dt*k3.y, s.z + dt*k3.z };
    State k4 = deriv(s4);

    return {
        s.x + (dt/6.0)*(k1.x + 2*k2.x + 2*k3.x + k4.x),
        s.y + (dt/6.0)*(k1.y + 2*k2.y + 2*k3.y + k4.y),
        s.z + (dt/6.0)*(k1.z + 2*k2.z + 2*k3.z + k4.z)
    };
}

The maximal Lyapunov exponent $\lambda \approx 0.9056$ implies prediction horizon $T_{\text{horizon}} \sim \frac{1}{\lambda} \ln \frac{\epsilon_{\text{tolerance}}}{\delta_0}$ .