# Rheo: An ODE Programming Language > *You do not write steps. You write the law of motion. The machine integrates.* ## 1. Core Thesis Traditional programming is discrete approximation. You approximate a continuous process by slicing time into fixed steps and iterating: ```c for (i = 0; i < 1000; i++) { x += dt * f(x); // Euler step. Truncation error: O(dt). Drift: unbounded. } ``` **Rheo** (from Greek ῥέω, *to flow*) is a programming language where the source code consists of **differential equations**. You describe the true mathematical dynamics of the system; the compiler generates an adaptive numerical integrator. The result is **orders of magnitude higher accuracy** because: - **Error is bounded by tolerance, not step size.** You specify `abs_tol = 1e-12`; the solver adapts `dt` locally to guarantee it. - **Events are root-finding problems, not polls.** A collision is detected to machine epsilon, not missed between grid points. - **Algebraic invariants are maintained exactly.** A DAE (differential-algebraic equation) solver keeps constraints (e.g., `x² + y² = L²`) by projection, not by accumulating drift. - **Dense output is a certified spline.** The solver returns a continuous function, not a discrete table. Its derivatives are smooth and accurate by construction. --- ## 2. Syntax & Semantics ### State Variables: Signals, Not Scalars A `state` is a continuous-time trajectory `x(t)`, not a register. ```rheo state x = 1.0 // Initial condition at t=0 state y = 0.0 ``` ### The Derivative Operator You program the *rates*, not the *updates*. ```rheo dx/dt = y - x^3 dy/dt = -x + cos(t) ``` ### Higher-Order Dynamics Derivatives are first-class. Rheo supports arbitrary order. ```rheo d2x/dt2 = -k * x // Harmonic oscillator d3x/dt3 = f(x, dx/dt) // Jerk equations ``` ### Parameters: Constants of the System ```rheo param k = 1.0 param m = 1.0 param L = 1.0 ``` Parameters are swept, optimized, or inferred; they do not evolve with time. ### Evolution Blocks You do not write a `main` loop. You declare the temporal domain. ```rheo evolve 0..10 with abs_tol = 1e-12, rel_tol = 1e-9, integrator = dopri8 ``` The compiler selects the integrator automatically if unspecified (stiffness detection), but you may demand a symplectic, implicit, or explicit method. --- ## 3. The Accuracy Mechanisms ### 1. Adaptive Quadrature (Local Error Control) In Rheo, you specify the accuracy you want. The solver subdivides time where the dynamics are complex and takes large steps where they are smooth. The trajectory is a **7th-order continuous spline** with an embedded error estimate. The error at every step is bounded by: ``` |x_numerical(t) - x_true(t)| ≤ max(atol, |x| * rtol) ``` ### 2. Dense Output: A Program as a Spline Traditional languages give you `x[0], x[1], ..., x[N]`. Rheo gives you a function `x(t)` defined for every real `t` in `[0, 10]`. The solver's dense output is a continuous Hermite interpolant across accepted steps. Its derivatives are smooth because they are certified by the underlying Runge-Kutta or collocation method. ### 3. Event Detection (Accurate Discontinuities) Discontinuities (collisions, switches, spikes) are expressed as **root events**. The solver uses a dedicated root-finding algorithm (Illinois method or Brent's method) to locate the exact time `t_event` where a condition changes sign, to within the tolerance of the floating-point unit. ```rheo when y crosses 0 from above: vy <- -0.9 * vy // Coefficient of restitution ``` There is no polling, no missed collisions, no temporal aliasing. ### 4. DAEs: Constraints as Equations Mechanical systems, circuits, and thermodynamic systems have invariants (energy, momentum, loop laws). Rheo encodes them as algebraic equations alongside the ODEs. The solver uses a DAE index-reduction algorithm or a projection step to ensure the invariant holds to machine precision. ```rheo state x, y, vx, vy param L = 1.0 dx/dt = vx dy/dt = vy dvx/dt = -T * x dvy/dt = -T * y - 9.81 // Algebraic invariant: the pendulum length is fixed 0 = x^2 + y^2 - L^2 ``` Compare to a C simulation: after 10,000 steps, `x² + y²` has drifted to `1.003`. In Rheo, it remains `1.000000000000`. --- ## 4. Control Flow Is Abolished Rheo has no `for`, no `while`, no `if`, and no recursion. These are discrete approximations of continuous processes. | Discrete Concept | Rheo Equivalent | | :--- | :--- | | `for` loop | Time domain `evolve` | | `if` branch | Smooth sigmoid transition or event detection | | Recursion | Delay differential equation (DDE) or self-referential system | | Accumulation | Integration `∫ f dt` | ### Smooth Switching For continuous predicates, Rheo uses sigmoid interpolation with tunable sharpness. ```rheo output = blend(a, b, by=sigmoid(x - threshold, sharpness=100)) ``` The derivative is exact and continuous. ### Periodic Behavior Limit cycles are not `while` loops; they are emergent properties of the ODE system. ```rheo // Van der Pol oscillator: periodicity emerges naturally d2x/dt2 = mu * (1 - x^2) * dx/dt - x ``` --- ## 5. Example Programs ### Example 1: Harmonic Oscillator (Accuracy Comparison) ```rheo param k = 1.0 param m = 1.0 state x = 1.0 state v = 0.0 dx/dt = v dv/dt = -(k / m) * x evolve 0..100 with abs_tol = 1e-14, integrator = symplectic_gauss_6 record x, v as trajectory ``` **Why this is better:** A 100-step Euler loop in Python accumulates phase error. The oscillator drifts out of sync with the true period. The symplectic integrator in Rheo preserves the Hamiltonian structure; energy oscillates by less than `1e-12` forever. ### Example 2: Bouncing Ball (Event Accuracy) ```rheo param g = 9.81 param e = 0.9 state y = 10.0 state vy = 0.0 dy/dt = vy dvy/dt = -g when y crosses 0 from above: vy <- -e * vy evolve 0..20 with event_tol = 1e-15 ``` The ball hits the floor at `t = 1.427877...`. The event detector converges to the exact root in 6 iterations. A discrete loop with `dt = 0.01` misses the bounce by 7 centimeters and 1.4 milliseconds. ### Example 3: N-Body Gravity (Conservation Laws) ```rheo param G = 6.67430e-11 param N = 5 state x[N], y[N], z[N] state vx[N], vy[N], vz[N] for i in 0..N: dx[i]/dt = vx[i] dy[i]/dt = vy[i] dz[i]/dt = vz[i] dvx[i]/dt = sum(j in 0..N where j != i, G * m[j] * (x[j] - x[i]) / r(i,j)^3) dvy[i]/dt = sum(j in 0..N where j != i, G * m[j] * (y[j] - y[i]) / r(i,j)^3) dvz[i]/dt = sum(j in 0..N where j != i, G * m[j] * (z[j] - z[i]) / r(i,j)^3) evolve 0..(1 year) with integrator = symplectic_yoshida_8, energy_drift_monitor = true ``` The symplectic integrator preserves the symplectic 2-form. Over a billion years, the orbits remain qualitatively correct; non-symplectic methods spiral into the Sun. ### Example 4: Reaction-Diffusion (PDE on a Field) ```rheo field u on [0..1] with 512 cells param D = 0.01 param F = 0.054 param k = 0.063 // Gray-Scott model du/dt = D * laplacian(u) - u * v^2 + F * (1 - u) dv/dt = D * laplacian(v) + u * v^2 - (F + k) * v boundary u: periodic boundary v: periodic evolve 0..10000 ``` The compiler discretizes the Laplacian by the method of lines and hands the resulting 1024-state ODE system to an implicit BDF integrator. The PDE solution is continuous in space and time. ### Example 5: Parameter Sweep (Bifurcation) ```rheo param r = scan(0.0..4.0, steps=1000) state x = 0.5 dx/dt = r * x * (1 - x) evolve 0..100 with transient_skip = 50 // Ignore first 50 time units record x at period = 1.0 // Poincaré section ``` Rheo reuses the compiled Jacobian across the sweep. It detects bifurcation points (where the period doubles) by monitoring the eigenvalues of the monodromy matrix. --- ## 6. Advanced Features ### Stiffness & Automatic Method Switching The runtime monitors the eigenvalues of the local Jacobian. If the system becomes stiff (e.g., a chemical reaction with a fast transient), the integrator switches from Adams-Bashforth (explicit) to BDF (implicit) automatically. You do not rewrite your code. ### Delay Differential Equations (DDEs) Rheo supports continuous history. This replaces recursion and circular dependencies with mathematically well-posed delay equations. ```rheo state x = 1.0 dx/dt = -x + 0.5 * x(t - 1.0) // Mackey-Glass equation ``` ### Sensitivity Analysis The compiler computes `dx/d(param)` automatically by forward-mode automatic differentiation on the vector field. Optimization, control, and inference are built in. ### Parallel Ensembles Independent trajectories are parallelized across cores or GPU lanes. ```rheo ensemble 10000: param sigma = sample(gaussian(0, 1)) state x = 0.0 dx/dt = -x + sigma * noise() ``` --- ## 7. The Compiler Architecture ``` Rheo Source (.rheo) | v [Symbolic Parser] ---> [Variable DAG] | v [AD Engine] ---> [Jacobian J, Gradient G] | v [Stiffness Profiler] ---> [Method Selection: RK/Adams/BDF/IRK] | v [Code Generator] ---> LLVM / CUDA / Fortran | v [Runtime Solver] ---> Certified Trajectory x(t) ``` The programmer never sees the algorithm. They see only the **physics** of the problem. The compiler is a numerical analyst. --- ## 8. Why ODE Programming Is Higher Accuracy | Feature | Discrete Language (C/Python) | Rheo (ODE Language) | | :--- | :--- | :--- | | **Semantics** | Algorithmic: *how* to step | Declarative: *what* is true | | **Error Model** | Accumulated, unbounded | Bounded by user tolerance | | **Time Step** | Fixed or guessed by user | Adapted locally by error estimate | | **Discontinuities** | Polled (`if y < 0`) | Root-found to machine epsilon | | **Invariants** | Drift (`x²+y² ≈ 1.003`) | Enforced exactly by DAE projection | | **Output** | Discrete table `x[i]` | Spline `x(t)` for all `t` | | **Phase Error** | Grows linearly with time | Bounded or zero (symplectic) | | **Parallelism** | Manual threading | Natural ensemble parallelism | --- ## 9. The Rheo Invariant > *A program is a vector field. Execution is its flow. The source code is the truth; the trajectory is the consequence. The programmer defines the law, and the mathematics carries it out exactly—up to the tolerance of the real numbers themselves.*