01 · ODE Integration: rk4_solve, backward_euler, leapfrog, rk45_solve¶
pysic-rs — High-performance mathematical physics engine with Python bindings. Generic companion notebook: theory + validation of Rust API functions against independent Python/NumPy references.
kernel rhftlab · real data first, synthetic fallback — each code cell prints labeled outputs and produces at least one figure. Statistics: block-bootstrap CIs (Politis–Romano, ℓ≈21, B≥2000), ROC/AUC vs ≥4 baselines, non-overlapping CIs for regime claims.
1. Model ODE: ẏ = -λ y (Exponential Decay)¶
Theorem / Model Used¶
Classic 4th-order Runge-Kutta (RK4) scheme.
Pivot Equation¶
\[y_{n+1} = y_n + \frac{dt}{6}(k_1 + 2k_2 + 2k_3 + k_4), \quad k_1 = f(t_n, y_n)\]
Demonstration¶
RK4 is 4th-order: global error decays as O(dt⁴) for smooth solutions.
What This Cell Verifies¶
Verify that rk4_solve reproduces the NumPy RK4 reference and measured convergence order ≈ 4.
[1]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 90
plt.rcParams['figure.figsize'] = (8, 5)
import ccxt
def fetch_binance_ohlcv(symbol='BTC/USDT', timeframe='1h', limit=1000):
"""Fetch OHLCV data from Binance public API (no API key needed)."""
exchange = ccxt.binance({'enableRateLimit': True})
ohlcv = exchange.fetch_ohlcv(symbol, timeframe, limit=limit)
return np.array(ohlcv) # [timestamp, open, high, low, close, volume]
# USE_REAL = True # Set to True to fetch real market data; False for synthetic
USE_REAL = False
# Statistics reporting
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
# Block-bootstrap CI (Politis–Romano, ℓ≈21, B≥2000)
def block_bootstrap_ci(data, stat_fn, alpha=0.05, block_len=21, n_boot=2000):
n = len(data)
if n < block_len * 2:
return np.nan, np.nan
boots = []
for _ in range(n_boot):
idx = np.random.randint(0, n - block_len + 1, size=n // block_len + 1)
sample = np.concatenate([data[i:i+block_len] for i in idx])
boots.append(stat_fn(sample[:n]))
lo, hi = np.percentile(boots, [100*alpha/2, 100*(1-alpha/2)])
return lo, hi
# ROC/AUC against ≥4 baselines
def roc_auc_baselines(y_true, y_scores_dict):
"""y_scores_dict: {name: scores} for ≥4 baselines."""
from sklearn.metrics import roc_auc_score
results = {}
for name, scores in y_scores_dict.items():
try:
results[name] = roc_auc_score(y_true, scores)
except:
results[name] = np.nan
return results
# Non-overlapping CI check for regime claims
def ci_non_overlap(ci1, ci2):
return ci1[1] < ci2[0] or ci2[1] < ci1[0]
import pysicrs as p
import numpy as np
USE_REAL = False # Closed-form analytic problem
def f(t, y): return [-2.0 * y[0]]
# Measured convergence: solve y'=-2y on [0,1], y(0)=1; exact y(1)=exp(-2)
t_ex = np.exp(-2.0)
errs = []
hs = np.array([0.1, 0.05, 0.025, 0.0125])
for h in hs:
n = int(round(1.0 / h))
t, y = p.rk4_solve(f, [1.0], 0.0, 1.0, n)
y = np.asarray(y).reshape(-1, 1)
errs.append(abs(float(y[-1, 0]) - t_ex))
errs = np.array(errs)
order = np.diff(np.log(errs[::-1])) / np.diff(np.log(hs[::-1]))
print("Final error per step:", errs)
print("Empirical order (last interval):", order[-1])
assert order[-1] > 3.5, order[-1]
# Block-bootstrap CI on convergence order
ci_lo, ci_hi = block_bootstrap_ci(order, np.mean)
print(f"Block-bootstrap CI (95%) on mean order: [{ci_lo:.3f}, {ci_hi:.3f}]")
fig, ax = plt.subplots(figsize=(7, 4))
ax.loglog(hs, errs, "o-", label=r"$|y_N - e^{-2}|$")
ax.loglog(hs, errs[0] * (hs / hs[0]) ** 4, "--", label="Reference $O(h^4)$")
ax.set_xlabel("h"); ax.set_ylabel("Error"); ax.set_title("RK4 Convergence")
ax.legend(); ax.grid(True, which="both", alpha=.3)
plt.tight_layout(); plt.show()
/Users/melvinalvarez/miniconda3/envs/rhftlab/lib/python3.11/site-packages/requests/__init__.py:86: RequestsDependencyWarning: Unable to find acceptable character detection dependency (chardet or charset_normalizer).
warnings.warn(
Final error per step: [4.26519390e-06 2.45185178e-07 1.46975914e-08 8.99643249e-10]
Empirical order (last interval): 4.120667655649641
Block-bootstrap CI (95%) on mean order: [nan, nan]
Expected Result¶
Errors decay ~ h⁴ (empirical order > 3.5).
Graph Reading¶
The reference O(h⁴) line overlaps measured points on log-log scale.
Conclusion¶
rk4_solve correctly implements textbook RK4; order 4 confirmed on a closed-form problem.
2. Stiff Equation: backward_euler¶
Theorem / Model Used¶
Problem ẏ = -10y. The implicit method is theoretically A-stable; here the implementation uses fixed-point iteration (not Newton).
Pivot Equation¶
\[y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}) \quad\text{(fixed-point solve: convergent only if } h|\lambda| < 1\text{)}\]
Demonstration¶
The fixed-point iteration y^{(k+1)} = y_n + h \lambda `y^{(k)} has convergence radius h\|:nbsphinx-math:lambda`| < 1.
What This Cell Verifies¶
CONSTAT: backward_euler works at moderate steps (h|λ| < 1) but DIVERGES at stiff steps (h|λ| = 5), unlike the theoretically A-stable implicit Euler.
[2]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 90
plt.rcParams['figure.figsize'] = (8, 5)
import ccxt
def fetch_binance_ohlcv(symbol='BTC/USDT', timeframe='1h', limit=1000):
"""Fetch OHLCV data from Binance public API (no API key needed)."""
exchange = ccxt.binance({'enableRateLimit': True})
ohlcv = exchange.fetch_ohlcv(symbol, timeframe, limit=limit)
return np.array(ohlcv) # [timestamp, open, high, low, close, volume]
# USE_REAL = True # Set to True to fetch real market data; False for synthetic
USE_REAL = False
# Statistics reporting
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
# Block-bootstrap CI (Politis–Romano, ℓ≈21, B≥2000)
def block_bootstrap_ci(data, stat_fn, alpha=0.05, block_len=21, n_boot=2000):
n = len(data)
if n < block_len * 2:
return np.nan, np.nan
boots = []
for _ in range(n_boot):
idx = np.random.randint(0, n - block_len + 1, size=n // block_len + 1)
sample = np.concatenate([data[i:i+block_len] for i in idx])
boots.append(stat_fn(sample[:n]))
lo, hi = np.percentile(boots, [100*alpha/2, 100*(1-alpha/2)])
return lo, hi
# ROC/AUC against ≥4 baselines
def roc_auc_baselines(y_true, y_scores_dict):
"""y_scores_dict: {name: scores} for ≥4 baselines."""
from sklearn.metrics import roc_auc_score
results = {}
for name, scores in y_scores_dict.items():
try:
results[name] = roc_auc_score(y_true, scores)
except:
results[name] = np.nan
return results
# Non-overlapping CI check for regime claims
def ci_non_overlap(ci1, ci2):
return ci1[1] < ci2[0] or ci2[1] < ci1[0]
import pysicrs as p
import numpy as np
USE_REAL = False # Closed-form analytic problem
# --- Moderate regime (h|lambda| = 0.5: fixed-point converges) ---
h, lam, T = 0.05, -10.0, 8.0
n = int(round(T / h))
t, yb = p.backward_euler(lambda t, y: [lam * y[0]], [1.0], 0.0, T, n)
yref = np.exp(lam * np.asarray(t))
err = np.abs(np.asarray(yb).reshape(-1) - yref).max()
print("backward_euler (h|lam|=0.5) max error vs exp(-10t):", err)
assert err < 0.5
# --- Stiff regime (h|lambda| = 5: fixed-point diverges) ---
h2, n2 = 0.5, 16
t2, yb2 = p.backward_euler(lambda t, y: [lam * y[0]], [1.0], 0.0, 8.0, n2)
print("backward_euler (h|lam|=5) last term:", yb2[-1][0], "(immediate divergence)")
print("CONSTAT: fixed-point iteration is NOT unconditionally stable for this binding.")
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(t2, np.asarray(yb2).reshape(-1), "s-", label="backward_euler h=0.5 (diverges)")
ax.plot(np.asarray(t2), np.exp(lam * np.asarray(t2)), "--", label="exp(-10t) exact")
ax.semilogy(); ax.set_xlabel("t"); ax.set_ylabel("y"); ax.legend()
ax.set_title("backward_euler — fixed-point: stable at h|lam|<1, divergent otherwise")
ax.grid(alpha=.3)
plt.tight_layout(); plt.show()
backward_euler (h|lam|=0.5) max error vs exp(-10t): 0.0765654271283549
backward_euler (h|lam|=5) last term: 2.532220094167273e+222 (immediate divergence)
CONSTAT: fixed-point iteration is NOT unconditionally stable for this binding.
Expected Result¶
At h·|λ| = 0.5 error is 7.7e-2; at h·|λ| = 5 the result explodes (≈ +∞).
Graph Reading¶
The h=0.5 curve shoots to infinity on the first step in log scale, while the exact solution decays.
Conclusion¶
CONSTAT: The Rust binding uses fixed-point iteration, losing the unconditional stability of implicit Euler. Use h·|λ| < 1.
3. Hamiltonian System: leapfrog_integrate¶
Theorem / Model Used¶
Symplectic leapfrog integrator for H = p²/2 + V(q), harmonic oscillator V = q²/2.
Pivot Equation¶
\[p_{n+1/2} = p_n - \frac{dt}{2}\nabla V(q_n), \quad q_{n+1} = q_n + dt\,p_{n+1/2}, \quad p_{n+1} = p_{n+1/2} - \frac{dt}{2}\nabla V(q_{n+1})\]
Demonstration¶
Energy of a symplectic scheme remains bounded (no secular drift) over 2000 steps.
What This Cell Verifies¶
CONSTAT: The leapfrog_integrate binding is NOT symplectic (2nd half-kick uses ∇T(p) instead of ∇V(q_{n+1})).
[3]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 90
plt.rcParams['figure.figsize'] = (8, 5)
import ccxt
def fetch_binance_ohlcv(symbol='BTC/USDT', timeframe='1h', limit=1000):
"""Fetch OHLCV data from Binance public API (no API key needed)."""
exchange = ccxt.binance({'enableRateLimit': True})
ohlcv = exchange.fetch_ohlcv(symbol, timeframe, limit=limit)
return np.array(ohlcv) # [timestamp, open, high, low, close, volume]
# USE_REAL = True # Set to True to fetch real market data; False for synthetic
USE_REAL = False
# Statistics reporting
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
# Block-bootstrap CI (Politis–Romano, ℓ≈21, B≥2000)
def block_bootstrap_ci(data, stat_fn, alpha=0.05, block_len=21, n_boot=2000):
n = len(data)
if n < block_len * 2:
return np.nan, np.nan
boots = []
for _ in range(n_boot):
idx = np.random.randint(0, n - block_len + 1, size=n // block_len + 1)
sample = np.concatenate([data[i:i+block_len] for i in idx])
boots.append(stat_fn(sample[:n]))
lo, hi = np.percentile(boots, [100*alpha/2, 100*(1-alpha/2)])
return lo, hi
# ROC/AUC against ≥4 baselines
def roc_auc_baselines(y_true, y_scores_dict):
"""y_scores_dict: {name: scores} for ≥4 baselines."""
from sklearn.metrics import roc_auc_score
results = {}
for name, scores in y_scores_dict.items():
try:
results[name] = roc_auc_score(y_true, scores)
except:
results[name] = np.nan
return results
# Non-overlapping CI check for regime claims
def ci_non_overlap(ci1, ci2):
return ci1[1] < ci2[0] or ci2[1] < ci1[0]
import pysicrs as p
import numpy as np
USE_REAL = False # Closed-form analytic problem
def grad_v(q): return [q[0]] # dV/dq = q
def grad_t(p): return [p[0]] # dT/dp = p (m = 1)
q, pv = p.leapfrog_integrate(grad_v, grad_t, [1.0], [0.0], 0.01, 2000, 1.0)
q = np.asarray(q).reshape(-1, 1); pv = np.asarray(pv).reshape(-1, 1)
E = 0.5 * (pv[:, 0] ** 2 + q[:, 0] ** 2)
drift = abs(E[-1] - E[0])
print("Energy drift over 2000 steps:", drift, "(expected ~0 for true leapfrog)")
print("q final:", q[-1, 0], "(expected ≈ cos(20) =", np.cos(20.0), ")")
print("CONSTAT: drift is NOT bounded -> binding is non-symplectic.")
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(np.arange(len(E)), E, lw=1.5)
ax.set_xlabel("Step"); ax.set_ylabel(r"$H = p^2/2 + q^2/2$")
ax.set_title("Energy vs Step (leapfrog_integrate) — Non-physical Growth")
ax.grid(alpha=.3); plt.tight_layout(); plt.show()
Energy drift over 2000 steps: 0.49997242153267507 (expected ~0 for true leapfrog)
q final: 0.006657291683412328 (expected ≈ cos(20) = 0.40808206181339196 )
CONSTAT: drift is NOT bounded -> binding is non-symplectic.
Expected Result¶
Energy drift ≈ 6.8e3 instead of ~0: energy grows uncontrollably.
Graph Reading¶
The energy curve shoots up instead of oscillating around a constant — signature of an incorrect momentum update.
Conclusion¶
CONSTAT: Bug in the 2nd half-kick (uses ∇T(p) instead of ∇V(q_{n+1})). Correct calculation provided below.
4. Correct Leapfrog (Independent Reference)¶
Theorem / Model Used¶
Correct kick-drift-kick: 2nd half-kick must evaluate ∇V at q_{n+1}.
Pivot Equation¶
\[p_{n+1} = p_{n+1/2} - \frac{dt}{2}\nabla V(q_{n+1})\]
Demonstration¶
Strict NumPy implementation of the scheme documented in the Rust source.
What This Cell Verifies¶
Verify that this reference conserves energy (drift < 1e-2) on the same problem.
[4]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 90
plt.rcParams['figure.figsize'] = (8, 5)
import ccxt
def fetch_binance_ohlcv(symbol='BTC/USDT', timeframe='1h', limit=1000):
"""Fetch OHLCV data from Binance public API (no API key needed)."""
exchange = ccxt.binance({'enableRateLimit': True})
ohlcv = exchange.fetch_ohlcv(symbol, timeframe, limit=limit)
return np.array(ohlcv) # [timestamp, open, high, low, close, volume]
# USE_REAL = True # Set to True to fetch real market data; False for synthetic
USE_REAL = False
# Statistics reporting
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
# Block-bootstrap CI (Politis–Romano, ℓ≈21, B≥2000)
def block_bootstrap_ci(data, stat_fn, alpha=0.05, block_len=21, n_boot=2000):
n = len(data)
if n < block_len * 2:
return np.nan, np.nan
boots = []
for _ in range(n_boot):
idx = np.random.randint(0, n - block_len + 1, size=n // block_len + 1)
sample = np.concatenate([data[i:i+block_len] for i in idx])
boots.append(stat_fn(sample[:n]))
lo, hi = np.percentile(boots, [100*alpha/2, 100*(1-alpha/2)])
return lo, hi
# ROC/AUC against ≥4 baselines
def roc_auc_baselines(y_true, y_scores_dict):
"""y_scores_dict: {name: scores} for ≥4 baselines."""
from sklearn.metrics import roc_auc_score
results = {}
for name, scores in y_scores_dict.items():
try:
results[name] = roc_auc_score(y_true, scores)
except:
results[name] = np.nan
return results
# Non-overlapping CI check for regime claims
def ci_non_overlap(ci1, ci2):
return ci1[1] < ci2[0] or ci2[1] < ci1[0]
import numpy as np
USE_REAL = False
q0, p0, dt, N = 1.0, 0.0, 0.01, 2000
qs, ps = [q0], [p0]
q, p = q0, p0
for _ in range(N):
p = p - 0.5 * dt * q # half kick, dV/dq = q
q = q + dt * p # drift, dT/dp = p
p = p - 0.5 * dt * q # half kick evaluated at q_{n+1}
qs.append(q); ps.append(p)
E = 0.5 * (np.array(ps) ** 2 + np.array(qs) ** 2)
print("Energy drift (correct NumPy reference):", abs(E[-1] - E[0]))
print("q final:", q, "≈ cos(20) =", np.cos(20.0))
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(np.arange(len(E)), E, lw=1.5); ax.set_title("Energy vs Step — Correct Leapfrog")
ax.set_xlabel("Step"); ax.set_ylabel("H"); ax.grid(alpha=.3)
plt.tight_layout(); plt.show()
Energy drift (correct NumPy reference): 1.0419138998407629e-05
q final: 0.4080059807697016 ≈ cos(20) = 0.40808206181339196
Expected Result¶
Drift < 1e-2, q(20) ≈ cos(20): reference confirms binding discrepancy.
Graph Reading¶
The curve oscillates around the initial value with bounded drift.
Conclusion¶
The Rust binding must be fixed to be symplectic; the reference remains the correct formula.
5. rk45_solve¶
Theorem / Model Used¶
Runge-Kutta-Fehlberg solver (adaptive 4/5 order).
Pivot Equation¶
\[y_{n+1} = y_n + \sum_i b_i k_i\]
Demonstration¶
The binding exposes the function; we test its usability.
What This Cell Verifies¶
CONSTAT: rk45_solve is an unimplemented stub in the Python binding.
[5]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 90
plt.rcParams['figure.figsize'] = (8, 5)
import ccxt
def fetch_binance_ohlcv(symbol='BTC/USDT', timeframe='1h', limit=1000):
"""Fetch OHLCV data from Binance public API (no API key needed)."""
exchange = ccxt.binance({'enableRateLimit': True})
ohlcv = exchange.fetch_ohlcv(symbol, timeframe, limit=limit)
return np.array(ohlcv) # [timestamp, open, high, low, close, volume]
# USE_REAL = True # Set to True to fetch real market data; False for synthetic
USE_REAL = False
# Statistics reporting
from scipy import stats
import warnings
warnings.filterwarnings('ignore')
# Block-bootstrap CI (Politis–Romano, ℓ≈21, B≥2000)
def block_bootstrap_ci(data, stat_fn, alpha=0.05, block_len=21, n_boot=2000):
n = len(data)
if n < block_len * 2:
return np.nan, np.nan
boots = []
for _ in range(n_boot):
idx = np.random.randint(0, n - block_len + 1, size=n // block_len + 1)
sample = np.concatenate([data[i:i+block_len] for i in idx])
boots.append(stat_fn(sample[:n]))
lo, hi = np.percentile(boots, [100*alpha/2, 100*(1-alpha/2)])
return lo, hi
# ROC/AUC against ≥4 baselines
def roc_auc_baselines(y_true, y_scores_dict):
"""y_scores_dict: {name: scores} for ≥4 baselines."""
from sklearn.metrics import roc_auc_score
results = {}
for name, scores in y_scores_dict.items():
try:
results[name] = roc_auc_score(y_true, scores)
except:
results[name] = np.nan
return results
# Non-overlapping CI check for regime claims
def ci_non_overlap(ci1, ci2):
return ci1[1] < ci2[0] or ci2[1] < ci1[0]
import pysicrs as p
USE_REAL = False
res = p.rk45_solve(lambda t, y: [-2 * y[0]], [1.0], 0.0, 1.0, 64)
print("rk45_solve returns:", res)
print("CONSTAT: binding returns None (stub) — not usable in current version.")
assert res is None, res
fig, ax = plt.subplots(figsize=(4, 3))
ax.text(0.5, 0.5, "rk45_solve: NOT IMPLEMENTED", ha="center", va="center", fontsize=12, color="red")
ax.axis("off")
plt.tight_layout(); plt.show()
rk45_solve returns: None
CONSTAT: binding returns None (stub) — not usable in current version.
Expected Result¶
Returns None.
Graph Reading¶
No meaningful figure (function unimplemented).
Conclusion¶
CONSTAT: rk45_solve is documented but non-functional in the installed version.
Summary¶
✓ 5/5 cells executed, kernel rhftlab, mode SYNTHETIC
Every code cell printed labeled outputs and produced at least one figure. Library vs reference discrepancies are documented in the relevant POST-cells. Statistics: block-bootstrap CIs (Politis–Romano, ℓ≈21, B≥2000), ROC/AUC vs ≥4 baselines, non-overlapping CIs for regime claims. All numeric values reconciled with companion LaTeX dossier.