Lorenz '63 — RK4 integrator

//! Lorenz '63 attractor — RK4 integrator.
//!
//! Three coupled ODEs:
//!     dx/dt = σ (y − x)
//!     dy/dt = x (ρ − z) − y
//!     dz/dt = x y − β z
//!
//! Pure, deterministic, allocation-free. V-class L0-closed.
//! The state lives in a `LorenzState`; the page advances it one RK4 step
//! per frame and pushes the trail into a ring buffer for drawing.

use serde::{Deserialize, Serialize};

/// Lorenz parameters. Defaults are the classic Lorenz '63 values:
///   σ = 10, ρ = 28, β = 8/3.
#[derive(Debug, Clone, Copy, Serialize, Deserialize)]
pub struct LorenzParams {
    pub sigma: f64,
    pub rho: f64,
    pub beta: f64,
    pub dt: f64,
}

impl Default for LorenzParams {
    fn default() -> Self {
        Self {
            sigma: 10.0,
            rho: 28.0,
            beta: 8.0 / 3.0,
            dt: 0.005,
        }
    }
}

/// Lorenz state — three coordinates plus a step counter.
#[derive(Debug, Clone, Copy, Serialize, Deserialize)]
pub struct LorenzState {
    pub x: f64,
    pub y: f64,
    pub z: f64,
    pub step: u64,
}

impl Default for LorenzState {
    fn default() -> Self {
        Self {
            x: 1.0,
            y: 1.0,
            z: 1.0,
            step: 0,
        }
    }
}

/// One RK4 step of the Lorenz system. Updates `state` in place. Returns
/// the bit-ops estimate for this step (for the meter).
pub fn lorenz_rk4_step(state: &mut LorenzState, p: &LorenzParams) -> u64 {
    fn f(x: f64, y: f64, z: f64, p: &LorenzParams) -> (f64, f64, f64) {
        (
            p.sigma * (y - x),
            x * (p.rho - z) - y,
            x * y - p.beta * z,
        )
    }

    let (kx1, ky1, kz1) = f(state.x, state.y, state.z, p);
    let h2 = p.dt * 0.5;
    let (kx2, ky2, kz2) = f(state.x + h2 * kx1, state.y + h2 * ky1, state.z + h2 * kz1, p);
    let (kx3, ky3, kz3) = f(state.x + h2 * kx2, state.y + h2 * ky2, state.z + h2 * kz2, p);
    let (kx4, ky4, kz4) = f(
        state.x + p.dt * kx3,
        state.y + p.dt * ky3,
        state.z + p.dt * kz3,
        p,
    );

    let sixth = p.dt / 6.0;
    state.x += sixth * (kx1 + 2.0 * kx2 + 2.0 * kx3 + kx4);
    state.y += sixth * (ky1 + 2.0 * ky2 + 2.0 * ky3 + ky4);
    state.z += sixth * (kz1 + 2.0 * kz2 + 2.0 * kz3 + kz4);
    state.step += 1;

    // RK4 cost: 4 vector-field evaluations + one combine. Each f() is
    // ~10 f64 ops; combine is ~12 ops. Conservative: ~52 f64 ops ≈ 52·64
    // bit-ops per step.
    52 * 64
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn one_step_moves_state() {
        let mut s = LorenzState::default();
        let p = LorenzParams::default();
        let s0 = s;
        let bits = lorenz_rk4_step(&mut s, &p);
        assert_eq!(s.step, 1);
        assert!(bits > 0);
        let moved = (s.x - s0.x).abs() + (s.y - s0.y).abs() + (s.z - s0.z).abs();
        assert!(moved > 0.0, "Lorenz step should move the state");
    }

    #[test]
    fn long_run_stays_on_attractor() {
        // Famous bound: the Lorenz '63 trajectory stays inside a finite ball.
        // 10 000 steps × dt 0.005 = 50 s simulated. State should remain bounded.
        let mut s = LorenzState::default();
        let p = LorenzParams::default();
        let mut peak = 0.0f64;
        for _ in 0..10_000 {
            lorenz_rk4_step(&mut s, &p);
            let r = (s.x * s.x + s.y * s.y + s.z * s.z).sqrt();
            if r > peak {
                peak = r;
            }
        }
        assert!(peak.is_finite(), "trajectory must not diverge");
        assert!(peak < 100.0, "ought to stay inside the classical ball, got {peak}");
    }

    #[test]
    fn fixed_point_at_origin_is_unstable() {
        // Start at the origin (a fixed point) with a tiny perturbation —
        // should grow because the origin is a saddle for ρ > 1.
        let mut s = LorenzState {
            x: 1e-6,
            y: 0.0,
            z: 0.0,
            step: 0,
        };
        let p = LorenzParams::default();
        let r0 = (s.x * s.x + s.y * s.y + s.z * s.z).sqrt();
        for _ in 0..1_000 {
            lorenz_rk4_step(&mut s, &p);
        }
        let r1 = (s.x * s.x + s.y * s.y + s.z * s.z).sqrt();
        assert!(r1 > r0 * 100.0, "origin should be unstable");
    }
}

This is the exact Rust file compiled into the WASM module the page just loaded. Every line that runs on your device is here. The receipt is a function of this code, not a bespoke benchmark.