Simulating Chaotic Dynamical Systems & N-Body Orbits with Runge-Kutta Methods
Nonlinear physical systems, from the Lorenz atmospheric attractor to gravitational -body celestial mechanics, exhibit sensitive dependence on initial conditions.
The Lorenz Attractor Equations
The continuous 3D phase-space system is governed by:
\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$: $$\mathbf{y}_{n+1} = \mathbf{y}_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$ Where: $$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)$$ $$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}$.