Incorporated the other features: genericization and autodiff mainly

readme.md

@@ -27,6 +27,8 @@ See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analy
 - **Dense output interpolation** at any time point
 - **Parameters** in derivative and callback functions
 - **Lazy computation** of extra interpolation stages (Vern7)
+- **PID step size control** as an alternative to the PI controller
+- **Automatic Differentiation**, the same problem can be solved with floats or also differentiated with dual numbers
 
 ### Future Improvements
 
@@ -68,15 +70,77 @@ let ode = ODE::new(&derivative, 0.0, 6.3, y0, params);
 let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6);
 let controller = PIController::default();
 
-let value_too_high = Callback {
-    event: &|t: f64, _y: Vector2<f64>, _p: &Params| { 5.0 - t },
+let value_too_high = FnCallback {
+    event: &|t: f64, _y: &Vector2<f64>, _p: &Params| { 5.0 - t },
     effect: &stop,
 };
 
 // Solve the problem
-let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high);
+let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high);
 let solution = problem.solve();
 
 // Can interpolate solutions to whatever you want
 let _interpolated_answer = solution.interpolate(4.4);
 ```
+
+## Automatic differentiation
+
+State, time and step size are generic over the `Scalar` trait. `f64` implements it by default,
+and so do the forward-mode dual numbers of [`num-dual`](https://crates.io/crates/num-dual) (the
+`autodiff` feature, on by default). Solving at a dual scalar propagates exact derivatives alongside
+the state.
+
+The only change needed is to the signature of the right-hand side. Here is the pendulum from above,
+before and after:
+
+```rust,ignore
+fn derivative(_t: f64, y: Vector2<f64>, p: &Params) -> Vector2<f64> {
+    let &(g, l) = p;
+    Vector2::new( y[1], -(g/l) * y[0].sin() )
+}
+
+fn derivative<T: Scalar>(_t: T, y: Vector2<T>, p: &Params) -> Vector2<T> {
+    let &(g, l) = p;
+    Vector2::new( y[1], -y[0].sin() * (g/l) )
+}
+```
+
+`Params` stays `(f64, f64)` — parameters only need to become `T` themselves if you want to
+differentiate *with respect to* them. Real constants mix into `T` arithmetic directly; `Scalar::lift`
+is only needed to scale a vector by a real, or to build a `T` out of nothing.
+
+With that one change, the same function can be solved or differentiated:
+
+```rust
+use nalgebra::{Matrix2, Vector2};
+use num_dual::DualSVec64;
+use ordinary_diffeq::prelude::*;
+use std::f64::consts::PI;
+
+type Params = (f64, f64);
+
+fn derivative<T: Scalar>(_t: T, y: Vector2<T>, p: &Params) -> Vector2<T> {
+    let &(g, l) = p;
+    Vector2::new(y[1], -y[0].sin() * (g / l))
+}
+
+// One propagation in `DualSVec64<2>` yields the whole 2x2 Jacobian of the
+// final state with respect to the initial one.
+let (y_end, jacobian): (Vector2<f64>, Matrix2<f64>) = num_dual::jacobian(
+    |y0: Vector2<DualSVec64<2>>| {
+        type T = DualSVec64<2>;
+        let ode = ODE::new(&derivative, T::lift(0.0), T::lift(1.0), y0, (9.81, 1.0));
+        let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
+        let mut problem = Problem::new(ode, vern7, PIController::default());
+        let solution = problem.solve();
+        *solution.states.last().unwrap()
+    },
+    &Vector2::new(0.0, PI / 2.0),
+);
+
+// The pendulum is Hamiltonian, so its flow is area-preserving
+assert!((jacobian.determinant() - 1.0).abs() < 1e-9);
+```
+
+These are exact derivatives of the scheme as it ran, at machine precision, with no finite
+differencing and no step size to tune.