9. Differential Equations, Runge-Kutta (RK4), Non-Equilibrium Thermodynamics, and Wave Superposition

Why the Planet’s Weather and Your Ear’s Cochlea Do Not Use Binary Long Division

Look at how a digital computer simulates a weight bouncing on a spring or a sound wave traveling through air:

1. It chops continuous time $t$ into discrete ticks $\Delta t$ and continuous voltage into 32-bit binary strings (010000010...).
2. If it uses a naïve first-order step ($\mathbf{x}_{n+1} = \mathbf{x}_n + \Delta t \cdot \mathbf{f}(\mathbf{x}_n)$, Forward Euler), something bizarre happens: every single step overshoots the curved trajectory tangentially outward, injecting fake energy into the spring until the simulation explodes to infinity!
3. Even when the digital computer uses a carefully balanced 4th-order integrator (RK4), it must shuttle millions of 64-bit numbers back and forth between DRAM and the ALU ($640\text{ pJ}$ per DRAM access) just to approximate what a physical capacitor and inductor ($L\ddot{q} + R\dot{q} + \frac{1}{C}q = V(t)$) or a continuous acoustic wave does instantly, in zero clock cycles, by obeying the laws of physics.

And look at a second paradox: if you try to train a standard neural network with piecewise-linear ReLU activations ($\mathrm{ReLU}(z) = \max(0, z)$) to solve a physical wave equation $\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$, its second derivative is identically zero everywhere ($\frac{d^2}{dz^2}\mathrm{ReLU}(z) = 0$)! A ReLU network is structurally blind to second-order curvature and wave interference.

Continuous physical and biological substrates—including Prigogine’s Brusselator chemical oscillator, Lovelock’s Daisyworld biosphere ODEs, RLC-mechanical analog computers, Liquid Time-Constant (LTC) neural ODEs, and SIREN / Wave-Physics PDEs—compute directly through coupled differential equations, non-equilibrium entropy export, and wave superposition.

9.1 Dynamical Systems and Ordinary Differential Equations (ODEs)

From Rates of Change to State-Space Vectors

Whenever a physical, biological, or economic state $\mathbf{x}(t) = (x_1(t), \dots, x_d(t))^\top \in \mathbb{R}^d$ evolves continuously in time, its instantaneous velocity vector is the time derivative $\dot{\mathbf{x}}(t) = \frac{d\mathbf{x}}{dt}$. An autonomous system of first-order Ordinary Differential Equations (ODEs) specifies this velocity as a function of the current state:

$$\frac{d\mathbf{x}}{dt} = \mathbf{f}(\mathbf{x})$$

How to Convert Any Second-Order Physical Law into a First-Order State-Space System

Many laws of physics—such as Newton’s Second Law $m\ddot{x} + c\dot{x} + kx = F(t)$ or Kirchhoff’s RLC loop equation $L\ddot{q} + R\dot{q} + \frac{1}{C}q = V(t)$—involve a second derivative (acceleration $\ddot{x} = \frac{d^2 x}{dt^2}$).

How do we solve a second-order ODE using a standard first-order solver? By introducing velocity $v = \dot{x}$ as a second state variable! Defining the 2D state vector $\mathbf{z}(t) = \begin{bmatrix} x(t) \ v(t) \end{bmatrix}$ transforms one second-order ODE into two coupled first-order ODEs:

$$\frac{d}{dt}\begin{bmatrix} x \ v \end{bmatrix} = \begin{bmatrix} v \ \frac{1}{m}\left(F(t) - c v - k x\right) \end{bmatrix}$$

Look at the exact correspondence between a mechanical mass-spring-damper and an electrical RLC circuit (used in Chapter 11):

Because both physical systems obey the exact same differential equation, measuring the voltage across the capacitor in an RLC circuit solves the mechanical differential equation in real time at the speed of electromagnetic propagation!

9.2 Fixed Points, Jacobian Stability, and Hopf Bifurcations (Limit Cycles)

Equilibrium (Fixed) Points and Linearization

A state $\mathbf{x}^\ast \in \mathbb{R}^d$ where all rates of change vanish simultaneously,

$$\mathbf{f}(\mathbf{x}^\ast) = \mathbf{0}$$

is called a fixed point (equilibrium) of the dynamical system. If the system is nudged slightly away from $\mathbf{x}^\ast$ by a small perturbation $\boldsymbol{\xi}(t) = \mathbf{x}(t) - \mathbf{x}^\ast$, does it return to $\mathbf{x}^\ast$ or fly away?

