Calculations with dual numbers for computing derivatives.
If I is the identity matrix and J is the matrix with 1's above the main diagonal and 0's everywhere else, then a dual number has the form
a = a0 I + a1 J + a2 J2 + ...
For sufficiently smooth functions
f(xI + J) = f(x) I + f'(x) J + f"(x)/2! J2 + ...
This allows us compute derivatives without taking limits of difference quotients.
(xI + J)3 = x3 I + 3 x2 J + 3 x J2 + J3.
We can read off the derivatives of f(x) = x3
: f'(x) = 3 x2
and f"(x)/2! = 3 x, so f"(x) = 6x and the third derivative is 3!.
This can be written in C++ as
double x = 0.5;
dual::number<double, 3> X(x, 1), Y;
Y = X*X*X;
ensure (Y._(0) == x*x*x);
ensure (Y._(1) == 3*x*x);
ensure (Y._(2) == 6*x);
Yes, those really are =='s. Dual numbers automatically calculate derivatives to bit-precision, not just machine precision.