Automatic Differentiation (aka AD) is a family of methods that are used to evaluate the derivative of a coded function. These methods are far more accurate than finite differences, since they are theoretically exact in the absence of floating point roundoff error.
AD is, however, subtly different than symbolic differentiation. The key difference here is that computer algebra systems (like Mathematica) will return a function that can be evaluated for a derivative. In many implementations of Automatic differentiation, however, the derivative is computed as a side-effect of evaluating the original function.
There are a number of different ways to implement AD, and they are outlined in a reasonable amount of detail on the wikipedia page.
https://en.wikipedia.org/wiki/Automatic_differentiation
Simple Implementation Method
Perhaps the most straightforward method to implement AD is using "Dual Numbers" and operator overloading. The key idea is to introduce a new type of scalar number, called a Dual Number. The dual number augments the algebra of real numbers by replacing every number $x$ with $x + x' \epsilon$, where $\epsilon$ has the property $\epsilon^2 = 0$.
Using this, you can define the algebra of dual numbers as you might expect:
\begin{align}
(u+u'\epsilon) + (v + v'\epsilon) &= (u+v) + (u' + v')\epsilon \\
(u+u'\epsilon) - (v + v'\epsilon) &= (u-v) + (u' - v')\epsilon \\
(u+u'\epsilon) * (v+v'\epsilon) &= (uv) + (uv' + vu')\epsilon \\
(u + u'\epsilon) / (v+v'\epsilon) &= (u/v) + \frac{u'v - uv'}{v^2}\epsilon
\end{align}
See the wikipedia page for extensions of this to other functions ($\sin$, $\cos$, $\exp$, ...)
The astute reader will notice that the $\epsilon$ component of the result of any operation between dual numbers implements the chain rule for that operation. Due to this, any sequence of operations (i.e. a coded function) will simultaneously compute the function value and its derivative. The only loss of accuracy that one incurs using this method is floating-point roundoff errors (so you still need to use numerically stable algorithms).
On the software side, the primary burden is the implementation of a Dual Number class and to overload any operators/functions that you want to use.
In addition, any functions that you wish to differentiate should be (in C++ parlance) templated on the input type. Eg, rather than:
double f(double x) { ... }
you would write
template<typename T>
T f(T x) {...}
Doing so will allow you to re-use the same code to do normal function calls or to do the AD function calls.
Example Code
Here is a link to my C++ implementation of the method that I use in my finite element library. The code is free to take/use/modify if you'd like.
You'll have to strip out the non-standard header and the namespace open/close lines, but it shouldn't have any dependencies to the rest of the library.
https://github.com/tjolsen/YAFEL/blob/master/include/utils/DualNumber.hpp
Final Thoughts
The ideas described here can be generalized to vector-valued inputs/outputs. By doing so, you can compute Jacobians without having to derive what each partial derivative should be. Depending on the dimension of the inputs and outputs, other methods than the Dual Number method described here may be more efficient, so it's worth reading about in more detail.
