Changelog.

Unreleased

0.1.5 - 2026-05-22

Added

• F-IC-SENS-MASS-SURFACE (additive, non-breaking) — numra-ode: ParametricOdeSystem gains the four-method mass-matrix surface mirroring OdeSystem — has_mass_matrix() -> bool (default false), mass_matrix(&self, mass: &mut [S]) (default fills row-major identity, length n²), is_singular_mass() -> bool (default false), algebraic_indices() -> Vec<usize> (default empty) — plus is_autonomous() -> bool (default false). Defaults reproduce the pre-0.1.5 implicit behavior exactly: every existing implementor (AugmentedSystem, ClosureSystem, ParametricMOLSystem2D/3D, IcAsParametric, the &T blanket impl, every user-defined external impl) compiles unchanged with no source-level changes. The &T blanket impl forwards the new methods so &Sys carries M into solve_forward_sensitivity without ownership transfer. Foundation Spec §3.8 added (was missing for ParametricOdeSystem — a documentation drift this entry closes); §3.3 “Mass matrix as a first-class composable” deferred bullet updated to note partial concretization; §6 #14 records the decision; §7 #3 updated (split into two independent sub-questions with separate triggers: identity-default removal vs DaeSystem peer trait — see §7 #3); §7 #9 added (IC manifold-projection follow-up, four-axis design space).

Fixed

