numbat/examples/numerical_diff.nbt

let eps = 1e-10

fn diff<X: Dim, Y: Dim>(f: Fn[(X) -> Y], x: X) -> Y / X =
  (f(x + eps · unit_of(x)) - f(x)) / (eps · unit_of(x))

assert_eq(diff(log, 2.0), 0.5, 1e-5)
assert_eq(diff(sin, 0.0), 1.0, 1e-5)

assert_eq(diff(sqrt, 1.0), 0.5, 1e-5)

fn f(x: Scalar) -> Scalar = x * x + 4 * x + 1

assert_eq(diff(f, 2.0), 8.0, 1e-5)

fn dist(t: Time) -> Length = 0.5 g0 t^2
fn velocity(t: Time) -> Velocity = diff(dist, t)

assert_eq(velocity(2.0 s), 2.0 s × g0, 1e-3 m/s)