11 — Dynamical Systems & Chaos Theory

A pure mathematics showcase — no downloaded dataset, the "real data" here is two of the canonical objects of dynamical systems theory: the logistic map and the Lorenz system. Every numerical claim below is cross-checked against known analytical facts (the logistic map's period-doubling cascade, the Lorenz system's critical bifurcation parameter) or against numpy directly.

This notebook also doubles as a real end-to-end regression check for two engine bugs found and fixed this session (engine v0.1.83, linaldb PyPI package 0.1.10), found via a deep audit of docs/DSL_REFERENCE.md against a live build of the engine:

1. Bool column predicates (WHERE <bool_col> = 1, and a bare WHERE <bool_col>) used to silently match zero rows instead of comparing correctly.
2. In-place TRANSFORM (no INTO) used to corrupt a dataset's schema when the projection changed the column set, breaking every later read.

Both get exercised directly below, in a real pipeline — not a synthetic repro.

import shutil
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

import linaldb

print("linaldb version:", linaldb.__version__)

DATA_DIR = Path("data_dynamical_systems")
shutil.rmtree(DATA_DIR, ignore_errors=True)
db = linaldb.Db(data_dir=str(DATA_DIR))

linaldb version: 0.1.10

Part 1 — The logistic map: bifurcation and the onset of chaos

$$x_{n+1} = r \, x_n (1 - x_n)$$

For $r$ swept from 2.5 to 4.0, this single line produces one of the most famous bifurcation diagrams in all of mathematics: a cascade of period-doublings (period 1 → 2 → 4 → 8 → ...) accumulating at the Feigenbaum point $r_\infty \approx 3.56995$, beyond which the map is chaotic almost everywhere.

