A primer on numerical differentiation

In order to numerically evaluate a derivative $y'(x)=dy/dx$ at point $x_0$, we approximate is by using finite differences: Therefore we find: $$\begin{darray}{rcl} && dx \approx \Delta x &=&x_1-x_0, \\ && dy \approx \Delta y &=&y_1-y_0 = y(x_1)-y(x_0) = y(x_0+\Delta_x)-y(x_0),\end{darray}$$

Then we re-write the derivative in terms of discrete differences as: $$\frac{dy}{dx} \approx \frac{\Delta y}{\Delta x}$$

Example

Let's look at the accuracy of this approximation in terms of the interval $\Delta x$. In our first example we will evaluate the derivative of $y=x^2$ at $x=1$.

Why is it that the sequence does not converge? This is due to the round-off errors in the representation of the floating point numbers. To see this, we can simply type:

Let's try using powers of 1/2

In addition, one could consider the midpoint difference, defined as: $$dy \approx \Delta y = y(x_0+\frac{\Delta_x}{2})-y(x_0-\frac{\Delta_x}{2}).$$

For a more complex function we need to import it from the math module. For instance, let's calculate the derivative of $sin(x)$ at $x=\pi/4$, including both the forward and midpoint differences.

A more in-depth discussion about round-off errors in numerical differentiation can be found here

Special functions in numpy

numpy provides a simple method diff() to calculate the numerical derivatives of a dataset stored in an array by forward differences. The function gradient() will calculate the derivatives by midpoint (or central) difference, that provides a more accurate result.

Notice above that gradient() uses forward and backward differences at the two ends.

More discussion about numerical differenciation, including higher order methods with error extrapolation can be found here.

The module scipy also includes methods to accurately calculate derivatives:

One way to improve the roundoff errors is by simply using the decimal package