Mpmath: Numerical differentiation

First derivative, using the complex-step derivative approximation

ctx.deriv1_c(f, x0, solver='secant', tol=None, verbose=False, verify=True, **kwargs)

where ctx is dec, mpm, or gmp.

See also: https://nhigham.com/2020/10/06/what-is-the-complex-step-approximation/

See Lai and Crassidis

The complex-step derivative approximation can be derived by approximating a nonlinear function with a complex variable using a Taylor’s series expansion:

\[f(x+ih) = f(x) + ihf'(x) - \frac{h^2}{2!}f''(x) - i \frac{h^3}{3!}f^{(3)}(x) + \frac{h^4}{4!}f^{(4)}(x) + \cdots\]

Taking only the imaginary parts of both sides, dividing by \(h\) and rearranging gives

\[f'(x) = \frac{\Im \left ( f(x+ih) \right )}{h} + O(h^2)\]

Terms with order \(h^2\) or higher can be ignored since the interval \(h\) can be chosen up to machine precision. Thus, to within first order the complex-step derivative approximation is given by

\[f'(x) = \frac{\Im \left ( f(x+ih) \right )}{h}, E_{\text{trunc}}(h) = \frac{h^2}{6}f^{(3)}(x)\]

where \(E_{\text{trunc}}(h)\) denotes the truncation error. Note that this solution is not a function of differences, which ultimately provides better roundoff characteristics than a standard finite difference.

Second derivative, using the complex-step derivative approximation

ctx.deriv2_c(f, x0, solver='secant', tol=None, verbose=False, verify=True, **kwargs)

where ctx is dec, mpm, or gmp.

See also: https://nhigham.com/2020/10/06/what-is-the-complex-step-approximation/

See Lai and Crassidis

Gradient, using the complex-step derivative approximation

ctx.gradient_c(f, x0, solver='secant', tol=None, verbose=False, verify=True, **kwargs)

where ctx is dec, mpm, or gmp.

See also: https://nhigham.com/2020/10/06/what-is-the-complex-step-approximation/

The gradient of a vector function is a simple extension of the scalar case.

Jacobi matrix, using the complex-step derivative approximation

ctx.jacobi_c(f, x0, solver='secant', tol=None, verbose=False, verify=True, **kwargs)

where ctx is dec, mpm, or gmp.

See also: https://nhigham.com/2020/10/06/what-is-the-complex-step-approximation/

The Jacobian of a vector function is a simple extension of the scalar case. This Jacobian is defined by

\[\begin{split}F_x = \begin{pmatrix} \frac{\partial f_1(\boldsymbol{x})}{\partial x_1} & \frac{\partial f_1(\boldsymbol{x})}{\partial x_2} & \cdots & \frac{\partial f_1(\boldsymbol{x})}{\partial x_p} & \cdots & \frac{\partial f_1(\boldsymbol{x})}{\partial x_n} \\ \frac{\partial f_2(\boldsymbol{x})}{\partial x_1} & \frac{\partial f_2(\boldsymbol{x})}{\partial x_2} & \cdots & \frac{\partial f_2(\boldsymbol{x})}{\partial x_p} & \cdots & \frac{\partial f_2(\boldsymbol{x})}{\partial x_n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{\partial f_q(\boldsymbol{x})}{\partial x_1} & \frac{\partial f_q(\boldsymbol{x})}{\partial x_2} & \cdots & \frac{\partial f_q(\boldsymbol{x})}{\partial x_p} & \cdots & \frac{\partial f_q(\boldsymbol{x})}{\partial x_n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \frac{\partial f_m(\boldsymbol{x})}{\partial x_1} & \frac{\partial f_m(\boldsymbol{x})}{\partial x_2} & \cdots & \frac{\partial f_m(\boldsymbol{x})}{\partial x_p} & \cdots & \frac{\partial f_m(\boldsymbol{x})}{\partial x_n} \\ \end{pmatrix}\end{split}\]

The complex-step approximation of the Jacobian is defined by

\[\begin{split}F_x = \frac{1}{h} \Im \begin{pmatrix} f_1(\boldsymbol{x} + ih\boldsymbol{e}_1) & f_1(\boldsymbol{x} + ih\boldsymbol{e}_2) & \cdots & f_1(\boldsymbol{x} + ih\boldsymbol{e}_q) & \cdots & f_1(\boldsymbol{x} + ih\boldsymbol{e}_n) \\ f_2(\boldsymbol{x} + ih\boldsymbol{e}_1) & f_2(\boldsymbol{x} + ih\boldsymbol{e}_2) & \cdots & f_2(\boldsymbol{x} + ih\boldsymbol{e}_q) & \cdots & f_2(\boldsymbol{x} + ih\boldsymbol{e}_n) \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ f_q(\boldsymbol{x} + ih\boldsymbol{e}_1) & f_q(\boldsymbol{x} + ih\boldsymbol{e}_2) & \cdots & f_q(\boldsymbol{x} + ih\boldsymbol{e}_q) & \cdots & f_q(\boldsymbol{x} + ih\boldsymbol{e}_n) \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ f_m(\boldsymbol{x} + ih\boldsymbol{e}_1) & f_m(\boldsymbol{x} + ih\boldsymbol{e}_2) & \cdots & f_m(\boldsymbol{x} + ih\boldsymbol{e}_q) & \cdots & f_m(\boldsymbol{x} + ih\boldsymbol{e}_n) \\ \end{pmatrix}\end{split}\]

where \(\boldsymbol{e}_p\) is the \(p^{\text{th}}\) column of an \(n^{\text{th}}\)-order identity matrix and \(f_q\) is the \(q^{\text{th}}\) equation of \(\boldsymbol{f(x)}\).

!!! CHECK CODE IN MPMATH !!!

Nth numerical (partial) derivative, using finite differences or numerical quadrature

ctx.diff(f, x, n=1, **options)

where ctx is dec, mpm, or gmp.

Numerically computes the derivative of \(f\), \(f'(x)\), or generally for an integer \(n \ge 0\), the \(n\)-th derivative \(f^{(n)}(x)\). A few basic examples are:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> diff(lambda x: x**2 + x, 1.0)
3.0
>>> diff(lambda x: x**2 + x, 1.0, 2)
2.0
>>> diff(lambda x: x**2 + x, 1.0, 3)
0.0
>>> nprint([diff(exp, 3, n) for n in range(5)])   # exp'(x) = exp(x)
[20.0855, 20.0855, 20.0855, 20.0855, 20.0855]

Even more generally, given a tuple of arguments \((x_1, \ldots, x_k)\) and order \((n_1, \ldots, n_k)\), the partial derivative \(f^{(n_1,\ldots,n_k)}(x_1,\ldots,x_k)\) is evaluated. For example:

>>> diff(lambda x,y: 3*x*y + 2*y - x, (0.25, 0.5), (0,1))
2.75
>>> diff(lambda x,y: 3*x*y + 2*y - x, (0.25, 0.5), (1,1))
3.0

Options

The following optional keyword arguments are recognized:

method

Supported methods are 'step' or 'quad': derivatives may be computed using either a finite difference with a small step size \(h\) (default), or numerical quadrature.

direction

Direction of finite difference: can be -1 for a left difference, 0 for a central difference (default), or +1 for a right difference; more generally can be any complex number.

addprec

Extra precision for \(h\) used to account for the function’s sensitivity to perturbations (default = 10).

relative

Choose \(h\) relative to the magnitude of \(x\), rather than an absolute value; useful for large or tiny \(x\) (default = False).

h

As an alternative to addprec and relative, manually select the step size \(h\).

singular

If True, evaluation exactly at the point \(x\) is avoided; this is useful for differentiating functions with removable singularities. Default = False.

radius

Radius of integration contour (with method = 'quad'). Default = 0.25. A larger radius typically is faster and more accurate, but it must be chosen so that \(f\) has no singularities within the radius from the evaluation point.

A finite difference requires \(n+1\) function evaluations and must be performed at \((n+1)\) times the target precision. Accordingly, \(f\) must support fast evaluation at high precision.

With integration, a larger number of function evaluations is required, but 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.

Further examples

The direction option is useful for computing left- or right-sided derivatives of nonsmooth functions:

>>> diff(abs, 0, direction=0)
0.0
>>> diff(abs, 0, direction=1)
1.0
>>> 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:

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

With integration, the result may have a small imaginary part even even if the result is purely real:

>>> diff(sqrt, 1, method='quad')
(0.5 - 4.59...e-26j)
>>> chop(_)
0.5

Adding precision to obtain an accurate value:

>>> diff(cos, 1e-30)
0.0
>>> diff(cos, 1e-30, h=0.0001)
-9.99999998328279e-31
>>> diff(cos, 1e-30, addprec=100)
-1.0e-30

Function object which evaluates the nth derivative of a given function

ctx.diffun(f, n=1, **options)

where ctx is dec, mpm, or gmp.

Given a function \(f\), returns a function \(g(x)\) that evaluates the nth derivative \(f^{(n)}(x)\):

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> cos2 = diffun(sin)
>>> sin2 = diffun(sin, 4)
>>> cos(1.3), cos2(1.3)
(0.267498828624587, 0.267498828624587)
>>> sin(1.3), sin2(1.3)
(0.963558185417193, 0.963558185417193)

The function \(f\) must support arbitrary precision evaluation. See diff() for additional details and supported keyword options.

Forward difference, based on a given sequence

ctx.difference(s, n)

where ctx is dec, mpm, or gmp.

Given a sequence \((s_k)\) containing at least \(n+1\) items, returns the \(n\)-th forward difference,

\[\Delta^n = \sum_{k=0}^{\infty} (-1)^{k+n} {n \choose k} s_k.\]

EXAMPLE !!!!!

Generating a sequence of derivatives

ctx.diffs(f, x, n=None, **options)

where ctx is dec, mpm, or gmp.

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.

Options are the same as for diff().

Examples

>>> from mpmath import *
>>> mp.dps = 15
>>> 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("%s %s" % (i, d))
...
0 0.54030230586814
1 -0.841470984807897
2 -0.54030230586814
3 0.841470984807897
4 0.54030230586814
5 -0.841470984807897

Composition of derivatives

ctx.diffs_prod(factors)

where ctx is dec, mpm, or gmp.

Given a list of \(N\) iterables or generators yielding \(f_k(x), f'_k(x), f''_k(x), \ldots\) for \(k = 1, \ldots, N\), generate \(g(x), g'(x), g''(x), \ldots\) where \(g(x) = f_1(x) f_2(x) \cdots f_N(x)\).

At high precision and for large orders, this is typically more efficient than numerical differentiation if the derivatives of each \(f_k(x)\) admit direct computation.

Note: This function does not increase the working precision internally, so guard digits may have to be added externally for full accuracy.

Examples

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> f = lambda x: exp(x)*cos(x)*sin(x)
>>> u = diffs(f, 1)
>>> v = mp.diffs_prod([diffs(exp,1), diffs(cos,1), diffs(sin,1)])
>>> next(u); next(v)
1.23586333600241
1.23586333600241
>>> next(u); next(v)
0.104658952245596
0.104658952245596
>>> next(u); next(v)
-5.96999877552086
-5.96999877552086
>>> next(u); next(v)
-12.4632923122697
-12.4632923122697

Composition of exponential of derivatives

ctx.diffs_exp(fdiffs)

where ctx is dec, mpm, or gmp.

Given an iterable or generator yielding \(f(x), f'(x), f''(x), \ldots\) generate \(g(x), g'(x), g''(x), \ldots\) where \(g(x) = \exp(f(x))\).

At high precision and for large orders, this is typically more efficient than numerical differentiation if the derivatives of \(f(x)\) admit direct computation.

Note: This function does not increase the working precision internally, so guard digits may have to be added externally for full accuracy.

Examples

The derivatives of the gamma function can be computed using logarithmic differentiation:

.. code-block:: pycon

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>>
>>> def diffs_loggamma(x):
...     yield loggamma(x)
...     i = 0
...     while 1:
...         yield psi(i,x)
...         i += 1
...
>>> u = diffs_exp(diffs_loggamma(3))
>>> v = diffs(gamma, 3)
>>> next(u); next(v)
2.0
2.0
>>> next(u); next(v)
1.84556867019693
1.84556867019693
>>> next(u); next(v)
2.49292999190269
2.49292999190269
>>> next(u); next(v)
3.44996501352367
3.44996501352367

Fractional derivatives / differintegration

ctx.differint(f, x, n=1, x0=0)

where ctx is dec, mpm, or gmp.

Calculates the Riemann-Liouville differintegral, or fractional derivative, defined by

\[\,_{x_0}{\mathbb{D}}^n_xf(x) = \frac{1}{\Gamma(m-n)} \frac{\mathrm{d}^m}{\mathrm{d} x^m} \int_{x_0}^{x}(x-t)^{m-n-1}f(t)\mathrm{d} t\]

where \(f\) is a given (presumably well-behaved) function, \(x\) is the evaluation point, \(n\) is the order, and \(x_0\) is the reference point of integration (\(m\) is an arbitrary parameter selected automatically).

With \(n = 1\), this is just the standard derivative \(f'(x)\); with \(n = 2\), the second derivative \(f''(x)\), etc. With \(n = -1\), it gives \(\int_{x_0}^x f(t) \mathrm{d} t\), with \(n = -2\) it gives \(\int_{x_0}^x \left( \int_{x_0}^t f(u) du \right) \mathrm{d} t\), etc.

As \(n\) is permitted to be any number, this operator generalizes iterated differentiation and iterated integration to a single operator with a continuous order parameter.

Examples

There is an exact formula for the fractional derivative of a monomial \(x^p\), which may be used as a reference. For example, the following gives a half-derivative (order 0.5):

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> x = mpf(3); p = 2; n = 0.5
>>> differint(lambda t: t**p, x, n)
7.81764019044672
>>> gamma(p+1)/gamma(p-n+1) * x**(p-n)
7.81764019044672

Another useful test function is the exponential function, whose integration / differentiation formula easy generalizes to arbitrary order. Here we first compute a third derivative, and then a triply nested integral. (The reference point \(x_0\) is set to \(-\infty\) to avoid nonzero endpoint terms.):

>>> differint(lambda x: exp(pi*x), -1.5, 3)
0.278538406900792
>>> exp(pi*-1.5) * pi**3
0.278538406900792
>>> differint(lambda x: exp(pi*x), 3.5, -3, -inf)
1922.50563031149
>>> exp(pi*3.5) / pi**3
1922.50563031149

However, for noninteger \(n\), the differentiation formula for the exponential function must be modified to give the same result as the Riemann-Liouville differintegral:

>>> x = mpf(3.5)
>>> c = pi
>>> n = 1+2*j
>>> differint(lambda x: exp(c*x), x, n)
(-123295.005390743 + 140955.117867654j)
>>> x**(-n) * exp(c)**x * (x*c)**n * gammainc(-n, 0, x*c) / gamma(-n)
(-123295.005390743 + 140955.117867654j)

Taylor series

ctx.taylor(f, x, n, **options)

where ctx is dec, mpm, or gmp.

Produces a degree-\(n\) Taylor polynomial around the point \(x\) of the given function \(f\). The coefficients are returned as a list.

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> nprint(chop(taylor(sin, 0, 5)))
[0.0, 1.0, 0.0, -0.166667, 0.0, 0.00833333]

The coefficients are computed using high-order numerical differentiation. The function must be possible to evaluate to arbitrary precision. See diff() for additional details and supported keyword options.

Note that to evaluate the Taylor polynomial as an approximation of \(f\), e.g. with polyval(), the coefficients must be reversed, and the point of the Taylor expansion must be subtracted from the argument:

>>> p = taylor(exp, 2.0, 10)
>>> polyval(p[::-1], 2.5 - 2.0)
12.1824939606092
>>> exp(2.5)
12.1824939607035

Solving an ODE using high-order Taylor series

ctx.odefun(F, x0, y0, tol=None, degree=None, method='taylor', verbose=False)

where ctx is dec, mpm, or gmp.

Returns a function \(y(x) = [y_0(x), y_1(x), \ldots, y_n(x)]\) that is a numerical solution of the \(n+1\)-dimensional first-order ordinary differential equation (ODE) system

\[ \begin{align}\begin{aligned}y_0'(x) = F_0(x, [y_0(x), y_1(x), \ldots, y_n(x)])\\y_1'(x) = F_1(x, [y_0(x), y_1(x), \ldots, y_n(x)])\\\vdots\\y_n'(x) = F_n(x, [y_0(x), y_1(x), \ldots, y_n(x)])\end{aligned}\end{align} \]

The derivatives are specified by the vector-valued function F that evaluates \([y_0', \ldots, y_n'] = F(x, [y_0, \ldots, y_n])\). The initial point \(x_0\) is specified by the scalar argument x0, and the initial value \(y(x_0) = [y_0(x_0), \ldots, y_n(x_0)]\) is specified by the vector argument y0.

For convenience, if the system is one-dimensional, you may optionally provide just a scalar value for y0. In this case, F should accept a scalar y argument and return a scalar. The solution function y will return scalar values instead of length-1 vectors.

Evaluation of the solution function \(y(x)\) is permitted for any \(x \ge x_0\).

A high-order ODE can be solved by transforming it into first-order vector form. This transformation is described in standard texts on ODEs. Examples will also be given below.

Options, speed and accuracy

By default, odefun() uses a high-order Taylor series method. For reasonably well-behaved problems, the solution will be fully accurate to within the working precision. Note that F must be possible to evaluate to very high precision for the generation of Taylor series to work.

To get a faster but less accurate solution, you can set a large value for tol (which defaults roughly to eps). If you just want to plot the solution or perform a basic simulation, tol = 0.01 is likely sufficient.

The degree argument controls the degree of the solver (with method=’taylor’, this is the degree of the Taylor series expansion). A higher degree means that a longer step can be taken before a new local solution must be generated from F, meaning that fewer steps are required to get from \(x_0\) to a given \(x_1\). On the other hand, a higher degree also means that each local solution becomes more expensive (i.e., more evaluations of F are required per step, and at higher precision).

The optimal setting therefore involves a tradeoff. Generally, decreasing the degree for Taylor series is likely to give faster solution at low precision, while increasing is likely to be better at higher precision.

The function object returned by odefun() caches the solutions at all step points and uses polynomial interpolation between step points. Therefore, once \(y(x_1)\) has been evaluated for some \(x_1\), \(y(x)\) can be evaluated very quickly for any \(x_0 \le x \le x_1\). and continuing the evaluation up to \(x_2 > x_1\) is also fast.

Examples of first-order ODEs

We will solve the standard test problem \(y'(x) = y(x), y(0) = 1\) which has explicit solution \(y(x) = \exp(x)\):

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> f = odefun(lambda x, y: y, 0, 1)
>>> for x in [0, 1, 2.5]:
...     print((f(x), exp(x)))
...
(1.0, 1.0)
(2.71828182845905, 2.71828182845905)
(12.1824939607035, 12.1824939607035)

The solution with high precision:

>>> mp.dps = 50
>>> f = odefun(lambda x, y: y, 0, 1)
>>> f(1)
2.7182818284590452353602874713526624977572470937
>>> exp(1)
2.7182818284590452353602874713526624977572470937

Using the more general vectorized form, the test problem can be input as (note that f returns a 1-element vector):

>>> mp.dps = 15
>>> f = odefun(lambda x, y: [y[0]], 0, [1])
>>> f(1)
[2.71828182845905]

odefun() can solve nonlinear ODEs, which are generally impossible (and at best difficult) to solve analytically. As an example of a nonlinear ODE, we will solve \(y'(x) = x \sin(y(x))\) for \(y(0) = \pi/2\). An exact solution happens to be known for this problem, and is given by \(y(x) = 2 \tan^{-1}\left(\exp\left(x^2/2\right)\right)\):

>>> f = odefun(lambda x, y: x*sin(y), 0, pi/2)
>>> for x in [2, 5, 10]:
...     print((f(x), 2*atan(exp(mpf(x)**2/2))))
...
(2.87255666284091, 2.87255666284091)
(3.14158520028345, 3.14158520028345)
(3.14159265358979, 3.14159265358979)

If \(F\) is independent of \(y\), an ODE can be solved using direct integration. We can therefore obtain a reference solution with quad():

>>> f = lambda x: (1+x**2)/(1+x**3)
>>> g = odefun(lambda x, y: f(x), pi, 0)
>>> g(2*pi)
0.72128263801696
>>> quad(f, [pi, 2*pi])
0.72128263801696

Examples of second-order ODEs

We will solve the harmonic oscillator equation \(y''(x) + y(x) = 0\). To do this, we introduce the helper functions \(y_0 = y, y_1 = y_0'\) whereby the original equation can be written as \(y_1' + y_0' = 0\). Put together, we get the first-order, two-dimensional vector ODE

\[\begin{split}\begin{cases} y_0' = y_1 \\ y_1' = -y_0 \end{cases}\end{split}\]

To get a well-defined IVP, we need two initial values. With \(y(0) = y_0(0) = 1\) and \(-y'(0) = y_1(0) = 0\), the problem will of course be solved by \(y(x) = y_0(x) = \cos(x)\) and \(-y'(x) = y_1(x) = \sin(x)\). We check this:

>>> f = odefun(lambda x, y: [-y[1], y[0]], 0, [1, 0])
>>> for x in [0, 1, 2.5, 10]:
...     nprint(f(x), 15)
...     nprint([cos(x), sin(x)], 15)
...     print("---")
...
[1.0, 0.0]
[1.0, 0.0]
---
[0.54030230586814, 0.841470984807897]
[0.54030230586814, 0.841470984807897]
---
[-0.801143615546934, 0.598472144103957]
[-0.801143615546934, 0.598472144103957]
---
[-0.839071529076452, -0.54402111088937]
[-0.839071529076452, -0.54402111088937]
---

Note that we get both the sine and the cosine solutions simultaneously.