Coverage for /usr/lib/python3/dist-packages/mpmath/calculus/differentiation.py: 12%
226 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from ..libmp.backend import xrange
2from .calculus import defun
4try:
5 iteritems = dict.iteritems
6except AttributeError:
7 iteritems = dict.items
9#----------------------------------------------------------------------------#
10# Differentiation #
11#----------------------------------------------------------------------------#
13@defun
14def difference(ctx, s, n):
15 r"""
16 Given a sequence `(s_k)` containing at least `n+1` items, returns the
17 `n`-th forward difference,
19 .. math ::
21 \Delta^n = \sum_{k=0}^{\infty} (-1)^{k+n} {n \choose k} s_k.
22 """
23 n = int(n)
24 d = ctx.zero
25 b = (-1) ** (n & 1)
26 for k in xrange(n+1):
27 d += b * s[k]
28 b = (b * (k-n)) // (k+1)
29 return d
31def hsteps(ctx, f, x, n, prec, **options):
32 singular = options.get('singular')
33 addprec = options.get('addprec', 10)
34 direction = options.get('direction', 0)
35 workprec = (prec+2*addprec) * (n+1)
36 orig = ctx.prec
37 try:
38 ctx.prec = workprec
39 h = options.get('h')
40 if h is None:
41 if options.get('relative'):
42 hextramag = int(ctx.mag(x))
43 else:
44 hextramag = 0
45 h = ctx.ldexp(1, -prec-addprec-hextramag)
46 else:
47 h = ctx.convert(h)
48 # Directed: steps x, x+h, ... x+n*h
49 direction = options.get('direction', 0)
50 if direction:
51 h *= ctx.sign(direction)
52 steps = xrange(n+1)
53 norm = h
54 # Central: steps x-n*h, x-(n-2)*h ..., x, ..., x+(n-2)*h, x+n*h
55 else:
56 steps = xrange(-n, n+1, 2)
57 norm = (2*h)
58 # Perturb
59 if singular:
60 x += 0.5*h
61 values = [f(x+k*h) for k in steps]
62 return values, norm, workprec
63 finally:
64 ctx.prec = orig
67@defun
68def diff(ctx, f, x, n=1, **options):
69 r"""
70 Numerically computes the derivative of `f`, `f'(x)`, or generally for
71 an integer `n \ge 0`, the `n`-th derivative `f^{(n)}(x)`.
72 A few basic examples are::
74 >>> from mpmath import *
75 >>> mp.dps = 15; mp.pretty = True
76 >>> diff(lambda x: x**2 + x, 1.0)
77 3.0
78 >>> diff(lambda x: x**2 + x, 1.0, 2)
79 2.0
80 >>> diff(lambda x: x**2 + x, 1.0, 3)
81 0.0
82 >>> nprint([diff(exp, 3, n) for n in range(5)]) # exp'(x) = exp(x)
83 [20.0855, 20.0855, 20.0855, 20.0855, 20.0855]
85 Even more generally, given a tuple of arguments `(x_1, \ldots, x_k)`
86 and order `(n_1, \ldots, n_k)`, the partial derivative
87 `f^{(n_1,\ldots,n_k)}(x_1,\ldots,x_k)` is evaluated. For example::
89 >>> diff(lambda x,y: 3*x*y + 2*y - x, (0.25, 0.5), (0,1))
90 2.75
91 >>> diff(lambda x,y: 3*x*y + 2*y - x, (0.25, 0.5), (1,1))
92 3.0
94 **Options**
96 The following optional keyword arguments are recognized:
98 ``method``
99 Supported methods are ``'step'`` or ``'quad'``: derivatives may be
100 computed using either a finite difference with a small step
101 size `h` (default), or numerical quadrature.
102 ``direction``
103 Direction of finite difference: can be -1 for a left
104 difference, 0 for a central difference (default), or +1
105 for a right difference; more generally can be any complex number.
106 ``addprec``
107 Extra precision for `h` used to account for the function's
108 sensitivity to perturbations (default = 10).
109 ``relative``
110 Choose `h` relative to the magnitude of `x`, rather than an
111 absolute value; useful for large or tiny `x` (default = False).
112 ``h``
113 As an alternative to ``addprec`` and ``relative``, manually
114 select the step size `h`.
115 ``singular``
116 If True, evaluation exactly at the point `x` is avoided; this is
117 useful for differentiating functions with removable singularities.
118 Default = False.
119 ``radius``
120 Radius of integration contour (with ``method = 'quad'``).
121 Default = 0.25. A larger radius typically is faster and more
122 accurate, but it must be chosen so that `f` has no
123 singularities within the radius from the evaluation point.
125 A finite difference requires `n+1` function evaluations and must be
126 performed at `(n+1)` times the target precision. Accordingly, `f` must
127 support fast evaluation at high precision.
129 With integration, a larger number of function evaluations is
130 required, but not much extra precision is required. For high order
131 derivatives, this method may thus be faster if f is very expensive to
132 evaluate at high precision.
134 **Further examples**
136 The direction option is useful for computing left- or right-sided
137 derivatives of nonsmooth functions::
139 >>> diff(abs, 0, direction=0)
140 0.0
141 >>> diff(abs, 0, direction=1)
142 1.0
143 >>> diff(abs, 0, direction=-1)
144 -1.0
146 More generally, if the direction is nonzero, a right difference
147 is computed where the step size is multiplied by sign(direction).
148 For example, with direction=+j, the derivative from the positive
149 imaginary direction will be computed::
151 >>> diff(abs, 0, direction=j)
152 (0.0 - 1.0j)
154 With integration, the result may have a small imaginary part
155 even even if the result is purely real::
157 >>> diff(sqrt, 1, method='quad') # doctest:+ELLIPSIS
158 (0.5 - 4.59...e-26j)
159 >>> chop(_)
160 0.5
162 Adding precision to obtain an accurate value::
164 >>> diff(cos, 1e-30)
165 0.0
166 >>> diff(cos, 1e-30, h=0.0001)
167 -9.99999998328279e-31
168 >>> diff(cos, 1e-30, addprec=100)
169 -1.0e-30
171 """
172 partial = False
173 try:
174 orders = list(n)
175 x = list(x)
176 partial = True
177 except TypeError:
178 pass
179 if partial:
180 x = [ctx.convert(_) for _ in x]
181 return _partial_diff(ctx, f, x, orders, options)
182 method = options.get('method', 'step')
183 if n == 0 and method != 'quad' and not options.get('singular'):
184 return f(ctx.convert(x))
185 prec = ctx.prec
186 try:
187 if method == 'step':
188 values, norm, workprec = hsteps(ctx, f, x, n, prec, **options)
189 ctx.prec = workprec
190 v = ctx.difference(values, n) / norm**n
191 elif method == 'quad':
192 ctx.prec += 10
193 radius = ctx.convert(options.get('radius', 0.25))
194 def g(t):
195 rei = radius*ctx.expj(t)
196 z = x + rei
197 return f(z) / rei**n
198 d = ctx.quadts(g, [0, 2*ctx.pi])
199 v = d * ctx.factorial(n) / (2*ctx.pi)
200 else:
201 raise ValueError("unknown method: %r" % method)
202 finally:
203 ctx.prec = prec
204 return +v
206def _partial_diff(ctx, f, xs, orders, options):
207 if not orders:
208 return f()
209 if not sum(orders):
210 return f(*xs)
211 i = 0
212 for i in range(len(orders)):
213 if orders[i]:
214 break
215 order = orders[i]
216 def fdiff_inner(*f_args):
217 def inner(t):
218 return f(*(f_args[:i] + (t,) + f_args[i+1:]))
219 return ctx.diff(inner, f_args[i], order, **options)
220 orders[i] = 0
221 return _partial_diff(ctx, fdiff_inner, xs, orders, options)
223@defun
224def diffs(ctx, f, x, n=None, **options):
225 r"""
226 Returns a generator that yields the sequence of derivatives
228 .. math ::
230 f(x), f'(x), f''(x), \ldots, f^{(k)}(x), \ldots
232 With ``method='step'``, :func:`~mpmath.diffs` uses only `O(k)`
233 function evaluations to generate the first `k` derivatives,
234 rather than the roughly `O(k^2)` evaluations
235 required if one calls :func:`~mpmath.diff` `k` separate times.
237 With `n < \infty`, the generator stops as soon as the
238 `n`-th derivative has been generated. If the exact number of
239 needed derivatives is known in advance, this is further
240 slightly more efficient.
242 Options are the same as for :func:`~mpmath.diff`.
244 **Examples**
246 >>> from mpmath import *
247 >>> mp.dps = 15
248 >>> nprint(list(diffs(cos, 1, 5)))
249 [0.540302, -0.841471, -0.540302, 0.841471, 0.540302, -0.841471]
250 >>> for i, d in zip(range(6), diffs(cos, 1)):
251 ... print("%s %s" % (i, d))
252 ...
253 0 0.54030230586814
254 1 -0.841470984807897
255 2 -0.54030230586814
256 3 0.841470984807897
257 4 0.54030230586814
258 5 -0.841470984807897
260 """
261 if n is None:
262 n = ctx.inf
263 else:
264 n = int(n)
265 if options.get('method', 'step') != 'step':
266 k = 0
267 while k < n + 1:
268 yield ctx.diff(f, x, k, **options)
269 k += 1
270 return
271 singular = options.get('singular')
272 if singular:
273 yield ctx.diff(f, x, 0, singular=True)
274 else:
275 yield f(ctx.convert(x))
276 if n < 1:
277 return
278 if n == ctx.inf:
279 A, B = 1, 2
280 else:
281 A, B = 1, n+1
282 while 1:
283 callprec = ctx.prec
284 y, norm, workprec = hsteps(ctx, f, x, B, callprec, **options)
285 for k in xrange(A, B):
286 try:
287 ctx.prec = workprec
288 d = ctx.difference(y, k) / norm**k
289 finally:
290 ctx.prec = callprec
291 yield +d
292 if k >= n:
293 return
294 A, B = B, int(A*1.4+1)
295 B = min(B, n)
297def iterable_to_function(gen):
298 gen = iter(gen)
299 data = []
300 def f(k):
301 for i in xrange(len(data), k+1):
302 data.append(next(gen))
303 return data[k]
304 return f
306@defun
307def diffs_prod(ctx, factors):
308 r"""
309 Given a list of `N` iterables or generators yielding
310 `f_k(x), f'_k(x), f''_k(x), \ldots` for `k = 1, \ldots, N`,
311 generate `g(x), g'(x), g''(x), \ldots` where
312 `g(x) = f_1(x) f_2(x) \cdots f_N(x)`.
314 At high precision and for large orders, this is typically more efficient
315 than numerical differentiation if the derivatives of each `f_k(x)`
316 admit direct computation.
318 Note: This function does not increase the working precision internally,
319 so guard digits may have to be added externally for full accuracy.
321 **Examples**
323 >>> from mpmath import *
324 >>> mp.dps = 15; mp.pretty = True
325 >>> f = lambda x: exp(x)*cos(x)*sin(x)
326 >>> u = diffs(f, 1)
327 >>> v = mp.diffs_prod([diffs(exp,1), diffs(cos,1), diffs(sin,1)])
328 >>> next(u); next(v)
329 1.23586333600241
330 1.23586333600241
331 >>> next(u); next(v)
332 0.104658952245596
333 0.104658952245596
334 >>> next(u); next(v)
335 -5.96999877552086
336 -5.96999877552086
337 >>> next(u); next(v)
338 -12.4632923122697
339 -12.4632923122697
341 """
342 N = len(factors)
343 if N == 1:
344 for c in factors[0]:
345 yield c
346 else:
347 u = iterable_to_function(ctx.diffs_prod(factors[:N//2]))
348 v = iterable_to_function(ctx.diffs_prod(factors[N//2:]))
349 n = 0
350 while 1:
351 #yield sum(binomial(n,k)*u(n-k)*v(k) for k in xrange(n+1))
352 s = u(n) * v(0)
353 a = 1
354 for k in xrange(1,n+1):
355 a = a * (n-k+1) // k
356 s += a * u(n-k) * v(k)
357 yield s
358 n += 1
360def dpoly(n, _cache={}):
361 """
362 nth differentiation polynomial for exp (Faa di Bruno's formula).
364 TODO: most exponents are zero, so maybe a sparse representation
365 would be better.
366 """
367 if n in _cache:
368 return _cache[n]
369 if not _cache:
370 _cache[0] = {(0,):1}
371 R = dpoly(n-1)
372 R = dict((c+(0,),v) for (c,v) in iteritems(R))
373 Ra = {}
374 for powers, count in iteritems(R):
375 powers1 = (powers[0]+1,) + powers[1:]
376 if powers1 in Ra:
377 Ra[powers1] += count
378 else:
379 Ra[powers1] = count
380 for powers, count in iteritems(R):
381 if not sum(powers):
382 continue
383 for k,p in enumerate(powers):
384 if p:
385 powers2 = powers[:k] + (p-1,powers[k+1]+1) + powers[k+2:]
386 if powers2 in Ra:
387 Ra[powers2] += p*count
388 else:
389 Ra[powers2] = p*count
390 _cache[n] = Ra
391 return _cache[n]
393@defun
394def diffs_exp(ctx, fdiffs):
395 r"""
396 Given an iterable or generator yielding `f(x), f'(x), f''(x), \ldots`
397 generate `g(x), g'(x), g''(x), \ldots` where `g(x) = \exp(f(x))`.
399 At high precision and for large orders, this is typically more efficient
400 than numerical differentiation if the derivatives of `f(x)`
401 admit direct computation.
403 Note: This function does not increase the working precision internally,
404 so guard digits may have to be added externally for full accuracy.
406 **Examples**
408 The derivatives of the gamma function can be computed using
409 logarithmic differentiation::
411 >>> from mpmath import *
412 >>> mp.dps = 15; mp.pretty = True
413 >>>
414 >>> def diffs_loggamma(x):
415 ... yield loggamma(x)
416 ... i = 0
417 ... while 1:
418 ... yield psi(i,x)
419 ... i += 1
420 ...
421 >>> u = diffs_exp(diffs_loggamma(3))
422 >>> v = diffs(gamma, 3)
423 >>> next(u); next(v)
424 2.0
425 2.0
426 >>> next(u); next(v)
427 1.84556867019693
428 1.84556867019693
429 >>> next(u); next(v)
430 2.49292999190269
431 2.49292999190269
432 >>> next(u); next(v)
433 3.44996501352367
434 3.44996501352367
436 """
437 fn = iterable_to_function(fdiffs)
438 f0 = ctx.exp(fn(0))
439 yield f0
440 i = 1
441 while 1:
442 s = ctx.mpf(0)
443 for powers, c in iteritems(dpoly(i)):
444 s += c*ctx.fprod(fn(k+1)**p for (k,p) in enumerate(powers) if p)
445 yield s * f0
446 i += 1
448@defun
449def differint(ctx, f, x, n=1, x0=0):
450 r"""
451 Calculates the Riemann-Liouville differintegral, or fractional
452 derivative, defined by
454 .. math ::
456 \,_{x_0}{\mathbb{D}}^n_xf(x) = \frac{1}{\Gamma(m-n)} \frac{d^m}{dx^m}
457 \int_{x_0}^{x}(x-t)^{m-n-1}f(t)dt
459 where `f` is a given (presumably well-behaved) function,
460 `x` is the evaluation point, `n` is the order, and `x_0` is
461 the reference point of integration (`m` is an arbitrary
462 parameter selected automatically).
464 With `n = 1`, this is just the standard derivative `f'(x)`; with `n = 2`,
465 the second derivative `f''(x)`, etc. With `n = -1`, it gives
466 `\int_{x_0}^x f(t) dt`, with `n = -2`
467 it gives `\int_{x_0}^x \left( \int_{x_0}^t f(u) du \right) dt`, etc.
469 As `n` is permitted to be any number, this operator generalizes
470 iterated differentiation and iterated integration to a single
471 operator with a continuous order parameter.
473 **Examples**
475 There is an exact formula for the fractional derivative of a
476 monomial `x^p`, which may be used as a reference. For example,
477 the following gives a half-derivative (order 0.5)::
479 >>> from mpmath import *
480 >>> mp.dps = 15; mp.pretty = True
481 >>> x = mpf(3); p = 2; n = 0.5
482 >>> differint(lambda t: t**p, x, n)
483 7.81764019044672
484 >>> gamma(p+1)/gamma(p-n+1) * x**(p-n)
485 7.81764019044672
487 Another useful test function is the exponential function, whose
488 integration / differentiation formula easy generalizes
489 to arbitrary order. Here we first compute a third derivative,
490 and then a triply nested integral. (The reference point `x_0`
491 is set to `-\infty` to avoid nonzero endpoint terms.)::
493 >>> differint(lambda x: exp(pi*x), -1.5, 3)
494 0.278538406900792
495 >>> exp(pi*-1.5) * pi**3
496 0.278538406900792
497 >>> differint(lambda x: exp(pi*x), 3.5, -3, -inf)
498 1922.50563031149
499 >>> exp(pi*3.5) / pi**3
500 1922.50563031149
502 However, for noninteger `n`, the differentiation formula for the
503 exponential function must be modified to give the same result as the
504 Riemann-Liouville differintegral::
506 >>> x = mpf(3.5)
507 >>> c = pi
508 >>> n = 1+2*j
509 >>> differint(lambda x: exp(c*x), x, n)
510 (-123295.005390743 + 140955.117867654j)
511 >>> x**(-n) * exp(c)**x * (x*c)**n * gammainc(-n, 0, x*c) / gamma(-n)
512 (-123295.005390743 + 140955.117867654j)
515 """
516 m = max(int(ctx.ceil(ctx.re(n)))+1, 1)
517 r = m-n-1
518 g = lambda x: ctx.quad(lambda t: (x-t)**r * f(t), [x0, x])
519 return ctx.diff(g, x, m) / ctx.gamma(m-n)
521@defun
522def diffun(ctx, f, n=1, **options):
523 r"""
524 Given a function `f`, returns a function `g(x)` that evaluates the nth
525 derivative `f^{(n)}(x)`::
527 >>> from mpmath import *
528 >>> mp.dps = 15; mp.pretty = True
529 >>> cos2 = diffun(sin)
530 >>> sin2 = diffun(sin, 4)
531 >>> cos(1.3), cos2(1.3)
532 (0.267498828624587, 0.267498828624587)
533 >>> sin(1.3), sin2(1.3)
534 (0.963558185417193, 0.963558185417193)
536 The function `f` must support arbitrary precision evaluation.
537 See :func:`~mpmath.diff` for additional details and supported
538 keyword options.
539 """
540 if n == 0:
541 return f
542 def g(x):
543 return ctx.diff(f, x, n, **options)
544 return g
546@defun
547def taylor(ctx, f, x, n, **options):
548 r"""
549 Produces a degree-`n` Taylor polynomial around the point `x` of the
550 given function `f`. The coefficients are returned as a list.
552 >>> from mpmath import *
553 >>> mp.dps = 15; mp.pretty = True
554 >>> nprint(chop(taylor(sin, 0, 5)))
555 [0.0, 1.0, 0.0, -0.166667, 0.0, 0.00833333]
557 The coefficients are computed using high-order numerical
558 differentiation. The function must be possible to evaluate
559 to arbitrary precision. See :func:`~mpmath.diff` for additional details
560 and supported keyword options.
562 Note that to evaluate the Taylor polynomial as an approximation
563 of `f`, e.g. with :func:`~mpmath.polyval`, the coefficients must be reversed,
564 and the point of the Taylor expansion must be subtracted from
565 the argument:
567 >>> p = taylor(exp, 2.0, 10)
568 >>> polyval(p[::-1], 2.5 - 2.0)
569 12.1824939606092
570 >>> exp(2.5)
571 12.1824939607035
573 """
574 gen = enumerate(ctx.diffs(f, x, n, **options))
575 if options.get("chop", True):
576 return [ctx.chop(d)/ctx.factorial(i) for i, d in gen]
577 else:
578 return [d/ctx.factorial(i) for i, d in gen]
580@defun
581def pade(ctx, a, L, M):
582 r"""
583 Computes a Pade approximation of degree `(L, M)` to a function.
584 Given at least `L+M+1` Taylor coefficients `a` approximating
585 a function `A(x)`, :func:`~mpmath.pade` returns coefficients of
586 polynomials `P, Q` satisfying
588 .. math ::
590 P = \sum_{k=0}^L p_k x^k
592 Q = \sum_{k=0}^M q_k x^k
594 Q_0 = 1
596 A(x) Q(x) = P(x) + O(x^{L+M+1})
598 `P(x)/Q(x)` can provide a good approximation to an analytic function
599 beyond the radius of convergence of its Taylor series (example
600 from G.A. Baker 'Essentials of Pade Approximants' Academic Press,
601 Ch.1A)::
603 >>> from mpmath import *
604 >>> mp.dps = 15; mp.pretty = True
605 >>> one = mpf(1)
606 >>> def f(x):
607 ... return sqrt((one + 2*x)/(one + x))
608 ...
609 >>> a = taylor(f, 0, 6)
610 >>> p, q = pade(a, 3, 3)
611 >>> x = 10
612 >>> polyval(p[::-1], x)/polyval(q[::-1], x)
613 1.38169105566806
614 >>> f(x)
615 1.38169855941551
617 """
618 # To determine L+1 coefficients of P and M coefficients of Q
619 # L+M+1 coefficients of A must be provided
620 if len(a) < L+M+1:
621 raise ValueError("L+M+1 Coefficients should be provided")
623 if M == 0:
624 if L == 0:
625 return [ctx.one], [ctx.one]
626 else:
627 return a[:L+1], [ctx.one]
629 # Solve first
630 # a[L]*q[1] + ... + a[L-M+1]*q[M] = -a[L+1]
631 # ...
632 # a[L+M-1]*q[1] + ... + a[L]*q[M] = -a[L+M]
633 A = ctx.matrix(M)
634 for j in range(M):
635 for i in range(min(M, L+j+1)):
636 A[j, i] = a[L+j-i]
637 v = -ctx.matrix(a[(L+1):(L+M+1)])
638 x = ctx.lu_solve(A, v)
639 q = [ctx.one] + list(x)
640 # compute p
641 p = [0]*(L+1)
642 for i in range(L+1):
643 s = a[i]
644 for j in range(1, min(M,i) + 1):
645 s += q[j]*a[i-j]
646 p[i] = s
647 return p, q