Expanding $\mathbf{f}(\mathbf{x}^\ast + \boldsymbol{\xi})$ in a first-order Taylor series around $\mathbf{x}^\ast$ gives the linearized system:

$$\frac{d\boldsymbol{\xi}}{dt} \approx J(\mathbf{x}^\ast),\boldsymbol{\xi}, \qquad \text{where } J_{ij}(\mathbf{x}^\ast) = \frac{\partial f_i}{\partial x_j}(\mathbf{x}^\ast)$$

is the $d \times d$ Jacobian Matrix at $\mathbf{x}^\ast$. For a 2D system ($d = 2$), the eigenvalues $\lambda_{1,2}$ of $J$ are determined completely by its trace $\tau = \mathrm{tr}(J) = J_{11} + J_{22}$ and determinant $\Delta = \det(J) = J_{11}J_{22} - J_{12}J_{21}$ via the characteristic quadratic $\lambda^2 - \tau\lambda + \Delta = 0$:

$$\lambda_{1,2} = \frac{\tau \pm \sqrt{\tau^2 - 4\Delta}}{2}$$

When $\Delta > 0$ and $\tau^2 < 4\Delta$, the eigenvalues are a complex-conjugate pair $\lambda_{1,2} = \alpha \pm i\omega$ with real part $\alpha = \frac{\tau}{2}$:

1. When $\tau = \mathrm{tr}(J) < 0$ ($\alpha < 0$): Perturbations decay as $e^{\alpha t}\cos(\omega t)$, spiraling inward to a stable equilibrium (like a damped pendulum coming to rest).
2. When $\tau = \mathrm{tr}(J) > 0$ ($\alpha > 0$): Perturbations grow as $e^{\alpha t}\cos(\omega t)$, spiraling outward away from the unstable fixed point! If nonlinear terms bound the system at large amplitudes, the outward spiral settles onto a self-sustained periodic orbit called a Limit Cycle.
3. Hopf Bifurcation ($\tau = \mathrm{tr}(J) = 0$): The critical threshold where $\mathrm{tr}(J)$ crosses from negative to positive—transforming a dead equilibrium into a living rhythmic clock (such as the heartbeat, circadian rhythm, or Prigogine’s Brusselator chemical oscillator at $B > 1 + A^2$ in Chapter 10)—is called a Hopf Bifurcation!

9.3 Numerical Integration: Why Forward Euler Fails and How RK4 Works

When an ODE $\dot{\mathbf{x}} = \mathbf{f}(t, \mathbf{x})$ has no closed-form paper solution, how do we step from $\mathbf{x}_n \approx \mathbf{x}(t_n)$ to $\mathbf{x}_{n+1} \approx \mathbf{x}(t_n + h)$ over a small step size $h = \Delta t$?

Why Forward Euler Artificially Gains Energy on Oscillators

The simplest scheme is Forward Euler, which follows the tangent slope at the start of the interval:

$$\mathbf{x}_{n+1}^{\text{Euler}} = \mathbf{x}_n + h,\mathbf{f}(t_n, \mathbf{x}_n)$$

Consider the pure harmonic oscillator $\dot{x} = v, ; \dot{v} = -x$, whose true trajectories are circles $x(t)^2 + v(t)^2 = R^2$ of constant energy. Look at what one step of Forward Euler does to the squared radius:

$$x_{n+1}^2 + v_{n+1}^2 = (x_n + h v_n)^2 + (v_n - h x_n)^2 = (1 + h^2)(x_n^2 + v_n^2) > x_n^2 + v_n^2$$

Every single Forward Euler step multiplies the energy by $(1 + h^2) > 1$, causing the trajectory to spiral outward (Panel 1 of Figure 9.1) with global error $\mathcal{O}(h)$!

The Classical 4th-Order Runge-Kutta (RK4) Algorithm

To cancel out first-, second-, and third-order curvature errors without needing higher derivatives of $\mathbf{f}$, Carl Runge and Wilhelm Kutta designed a 4-stage probe across the interval $[t_n, t_n + h]$:

1. Stage 1 ($k_1$, slope at the start): $$\mathbf{k}_1 = \mathbf{f}(t_n, ; \mathbf{x}_n)$$
2. Stage 2 ($k_2$, slope at the midpoint predicted by $k_1$): $$\mathbf{k}_2 = \mathbf{f}!\left(t_n + \frac{h}{2}, ; \mathbf{x}_n + \frac{h}{2}\mathbf{k}_1\right)$$
3. Stage 3 ($k_3$, corrected slope at the midpoint predicted by $k_2$): $$\mathbf{k}_3 = \mathbf{f}!\left(t_n + \frac{h}{2}, ; \mathbf{x}_n + \frac{h}{2}\mathbf{k}_2\right)$$
4. Stage 4 ($k_4$, slope at the far endpoint predicted by $k_3$): $$\mathbf{k}_4 = \mathbf{f}(t_n + h, ; \mathbf{x}_n + h\mathbf{k}_3)$$