• F-IC-SENS-MASS-SURFACE (correctness fix) — numra-ode: solve_initial_condition_sensitivity{,_with} now honors the host OdeSystem’s mass matrix. Closes a correctness gap in 0.1.4’s F-IC-SENS: the IcAsParametric wrapper introduced in 0.1.4 (§6 #13) silently discarded the host’s has_mass_matrix / mass_matrix / is_singular_mass / algebraic_indices / is_autonomous declarations when re-presenting the system as ParametricOdeSystem, so the augmented variational integration ran against an identity mass matrix regardless of what the user’s OdeSystem declared. Fixed at two boundaries: (i) IcAsParametric now forwards all five methods from the wrapped OdeSystem into the now-richer ParametricOdeSystem; (ii) AugmentedSystem’s OdeSystem impl lifts the host’s N × N mass matrix block-diagonally onto the augmented (N · (1 + N_s)) × (N · (1 + N_s)) state, with (N_s + 1) copies of M on the diagonal (one for the state block, one per sensitivity-column block — the variational equation M · (∂y/∂p_k)' = J_y · (∂y/∂p_k) + J_p_{:,k} places each sensitivity column in the same M-weighted space as the state). Host algebraic indices lift correspondingly: host index i maps to augmented indices {i, N+i, 2N+i, …, N_s·N+i}. The has_mass_matrix() = false branch of the lifted mass_matrix preserves the semantic distinction between “host has M, M happens to be identity” and “host has no M” — downstream solvers (Radau5, BDF) dispatch on it.

  Precisely-supported affected-scope claim for 0.1.4 (verified empirically against HEAD:numra-ode/src/sensitivity.rs:152-281 at commit 98ee6ca, 0.1.4 head: the pre-extension trait body contains exactly n_states, n_params, params, rhs_with_params, rhs, jacobian_y, jacobian_p, initial_sensitivity, has_analytical_jacobian_y, has_analytical_jacobian_p and none of the M-surface methods):

  • Users of solve_forward_sensitivity{,_with} with ParametricOdeSystem implementors: unaffected — the pre-0.1.5 trait had no mass-matrix surface, so users had no way to declare M through this path; the missing delegation in AugmentedSystem was a compile-time fact (the problem couldn’t be expressed), not a runtime numerical failure.
  • Users of solve_initial_condition_sensitivity{,_with} with identity-M OdeSystem (the common ODE case): unaffected — the silently-dropped M was identity anyway.
  • Users of solve_initial_condition_sensitivity{,_with} with non-identity M (e.g. finite-element-discretized PDEs): received Φ off by an M⁻¹ factor in the matrix-exponential’s argument — Φ_bug(t) = expm(A·t) instead of the analytical Φ_truth(t) = expm(M⁻¹·A·t). Pinned by the pre-fix-bug numerical value exp(-2) ≈ 0.13534 in numra-ode/tests/ic_sensitivity_mass_matrix.rs::path_b_non_identity_mass_ic_sensitivity, derived in closed form: the test uses the scalar DAE 2y' = -y (M = 2, A = -1); M-blind integration drops M and runs y' = -y, giving y_bug(t) = exp(-t) and Φ_bug(t) = exp(-t), versus truth Φ_truth(t) = exp(-t/2). At t = 2: Φ_bug = exp(-2) ≈ 0.13534, Φ_truth = exp(-1) ≈ 0.36788. The scalar choice makes the closed-form derivation visible; higher-dimensional non-identity-M systems exhibit the same shape (matrix-exponential argument off by M⁻¹) but produce no scalar reduction. Revert-confirm protocol verified: revert the wiring → got 0.13534, restore → got 0.36788 = exp(-1).
  • Users with singular M (semi-explicit index-1 DAEs): received variational integrations that treated algebraic rows as differential — wrong dynamics, not just wrong scaling — pinned by the pre-fix-bug numerical value y₂(2) ≈ 0.73039 (= (1/2)·(sin(2) − cos(2) + e⁻²), the closed-form solution to y₂' + y₂ = sin(t) that the M-blind integration actually solves) in path_b_singular_mass_dae_ic_sensitivity.

  Upgrade guidance, severity-split:

  • Upgrade recommended for any 0.1.4 user of solve_initial_condition_sensitivity{,_with} with non-identity M: the pre-fix result is quantitatively wrong (an M⁻¹-factor distortion of Φ), not qualitatively. Existing analyses may be approximately correct depending on how the user used Φ.
  • Upgrade required for any 0.1.4 user of solve_initial_condition_sensitivity{,_with} with singular M (DAE): the pre-fix result is qualitatively wrong — algebraic constraints were treated as differential equations, producing trajectories that satisfy the wrong dynamical system entirely. Any downstream analysis using the pre-fix Φ or state trajectory on a singular-M system should be re-run on 0.1.5.

  Yank decision: 0.1.4 is not yanked. Bounded scope (specific feature path × specific user category); academic citation record references 0.1.4 immutably (Zenodo DOI back-port landed in 47d25a9); 0.1.5 supersedes with a clear CHANGELOG connection.

  Out-of-scope, tracked separately: the finite-difference IC perturbation in fd_jacobian_y_inline (sensitivity.rs:394) does not project the perturbed initial condition onto the constraint manifold for singular-M systems — post-fix DAE sensitivity is “wrong sensitivity-w.r.t.-perturbed-IC because the perturbed IC may drift off-manifold, correctly integrated because the variational integration is now M-aware” — strictly improved over 0.1.4’s both-wrong state, but not manifold-projected. Tracked as Foundation Spec §7 #9 (four-axis design space: when to project, in which direction, via which API, under which correctness criterion — linear vs nonlinear constraints have different criteria).

  Six tests pin the fix end-to-end at numra-ode/tests/ic_sensitivity_mass_matrix.rs plus the symmetric structural test at numra-ode/src/sensitivity.rs::tests::ic_as_parametric_forwards_mass_surface (lives in the lib because IcAsParametric is private): two Path A analytic-oracle tests (post-fix only — compile-time fact, would not compile against 0.1.4), two Path B analytic-oracle tests with falsifiability-gate assertions on the pre-fix-bug values (runtime fact, exact numerical match verified), and two structural unit tests completing the boundary-symmetry pair. Revert-confirm-restore protocol applied: temporarily reverting IcAsParametric::has_mass_matrix / mass_matrix / is_singular_mass / algebraic_indices to defaults → Path B tests fail with the exact predicted pre-fix-bug numerical values (exp(-2) for non-identity-M, 0.73039 for singular-M DAE) → restore wiring → tests pass. A seventh sanity test (augmented_system_identity_m_default_path_unchanged) pins the identity-M path so the semantic distinction between “host has no M” and “host has identity M” survives future refactors.

  Saved memory (feedback_wrapper_boundary_symmetry): wrappers between foundation types verify boundary discipline in both directions (outward: no leakage; inward: no loss of foundation capabilities). F-IC-SENS’s 0.1.4 boundary check was outward-only; this entry closes the inward gap and records the lesson for future wrappers.

Changed

• Workspace version 0.1.4 → 0.1.5 prep (additive correctness release): workspace [workspace.package].version and all 20 internal-dependency pins in [workspace.dependencies] bumped in lockstep (all 21 publishable crates inherit via version.workspace = true); CITATION.cff version: field 0.1.4 → 0.1.5, date-released 2026-05-18 → 2026-05-22. The version-DOI block in identifiers is intentionally left at 0.1.4 values — Zenodo mints the 0.1.5 DOI on GitHub Release, after which a follow-up docs(citation): back-port Zenodo DOIs for v0.1.5 commit fills it in. Same pattern as 5aa85d0 (0.1.2) / 8b3ec0a (0.1.3) / 47d25a9 (0.1.4). No existing public signature changes meaning: every pre-0.1.5 ParametricOdeSystem implementor remains source-compatible with the new trait default-impl additions; solve_forward_sensitivity{,_with} users are unaffected at runtime; solve_initial_condition_sensitivity{,_with} users see only the corrected M-aware variational integration (no change for identity-M; correctness restoration for non-identity-M and singular-M).

0.1.4 - 2026-05-18

Added

• F-IC-SENS (additive) — numra-ode: first-class initial-condition sensitivity. New public functions solve_initial_condition_sensitivity and solve_initial_condition_sensitivity_with, plus the result type StateTransitionResult<S>, compute the state-transition matrix Φ(t) = ∂y(t)/∂y₀ of the linearised flow dy/dt = f(t, y) over any [Solver]. Implemented as a private IcAsParametric<'a, S, Sys> wrapper that reinterprets a plain [OdeSystem<S>] as a [ParametricOdeSystem<S>] with n_params := n_states, J_p ≡ 0, and S(t₀) = I_N, then delegates to the existing solve_forward_sensitivity / AugmentedSystem machinery — no fork of the variational integrator. The Phase-0 discovery established two architecturally available routes (variational seeding vs AD-through-y₀); the AD-through-y₀ route would require either making Dual<f64>: Scalar with dual-aware step-control norms across every solver (the existing Dual<S>::PartialOrd is primal-only, so the tangent — the STM itself — would integrate with no adaptive error control) or a fork of the integrator. Both are forbidden by the foundation/composability story for this scope; the variational route is the architecturally correct choice given the actual codebase, not a lesser-route-for-task-scope compromise. The would-be unification of Features 1+2 via AD-through-y₀ remains a sensible future foundation-design pass (Foundation Spec §7 #1’s per-call/per-system/global question) but is not preempted by this PR. Constraint-1 seam: the ParametricOdeSystem length assertion at sensitivity.rs:852-858 requires a params() slice of length n_params; the IC wrapper holds a zeroed dummy Vec<S> of length N for this purpose, fully encapsulated — the user never constructs, passes, or reasons about a parameter vector, and no parameter concept appears in solve_initial_condition_sensitivity / StateTransitionResult’s signatures or rustdoc examples. Invariant IC-NO-PARAM-READ is stated as a mathematical fact: IcAsParametric::rhs_with_params(t, y, p, dy) provably never reads p because the wrapped OdeSystem::rhs(t, y, dy) has no p argument; therefore ∂f/∂p ≡ 0 is mathematically exact, the wrapper’s jacobian_p writing zeros with has_analytical_jacobian_p() = true is a truthful analytical declaration, and the AugmentedSystem debug consistency check at sensitivity.rs:373 is correctly suppressed by the flag (not silenced). The invariant is pinned by ic_dummy_params_are_unread (numra-ode/tests/ic_sensitivity_regression.rs), which runs the variational integration three times against the same wrapper shape with dummy params() = [0.0], [f64::NAN], [f64::INFINITY] — bit-for-bit identical results across all three prove the RHS does not consume p (if it did, NaN/∞ would propagate through the augmented RHS into y(t) and Φ(t)). Same revert-confirm-restore falsifiability shape as F-SOLVER-FIELDS’s pinning test. Specialization-vs-core regression pinned by ic_matches_hand_seeded_variational: asserts bit-for-bit equality between the IC API result and a hand-built ParametricOdeSystem configured with N_s = N, J_p ≡ 0, S₀ = I (the same shape the wrapper produces). This test passes before the analytic-oracle tests (sequencing established during Phase 2 implementation) so wrapper-correctness is isolated from variational-math correctness — a Phase-2 failure would be diagnostic rather than ambiguous. Analytic-oracle tests scalar_linear_matches_exp_minus_kt and linear_2x2_matches_expm_a_t verify Φ(t) = exp(-k t) on dy/dt = -k y and Φ(t) = expm(A·t) on a 2×2 system with closed-form solution. Boundary check (deliverable E): the public API surface mentions only t, y, phi, phi_ij, phi_column, final_phi, final_y, dim, len, is_empty, success, message, stats, with parameter names i (time index), row, col (state-transition matrix indices in [0, N)); no name in the public API mentions parameter, sensitivity, p, sigma, uncertainty, measurement, or any application-domain concept; StateTransitionResult<S> holds the underlying SensitivityResult<S> in a private field with no Deref impl and no re-export, so no *_param-style accessor leaks through. The rustdoc worked example is dy/dt = -k y ⇒ Φ(t) = exp(-k t), scalar, no parameters present.

• F-AUTODIFF-JY (additive, behind cargo feature autodiff) — numra-ode: an [OdeSystem<S>] adapter AutodiffJacobianSystem<S, F> whose jacobian() is computed by forward-mode automatic differentiation through numra-autodiff::Dual<S>. The adapter owns a Dual<S>-typed RHS closure and a pair of pre-allocated dual buffers (RefCell<Vec<Dual<S>>> × 2), so per-step cost is zero allocation plus N + 1 Dual<S>-RHS evaluations per integrator step (1 for the primal RHS, N for the N-direction Jacobian sweep). The existing FD-default Jacobian path performs N + 1 S-RHS evaluations; the AD adapter is therefore in the same complexity class (O(N) RHS evals per Jacobian), at a ~2× per-eval cost (each Dual<S> op carries one extra fused-multiply-add for the tangent), in exchange for a Jacobian accurate to round-off rather than ~sqrt(S::EPSILON) ≈ 1.5e-8. The accuracy delta materially benefits stiff implicit solvers (Radau5, BDF, ESDIRK*) whose Newton convergence depends on Jacobian fidelity. Behind cargo feature autodiff (default off): callers who do not opt in do not compile or pull numra-autodiff and pay zero in compile time and binary size. The autodiff path is never silently substituted — selection is explicit by constructor (AutodiffJacobianSystem::new(rhs_dual, dim)); a user supplying analytical J_y via OdeSystem::jacobian override is entirely unaffected in behavior and performance. The adapter composes orthogonally with F-IC-SENS: passing an AutodiffJacobianSystem to solve_initial_condition_sensitivity obtains Φ(t) with AD-exact J_y. Correctness pinned by autodiff_jy_matches_analytic_jy_on_lotka_volterra (numra-ode/tests/ic_sensitivity_autodiff.rs, #[cfg(feature = "autodiff")]): runs the IC API twice on Lotka-Volterra — once with a hand-written analytical Jacobian, once with AutodiffJacobianSystem wrapping the same RHS — and asserts Φ(t_f) agrees to < 1e-8 and y(t_f) to < 1e-9. Three additional unit tests (ad_jacobian_2x2_nonlinear, ad_rhs_matches_f64_primal, ad_jacobian_reentrant) cover the adapter in isolation. Boundary check (deliverable E): AutodiffJacobianSystem<S, F>::new(rhs_dual, dim) and the OdeSystem<S> methods name only t, y, dy, jy, dim — standard ODE / Jacobian vocabulary; the rustdoc worked example is the Lotka-Volterra RHS typed against Dual<f64>, no parameter / sigma / uncertainty / domain concept in the signature or example. Resolves Foundation Specification §7 #1 (open question “Autodiff integration shape — per-call adapter, per-system constructor, or global feature flag”) for the Jacobian-source axis: the answer this PR commits to is per-system adapter (OdeSystem<S> impl, bound to the system once at construction) made opt-in via cargo feature gate. The three options in §7 #1 collapse: adapter-as-OdeSystem (option a) is the API mechanism; per-system constructor (option b) is the entry-point shape; cargo-feature gate (option c) is the compile-time discoverability. §3.3’s deferred “Jacobian enum (AutoDiff | Analytical | FD) would change the trait” foundational change is not preempted by this PR — the adapter composes below the OdeSystem trait shape; if a future foundation-design pass adopts the trait-shape change, this adapter may be refactored or absorbed, but the public API surface today is not load-bearing on that future decision. The §7 #1 narrower scope (“every Jacobian/gradient call site”) is not fully closed here — gradient-side wiring (numra-optim, numra-fit) remains open; only the OdeSystem::jacobian axis concretises.

• Foundation Specification §6 decision-log entry (13) — pending SHA: 0.1.4 prep PR (will substitute the actual commit SHA at merge time, per the two-commit <pending> → SHA pattern established by §6 #12). Records the F-AUTODIFF-JY resolution above as the per-system-adapter-plus-cargo-feature answer to §7 #1’s Jacobian-source axis; §7 #1 is rewritten in the same PR to reflect the narrower remaining scope (gradient-side wiring open). The Foundation Spec §3.3 deferred Jacobian-enum question remains untouched by this entry.

• Interop workflow 7 (numra/tests/interop_workflows.rs): workflow_ic_sensitivity_interp_integrate — solves the IC sensitivity equations for dy/dt = -k y, extracts Φ_{0,0}(t) across the trajectory, fits a numra-interp::CubicSpline, and integrates the spline with numra-integrate::quad. Verifies StateTransitionResult flows into both numra-interp and numra-integrate without adapter code, returns Result<(), NumraError> via ? across three crate boundaries — closing composability contract item 5 (foundation-flowable output) + item 7 (interop test) for the IC API. Analytical oracle: ∫_0^3 exp(-k t) dt = (1 - exp(-3k)) / k; assertion tolerance < 1e-3.

Changed

• Workspace version 0.1.3 → 0.1.4 prep (additive minor release): workspace [workspace.package].version and all 21 internal-dependency pins in the workspace Cargo.toml bumped; CITATION.cff version: field bumped (date stamp deferred to the actual release-prep PR, per the established convention from 514fdf7). No existing public signature changes meaning: ParametricOdeSystem, AugmentedSystem, SensitivityResult, solve_forward_sensitivity{,_with}, ClosureSystem, OdeSystem — all unchanged in shape and semantics; the 20 pre-existing sensitivity tests (10 unit + 10 regression) remain green (D5 baseline re-run during Phase 2 sign-off). numra-ode adds the optional numra-autodiff workspace-dep and a new autodiff cargo feature (default off). The numra facade (numra/src/lib.rs) re-exports solve_initial_condition_sensitivity, solve_initial_condition_sensitivity_with, and StateTransitionResult at the top level, symmetric with the existing solve_forward_sensitivity promotion.

0.1.3 - 2026-05-17

Changed

• F-SOLVER-FIELDS (additive) — numra-ode: SolverOptions gains two top-level Option<usize> fields with corresponding builders, max_order and min_order, which control the BDF adaptive-order window (clamped to [1, MAX_ORDER=5]). Default None preserves prior behavior exactly (BDF cap 5, floor 1) — every existing Bdf::solve call site produces byte-identical results. The new fields are consumed by BDF and ignored by every other solver, matching the established “applicable to some solvers, ignored by the rest” pattern of dense_output and events. Closes the F-SOLVER-FIELDS functional gap: Bdf::with_max_order / Bdf::fixed_order were silently non-functional at the static Solver::solve trait surface (the trait method discarded self and re-instantiated Bdf::new()); the wiring now flows through SolverOptions. Also closes a documented Foundation Specification §3.7 spec-vs-reality drift: the spec previously claimed BDF’s max order lived in SolverOptions while reality had it on the Bdf struct as a non-functional builder; the §3.7 paragraph is rewritten in the same PR to describe the now-real mechanism. Per Foundation Spec decision-log §6 #12, which records the rationale (option (b) of the original tension, top-level shape per §2.5’s centralized-configuration principle, sub-struct path explicitly deferred to a future foundation-design pass triggered by multi-family namespace pressure). Migration: SolverOptions::default().max_order(2) caps BDF order; .max_order(n).min_order(n) pins it once adaptive selection reaches n (BDF always starts at order 1). Regression pinned by numra-ode/tests/solver_options_order_regression.rs with revert-confirm-restore protocol applied — temporarily reverting the wiring in bdf.rs::solve_internal makes the pinning test fail (pinned == default n_accept proves the knob has no effect); restoring it makes the test pass (pinned ≫ 5 × default proves the wiring is real).

• F-OPTS (rustdoc only, no behavior change): the five divergent solver-family options structs (numra_sde::SdeOptions, numra_fde::FdeOptions, numra_ide::IdeOptions, numra_dde::DdeOptions, numra_spde::SpdeOptions) now carry paragraph-form rustdoc rationales explaining why each diverges from numra_ode::SolverOptions, per Foundation Spec §2.5’s “Configuration objects are shared across solver families” principle (which explicitly names F-OPTS as this documentation follow-up at §2.5’s audit baseline). Per-family rationales: SDE/SPDE carry stochastic-noise time control plus seed for reproducibility; FDE/IDE are inherently fixed-step (Caputo-derivative convolution memory; integral-kernel quadrature) so adaptive rtol/atol would be dead knobs; DDE adds discontinuity-propagation fields with no ODE analog and defaults dense_output to true because Method-of-Steps history interpolation requires it. Audit pass surfaced 2 additional sites beyond the entry’s named 3 (DdeOptions, SpdeOptions) — same audit-surfaces-more-than-named pattern as F-FD-STEP / F-CI-NODE20 / F-FD-CROSSCRATE / F-ERR (5-of-6 follow-ups now). Per-family claims pre-commit-verified against consuming solver code: SRA-family (numra-sde/src/sra.rs:128,296) confirmed reading options.rtol/atol for adaptive error scaling; DDE dense_output: true default confirmed at numra-dde/src/system.rs:87; SPDE adaptive: bool gates real branching at numra-spde/src/solver.rs:395 between fixed-step EM/Milstein dispatch and adaptive-variant code. The optional “consider whether the divergence is still justified” retrofit question in the original entry is not rolled into this close — DdeOptions is the most plausible retrofit candidate (closest shape to SolverOptions); if pursued it would open as a separate follow-up rather than being absorbed speculatively. Closes F-OPTS. Remaining trio items F-SOLVER-FIELDS and F-MATRIX-SHAPE are deferred foundation work tracked in docs/internal-followups.md, not on the 0.1.3 path.

Fixed

• F-WEBSITE-AUDIT-GATES (partial close — CI/infrastructure, not a crate change): the four website audit gates (Lighthouse (marketing site), Lighthouse (book), Accessibility (pa11y-ci), Playwright (dark-mode regression)) had never executed in the repo’s history before F-CI-NODE20’s PR #5 — every prior website-touching commit landed via direct-push to main, where their if: github.event_name == 'pull_request' guards keep them dormant. The investigative audit pass revealed three distinct failure modes, not the single “stale Lighthouse config” mode the follow-up entry pre-diagnosed. Fixed in this PR:

  • Config staleness on both Lighthouse configs (website/ci/lighthouserc.json, lighthouserc-book.json): the _comment_pwa, _comment_third_party_summary, and _comment_book_specific keys sitting inside assertions{} removed (LHCI parses any key in that object as an audit ID and reports “is not a known audit”). Three deprecated audit assertions removed (no-unload-listeners, no-vulnerable-libraries, uses-https — all removed upstream in Lighthouse 12; is-on-https is the surviving HTTPS-validation audit and remains). The render-blocking-resources maxNumericValue override removed (no-op in Lighthouse 12, which reshaped the audit to return a list of items — the preset’s maxLength is what actually fires). Top-level _comment fields in both configs expanded to document the audits that were removed upstream and to warn future readers that assertions{} keys are parsed as audit IDs.
  • Workflow URL list mismatch on the book Lighthouse job (.github/workflows/website.yml): the hardcoded chapter URL /ch01-fundamentals/numerical-stability/ doesn’t exist in the deployed book — the actual directory is ch01-introduction/ with no numerical-stability page, likely a holdover from an early outline. The 404 short-circuited Lighthouse-book’s run before any assertions could fire. Replaced with /ch01-introduction/installation/ (the hero-linked entry from the book homepage); a short inline comment documents the trailingSlash: 'always' convention from website/book/astro.config.mjs so the next URL-list edit doesn’t drop the trailing slashes.
  • Playwright wrong-expectation correction (website/tests/specs/dark-mode.spec.ts): three marketing dark-mode tests removed — two system-dark paint assertions and one stored-preference data-theme="dark" assertion. The marketing site is deliberately light-only by design; the removed assertions encoded a non-requirement that had been silently failing since the gate first ran on a pull_request event. The spec file’s docstring rewritten to state the intentional asymmetry with the book (which DOES support dark mode via Starlight and whose dark-mode tests are kept), and to forbid future contributors from re-introducing the wrong assertions for “consistency”. This is gate-correction-not-gate-silencing, same shape as removing the stale Lighthouse audit IDs above — the gate now honestly asserts the actual design intent.
  • Partial close, with four genuine-site-issue follow-ups opened: the audit surfaced that three of the four gates are catching real deployed-site regressions, not config drift. The marketing site has 78 WCAG2AA color-contrast violations across two pages (one Astro component instanced many times), is-crawlable: 0 site-wide on both the marketing site and the book (the entire Numra web presence is currently un-indexable by search engines), real CLS / render-blocking / image-delivery issues on marketing, and book-specific Lighthouse failures (label-content-name-mismatch, network-dependency-tree-insight, font-display-insight, plus lcp-discovery-insight + lcp-lazy-loaded on /ch13-performance/comparisons/) unmasked by this PR’s URL fix. Tracked as four new follow-ups in docs/internal-followups.md: F-WEBSITE-SEO (highest priority in the entire backlog — fixes site-wide discoverability across both subdomains; entry instructs that if the root cause is a trivial _headers / robots.txt config fix, ship as an expedited tiny PR rather than waiting for normal scheduling), F-WEBSITE-MARKETING-A11Y, F-WEBSITE-MARKETING-PERF, F-WEBSITE-BOOK-LHC-FIXES. label-content-name-mismatch appears in both the marketing-A11Y and book-LHC-FIXES entries; the two share an audit name but not a fix surface (custom Astro components vs. Starlight upstream templates) — both entries explicitly note the same-symptom, different-source relationship.
  • Release-independence note: this is CI configuration plus a website-workflow URL-list fix plus a Playwright test removal; it does not change any crate code and does not ship through cargo publish. The four split-out follow-ups are also website-quality work that lands on main and redeploys via Cloudflare Pages — release-independent in the same way. This entry sits in Unreleased because the convention is to record CHANGELOG-relevant main-branch changes here, not because it gates a crate release.
  • Audit’s evolution and the hard-stop discipline: the audit accurately diagnosed each gate’s failure from the PR #5 CI logs, but two corrections landed during pre-merge hard-stop review: (1) the URL-fix hypothesis (Lighthouse-book greens under the URL fix alone) was falsified by this PR’s own gate run, which unmasked book-specific Lighthouse failures previously hidden by the 404 short-circuit — handled by adding F-WEBSITE-BOOK-LHC-FIXES and rescoping F-WEBSITE-MARKETING-SEO → F-WEBSITE-SEO once is-crawlable: 0 was observed on the book as well; (2) the Playwright marketing dark-mode failure was initially audited as a “real regression” because the audit had no access to the design-intent input that the marketing site is deliberately light-only — surfaced by the user during pre-merge review, corrected here by removing the wrong assertions in-scope rather than opening a “regression” follow-up that didn’t exist. The investigative-audit framing was load-bearing: had this been treated as enumerable config-fixing work per the entry’s pre-diagnosis, the resulting PR would have “fixed” the gates by silencing them while leaving 78 a11y violations and a site-wide search-indexing block live in production. The hard-stop discipline was equally load-bearing: had the URL-fix and Playwright corrections not been gated through pre-merge review, falsifiable claims (a “Lighthouse-book expected to pass” prediction; a non-existent dark-mode regression) would have landed in PR history.
  • Gate-greenness on this PR’s own run: the Lighthouse-book gate is expected to remain failing under the in-scope fixes — the URL fix unmasks genuine book-side Lighthouse failures (tracked as F-WEBSITE-BOOK-LHC-FIXES). The Playwright gate is expected to pass under the test removal (the remaining tests — book dark, book light, marketing light — already passed in the prior run). The Lighthouse-marketing and pa11y gates still fail because they correctly detect the genuine site regressions; this PR is admin-merged over those three gates with per-gate documentation citing each split-out follow-up, same documented-exception discipline as F-CI-NODE20’s PR #5.

• F-WEBSITE-SEO (closed via reframe — CI/configuration, no crate change, no site fix): the follow-up entry’s “the entire Numra web presence is currently un-indexable by search engines” framing (cited in F-WEBSITE-AUDIT-GATES above and in the originally-opened F-WEBSITE-SEO entry) was based on a misread: the is-crawlable: 0 Lighthouse audit ran against PR preview deployments, which Cloudflare Pages deliberately injects X-Robots-Tag: noindex on to prevent duplicate-content with production. Direct checks against production on 2026-05-16 disconfirm the config-block hypothesis: curl -I against numra-rs.org/ and book.numra-rs.org/ returned HTTP 200 with no x-robots-tag header; zero <meta name="robots"> in either site’s built HTML; neither _headers file declares an X-Robots-Tag; website/site/public/robots.txt is Allow: /; Cloudflare’s documented preview-only noindex behavior explicitly excludes production custom domains. Production custom domains are not config-blocked from indexing as of those checks. Whether the production sites are actually indexed in practice (sitemap pickup, Search Console coverage, organic crawler discovery) is unverified and separately tracked — “config isn’t blocking” is necessary-not-sufficient, and that narrower question is held open as F-WEBSITE-SEO-VERIFY rather than collapsed into this retraction. The audit caught the misread at the hard-stop before any edit. In-scope fix: disabled is-crawlable in both Lighthouse preset configs (website/ci/lighthouserc.json, website/ci/lighthouserc-book.json) with a foreclosing _comment documenting the preview-noindex semantics and forbidding re-enablement without first implementing a production-URL audit path. Gate-correction-not-gate-silencing — same preview-only-false-signal shape as the marketing dark-mode Playwright test removal in F-WEBSITE-AUDIT-GATES above, though structurally narrower: the dark-mode case was a wrong property (marketing is light-only by design); the SEO case is the right property against the wrong URL set. Re-introducing the SEO assertion against production URLs would be correct (tracked as F-WEBSITE-SEO-PROD-PROBE); re-introducing dark-mode marketing assertions against any URL would not. The “highest priority in the entire backlog” / “every day of delay” / expedite-if-trivial framing is fully retracted: it was built on the demonstrated misread of preview-noindex as production-block, and the direct checks disconfirm that misread. Two narrower follow-ups forked from the audit’s residual observations and tracked in docs/internal-followups.md: F-WEBSITE-SEO-PROD-PROBE (medium — the unimplemented production-URL audit path the workflow comment at .github/workflows/website.yml:213-215 aspires to but doesn’t ship), F-WEBSITE-SEO-VERIFY (low — non-engineering go-look-at-the-world sanity check that production is actually indexed in practice, since “config isn’t blocking” is necessary-not-sufficient). Neither fork inherits the highest-priority framing. Honest framing of this correction: the F-WEBSITE-AUDIT-GATES entry above and its companion docs/internal-followups.md entry both amplified the un-indexable-in-production claim; both records are amended in the same pass as this retirement (docs/internal-followups.md gains an explicit Amendment 3 to F-WEBSITE-AUDIT-GATES retracting the production-block framing in the same record that made the claim, and F-WEBSITE-PR-FLOW’s worked-example list is recategorized so is-crawlable joins the dark-mode tests as a “preview-only false signal” example rather than a “gate catches real bug” example, with the structural narrower-distinction noted above). The original F-WEBSITE-AUDIT-GATES bullet’s “site-wide search-indexing block” phrasing earlier in this CHANGELOG is left verbatim with this cross-referencing retraction, rather than amended in place — preserves the audit trail of the original wrong claim, consistent with how Amendments 1 and 2 handle their own corrections. Pattern lesson: gates that run against preview URLs cannot assert properties whose semantics differ between preview and production — third worked example after F-WEBSITE-AUDIT-GATES’s stale-config and dark-mode-wrong-expectation corrections. The hard-stop audit discipline caught the misread before any code edit, any tiny-PR expedite, or any CHANGELOG claim of “site-wide search-indexing block fixed.” Symmetric-overclaim guard: this retraction was further reviewed against the inverse risk — replacing “production is blocked” with “production was always correctly indexable” would have been one unverified absolute swapped for another, an emotionally-satisfying shape a fourth correction in this arc was specifically vulnerable to. The replacement claim is stated at the precision the direct checks support (“not config-blocked from indexing, as of 2026-05-16 checks”); the in-practice question is held open via F-WEBSITE-SEO-VERIFY rather than closed by inference. Same necessary-not-sufficient discipline that caught the original preview-vs-production misread, applied one level deeper, to the correction itself. Release-independence note: this is CI configuration plus docs/internal-followups.md record corrections; no crate code changed, no cargo publish involved. The CI-config change’s practical effect is “the gate stops asserting a deterministically-failing thing on future PR previews,” not a production site change (production was not-config-blocked as of the direct checks on 2026-05-16; in-practice indexing separately tracked).

• F-FD-NOSCALE-BUG: Four FD utilities used a hardcoded h = 1e-8 step without scaling by |x|, which silently degraded gradient / Jacobian outputs for callers with |x| > ~5e7 (the precision floor where x + 1e-8 rounds back to x in f64). All four now use canonical precision-aware steps with additive scaling: central-FD sites get cbrt(S::EPSILON) * (1 + |x|); the forward-FD diffusion_derivative gets sqrt(S::EPSILON) * (1 + |x|).

  • Affected utilities: numra-optim::finite_diff_gradient (public free function), numra-optim::finite_diff_jacobian (public free function), numra-dde::History::evaluate_derivative (initial-history branch), numra-sde::SdeSystem::diffusion_derivative (trait default for the Milstein method). The first three are central FD; the fourth is forward FD (different canonical step).
  • Affected users: OptimProblem::solve() callers with large-magnitude parameters and no analytical gradient — the bug propagates through numra-optim::auto’s solver dispatcher (4 finite_diff_gradient / 1 finite_diff_jacobian call sites) and numra-optim::augmented_lagrangian’s inner loops (5 call sites); DDE callers querying derivatives in the initial-history region with large-magnitude time arguments; SDE Milstein-method users who don’t override diffusion_derivative.
  • Failure-mode shape: at |x| = 1e8, the previous unscaled formula returns either 0 (clean catastrophic case for the central-FD sites — both x+h and x-h round back to x) or a ~50%-distorted answer (the forward-FD site, where x+h rounds to x + ulp rather than x). Both regimes pass through the optimiser as silently-wrong gradient information. The fix recovers analytical accuracy in both regimes.
  • Pinned with four regression tests at |x| = 1e8: test_finite_diff_gradient_large_x_no_scaling_bug and test_finite_diff_jacobian_large_x_no_scaling_bug in numra-optim/src/problem.rs; test_evaluate_derivative_large_t_no_scaling_bug in numra-dde/src/history.rs; test_diffusion_derivative_default_large_x_no_scaling_bug in numra-sde/src/system.rs. Each test asserts proximity to the analytical value within 1e-3 relative tolerance — wide enough to be CI-stable across f64 platforms, narrow enough that the previous formula’s outputs (0 or ~50% off) fail by orders of magnitude. Structural-correctness check verified for the forward-FD site (revert formula → confirm test fails → restore).
  • Public-API rustdoc: finite_diff_gradient and finite_diff_jacobian rustdoc updated to document the canonical step formula, the additive-scaling rationale, and the |x| > ~5e7 precision floor — durable contract documentation for the next reader of the public API. The two trait/method sites (History::evaluate_derivative, SdeSystem::diffusion_derivative) were not annotated since they don’t have established rustdoc conventions for the FD-step choice.

This closes the no-scaling correctness portion of the FD-step audit. The 0.1.2 CHANGELOG’s F-FD-NOSCALE-BUG follow-up note over-listed numra-optim::robust as a no-scaling site; the F-FD-CROSSCRATE audit subsequently surfaced that those sites already have proper additive scaling and are F-FD-CROSSCRATE-class hygiene (hardcoded 1e-8, additive-scaled), not F-FD-NOSCALE-BUG-class correctness. The robust.rs sites remain a small leftover hygiene item, deferred.

Removed

• F-SOLVER-FIELDS (breaking removal) — numra-ode: deleted vestigial public surface that was silently non-functional under the pre-fix code path. All removed items had no correct behavior reachable through them prior to this PR; the breaking-ness is documented honestly but no correct program depended on them because no correct behavior was reachable. Removed items:
  • Bdf::new(), Bdf::with_max_order(n), Bdf::fixed_order(n) — the silently-non-functional instance builders; their max_order / min_order fields on the Bdf struct were unreadable from the static Solver::solve trait method that discarded self. Every caller of these builders got their order-control request silently ignored at solve time. The Bdf struct is now a unit struct; order control flows through SolverOptions::max_order / SolverOptions::min_order (see the additive ### Changed entry above).
  • Auto struct, Auto::new, Auto::with_hints(hints), impl Solver<S> for Auto, the Auto::hints field — same trait-discards-self pattern as the BDF builders. Auto::solve(...) and <Auto as Solver<S>>::solve(...) ignored the configured SolverHints and called SolverHints::new() internally. The Accuracy, Stiffness, and SolverHints types and the free-function user-facing API (auto_solve, auto_solve_with_hints) are unaffected — auto_solve and auto_solve_with_hints are now the only entry points to automatic solver selection. The previously-exported Auto symbol is removed from the public re-export list at numra-ode/src/lib.rs.
  • Justification — silent hint-loss producing potentially-inaccurate uncertainty quantification, not just vestigial surface: the deleted impl Solver<S> for Auto silently discarded any SolverHints configured by a caller — <Auto as Solver<S>>::solve(...) constructed SolverHints::new() internally and ignored anything Auto::with_hints(...) had been passed. The documented solve_with_uncertainty::<Auto, f64>(...) example (formerly in numra-ode/src/uncertainty.rs) therefore used auto-detected dispatch regardless of user intent. This can produce inaccurate uncertainty-quantification results when the user’s hints encoded problem structure that auto-detection’s heuristics don’t reconstruct (e.g., stiffness profiles the simple Jacobian-ratio detector at detect_stiffness misses, or accuracy targets the rtol-classifier categorizes differently than the user) — not unconditionally, only when hints carry auto-undetectable information. That conditional is the precisely-supported scope of the claim: a user who never built SolverHints and called auto_solve directly was unaffected; a user who built hints and passed Auto as a Solver<S> type parameter could be silently degraded. This is still a stronger removal justification than “vestigial”: the surface was structurally lossy in a way that could silently degrade scientific output for the documented use case, not merely unused. The doc example is updated to use ::<Tsit5> / ::<DoPri5> / ::<Radau5> (working solvers); users who want auto-selection should call auto_solve / auto_solve_with_hints directly on a single problem.
  • Migration:
    • Bdf::with_max_order(n) → SolverOptions::default().max_order(n)
    • Bdf::fixed_order(n) → SolverOptions::default().max_order(n).min_order(n)
    • Bdf::new() → Bdf (the struct is now a zero-size unit type; just construct it directly, or use the Bdf::solve static method as every workspace call site already does)
    • Auto::solve(problem, t0, tf, y0, options) → auto_solve(problem, t0, tf, y0, options)
    • Auto::with_hints(hints).solve(problem, t0, tf, y0, options) → auto_solve_with_hints(problem, t0, tf, y0, options, &hints)
    • solve_with_uncertainty::<Auto, f64>(...) → choose an explicit solver type parameter (Tsit5, DoPri5, Radau5, or any Solver<S> implementor); automatic selection within an uncertainty quantification call is not supported.
  • Foundation Spec §3.4 / §3.7 / §6 / §7 updated in the same PR: §3.4’s (a)/(b)/(c) design-tension paragraph rewritten to reflect the resolution; §3.7 drift closed (the spec’s claim that BDF max order lives in SolverOptions is now true); §6 decision-log entry #12 records the change with full rationale; §7 #8 records option (c) (reshape Solver::solve to &self) as a standing foundation open question — deferred, not foreclosed.

0.1.2 - 2026-05-15

Changed

• numra-ode: OdeSystem::jacobian default’s finite-difference step changed from a hardcoded 1e-8 to the textbook precision-aware sqrt(S::EPSILON) * (1 + |y_j|). The previous step was below f32::EPSILON ≈ 1.19e-7, so the FD perturbation quantised to zero on f32 and the default Jacobian came back as all-zeros — silently. The new form is correct on every Scalar precision (f64 lands at ≈1.49e-8, within ~50% of the prior value; f32 at ≈3.45e-4, no longer quantised). Behavioural impact on f64: Jacobian columns shift in the last 1-2 digits at states with |y_j| ≈ 1. No regression-suite test asserted on the prior values; the numra-ode test suite passes unchanged. Pinned with a new f32 regression test.
• numra-core: Signal::eval_derivative default’s central-difference step changed from a hardcoded 1e-8 (no scaling) to the textbook precision-aware cbrt(S::EPSILON) * (1 + |t|) — optimal for central FD by the truncation/round-off balance, and useful at every Scalar precision (f64 at ≈6.06e-6, f32 at ≈4.92e-3). Same f32 quantisation-to-zero failure mode as OdeSystem::jacobian; same fix shape. Most built-in Signal impls (Harmonic, Ramp, Sum, Product, Scaled, Constant, Zero) override with closed-form derivatives, so the FD path mainly affects Tabulated, FromFile, and user-supplied signals without an analytical override. Pinned with a new f32 regression test.
• numra-pde: six MOL hybrid analytical-Jacobian implementations updated in lockstep with the trait-default formula change (MOLSystem2D::jacobian, MOLSystem3D::jacobian, and the four ParametricMOLSystem{2,3}D::jacobian_y/_p). The reaction-FD diagonal in each cited the trait default in code comments — keeping the comment-vs-code consistency the comment claims required moving the implementations together. The 2D / 3D / parametric MOL analytical-vs-FD-agreement regression tests pass unchanged at their existing tolerances.
• numra-ode: Radau5 module-level rustdoc updated to cite the new FD step formula. The AugmentedSystem inline FD paths in numra-ode/src/sensitivity.rs were already on sqrt(S::EPSILON) from the original forward-sensitivity work; no change.
• numra-bench: pde_mol::FdJacobianMol2D wrapper takes no jacobian override, so it automatically tracks whatever the current trait default is. bench_mol2d_radau5_jacobian_path’s “FD baseline” therefore reports the actual FD path a user without an analytical override would hit today; the wrapper rustdoc is updated to make this explicit.
• Cross-crate finite-difference step formulas reconciled to the canonical precision-aware forms across 18 sites in numra-nonlinear::newton, numra-ode::{auto, dae_init, esdirk, index_reduction}, numra-ocp::adjoint, numra-fit::curve_fit, numra-core::uncertainty, and numra-optim::optim_sensitivity. Forward-FD sites adopt sqrt(S::EPSILON) * (1 + |x_j|); central-FD sites adopt cbrt(S::EPSILON) * (1 + |x_j|). Previously these used a mix of hardcoded 1e-7, 1e-8, or 1e-5 constants — three families of recipes scattered across the workspace. The audit pass surfaced 3 sites beyond the entry’s named 15 (two ReducedDaeSystem.fd_eps field initializers in numra-ode/src/index_reduction.rs and the compute_param_sensitivity default in numra-optim/src/optim_sensitivity.rs); same audit-surfaces-more-than-named pattern as F-FD-STEP and F-CI-NODE20.
  • Behavioural delta to flag for users: adjoint and DAE-initialisation paths previously used 1e-7 and now use sqrt(f64::EPSILON) ≈ 1.49e-8 — a 6.6× smaller forward-FD step. Adjoint gradients (numra-ocp::adjoint::adjoint_gradient and friends) and dae_init/index_reduction structure-detection paths shift in the last few digits at states with |x| ≈ 1. Numerically this is a roundoff/discretization improvement (canonical step balances both); the existing regression suites pass unchanged at their existing tolerances. Users depending on bit-exact reproducibility of adjoint gradients should rerun their reference values.
  • Behavioural delta in numra-fit::curve_fit: the multiplicative central-FD recipe (h = if |orig| > 1e-10 { orig * 1e-7 } else { 1e-7 }) is replaced by the canonical additive form (h = cbrt(S::EPSILON) * (1 + |orig|)). At |orig| ≈ 1 the new step is ~120× larger (~1.21e-5 vs the prior ~1e-7), with better roundoff/discretization balance. Levenberg–Marquardt convergence paths may shift in the first few iterations on problems sensitive to Jacobian noise; final fits remain within tolerance and the existing curve-fit regression tests pass unchanged.
  • Pinned with a new regression test (numra-optim::optim_sensitivity::tests::test_sensitivity_canonical_default_eps_at_large_param) exercising compute_param_sensitivity with eps = None at |p| = 100 — asserts the canonical default reproduces the analytical sensitivity of dx*/dp = 1 for min (x - p)^2 to within 1e-2.

This closes the cross-crate-consistency portion of the FD-step audit. One related follow-up remains open and is not addressed here:

• F-FD-NOSCALE-BUG — the no-scaling correctness bug in numra-optim::finite_diff_gradient, numra-optim::finite_diff_jacobian, numra-optim::robust, numra-dde::history, numra-sde::system::diffusion_derivative. These public-API sites use h = eps with no (1 + |x|) scaling, so f64 callers with large-magnitude states see the perturbation quantised. Independent of the f32 framing in F-FD-STEP and the consistency framing in F-FD-CROSSCRATE; warrants its own correctness analysis and PR.

Closes F-FD-STEP and F-FD-CROSSCRATE. Remaining follow-up tracked in docs/internal-followups.md.

Tooling

• CI workflow actions upgraded to Node.js 24-runtime majors ahead of GitHub’s 2026-06-02 deprecation deadline: actions/checkout@v4 → @v6, actions/setup-node@v4 → @v6, actions/upload-artifact@v4 → @v7, pnpm/action-setup@v3 → @v6, cloudflare/wrangler-action@v3 → @v4. pnpm/action-setup blocks also drop the redundant with: version: 9 in favour of the packageManager: [email protected] field in package.json as the single source of truth (required for the v4+ strict version-disagreement check). cloudflare/wrangler-action@v4 defaults wrangler to v4; pages deploy syntax is stable across the bump. Closes F-CI-NODE20.

0.1.1 - 2026-05-14

Changed (breaking)

• numra-core: NumraError is now #[non_exhaustive]. Exhaustive match arms over NumraError outside of numra-core must add a _ => … catch-all. Rationale: future-proofs additive variant changes, so subsequent cross-crate From impls can extend the workspace error story without forcing a semver-major bump. Workspace convention going forward: every public enum should be #[non_exhaustive] by default unless there’s a specific reason not to.
• numra-core: NumraError::Optimization variant renamed to NumraError::NumericalOptim. The wrapped type stays numra_core::OptimizationError; only the NumraError variant tag changes. Rationale: disambiguates from the new NumraError::Optim variant (which carries numra-optim::OptimError, the optimization crate’s API-level errors). The renamed variant clarifies its scope — low-level numerics-failure modes from inside algorithms (line search, descent direction, convergence) — distinct from crate-level optimization errors.

Added

• numra-core: ten new NumraError variants for cross-crate error propagation — Ode(String), Optim(String), Ocp(String), Fit(String), Signal(String), LineSearch(String), Interp(String), Integrate(String), Special(String), Stats(String). Each variant tags the source crate so callers can match on which subsystem failed; the per-crate error’s detail flattens into the carried String via Display. Resolves the workspace composability contract item 3 (errors compose into the workspace error type).
• numra-ode, numra-optim, numra-ocp, numra-fit, numra-signal, numra-nonlinear: From<CrateError> for NumraError impls covering SolverError, OptimError, OcpError, FitError, SignalError, and LineSearchError respectively. Six previously-missing bridges now land directly into NumraError, so ? propagation across these crate boundaries works without manual .map_err.
• numra/tests/interop_workflows.rs: workflow_ode_interp_integrate rewritten to return Result<(), NumraError> and use ? across three crate boundaries (numra-ode → numra-interp → numra-integrate). Canonical structural CI signal that the workspace ?-propagation property is real — the test compiles only if all three From impls are present. Closes the gap surfaced in the foundation-pass verification (finding D4): no interop test previously exercised cross-crate ?-propagation.

Changed

• numra-interp, numra-integrate, numra-special, numra-stats: existing From<…> for NumraError impls migrated from collapsing into NumraError::InvalidInput(String) to the new tagged-variant pattern (NumraError::Interp, NumraError::Integrate, NumraError::Special, NumraError::Stats). Workspace-uniform error story: every external-crate error now lands in a tagged variant rather than half-tagged-half-flattened. Display output for these errors changes (e.g. "Invalid input: …" → "Interpolation error: …"); Display is not part of any stability guarantee in 0.1.x.

0.1.0 - 2026-05-13

First public release. Archived on Zenodo: concept DOI 10.5281/zenodo.20159709 (all versions; preferred for citation); version DOI 10.5281/zenodo.20159710 (this 0.1.0 release; for reproducibility).

Added

• Workspace facade crate covering ODE, SDE, DDE, FDE, IDE, PDE, SPDE, optimization, optimal control, linear algebra, quadrature, interpolation, special functions, FFT, statistics, fitting, signal processing, and autodiff.
• Public release audit artifacts under docs/audit/, including API inventory, book coverage matrix, correctness test map, and release checklist.
• CI and local audit gates for formatting, clippy, tests, docs, MSRV, and supply-chain checks.
• SolverResult.dense_output: Option<DenseOutput<S>> field so callers can actually use the dense interpolant they requested via SolverOptions::dense().
• numra facade: top-level promotion of the canonical forward-sensitivity surface — solve_forward_sensitivity, solve_forward_sensitivity_with, ParametricOdeSystem, SensitivityResult — symmetric with the existing compute_sensitivities parameter-uncertainty API. Advanced building blocks (AugmentedSystem, ClosureSystem) remain reachable via the numra::ode submodule. The crate-level docstring distinguishes the two distinct sensitivity concepts (parameter-uncertainty vs ODE forward sensitivity) so users land on the right one from the doc index.
• numra/examples/robertson_sensitivity.rs: complete worked forward-sensitivity example on the canonical Robertson stiff kinetics benchmark. Demonstrates analytical state and parameter Jacobians (with the matching has_analytical_jacobian_* flags), Radau5-driven augmented integration, the full accessor surface (final_state, dyi_dpj, sensitivity_for_param, normalized_sensitivity_at), and dimensionless influence ranking. Internal sanity checks pass on output: state mass conservation to 12 digits and ∑ ∂y/∂k_i = 0 per parameter.
• numra-bench/benches/sensitivity.rs: publication-grade Criterion harness for forward sensitivity. Three groups: sensitivity_param_scaling (cost vs N_s on a 2-state oscillator), sensitivity_jacobian_mode (FD vs analytical-J_y vs analytical-J_y+J_p on Lotka–Volterra under Radau5), sensitivity_solver_choice (Radau5/BDF/DoPri5/Tsit5 on Robertson). 20s measurement / 3s warmup per point. Flattened JSON in numra-bench/results/sensitivity_*.json.
• website/figures/perf/sensitivity_{param_scaling,jacobian_mode,solver_choice}.py: three SPEC §5.2-compliant figure scripts. Each consumes its respective Criterion-flattened JSON and renders an SVG (with provenance sidecar) into website/book/src/assets/figures/perf_sensitivity_*.svg. make perf auto-includes them via the existing perf/*.py glob.
• Book: comprehensive forward-sensitivity chapter at ch11-uncertainty/sensitivity-analysis covering the math (sensitivity equation, simultaneous-corrector augmented system, initial sensitivity), the API (solve_forward_sensitivity{,_with}, ParametricOdeSystem trait with default-FD vs analytical Jacobians and the has_analytical_jacobian_* flag contract + debug-build safety net), the result layout (column-major over parameters), a Robertson worked example tying back to numra/examples/robertson_sensitivity.rs, solver-choice guidance with the bench chart embedded, performance characterisation with the param-scaling and Jacobian-mode bench charts embedded, and edge-case discussion (parameter scaling, identifiability, tolerance inheritance, non-trivial IC, closure-form perf). Also includes a forward-vs-adjoint trade-off table pointing at the adjoint chapter for large parameter counts.
• Book: ch11-uncertainty/parameter-importance sibling page holding the relocated compute_sensitivities (numra-core) content — local sensitivity coefficients of a scalar function via central FD. The two pages now cleanly separate scalar-function parameter importance from full ODE-trajectory forward sensitivity, with each page pointing at the other for the complementary use case. Sidebar order in ch11-uncertainty group: Error Propagation → Parameter Importance → Sensitivity Analysis → Interval Arithmetic → Monte Carlo with ODEs.
• numra-pde: MOLSystem3D wrapper for 3D method-of-lines on uniform tensor-product grids (numra-pde/src/mol3d.rs). Mirrors the MOLSystem2D API shape — heat, laplacian, with_operator, with_reaction, OdeSystem<S> impl, build_full_solution over column-major interior layout — backed by the new Operator3DCoefficients (a*u_xx + b*u_yy + c*u_zz + d*u_x + e*u_y + f*u_z + g*u) and assemble_operator_3d. assemble_laplacian_3d now delegates to the general assembler, removing duplicated stencil logic and matching how the 2D file is organised. Closes the audited 3D-MOL gap (the foundation pieces — Grid3D, BoundaryConditions3D, sparse 7-point Laplacian — were already shipping; only the wrapper was missing).
• numra-pde: 3D convenience builders HeatEquation3D::build, AdvectionDiffusion3D::build, ReactionDiffusion3D::build, and ReactionDiffusion3D::fisher (numra-pde/src/equations3d.rs), mirroring equations2d.rs exactly with a z-axis added.
• numra-pde: MOLSystem2D and MOLSystem3D both override OdeSystem::jacobian. The override is a direct CSC-to-row-major-dense copy of the assembled spatial operator (which already equals ∂(L[u])/∂u by construction) plus a diagonal-FD reaction term. The diagonal property holds rigorously because pointwise reactions cannot populate off-diagonal entries — R(t, x_i, ...; u_i) depends only on local state, so ∂R_i/∂u_m = 0 for m ≠ i. Implicit solvers (Radau5, Bdf) consume the override automatically with no opt-in. A future nonlocal / integro-PDE / multi-component reaction would not satisfy the diagonal property and falls through to the trait-default full FD path; that’s tracked as a named follow-up.
• numra-bench: new bench_jacobian_unification group in solvers.rs covering Robertson at rtol=1e-8, Van der Pol stiff at rtol=1e-6 for μ ∈ {10, 1000}, and a 2-state smooth-linear case — locking in regression coverage on the workload shapes most likely to surface alloc / extra-rhs overhead from the Jacobian unification work. New bench_mol2d_radau5_jacobian_path group in pde_mol.rs directly compares the analytical-Jacobian path against an FdJacobianMol2D wrapper that delegates rhs to MOLSystem2D but does not override jacobian, exercising both code paths in the same Criterion run.
• Book: 3D Problems section in ch04-beyond-odes/partial-des.md documenting MOLSystem3D, Operator3DCoefficients, the convenience builders, and the analytical-Jacobian path that ships with the MOL systems. Three worked examples (3D heat with the analytic exp(-3π²αt) decay, 3D Fisher–KPP via with_reaction, advection–diffusion via Operator3DCoefficients) and a memory/runtime scaling table for grids up to 65³.

Changed

• numra-ode: Radau5 and Bdf now route their Jacobian computation through OdeSystem::jacobian instead of a solver-inlined finite-difference path. Systems that override the trait method (e.g. MOLSystem2D, MOLSystem3D) get an analytical Jacobian for free; systems that don’t fall through to the canonical forward-FD default in numra-ode/src/problem.rs (eps = 1e-8, step = eps * (1 + |y_j|), row-major dense output). The unified FD step formula resolves prior drift between Radau5’s inlined eps * max(1, |y_j|) and the canonical (1 + |y_j|) form used by Bdf, the trait default, and the Hairer–Wanner reference; the smooth (1 + |y_j|) version has no discontinuity at |y_j| = 1 and is what every textbook / Hairer-style implementation specifies. Behavioural impact: Radau5 Newton residuals and Jacobian-rebuild step counts shift in the last 2-3 digits at states with |y_j| ≈ 1. No test in the regression suite asserted on those exact values; the entire numra-ode test suite (138 unit tests + integration tests) passes unchanged.

• numra-core: renamed Sensitivity → ParameterSensitivity and SensitivityResult → ParameterSensitivityResult (and the corresponding numra facade re-exports) to disambiguate from the new ODE forward-sensitivity types. The compute_sensitivities function name is preserved; only the type names changed.

• numra-ode: renamed the SensitivityEquations trait to ParametricOdeSystem and reshaped it to take params() and rhs_with_params(t, y, p, dydt) (which lets the trait drive its own FD on J_p); added FD-default jacobian_y and jacobian_p overrides so a minimal implementation supplies only n_states, n_params, params, and rhs_with_params. Removed the unused SensitivityState helper. AugmentedSystem now implements OdeSystem directly and can be solved by any Solver. Sensitivity flattening is column-major (parameter-major), matching CVODES indexing.

• numra-ode: SensitivityResult reshaped to flat row-major over time and column-major over parameters within each time block, with y_at, sensitivity_at, sensitivity_for_param, dyi_dpj, final_state, final_sensitivity, and normalized_sensitivity_at accessors.

• numra-ode: solve_trajectory (uncertainty.rs) reimplemented on top of the canonical solve_forward_sensitivity_with primitive. The private AugmentedOdeSystem is gone — single source of truth. Public API surface unchanged; the closure bound relaxed from Fn(...) + Send + Sync to Fn(...) (loosening, non-breaking). Added ParametricOdeSystem::has_analytical_jacobian_y / has_analytical_jacobian_p flag methods (default false) so AugmentedSystem’s hot path can inline FD with reused scratch buffers when the user hasn’t supplied analytical overrides — and call the override directly when they have. Document contract: implementors of jacobian_y / jacobian_p should also override the corresponding flag method to return true, otherwise their override is silently bypassed when running through solve_forward_sensitivity.

• numra-ocp: forward_sensitivity reimplemented as a thin wrapper over the canonical numra_ode::sensitivity::solve_forward_sensitivity_with primitive; the duplicated private augmented_rhs helper is gone. numra_ocp::SensitivityResult is now a re-export of numra_ode::SensitivityResult — single shape, single accessor surface, single test suite as the source of truth. The sensitivity layout therefore changes from row-major over states (s[i*N*N_s + state*N_s + param]) to column-major over parameters (s[i*N*N_s + param*N + state]), matching CVODES indexing. Callers using the typed accessors (sensitivity_at, y_at) are unaffected at the call site for state-1 / param-1 problems but should switch to dyi_dpj(time, state, param) and sensitivity_for_param(time, param) for forward-compatible multi-state usage. The t_eval parameter was renamed output_times in the function signature; positional callers are unaffected (Rust has no named-argument calls).

Performance

• numra-ode: solve_trajectory (uncertainty.rs) is substantially faster (~40% on the LV benchmark with 4 parameters at rtol=1e-9) after unification onto the canonical sensitivity primitive eliminated per-call heap allocations in the augmented RHS hot path. The previous private AugmentedOdeSystem allocated six fresh Vecs per RHS evaluation; the new path reuses scratch via RefCell buffers on AugmentedSystem. Magnitudes will vary with problem size, parameter count, and tolerance — the underlying mechanism (no per-call allocation) applies to all solve_trajectory workloads.
• numra-pde: stiff PDE workloads on Radau5 / Bdf are faster because MOLSystem2D and MOLSystem3D now supply analytical state Jacobians instead of forcing the solver into trait-default forward FD. Measured ~2.8% faster on the canonical stiff 2D heat-with-reaction workload (Allen-Cahn-like u_t = α∇²u + u - u³, α=0.5, 19²=361 interior points, rtol=1e-6); the win comes from skipping the N+1 rhs evaluations FD needs per Jacobian rebuild, but its size is bounded by how often the solver rebuilds the Jacobian (Radau5 caches across steps when Newton converges quickly). Magnitudes will vary with problem size, stiffness, grid resolution, and Jacobian-rebuild frequency — workloads with more rebuilds (looser tolerances triggering more rejected steps, larger N where the FD sweep cost grows, stiffer problems where Newton struggles more often) see proportionally larger speedups. Tracked permanently by bench_mol2d_radau5_jacobian_path in numra-bench/benches/pde_mol.rs.
• numra-ode: the OdeSystem::jacobian unification is itself perf-neutral on workloads that don’t override the trait method. Across bench_jacobian_unification (Robertson at rtol=1e-8, Van der Pol stiff at rtol=1e-6 for μ ∈ {10, 1000}, smooth 2D linear at rtol=1e-6), Radau5 deltas span −0.58% to +0.58%, all within Criterion’s noise threshold. BDF on Van der Pol stiff at the legacy rtol=1e-4 shows +3.45% / +1.70% / +0.74% for μ=10/100/1000 — well within the ≤5% regression bound and explained by one extra rhs evaluation per Jacobian rebuild (BDF previously passed an externally-computed f_0 into its inlined FD, while the trait default recomputes f_0 internally).

Forward-sensitivity v1 scope. What shipped this release: the ParametricOdeSystem trait with FD-default and analytical-Jacobian paths plus the has_analytical_jacobian_* flag contract and a debug-build consistency-check safety net; AugmentedSystem and ClosureSystem impls of OdeSystem; solve_forward_sensitivity and solve_forward_sensitivity_with entry points usable with any Solver (including Auto); flat row-major-over-time, column-major-over-parameters SensitivityResult with the typed accessor surface; ports of numra-ocp::forward_sensitivity and numra-ode::uncertainty::solve_trajectory onto the canonical primitive (no more parallel implementations); a 10-test regression suite at numra-ode/tests/sensitivity_regression.rs (eight numerical, two API-contract); the numra/examples/robertson_sensitivity.rs worked example; a publication-grade Criterion harness at numra-bench/benches/sensitivity.rs with three perf figures; and the rewritten ch11-uncertainty/sensitivity-analysis book chapter plus the parameter-importance sibling page. What did not ship and is tracked as named follow-ups: block-diagonal-aware factorisation (the augmented Jacobian is block_diag(J_y, …, J_y) and Radau5/BDF currently dense-LU the full N(N_s+1)² matrix); a JVP variant of the trait for matrix-free / sparse paths; an AD-based ParametricOdeSystem impl using numra-autodiff; staggered (CV_STAGGERED) and per-parameter staggered (CV_STAGGERED1) correction strategies; separate sensitivity tolerances; ergonomic alternatives to the has_analytical_jacobian_* flag (proc macro / typestate / trait redesign).

Fixed

• numra-ode: the blanket impl ParametricOdeSystem for &T forwards n_states, n_params, params, rhs_with_params, jacobian_y, jacobian_p, and initial_sensitivity, but did not forward has_analytical_jacobian_y / has_analytical_jacobian_p. Since those flags have a default false impl, every caller using the documented solve_forward_sensitivity::<_, _, _>(&system, ...) pattern was silently dispatching through the FD path even when the underlying system had analytical overrides — measurable slowdown on stiff problems where analytical Jacobians materially change runtime. The blanket impl now forwards the flag methods. Pinned by blanket_ref_impl_forwards_analytical_flags in the regression suite.
• numra-ode: AugmentedSystem now ships a debug-build consistency check that catches the related “user implemented jacobian_* analytically but forgot to override has_analytical_jacobian_*” misconfiguration on the first RHS call. The check evaluates the user’s jacobian_* and inline FD at the initial state and panics in debug builds (zero overhead in release) if the Frobenius-relative deviation exceeds 1e-3. Pinned by analytical_jacobian_y_flag_forgotten_panics. The required-flag-override contract is itself a known footgun; a follow-up tracks ergonomic alternatives (proc macro, typestate, trait redesign) for post-1.0.
• numra-ode: DoPri5 was building a DenseOutput when SolverOptions::dense() was set but never returning it, so the interpolant was silently dropped at the end of integration. The new SolverResult.dense_output field is now populated on both the normal exit and the early-termination event path.
• numra-ode: Radau5 step controller rewritten against Hairer–Wanner ODE II §IV.8 and the SciPy Radau reference. Eight bugs fixed (most importantly: the error_estimate forcing term used y instead of f(t,y); Newton’s initial guess now extrapolates the previous step’s collocation polynomial; LU is reused unless the step ratio leaves [1.0, 1.2]; Gustafsson predictive controller). Step counts on the standard reference suite are now within ~1.5–2× of SciPy’s, and Van der Pol μ=10 at rtol=1e-4 runs in ~0.66 ms (vs ~200 ms before, and ~12 ms for SciPy).
• numra-ode: SolverOptions::t_eval is now honored by every solver. Previously it was parsed, stored, cloned, and ignored — every solver returned its natural step grid regardless of what the caller requested, including solve_forward_sensitivity. Now each solver emits (t, y) pairs at exactly the requested times via Hermite cubic interpolation between accepted step endpoints (closes GitHub #1). Endpoints (t_eval[0] == t0, t_eval[-1] == tf) come back bit-exact; interior points are 3rd-order accurate locally, which dominates the error budget only for very-high-order solvers (Vern8) at large step sizes. Backward integration and out-of-range / unsorted grids are validated up front. Pinned by numra-ode/tests/t_eval_regression.rs (covers DoPri5, Tsit5, Vern6/7/8, Radau5, BDF, Esdirk32/43/54, plus solve_forward_sensitivity and edge cases: bit-exact endpoints, dense grids that span multiple steps, backward integration).
• numra-ode: BDF rewritten as a Rust port of SciPy’s BDF (the canonical Shampine–Reichelt NDF / ode15s algorithm). Three independent algorithmic bugs fixed: (1) the Newton convergence test now compares the correction norm against tolerance rather than the residual norm — the previous test silently declared convergence at iteration 0 with y_new == y_predict, producing the reported “success=true with y_final == y_initial” mode; (2) the local-truncation-error estimate now uses the (order+1)-th modified divided difference (the cumulative Newton correction d), not a stale history[0] reference that always evaluated to ~0; (3) BDF coefficients are applied to a modified-divided-differences table that is rescaled in place when h changes, restoring formal order accuracy under variable steps. Newton failure with a stale Jacobian now retries at the same h after refreshing J (previously it shrank h on every Newton failure, which combined with bug 1 drove h to zero on benign problems). Step counts now match SciPy within ~1× on smooth problems and ~2× on stiff ones; the regression test test_bdf_regression_silent_wrong_answer pins the silent-corruption mode out.

Policy

• Root Cargo.lock is committed for reproducible public CI.
• The user-facing book is built from website/book/ (Astro + Starlight) and deployed to book.numra-rs.org by the Website workflow.