Numerical derivatives, integrals

[sharkdp]

It should actually be fairly easy to implement a simple version of this (if we disregard all of the real world numerical problems that can arise). The major blocker for this is the fact that we don't have function types and corresponding function objects (#159), so we cannot simply write something like diff(sqrt, 1.0)

With a hypothetical syntax for function types Fn[(X1, X2, X3, …) -> X_return], we could even implement numerical diff in Numbat itself using the new builtin unit_of function:

let eps = 1e-10

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

Obviously, we would like to have something a little bit more sophisticated for production though.

For a concrete function, this already works:

fn diff_sqrt<X>(x: X) -> X^(-1/2) = (sqrt(x + eps · unit_of(x)) - sqrt(x)) / (eps · unit_of(x))

assert_eq(diff_sqrt(1), 1/2, 1e-5)
assert_eq(diff_sqrt(2), 1/2 × 2^-0.5, 1e-5)

[sharkdp]

With #159 merged, we can now do this:

fn diff(f: Fn[(Scalar) -> Scalar], x: Scalar) -> Scalar =
    (f(x + eps) - f(x)) / eps

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

We still can't do the generic version though, because our handling of generics in function types is not yet powerful enough. We can't even pass a generic function like sqrt to the non-generic diff function.

[sharkdp]

The generic version of this is now also supported (see #451), so we can do things like

use numerics::diff

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)