The DSL has no loop construct, so the iteration itself runs in numpy (1000+ steps per $r$ value, exactly the kind of tight numerical loop it's built for) — linaldb's job starts once we have real per-$r$ results to store, query, and classify.

r_values = np.linspace(2.5, 4.0, 300)
n_transient, n_sample = 500, 200

lyapunov = np.zeros_like(r_values)
final_x = np.zeros((len(r_values), n_sample))

for i, r in enumerate(r_values):
    x = 0.5
    lyap_sum = 0.0
    for n in range(n_transient + n_sample):
        deriv = abs(r - 2 * r * x)
        if n >= n_transient:
            lyap_sum += np.log(deriv) if deriv > 0 else -50.0
            final_x[i, n - n_transient] = x
        x = r * x * (1 - x)
    lyapunov[i] = lyap_sum / n_sample

is_chaotic = lyapunov > 0
print(f"{is_chaotic.sum()} / {len(r_values)} sampled r values are chaotic "
      f"(positive Lyapunov exponent)")

77 / 300 sampled r values are chaotic (positive Lyapunov exponent)

Sanity check against known analytical facts, before trusting any of this further:

idx_period2 = np.argmin(np.abs(r_values - 3.2))
idx_chaotic = np.argmin(np.abs(r_values - 3.9))
onset_r = r_values[np.where(lyapunov > 0)[0][0]]

print(f"r=3.2 (known period-2 window): lyapunov = {lyapunov[idx_period2]:+.4f}  (expect < 0)")
print(f"r=3.9 (known chaotic region):  lyapunov = {lyapunov[idx_chaotic]:+.4f}  (expect > 0)")
print(f"first r with positive lyapunov: {onset_r:.4f}  (expect close to the Feigenbaum point ~3.5700)")

assert lyapunov[idx_period2] < 0
assert lyapunov[idx_chaotic] > 0
assert abs(onset_r - 3.56995) < 0.02

r=3.2 (known period-2 window): lyapunov = -0.9496  (expect < 0)
r=3.9 (known chaotic region):  lyapunov = +0.4755  (expect > 0)
first r with positive lyapunov: 3.5736  (expect close to the Feigenbaum point ~3.5700)

Into a real hybrid dataset

Scalar columns (r, lyapunov_exponent, is_chaotic) alongside a Vector(200) column holding each run's sampled long-term iterates — exactly the kind of hybrid table linaldb is built around.

db.execute(
    "DATASET logistic_runs COLUMNS (r: Double, lyapunov_exponent: Double, "
    "is_chaotic: Bool, sample_points: Vector(200))"
)
for i, r in enumerate(r_values):
    flag = "true" if lyapunov[i] > 0 else "false"
    vec = "[" + ", ".join(f"{v:.6f}" for v in final_x[i]) + "]"
    db.execute(
        f"INSERT INTO logistic_runs VALUES ({r:.6f}, {lyapunov[i]:.6f}, {flag}, {vec})"
    )

print(f"inserted {len(r_values)} logistic-map runs")

inserted 300 logistic-map runs

Regression check 1 — Bool column predicates

WHERE is_chaotic = 1, a bare WHERE is_chaotic, and WHERE is_chaotic = true must all agree — before the fix, the first two silently returned zero rows.

r_eq1 = db.execute("SELECT COUNT(*) AS n FROM logistic_runs WHERE is_chaotic = 1")
r_bare = db.execute("SELECT COUNT(*) AS n FROM logistic_runs WHERE is_chaotic")
r_true = db.execute("SELECT COUNT(*) AS n FROM logistic_runs WHERE is_chaotic = true")

print("WHERE is_chaotic = 1   ->", r_eq1.rows)
print("WHERE is_chaotic       ->", r_bare.rows)
print("WHERE is_chaotic = true ->", r_true.rows)

assert r_eq1.rows == r_bare.rows == r_true.rows == [[int(is_chaotic.sum())]]
print("\nAll three forms agree, and match the real chaotic-run count.")

WHERE is_chaotic = 1   -> [[77]]
WHERE is_chaotic       -> [[77]]
WHERE is_chaotic = true -> [[77]]

All three forms agree, and match the real chaotic-run count.

grouped = db.execute(
    "SELECT is_chaotic, AVG(lyapunov_exponent) AS avg_lyapunov, COUNT(*) AS n "
    "FROM logistic_runs GROUP BY is_chaotic"
)
for row in grouped.rows:
    label = "chaotic" if row[0] else "stable/periodic"
    print(f"{label:16s}  avg lyapunov = {row[1]:+.4f}   n = {row[2]}")

stable/periodic   avg lyapunov = -0.3718   n = 223
chaotic           avg lyapunov = +0.4018   n = 77

The bifurcation diagram

Every sampled $r$ contributes its 200 long-term iterates as one point cloud column — the classic period-doubling cascade, colored by whether linaldb's own Lyapunov-sign classification (computed above, stored in is_chaotic) calls it chaotic.

fig, ax = plt.subplots(figsize=(10, 6))
colors = np.where(is_chaotic, "#d62728", "#1f77b4")
for i, r in enumerate(r_values):
    ax.plot([r] * n_sample, final_x[i], ",", color=colors[i], alpha=0.35)
ax.axvline(3.56995, color="black", linestyle="--", linewidth=1, label="Feigenbaum point")
ax.set_xlabel("r")
ax.set_ylabel("long-term x")
ax.set_title("Logistic map bifurcation diagram (red = chaotic, blue = stable/periodic)")
ax.legend()
plt.tight_layout()
plt.show()

Regression check 2 — in-place TRANSFORM

Drop the heavy Vector(200) column now that the diagram is drawn — a real "clean up after computing" step. Without INTO, TRANSFORM overwrites the dataset in place; before the fix, changing the column set here left the dataset's schema stale while its rows moved on, corrupting it (any later read failed with a value-count mismatch).

print("schema before:")
print(db.execute("SHOW SCHEMA logistic_runs"))

db.execute(
    "TRANSFORM logistic_runs SELECT r, lyapunov_exponent, is_chaotic "
    "WHERE lyapunov_exponent > -2.0"
)

print("\nschema after in-place TRANSFORM:")
print(db.execute("SHOW SCHEMA logistic_runs"))

post = db.execute("SELECT * FROM logistic_runs WHERE is_chaotic ORDER BY r LIMIT 3")
print("\nreadable immediately afterward:", post.columns)
for row in post.rows:
    print(" ", row)
assert post.columns == ["r", "lyapunov_exponent", "is_chaotic"]

schema before:
Schema for dataset 'logistic_runs' (Legacy):
Field                Type                 Nullable  
----------------------------------------------------
r                    Float64              false     
lyapunov_exponent    Float64              false     
is_chaotic           Bool                 false     
sample_points        Vector(200)          false     


schema after in-place TRANSFORM:
Schema for dataset 'logistic_runs' (Legacy):
Field                Type                 Nullable  
----------------------------------------------------
r                    Float64              false     
lyapunov_exponent    Float64              false     
is_chaotic           Bool                 false     


readable immediately afterward: ['r', 'lyapunov_exponent', 'is_chaotic']
  [3.573579, 0.077525, True]
  [3.578595, 0.10068, True]
  [3.583612, 0.036414, True]

Part 2 — The Lorenz system: sensitive dependence on initial conditions

$$\dot{x} = \sigma(y - x), \quad \dot{y} = x(\rho - z) - y, \quad \dot{z} = xy - \beta z$$

The textbook chaotic parameters are $\sigma=10$, $\rho=28$, $\beta=8/3$. At $\rho=10$ the system instead settles onto a stable fixed point. We integrate both regimes (RK4, numpy again — the DSL has no ODE solver) from two initial conditions separated by only $10^{-4}$, and let linaldb's own DISTANCE function track how that tiny gap evolves.

sigma, beta = 10.0, 8.0 / 3.0

def lorenz_deriv(state, rho):
    x, y, z = state
    return np.array([sigma * (y - x), x * (rho - z) - y, x * y - beta * z])

def rk4_step(state, dt, rho):
    k1 = lorenz_deriv(state, rho)
    k2 = lorenz_deriv(state + dt / 2 * k1, rho)
    k3 = lorenz_deriv(state + dt / 2 * k2, rho)
    k4 = lorenz_deriv(state + dt * k3, rho)
    return state + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

def integrate(state0, steps, dt, rho):
    traj = np.zeros((steps, 3))
    state = np.array(state0, dtype=float)
    for i in range(steps):
        state = rk4_step(state, dt, rho)
        traj[i] = state
    return traj

steps, dt, sample_every = 2000, 0.01, 100
# Vector/Matrix/Tensor are f32-only by design in linaldb -- a separation near
# f32 epsilon (~1.19e-7) between ~O(1) coordinates would catastrophically
# cancel to exactly 0.0 in DISTANCE. 1e-4 stays well within real f32
# precision while the exponential-divergence story is unchanged.
state_a, state_b = [1.0, 1.0, 1.0], [1.0 + 1e-4, 1.0, 1.0]

db.execute("DATASET trajectory_pairs COLUMNS (rho: Double, step: Int, a: Vector(3), b: Vector(3))")
trajectories = {}
for rho in (10.0, 28.0):
    traj_a = integrate(state_a, steps, dt, rho)
    traj_b = integrate(state_b, steps, dt, rho)
    trajectories[rho] = (traj_a, traj_b)
    for step in range(0, steps, sample_every):
        av = "[" + ", ".join(f"{v:.8f}" for v in traj_a[step]) + "]"
        bv = "[" + ", ".join(f"{v:.8f}" for v in traj_b[step]) + "]"
        db.execute(f"INSERT INTO trajectory_pairs VALUES ({rho}, {step}, {av}, {bv})")

print("inserted trajectory samples for rho=10 (stable) and rho=28 (chaotic)")

inserted trajectory samples for rho=10 (stable) and rho=28 (chaotic)

sep_chaotic = db.execute(
    "SELECT step, DISTANCE(a, b) AS sep FROM trajectory_pairs WHERE rho = 28.0 ORDER BY step"
)
sep_stable = db.execute(
    "SELECT step, DISTANCE(a, b) AS sep FROM trajectory_pairs WHERE rho = 10.0 ORDER BY step"
)

print("rho=28 (chaotic)  -- separation grows:")
for row in sep_chaotic.rows[::4]:
    print(f"  step {row[0]:5d}  sep = {row[1]:.6e}")

print("\nrho=10 (stable)   -- separation shrinks:")
for row in sep_stable.rows[::4]:
    print(f"  step {row[0]:5d}  sep = {row[1]:.6e}")

# cross-check linaldb's DISTANCE against numpy at the same step
last_step, last_sep = sep_chaotic.rows[-1]
traj_a28, traj_b28 = trajectories[28.0]
np_sep = np.linalg.norm(traj_a28[last_step] - traj_b28[last_step])
print(f"\ncross-check at step {last_step}: linaldb DISTANCE = {last_sep:.6f}, "
      f"numpy norm = {np_sep:.6f}")
assert abs(last_sep - np_sep) < 1e-3

rho=28 (chaotic)  -- separation grows:
  step     0  sep = 9.530968e-05
  step   400  sep = 1.734377e-04
  step   800  sep = 2.018238e-04
  step  1200  sep = 1.090139e-04
  step  1600  sep = 5.795173e-03

rho=10 (stable)   -- separation shrinks:
  step     0  sep = 9.123656e-05
  step   400  sep = 4.354649e-05
  step   800  sep = 3.568323e-06
  step  1200  sep = 6.743496e-07
  step  1600  sep = 0.000000e+00

cross-check at step 1900: linaldb DISTANCE = 0.139688, numpy norm = 0.139687

fig, ax = plt.subplots(figsize=(9, 5))
steps_c = [row[0] for row in sep_chaotic.rows]
sep_c = [row[1] for row in sep_chaotic.rows]
steps_s = [row[0] for row in sep_stable.rows]
sep_s = [row[1] for row in sep_stable.rows]
ax.semilogy(steps_c, sep_c, label="rho=28 (chaotic)", color="#d62728")
ax.semilogy(steps_s, [max(v, 1e-10) for v in sep_s], label="rho=10 (stable)", color="#1f77b4")
ax.set_xlabel("integration step")
ax.set_ylabel("separation (log scale), via linaldb DISTANCE")
ax.set_title("Sensitive dependence on initial conditions")
ax.legend()
plt.tight_layout()
plt.show()

Part 3 — Attractor shape via EIGENVALUES

EIGENVALUES/EIGEN in linaldb only accept symmetric matrices (a real, documented constraint — the engine has no complex-number Value support, and a general matrix's eigenvalues can be complex). The Lorenz Jacobian is not symmetric, so we reach for a matrix that legitimately is: the covariance matrix of points sampled from the chaotic attractor itself. Its eigenvalues are the real variances along the attractor's three principal axes — a standard technique for describing an attractor's shape.

rho = 28.0
long_traj = integrate([1.0, 1.0, 1.0], 5000, dt, rho)
burn_in = long_traj[1000:]  # drop the initial transient
cov = np.cov(burn_in.T)

mat_literal = "[[" + "], [".join(
    ", ".join(f"{v:.6f}" for v in row) for row in cov
) + "]]"
db.execute(f"MATRIX attractor_cov = {mat_literal}")
db.execute("LET cov_eigvals = EIGENVALUES attractor_cov")
linaldb_eigvals = sorted(db.execute("SHOW cov_eigvals").data)
numpy_eigvals = sorted(np.linalg.eigvalsh(cov))

print("linaldb EIGENVALUES(covariance):", [f"{v:.4f}" for v in linaldb_eigvals])
print("numpy   eigvalsh(covariance):   ", [f"{v:.4f}" for v in numpy_eigvals])
for a, b in zip(linaldb_eigvals, numpy_eigvals):
    assert abs(a - b) / abs(b) < 1e-3

linaldb EIGENVALUES(covariance):

 ['8.4261', '73.2233', '135.2743']
numpy   eigvalsh(covariance):    ['8.4261', '73.2233', '135.2743']

Honest limitation — the raw Jacobian

For completeness: TRACE works on any square matrix (it's just a diagonal sum), so we can still get one real, correct number from the actual (non-symmetric) Jacobian at a fixed point — but not its full spectrum. EIGENVALUES on the same matrix correctly refuses, rather than silently returning a wrong or partial answer.

def lorenz_jacobian(x, y, z, rho):
    return np.array([[-sigma, sigma, 0.0], [rho - z, -1.0, -x], [y, x, -beta]])

xs = np.sqrt(beta * (rho - 1))
J = lorenz_jacobian(xs, xs, rho - 1, rho)
j_literal = "[[" + "], [".join(", ".join(f"{v:.6f}" for v in row) for row in J) + "]]"
db.execute(f"MATRIX lorenz_jacobian = {j_literal}")
db.execute("LET jac_trace = TRACE lorenz_jacobian")
linaldb_trace = db.execute("SHOW jac_trace").data[0]

print(f"linaldb TRACE(Jacobian) = {linaldb_trace:.4f}   numpy trace = {np.trace(J):.4f}")
assert abs(linaldb_trace - np.trace(J)) < 1e-2

try:
    db.execute("LET jac_eigvals = EIGENVALUES lorenz_jacobian")
    raise AssertionError("expected EIGENVALUES to reject a non-symmetric matrix")
except linaldb.LinalError as e:
    print(f"\nEIGENVALUES on the raw (non-symmetric) Jacobian correctly errors:\n  {e}")

linaldb TRACE(Jacobian) = -13.6667   numpy trace = -13.6667

EIGENVALUES on the raw (non-symmetric) Jacobian correctly errors:
  [line 359] Engine error: Invalid operation: EIGENVALUES: matrix is not symmetric (entries [0][1]=10 vs [1][0]=1 differ) -- only symmetric matrices are supported today

The attractor itself

For the record — the shape all of this has been describing:

fig = plt.figure(figsize=(9, 7))
ax = fig.add_subplot(projection="3d")
ax.plot(long_traj[:, 0], long_traj[:, 1], long_traj[:, 2], linewidth=0.4, color="#2ca02c")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_title("The Lorenz attractor (rho=28)")
plt.tight_layout()
plt.show()

Summary

• Logistic map: 300 $r$ values swept, classified chaotic/stable by Lyapunov-exponent sign, cross-checked against the known period-doubling cascade and Feigenbaum point.
• Bool predicate regression check: WHERE is_chaotic = 1, bare WHERE is_chaotic, and WHERE is_chaotic = true all agree — the exact bug pattern fixed in engine v0.1.83 / linaldb 0.1.10, now verified in a real pipeline, not just a synthetic repro.
• In-place TRANSFORM regression check: dropping a column in place left the dataset immediately, correctly readable — the exact corruption pattern from the same fix, also verified for real.
• Lorenz system: sensitive dependence on initial conditions demonstrated directly through linaldb's own DISTANCE function (cross-checked against numpy), and the attractor's shape described honestly through EIGENVALUES on a real symmetric matrix (the covariance) — with the raw, non-symmetric Jacobian's EIGENVALUES constraint documented rather than worked around.

Every numerical result in this notebook was cross-checked against either a known analytical fact or a direct numpy computation — no cherry-picking.