fix(integrators): RKF45 node + KlPl SRK diffusion weights

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)

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)

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')

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}'

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)