Combining these four slopes with Simpson’s $1\text{–}2\text{–}2\text{–}1$ quadrature weights gives the RK4 update:

$$\mathbf{x}_{n+1} = \mathbf{x}_n + \frac{h}{6}\left(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4\right)$$

Because the midpoint slopes $\mathbf{k}_2, \mathbf{k}_3$ receive double weight and cancel the tangential overshoot of $\mathbf{k}_1$, RK4 achieves local truncation error $\mathcal{O}(h^5)$ and global error $\mathcal{O}(h^4)$—cutting the step size $h$ by $10\times$ shrinks the integration error by $10\text{,}000\times$!

9.4 Open Non-Equilibrium Thermodynamics ($dS = d_e S + d_i S$)

How can a living cell, a biosphere, or a chemical clock spontaneously organize itself into rhythmic patterns without violating the Second Law of Thermodynamics?

Ilya Prigogine (1977 Nobel Prize) resolved this apparent paradox by dividing the total entropy change $dS$ of an open system (a system that exchanges energy and matter with its surroundings) into two distinct parts:

$$dS = d_e S + d_i S$$

1. Internal Entropy Production ($d_i S \ge 0$): Inside the system, irreversible chemical reactions, diffusion, and friction always produce non-negative entropy ($d_i S / dt \ge 0$, which is strictly required by the Second Law). For a chemical reaction network with forward flux $J_\rho^+$ and reverse flux $J_\rho^-$, Schnakenberg’s formula gives: $$\frac{d_i S}{dt} = k_B \sum_\rho \left(J_\rho^+ - J_\rho^-\right) \ln!\left(\frac{J_\rho^+}{J_\rho^-}\right) \ge 0$$ (Since $(a - b)\ln(a/b) \ge 0$ for all $a, b > 0$, this is zero if and only if every reaction is in dead detailed balance $J_\rho^+ = J_\rho^-$!)
2. External Entropy Exchange ($d_e S$, unrestricted in sign): An open system absorbs high-grade free energy (such as sunlight or chemical feed $A, B$) and exports low-grade heat and waste to its environment, making $d_e S < 0$.

Whenever the outward export of entropy exceeds internal production ($d_e S < -d_i S \le 0$), the net entropy of the system decreases ($dS = d_e S + d_i S < 0$) or settles into a non-equilibrium steady state ($dS = 0$ with $d_i S = -d_e S > 0$), sustaining Dissipative Structures far from thermal death!

9.5 Wave Superposition, Interference, and the Wave / Helmholtz PDEs

Amplitude, Frequency, Phase, and Superposition

A continuous harmonic wave in space $x$ and time $t$ is written as:

$$u(x, t) = A \sin(k x - \omega t + \phi)$$

where $A \ge 0$ is the amplitude (peak height; energy scales as $A^2$), $k = 2\pi / \lambda$ is the spatial wavenumber, $\omega = 2\pi f$ is the angular frequency, and $\phi \in [0, 2\pi)$ is the phase angle.

When two waves $u_1 = A\sin(\theta)$ and $u_2 = A\sin(\theta + \Delta\phi)$ of the same frequency meet at a point in space, the trigonometric identity $\sin(\theta) + \sin(\theta + \Delta\phi) = 2\cos(\Delta\phi / 2)\sin(\theta + \Delta\phi / 2)$ produces wave interference:

• Constructive Interference ($\Delta\phi = 0$, in phase): $\cos(0) = 1$, so the combined amplitude doubles to $2A$ (and intensity quadruples to $4A^2$).
• Destructive Interference ($\Delta\phi = \pi = 180^\circ$, out of phase): $\cos(\pi/2) = 0$, so the two waves annihilate each other completely ($u_1 + u_2 = 0$)!

Notice that a single point in a wave medium simultaneously processes both magnitude ($A$) and relative timing/sign ($\phi$) without a single digital gate.

Why Sinusoids (SIREN) Solve Second-Order Wave and Helmholtz PDEs While ReLU Fails

In physics, wave propagation through a medium with spatially varying wave speed $c(\mathbf{x})$ and damping $\gamma(\mathbf{x})$ is governed by the Second-Order Wave Partial Differential Equation (PDE):

