From 34d52dffdaa75a550874ff109d1d68438e9a20ba Mon Sep 17 00:00:00 2001 From: chaoming Date: Fri, 19 Jun 2026 02:08:21 +0800 Subject: [PATCH] fix(integrators): correct RKF45 node and KlPl SRK diffusion weights - RKF45 6th-stage node was c6=1/3 instead of 1/2, silently collapsing the method to order 1 for any time-dependent f (verified order 1.0 -> 5.0) (Critical) - KlPl (srk_scalar) final-stage diffusion weights were wrong so g1+g2 != I1 and the scheme did not converge; corrected to g1=I1-I11/sqrt(h), g2=I11/sqrt(h) (High) Findings recorded in docs/issues-found-20260619-integrators-ode-sde.md --- brainpy/integrators/ode/adaptive_rk.py | 5 +- .../ode/ode_method_adaptive_rk_test.py | 47 +++++++ brainpy/integrators/sde/srk_scalar.py | 10 +- brainpy/integrators/sde/srk_scalar_test.py | 98 +++++++++++++ ...sues-found-20260619-integrators-ode-sde.md | 133 ++++++++++++++++++ 5 files changed, 289 insertions(+), 4 deletions(-) create mode 100644 docs/issues-found-20260619-integrators-ode-sde.md diff --git a/brainpy/integrators/ode/adaptive_rk.py b/brainpy/integrators/ode/adaptive_rk.py index 85826c2ab..77f3f2f1a 100644 --- a/brainpy/integrators/ode/adaptive_rk.py +++ b/brainpy/integrators/ode/adaptive_rk.py @@ -319,7 +319,10 @@ class RKF45(AdaptiveRKIntegrator): ('-8/27', 2, '-3544/2565', '1859/4104', -0.275)] B1 = ['16/135', 0, '6656/12825', '28561/56430', -0.18, '2/55'] B2 = ['25/216', 0, '1408/2565', '2197/4104', -0.2, 0] - C = [0, 0.25, 0.375, '12/13', 1, '1/3'] + # The 6th node is c6 = 1/2 (its A-row sums to 1/2); using 1/3 violates the + # consistency condition sum_j a_{6j} = c_6 and collapses the method to order 1 + # for time-dependent f. + C = [0, 0.25, 0.375, '12/13', 1, 0.5] register_ode_integrator('rkf45', RKF45) diff --git a/brainpy/integrators/ode/ode_method_adaptive_rk_test.py b/brainpy/integrators/ode/ode_method_adaptive_rk_test.py index 90f5ba2e4..786521537 100644 --- a/brainpy/integrators/ode/ode_method_adaptive_rk_test.py +++ b/brainpy/integrators/ode/ode_method_adaptive_rk_test.py @@ -89,3 +89,50 @@ def test_all_methods(self): adaptive_rk.HeunEuler]: bm.random.seed() run_integrator(method, show=False) + + +def _convergence_order(method_cls, f, exact, y0, t_end=2.0, ns=(20, 40, 80, 160)): + """Measure the empirical convergence order of the higher-order (B1) solution. + + The integrator is run in *fixed-step* mode (``adaptive=False``) on a + time-dependent scalar ODE with a known exact solution; the average + log2 error-ratio over successive step halvings gives the order ``p``. + """ + errs = [] + for n in ns: + dt = t_end / n + intg = method_cls(f, adaptive=False, dt=dt) + y = y0 + for i in range(n): + y = intg(y, i * dt, dt=dt) + errs.append(abs(float(bm.as_jax(y)) - exact(t_end))) + rates = [np.log2(errs[i] / errs[i + 1]) for i in range(len(errs) - 1)] + return float(np.mean(rates)), errs + + +class TestRKF45NodeFix(unittest.TestCase): + """Regression for P6-C1: RKF45 6th-stage node ``c6`` must be 1/2, not 1/3. + + With the wrong node the consistency condition ``sum_j a_{6j} = c_6`` is + violated and the order conditions break for any time-dependent ``f``, + collapsing the advertised order-5 solution to order ~1. This is only + visible when ``f`` depends on ``t`` (autonomous smoke tests miss it). + """ + + def test_rkf45_is_order5_on_time_dependent_ode(self): + # dy/dt = cos(t), y(0) = 0 -> y(t) = sin(t). + # Measure under x64 so the truncation error (not float32 round-off) + # governs the rate. The genuine order-5 method exceeds order 4; the + # buggy c6=1/3 variant measures order ~1. + f = lambda y, t: bm.cos(t) + exact = lambda t: np.sin(t) + bm.enable_x64() + try: + order, errs = _convergence_order(adaptive_rk.RKF45, f, exact, y0=0.0, + t_end=2.0, ns=(8, 16, 32, 64)) + finally: + bm.disable_x64() + self.assertGreater(order, 4.0, msg=f'RKF45 order={order:.2f}, errs={errs}') + + def test_rkf45_node_value(self): + self.assertAlmostEqual(float(eval(str(adaptive_rk.RKF45.C[-1]))), 0.5) diff --git a/brainpy/integrators/sde/srk_scalar.py b/brainpy/integrators/sde/srk_scalar.py index 68d852451..c21dc1a2a 100644 --- a/brainpy/integrators/sde/srk_scalar.py +++ b/brainpy/integrators/sde/srk_scalar.py @@ -366,11 +366,15 @@ def build(self): # 2.4 final stage # ---- - # g1 = (I1 - I11 / dt_sqrt + I10 / dt) - # g2 = I11 / dt_sqrt + # The strong order 1.0 SRK diffusion update is + # g1 = I1 - I11 / dt_sqrt + # g2 = I11 / dt_sqrt + # so that g1 + g2 = I1 reproduces the leading g * dW term for additive + # noise. The previous weights (g1 = -I1 + I11/dt_sqrt + I10/dt) summed to + # -I1 + 2 I11/dt_sqrt + I10/dt != I1, so the scheme did not converge. # y1 = x + dt * f_H0s1 + g1 * g_H1s1 + g2 * g_H1s2 for var in self.variables: - self.code_lines.append(f' {var}_g1 = -{var}_I1 + {var}_I11/dt_sqrt + {var}_I10/{constants.DT}') + self.code_lines.append(f' {var}_g1 = {var}_I1 - {var}_I11/dt_sqrt') self.code_lines.append(f' {var}_g2 = {var}_I11 / dt_sqrt') self.code_lines.append(f' {var}_new = {var} + {constants.DT} * {var}_f_H0s1 + ' f'{var}_g1 * {var}_g_H1s1 + {var}_g2 * {var}_g_H1s2') diff --git a/brainpy/integrators/sde/srk_scalar_test.py b/brainpy/integrators/sde/srk_scalar_test.py index 86785f1a3..ee0b059f0 100644 --- a/brainpy/integrators/sde/srk_scalar_test.py +++ b/brainpy/integrators/sde/srk_scalar_test.py @@ -79,3 +79,101 @@ def test_sdeint_dispatch(self, method): for i in range(10): x = intg(x, i * 0.01) assert np.all(np.isfinite(np.asarray(bm.as_jax(x)))) + + +# --------------------------------------------------------------------------- +# Regression for P6-H1: the KlPl final-stage diffusion weights were wrong +# (g1 = -I1 + I11/dt_sqrt + I10/dt, g2 = I11/dt_sqrt), so g1 + g2 != I1 and +# the scheme did not converge. The corrected weights are +# g1 = I1 - I11/dt_sqrt, g2 = I11/dt_sqrt (=> g1 + g2 = I1). +# --------------------------------------------------------------------------- + +_A, _B, _X0 = -0.5, 0.3, 1.0 # dX = a X dt + b X dW (Ito) + + +def _gbm_exact(w_total, t): + return _X0 * np.exp((_A - _B * _B / 2) * t + _B * w_total) + + +def _ref_klpl_terminal(n, dW, dI0, t_end=1.0): + """NumPy reference of the *corrected* KlPl step, driven by a fixed path.""" + dt = t_end / n + ds = np.sqrt(dt) + x = _X0 + for i in range(n): + I1 = dW[i] + I11 = 0.5 * (I1 ** 2 - dt) + f1 = _A * x + g1s1 = _B * x + h1s2 = x + dt * f1 + ds * g1s1 + g1s2 = _B * h1s2 + gg1 = I1 - I11 / ds + gg2 = I11 / ds + x = x + dt * f1 + gg1 * g1s1 + gg2 * g1s2 + return x + + +class TestKlPlConvergence: + """P6-H1: KlPl must converge (strong order ~1.0) on a linear SDE.""" + + def test_generated_weights_corrected(self): + # The generated step code must encode g1 = I1 - I11/dt_sqrt and + # g2 = I11/dt_sqrt (so g1 + g2 = I1). The buggy version had + # g1 = -x_I1 + x_I11/dt_sqrt + x_I10/dt. + import io + from contextlib import redirect_stdout + buf = io.StringIO() + with redirect_stdout(buf): + sde.KlPl(f=f, g=g, dt=0.01, show_code=True) + code = buf.getvalue() + assert 'x_g1 = x_I1 - x_I11/dt_sqrt' in code, code + assert 'x_g2 = x_I11 / dt_sqrt' in code, code + # the final-stage diffusion weight g1 must no longer contain the + # spurious mixed-integral term I10/dt + g1_line = [ln for ln in code.splitlines() if 'x_g1 =' in ln][0] + assert 'I10' not in g1_line, g1_line + + def test_reference_strong_convergence(self): + # The corrected scheme (NumPy reference) converges to the exact GBM + # solution as dt -> 0; the buggy scheme had a flat ~0.3 error. + errs = [] + for n in (50, 100, 200, 400): + dt = 1.0 / n + rng = np.random.default_rng(7) + path_errs = [] + for _ in range(300): + dW = rng.normal(0, np.sqrt(dt), n) + dI0 = rng.normal(0, np.sqrt(dt), n) + x = _ref_klpl_terminal(n, dW, dI0) + path_errs.append(abs(x - _gbm_exact(dW.sum(), 1.0))) + errs.append(float(np.mean(path_errs))) + rates = [np.log2(errs[i] / errs[i + 1]) for i in range(len(errs) - 1)] + assert np.mean(rates) > 0.7, f'KlPl not converging: errs={errs}, rates={rates}' + assert errs[-1] < 1e-3, f'KlPl terminal error too large: {errs}' + + def test_integrator_matches_reference(self): + # Run the *actual* compiled KlPl integrator and a NumPy reference on + # the same Wiener path. The generated code draws, per step, + # I1 = dt_sqrt * randn(), I0 = dt_sqrt * randn() + # so seeding bm.random lets us mirror the exact draws. + n, dt = 8, 0.01 + ds = np.sqrt(dt) + bm.random.seed(2024) + # mirror the generated draw order (I1 then I0 per step) with the SAME + # RandomState by drawing from it directly first to fix the sequence. + rng_state = bm.random.RandomState(2024) + dW = np.empty(n) + dI0 = np.empty(n) + for i in range(n): + dW[i] = ds * float(rng_state.randn()) + dI0[i] = ds * float(rng_state.randn()) + + bm.random.seed(2024) + intg = sde.KlPl(f=lambda x, t: _A * x, g=lambda x, t: _B * x, dt=dt) + x = _X0 + for i in range(n): + x = intg(x, i * dt) + got = float(bm.as_jax(x)) + + ref = _ref_klpl_terminal(n, dW, dI0, t_end=n * dt) + assert np.isclose(got, ref, rtol=1e-4, atol=1e-5), f'got={got}, ref={ref}' diff --git a/docs/issues-found-20260619-integrators-ode-sde.md b/docs/issues-found-20260619-integrators-ode-sde.md new file mode 100644 index 000000000..588d7268f --- /dev/null +++ b/docs/issues-found-20260619-integrators-ode-sde.md @@ -0,0 +1,133 @@ +# Integrators (ODE/SDE) Audit — Issues Found (2026-06-19) + +**Reviewer role:** Senior numerical-methods + JAX expert (ODE/SDE solvers) +**Scope:** `brainpy/integrators/{base,constants,joint_eq,runner,utils}.py`, +`brainpy/integrators/ode/{adaptive_rk,base,common,explicit_rk,exponential,generic}.py`, +`brainpy/integrators/sde/{base,generic,normal,srk_scalar}.py` (+ co-located `*_test.py`). +**Environment (verified):** Python 3.13 · jax 0.10.2 · brainpy 2.7.8 (CPU). `import brainpy` works; high-severity +findings are runtime-reproduced via empirical convergence-order tests. + +--- + +## Summary + +| Severity | Count | +|----------|-------| +| Critical | 1 | +| High | 1 | +| Medium | 0 | +| Low | 5 | + +The two correctness bugs are both **silent**: the integrators run and produce finite numbers, but the +numerical *order of accuracy* collapses. They were found by measuring strong/weak convergence rates against +SDEs/ODEs with known exact solutions, not by inspection alone. + +Prior-audit cross-check (`dev/issues-found-20260618.md`): C-12 (adaptive-RK `tol` default) and C-13 +(SDE `NameError` on bad `intg_type`/Heun guard) are **already fixed** in this worktree — verified by repro +(see P6-L5). + +--- + +### P6-C1 — RKF45 last node `c6 = 1/3` instead of `1/2`: method degrades to order 1 [Critical] +- File: brainpy/integrators/ode/adaptive_rk.py:322 +- Category: numerics +- What: `RKF45.C = [0, 0.25, 0.375, '12/13', 1, '1/3']`. The 6th stage node must be `c6 = 1/2` + (its `A`-row `(-8/27, 2, -3544/2565, 1859/4104, -11/40)` sums to `1/2`, and the rendered docstring + tableau line shows `1/2`). The code evaluates stage `k6` at `t + dt*(1/3)` instead of `t + dt*(1/2)`. +- Why it's a bug: The Runge–Kutta consistency condition `sum_j a_{ij} = c_i` is violated for row 6. + Whenever the derivative `f` depends on time `t`, the order conditions break and the 5th-order + `B1` solution collapses to **first order**. The companion error estimate is also wrong, so adaptive + step control is misled. For autonomous `f` (no `t`) the node is unused, which is why existing + smoke tests (autonomous / `f=-x`, Lorenz) never caught it. This is the default behavior of + `odeint(..., method='rkf45')`. +- Repro (measured, `dy/dt = cos(t)`, `y(0)=0`, exact `sin(t)`): + ``` + c6=1/3 (current): errors per halving 4.2e-3 -> 2.1e-3 -> 1.1e-3 -> 5.3e-4 order ~ 1.0 + c6=1/2 (correct): errors per halving 9.3e-7 -> 2.9e-8 -> 9.0e-10 -> 2.8e-11 order ~ 5.0 + ``` +- Fix: `C = [0, 0.25, 0.375, '12/13', 1, 0.5]`. +- Tests: `ode/ode_method_adaptive_rk_test.py::TestRKF45NodeFix` (order-5 convergence on time-dependent ODE). +- Status: fixed + +### P6-H1 — `KlPl` SDE final-stage diffusion weights wrong: method does not converge [High] +- File: brainpy/integrators/sde/srk_scalar.py:373-374 +- Category: numerics +- What: The KlPl (Kloeden–Platen, strong order 1.0, scalar Wiener) final stage uses + `g1 = -I1 + I11/dt_sqrt + I10/dt`, `g2 = I11/dt_sqrt`. The correct order-1.0 SRK diffusion update is + `g1 = I1 - I11/dt_sqrt`, `g2 = I11/dt_sqrt` (so the total diffusion weight `g1+g2 = I1`, reproducing + `g·ΔW` for additive noise). The current code has the sign of `I1` flipped, the sign of `I11/dt_sqrt` + flipped, and a spurious `+ I10/dt` term. The in-code comment at line 369 (`g1 = (I1 - I11/dt_sqrt + I10/dt)`) + also disagrees with the code and is itself only partly right (the `I10/dt` should not be there). +- Why it's a bug: The lowest-order consistency requirement for the diffusion is `sum_i (weight_i) = I1` + when `g` is constant. Current weights sum to `-I1 + 2·I11/dt_sqrt + I10/dt ≠ I1`, so the leading + `g·ΔW` term is wrong and the scheme does **not** converge to the true solution. +- Repro (measured, geometric Brownian motion `dX = aX dt + bX dW`, exact `X0 exp((a-b²/2)T + bW)`, + same driving path, mean abs error over 300 paths): + ``` + n= 50 KlPl_current = 2.7e-1 KlPl_fixed = 2.1e-3 + n=100 KlPl_current = 2.6e-1 KlPl_fixed = 1.1e-3 + n=200 KlPl_current = 3.1e-1 KlPl_fixed = 5.4e-4 (fixed halves each doubling -> order ~1.0) + n=400 KlPl_current = 2.9e-1 KlPl_fixed = 2.7e-4 + ``` + Current error is flat (~0.3) regardless of step size — no convergence. +- Fix: `g1 = {var}_I1 - {var}_I11/dt_sqrt`, `g2 = {var}_I11 / dt_sqrt`. Comment updated to match. +- Tests: `sde/srk_scalar_test.py::TestKlPlConvergence` (strong-order-1.0 convergence on linear SDE). +- Status: fixed + +--- + +## Low (recorded only — not fixed per instructions) + +### P6-L1 — Dead helper functions in `sde/normal.py` [Low] +- File: brainpy/integrators/sde/normal.py:37-60 +- Category: style +- What: `df_and_dg`, `dfdt`, `noise_terms` are module-level helpers that are never called (the SDE + integrator classes build their steps inline with `bm.random.randn`). `noise_terms` even references a + `random.normal(...).value` API and a code-generation style no longer used. +- Why it's a bug: Dead code; misleading (`noise_terms` looks like the noise generator but is unused). +- Fix: recorded only. +- Status: recorded-only + +### P6-L2 — `SRK2W1` stage-4 drift uses `f_H0s1` where tableau/comment imply `f_H0s3` [Low] +- File: brainpy/integrators/sde/srk_scalar.py:293 (comment at :289 says `0.25 * f_H0s3`) +- Category: numerics +- What: `H1s4 = x + 0.25*dt*f_H0s1 + dt_sqrt*(2*g_H1s1 - g_H1s2 + 0.5*g_H1s3)`. The drift uses + `f_H0s1`, but the comment (`H1s4 = x + dt * 0.25 * f_H0s3 + ...`) and the `A^(0)` tableau row + (nonzero entry in the third drift column) call for `f_H0s3`. +- Why it's *Low*: Empirically both variants retain strong order ≈1.5 on a linear SDE (errors + 4.1e-5/3.5e-5 at n=50, both ~order 1.5); the difference is within the method's own error constant + because the `H1s4` stage only feeds the `g4 = I111/dt` weight, an `O(h^1.5)` term. No order loss in + the tested regime, so recorded only (Low, not fixed). +- Fix: recorded only (would change `f_H0s1` -> `f_H0s3` at line 293 to match the published tableau). +- Status: recorded-only + +### P6-L3 — Adaptive step controller exponent fixed at 0.2 regardless of method order [Low] +- File: brainpy/integrators/ode/adaptive_rk.py:226 +- Category: numerics +- What: `factor = 0.9 * (tol / (error + 1e-12)) ** 0.2`. The exponent `0.2 = 1/5` is appropriate for an + order-4/5 pair (RKF45, DOPRI, CashKarp) but is used unchanged for `RKF12`/`HeunEuler` (order 1/2), + where `1/(p+1) = 1/2` would be the standard PI-free controller exponent. +- Why it's *Low*: Step-size control efficiency only; results stay correct (controller still + shrinks/grows toward tolerance), just sub-optimal step selection for the low-order pairs. +- Fix: recorded only. +- Status: recorded-only + +### P6-L4 — Docstring tableau typo in `BogackiShampine` (`4/90` should be `4/9`) [Low] +- File: brainpy/integrators/ode/adaptive_rk.py:453 +- Category: style +- What: Rendered docstring tableau shows `2/9 & 1/3 & 4/90` for the B1 row; the code `B1 = ['2/9','1/3','4/9',0]` + is correct. Pure doc typo. +- Why it's *Low*: Documentation only; code is right. +- Fix: recorded only. +- Status: recorded-only + +### P6-L5 — Prior-audit C-12 / C-13 already fixed (verification note) [Low / informational] +- Files: brainpy/integrators/ode/adaptive_rk.py:70,187 ; brainpy/integrators/sde/base.py:20 ; sde/normal.py:20,225 +- Category: correctness (informational) +- What: The 2026-06-18 audit reported `TypeError` on `odeint(method='rkf45', adaptive=True)` with no `tol` + (C-12) and `NameError: errors` on bad `intg_type` / the `Heun` Ito guard (C-13). In this worktree + `adaptive_rk.py` already imports `from brainpy import _errors as errors` and sets + `code_scope['tol'] = self.tol`, and `sde/base.py` / `sde/normal.py` already import `errors`. + Verified by repro: `odeint(f, method='rkf45', adaptive=True, dt=0.1)(1.,0.)` returns a value + new dt + (no crash); the still-broken numeric is the unrelated `c6` node (P6-C1). +- Status: recorded-only (already-fixed elsewhere; no action)