Differentiation

Numerical derivatives (diff, diffs)

mpmath.calculus.diff(f, x, n=1, method='step', scale=1, direction=0)

Numerically computes the derivative of f, f'(x). Optionally, computes the n-th derivative, f^{(n)}(x), for any order n.

Basic examples

Derivatives of a simple function:

>>> from mpmath import *
>>> mp.dps = 15
>>> print diff(lambda x: x**2 + x, 1.0)
3.0
>>> print diff(lambda x: x**2 + x, 1.0, 2)
2.0
>>> print diff(lambda x: x**2 + x, 1.0, 3)
0.0

The exponential function is invariant under differentiation:

>>> nprint([diff(exp, 3, n) for n in range(5)])
[20.0855, 20.0855, 20.0855, 20.0855, 20.0855]

Method

One of two differentiation algorithms can be chosen with the method keyword argument. The two options are 'step', and 'quad'. The default method is 'step'.

'step':

The derivative is computed using a finite difference approximation, with a small step h. This requires n+1 function evaluations and must be performed at (n+1) times the target precison. Accordingly, f must support fast evaluation at high precision.

'quad':

The derivative is computed using complex numerical integration. This requires a larger number of function evaluations, but the advantage is that not much extra precision is required. For high order derivatives, this method may thus be faster if f is very expensive to evaluate at high precision.

With 'quad' the result is likely to have a small imaginary component even if the derivative is actually real:

>>> print diff(sqrt, 1, method='quad')
(0.5 - 9.44048454290863e-27j)

Scale

The scale option specifies the scale of variation of f. The step size in the finite difference is taken to be approximately eps*scale. Thus, for example if f(x) = \cos(1000 x), the scale should be set to 1/1000 and if f(x) = \cos(x/1000), the scale should be 1000. By default, scale = 1.

(In practice, the default scale will work even for \cos(1000 x) or \cos(x/1000). Changing this parameter is a good idea if the scale is something preposterous.)

If numerical integration is used, the radius of integration is taken to be equal to scale/2. Note that f must not have any singularities within the circle of radius scale/2 centered around x. If possible, a larger scale value is preferable because it typically makes the integration faster and more accurate.

Direction

By default, diff() uses a central difference approximation. This corresponds to direction=0. Alternatively, it can compute a left difference (direction=-1) or right difference (direction=1). This is useful for computing left- or right-sided derivatives of nonsmooth functions:

>>> print diff(abs, 0, direction=0)
0.0
>>> print diff(abs, 0, direction=1)
1.0
>>> print diff(abs, 0, direction=-1)
-1.0

More generally, if the direction is nonzero, a right difference is computed where the step size is multiplied by sign(direction). For example, with direction=+j, the derivative from the positive imaginary direction will be computed.

This option only makes sense with method=’step’. If integration is used, it is assumed that f is analytic, implying that the derivative is the same in all directions.

mpmath.calculus.diffs(f, x, n=mpf('+inf'), method='step', scale=1, direction=0)

Returns a generator that yields the sequence of derivatives

f(x), f'(x), f''(x), \ldots, f^{(k)}(x), \ldots

With method='step', diffs() uses only O(k) function evaluations to generate the first k derivatives, rather than the roughly O(k^2) evaluations required if one calls diff() k separate times.

With n < \infty, the generator stops as soon as the n-th derivative has been generated. If the exact number of needed derivatives is known in advance, this is further slightly more efficient.

Examples

>>> nprint(list(diffs(cos, 1, 5)))
[0.540302, -0.841471, -0.540302, 0.841471, 0.540302, -0.841471]
>>> for i, d in zip(range(6), diffs(cos, 1)): print i, d
...
0 0.54030230586814
1 -0.841470984807897
2 -0.54030230586814
3 0.841470984807897
4 0.54030230586814
5 -0.841470984807897