$$\frac{\partial^2 u}{\partial t^2} = c(\mathbf{x})^2 \nabla^2 u - \gamma(\mathbf{x})\frac{\partial u}{\partial t} + s(\mathbf{x}, t)$$

For a steady-state oscillation $u(x, t) = \psi(x)e^{-i\omega t}$, this reduces to the spatial Helmholtz Equation:

$$\frac{d^2 \psi(x)}{dx^2} + \omega_0^2 \psi(x) = s(x)$$

Now compare what happens when we differentiate two neural activation functions twice:

1. Piecewise-Linear ReLU ($\sigma(x) = \max(0, x)$): $$\sigma’(x) = \mathbf{1}{x > 0}, \qquad \sigma’‘(x) = \mathbf{0} \quad (\forall x \ne 0)$$ Because $\sigma’’(x) = 0$ everywhere, a ReLU network has zero second derivative and cannot represent the $\frac{d^2\psi}{dx^2}$ Laplacian term of the Helmholtz/Wave PDE!
2. Sinusoidal Activation (SIREN, Sitzmann et al., $\sigma(x) = \sin(\omega_0 x + b)$): $$\sigma’(x) = \omega_0 \cos(\omega_0 x + b), \qquad \sigma’’(x) = -\omega_0^2 \sin(\omega_0 x + b) = -\omega_0^2 \sigma(x)$$ The derivative of a SIREN network is another SIREN network (phase-shifted by $\pi/2$), and its second derivative is proportional to itself, matching the Helmholtz operator $\frac{d^2\psi}{dx^2} + \omega_0^2 \psi$ natively!

Computational Audit in Rust

All dynamical, thermodynamic, and wave-equation solvers built on these foundations live in crates/stochology-lenses/src/relational.rs and crates/stochology-lenses/src/analog.rs:

• PrigogineBrusselatorAudit (relational.rs): Integrates the 2D Brusselator ODE system via RK4 and evaluates Schnakenberg’s internal entropy production $d_i S / dt > 0$.
• LovelockDaisyworldAudit (relational.rs): Integrates the coupled Daisyworld biosphere-albedo ODEs via RK4.
• ContinuousOdeLtcIsomorphismAudit (analog.rs): Integrates the 2nd-order RLC / Mass-Spring-Damper ODE and Liquid Time-Constant (LTC) neural ODE via RK4.
• WavePhysicsSirenRnnAudit (analog.rs): Compares SIREN $\sin(\omega_0 x)$ against ReLU on the 2nd-order Helmholtz PDE and integrates the 2nd-order damped nonlinear Wave PDE.

Run Gate G8 and Lab 3 to inspect these continuous systems:

cargo run --release --bin labs -- --lenses
cargo run --release --bin gates

Exercises

1. One Step of RK4 vs. Forward Euler by Hand: Consider the exponential decay ODE $\dot{x} = -x$ with initial condition $x(0) = 1.0$ and step size $h = 0.5$ (exact solution at $t = 0.5$ is $e^{-0.5} \approx 0.60653$).
   • Compute the Forward Euler approximation $x_1^{\text{Euler}} = x_0 + h(-x_0)$.
   • Compute the four RK4 slopes $k_1, k_2, k_3, k_4$ and verify that $x_1^{\text{RK4}} = 1 + \frac{0.5}{6}(-1 + 2(-0.75) + 2(-0.8125) - 0.59375) \approx 0.60677$, reducing the error of Forward Euler ($0.1065$) by more than $400\times$ in a single step!
2. Hopf Bifurcation of the Brusselator: For the Brusselator reactions $\dot{X} = A - (B + 1)X + X^2 Y$ and $\dot{Y} = BX - X^2 Y$, verify that $(X^\ast, Y^\ast) = (A, B/A)$ is the unique fixed point ($\dot{X} = \dot{Y} = 0$). Compute the $2 \times 2$ Jacobian matrix $J(X^\ast, Y^\ast) = \begin{bmatrix} B - 1 & A^2 \ -B & -A^2 \end{bmatrix}$ and show that $\det(J) = A^2 > 0$ and $\mathrm{tr}(J) = B - 1 - A^2$, proving that the fixed point loses stability and enters a limit cycle precisely when $B > 1 + A^2$.
3. Constructive vs. Destructive Interference: Compute the superposition $u(t) = 3\sin(\omega t) + 3\sin(\omega t + \Delta\phi)$ at $\omega t = \pi/2$ for (a) $\Delta\phi = 0$, (b) $\Delta\phi = \pi/2$, and (c) $\Delta\phi = \pi$.