Genericization and Autodiff #4

2026-09-29T18:23:55-04:00

connor commented

2026-09-29 18:23:55 -04:00

Generic scalar type and forward-mode autodiff

Makes the solver generic over its scalar so a problem can be solved in f64 or
differentiated by solving it in dual numbers. Also fixes a step-controller bug
found on the way.

Autodiff

New scalar::Scalar trait (RealField + Copy + NumOps<f64, Self>, plus re()
and lift()). Implemented for f64, and for num-dual's forward-mode dual
numbers under the autodiff feature (on by default). State, time and step size
are all generic; all four integrators and both controllers are ported.

Solving at a dual scalar gives exact derivatives through the whole solve — the
initial state, t0, t_end, and any parameter the right-hand side reads. One
propagation in DualSVec64<6> yields a full 6x6 state transition matrix.

The only change a user needs is <T: Scalar> on the right-hand side signature;
the body is usually untouched, and Params can stay f64 unless you want to
differentiate with respect to it. A concrete f64 right-hand side still works
unchanged — the benches are the proof, they compile as they were.

Two things are deliberately not differentiated, because differentiating them
would be wrong rather than merely expensive: the adaptive step sequence (a
discrete accept/reject decision, as in DifferentialEquations.jl) and event
location (the event time is defined implicitly, so its derivative comes from the
implicit function theorem, not from differentiating Brent's iterations).

Step size controller fix

The rejection loop in Problem::solve judged a retried step by the previous
attempt's error, so a step could be accepted on the strength of a different
step's error. Its size depended on how often a method was rejected, so it hit
methods unequally: on the 6D orbit problem at 1e-10, DP45 took 2,023 steps
instead of 440, while Vern7 was mildly flattered (152 against 184).

Error now tracks tolerance cleanly: the harmonic oscillator gives 1.79e-6 at
1e-6, 1.85e-8 at 1e-8, 1.86e-10 at 1e-10. Both benchmark reports are
regenerated, which changes some headline ratios — Vern7 vs DP5 on the 6D orbit
is 1.1x, not 8.8x.

Breaking changes

Integrator::step returns StepResult { y, err, dense } instead of a
(y, err, dense) tuple
Callback is now a trait; FnCallback holds the closure pair
Problem is generic over C: Controller<T>, which is what makes
PIDController usable — it previously could not be passed to a Problem
PIDController's beta1/beta2/beta3 are now p_term/i_term/d_term
a_tol, a_tol_full, r_tol and with_callback are consuming builders
nalgebra 0.34 → 0.35 (required by num-dual)

Other

cargo clippy --all-targets -- -D warnings is clean (was 25 warnings)
readme is included as crate docs, so every example in it is compiled and run
DormandPrince45 gained a docstring; it notes that FSAL is a property of the
method that neither it nor BS3 currently exploits
roadmap: new feature 39 (per-step allocation), note added to feature 19
(Jacobian computation) on how autodiff arrived by a different route

Verification

43 unit tests + 7 doctests; 38 + 5 with --no-default-features. AD checked
three ways: machine precision against an exactly-integrable quadrature
(y' = p·t^4, matching to 1e-14), against the analytic sensitivity of a linear
problem, and on the orbit problem via det(STM) = 1 (Liouville) plus a
central-difference cross-check. Vern7's lazy extra stages are exercised in dual
arithmetic.

f64 numerics are bit-identical to v0.3.1 for the genericisation itself —
verified by step count and final state on three integrators before the
controller fix was applied.

## Generic scalar type and forward-mode autodiff Makes the solver generic over its scalar so a problem can be solved in `f64` or differentiated by solving it in dual numbers. Also fixes a step-controller bug found on the way. ### Autodiff New `scalar::Scalar` trait (`RealField + Copy + NumOps<f64, Self>`, plus `re()` and `lift()`). Implemented for `f64`, and for `num-dual`'s forward-mode dual numbers under the `autodiff` feature (on by default). State, time and step size are all generic; all four integrators and both controllers are ported. Solving at a dual scalar gives exact derivatives through the whole solve — the initial state, `t0`, `t_end`, and any parameter the right-hand side reads. One propagation in `DualSVec64<6>` yields a full 6x6 state transition matrix. The only change a user needs is `<T: Scalar>` on the right-hand side signature; the body is usually untouched, and `Params` can stay `f64` unless you want to differentiate with respect to it. A concrete `f64` right-hand side still works unchanged — the benches are the proof, they compile as they were. Two things are deliberately not differentiated, because differentiating them would be wrong rather than merely expensive: the adaptive step sequence (a discrete accept/reject decision, as in DifferentialEquations.jl) and event location (the event time is defined implicitly, so its derivative comes from the implicit function theorem, not from differentiating Brent's iterations). ### Step size controller fix The rejection loop in `Problem::solve` judged a retried step by the *previous* attempt's error, so a step could be accepted on the strength of a different step's error. Its size depended on how often a method was rejected, so it hit methods unequally: on the 6D orbit problem at 1e-10, DP45 took 2,023 steps instead of 440, while Vern7 was mildly flattered (152 against 184). Error now tracks tolerance cleanly: the harmonic oscillator gives 1.79e-6 at 1e-6, 1.85e-8 at 1e-8, 1.86e-10 at 1e-10. Both benchmark reports are regenerated, which changes some headline ratios — Vern7 vs DP5 on the 6D orbit is 1.1x, not 8.8x. ### Breaking changes - `Integrator::step` returns `StepResult { y, err, dense }` instead of a `(y, err, dense)` tuple - `Callback` is now a trait; `FnCallback` holds the closure pair - `Problem` is generic over `C: Controller<T>`, which is what makes `PIDController` usable — it previously could not be passed to a `Problem` - `PIDController`'s `beta1`/`beta2`/`beta3` are now `p_term`/`i_term`/`d_term` - `a_tol`, `a_tol_full`, `r_tol` and `with_callback` are consuming builders - nalgebra 0.34 → 0.35 (required by `num-dual`) ### Other - `cargo clippy --all-targets -- -D warnings` is clean (was 25 warnings) - readme is included as crate docs, so every example in it is compiled and run - `DormandPrince45` gained a docstring; it notes that FSAL is a property of the method that neither it nor BS3 currently exploits - roadmap: new feature 39 (per-step allocation), note added to feature 19 (Jacobian computation) on how autodiff arrived by a different route ### Verification 43 unit tests + 7 doctests; 38 + 5 with `--no-default-features`. AD checked three ways: machine precision against an exactly-integrable quadrature (`y' = p·t^4`, matching to 1e-14), against the analytic sensitivity of a linear problem, and on the orbit problem via `det(STM) = 1` (Liouville) plus a central-difference cross-check. Vern7's lazy extra stages are exercised in dual arithmetic. `f64` numerics are bit-identical to v0.3.1 for the genericisation itself — verified by step count and final state on three integrators before the controller fix was applied.

connor added 4 commits 2026-09-29 18:23:55 -04:00

Incorporated the other features: genericization and autodiff mainly f8a47d0126

Updated docs on benchmarks db14a02df4

Fixed a small bug on the step sizing c20a090036

Updated reports after bug fix 2053aaa68e

connor merged commit 788f93cccb into main

2026-09-29 18:24:02 -04:00

Sign in to join this conversation.

1 Participants

Notifications

Due Date

No due date set.

Dependencies

No dependencies set.

Reference: connor/ordinary-diffeq#4