Coverage for /usr/lib/python3/dist-packages/mpmath/calculus/quadrature.py: 11%
266 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import math
3from ..libmp.backend import xrange
5class QuadratureRule(object):
6 """
7 Quadrature rules are implemented using this class, in order to
8 simplify the code and provide a common infrastructure
9 for tasks such as error estimation and node caching.
11 You can implement a custom quadrature rule by subclassing
12 :class:`QuadratureRule` and implementing the appropriate
13 methods. The subclass can then be used by :func:`~mpmath.quad` by
14 passing it as the *method* argument.
16 :class:`QuadratureRule` instances are supposed to be singletons.
17 :class:`QuadratureRule` therefore implements instance caching
18 in :func:`~mpmath.__new__`.
19 """
21 def __init__(self, ctx):
22 self.ctx = ctx
23 self.standard_cache = {}
24 self.transformed_cache = {}
25 self.interval_count = {}
27 def clear(self):
28 """
29 Delete cached node data.
30 """
31 self.standard_cache = {}
32 self.transformed_cache = {}
33 self.interval_count = {}
35 def calc_nodes(self, degree, prec, verbose=False):
36 r"""
37 Compute nodes for the standard interval `[-1, 1]`. Subclasses
38 should probably implement only this method, and use
39 :func:`~mpmath.get_nodes` method to retrieve the nodes.
40 """
41 raise NotImplementedError
43 def get_nodes(self, a, b, degree, prec, verbose=False):
44 """
45 Return nodes for given interval, degree and precision. The
46 nodes are retrieved from a cache if already computed;
47 otherwise they are computed by calling :func:`~mpmath.calc_nodes`
48 and are then cached.
50 Subclasses should probably not implement this method,
51 but just implement :func:`~mpmath.calc_nodes` for the actual
52 node computation.
53 """
54 key = (a, b, degree, prec)
55 if key in self.transformed_cache:
56 return self.transformed_cache[key]
57 orig = self.ctx.prec
58 try:
59 self.ctx.prec = prec+20
60 # Get nodes on standard interval
61 if (degree, prec) in self.standard_cache:
62 nodes = self.standard_cache[degree, prec]
63 else:
64 nodes = self.calc_nodes(degree, prec, verbose)
65 self.standard_cache[degree, prec] = nodes
66 # Transform to general interval
67 nodes = self.transform_nodes(nodes, a, b, verbose)
68 if key in self.interval_count:
69 self.transformed_cache[key] = nodes
70 else:
71 self.interval_count[key] = True
72 finally:
73 self.ctx.prec = orig
74 return nodes
76 def transform_nodes(self, nodes, a, b, verbose=False):
77 r"""
78 Rescale standardized nodes (for `[-1, 1]`) to a general
79 interval `[a, b]`. For a finite interval, a simple linear
80 change of variables is used. Otherwise, the following
81 transformations are used:
83 .. math ::
85 \lbrack a, \infty \rbrack : t = \frac{1}{x} + (a-1)
87 \lbrack -\infty, b \rbrack : t = (b+1) - \frac{1}{x}
89 \lbrack -\infty, \infty \rbrack : t = \frac{x}{\sqrt{1-x^2}}
91 """
92 ctx = self.ctx
93 a = ctx.convert(a)
94 b = ctx.convert(b)
95 one = ctx.one
96 if (a, b) == (-one, one):
97 return nodes
98 half = ctx.mpf(0.5)
99 new_nodes = []
100 if ctx.isinf(a) or ctx.isinf(b):
101 if (a, b) == (ctx.ninf, ctx.inf):
102 p05 = -half
103 for x, w in nodes:
104 x2 = x*x
105 px1 = one-x2
106 spx1 = px1**p05
107 x = x*spx1
108 w *= spx1/px1
109 new_nodes.append((x, w))
110 elif a == ctx.ninf:
111 b1 = b+1
112 for x, w in nodes:
113 u = 2/(x+one)
114 x = b1-u
115 w *= half*u**2
116 new_nodes.append((x, w))
117 elif b == ctx.inf:
118 a1 = a-1
119 for x, w in nodes:
120 u = 2/(x+one)
121 x = a1+u
122 w *= half*u**2
123 new_nodes.append((x, w))
124 elif a == ctx.inf or b == ctx.ninf:
125 return [(x,-w) for (x,w) in self.transform_nodes(nodes, b, a, verbose)]
126 else:
127 raise NotImplementedError
128 else:
129 # Simple linear change of variables
130 C = (b-a)/2
131 D = (b+a)/2
132 for x, w in nodes:
133 new_nodes.append((D+C*x, C*w))
134 return new_nodes
136 def guess_degree(self, prec):
137 """
138 Given a desired precision `p` in bits, estimate the degree `m`
139 of the quadrature required to accomplish full accuracy for
140 typical integrals. By default, :func:`~mpmath.quad` will perform up
141 to `m` iterations. The value of `m` should be a slight
142 overestimate, so that "slightly bad" integrals can be dealt
143 with automatically using a few extra iterations. On the
144 other hand, it should not be too big, so :func:`~mpmath.quad` can
145 quit within a reasonable amount of time when it is given
146 an "unsolvable" integral.
148 The default formula used by :func:`~mpmath.guess_degree` is tuned
149 for both :class:`TanhSinh` and :class:`GaussLegendre`.
150 The output is roughly as follows:
152 +---------+---------+
153 | `p` | `m` |
154 +=========+=========+
155 | 50 | 6 |
156 +---------+---------+
157 | 100 | 7 |
158 +---------+---------+
159 | 500 | 10 |
160 +---------+---------+
161 | 3000 | 12 |
162 +---------+---------+
164 This formula is based purely on a limited amount of
165 experimentation and will sometimes be wrong.
166 """
167 # Expected degree
168 # XXX: use mag
169 g = int(4 + max(0, self.ctx.log(prec/30.0, 2)))
170 # Reasonable "worst case"
171 g += 2
172 return g
174 def estimate_error(self, results, prec, epsilon):
175 r"""
176 Given results from integrations `[I_1, I_2, \ldots, I_k]` done
177 with a quadrature of rule of degree `1, 2, \ldots, k`, estimate
178 the error of `I_k`.
180 For `k = 2`, we estimate `|I_{\infty}-I_2|` as `|I_2-I_1|`.
182 For `k > 2`, we extrapolate `|I_{\infty}-I_k| \approx |I_{k+1}-I_k|`
183 from `|I_k-I_{k-1}|` and `|I_k-I_{k-2}|` under the assumption
184 that each degree increment roughly doubles the accuracy of
185 the quadrature rule (this is true for both :class:`TanhSinh`
186 and :class:`GaussLegendre`). The extrapolation formula is given
187 by Borwein, Bailey & Girgensohn. Although not very conservative,
188 this method seems to be very robust in practice.
189 """
190 if len(results) == 2:
191 return abs(results[0]-results[1])
192 try:
193 if results[-1] == results[-2] == results[-3]:
194 return self.ctx.zero
195 D1 = self.ctx.log(abs(results[-1]-results[-2]), 10)
196 D2 = self.ctx.log(abs(results[-1]-results[-3]), 10)
197 except ValueError:
198 return epsilon
199 D3 = -prec
200 D4 = min(0, max(D1**2/D2, 2*D1, D3))
201 return self.ctx.mpf(10) ** int(D4)
203 def summation(self, f, points, prec, epsilon, max_degree, verbose=False):
204 """
205 Main integration function. Computes the 1D integral over
206 the interval specified by *points*. For each subinterval,
207 performs quadrature of degree from 1 up to *max_degree*
208 until :func:`~mpmath.estimate_error` signals convergence.
210 :func:`~mpmath.summation` transforms each subintegration to
211 the standard interval and then calls :func:`~mpmath.sum_next`.
212 """
213 ctx = self.ctx
214 I = total_err = ctx.zero
215 for i in xrange(len(points)-1):
216 a, b = points[i], points[i+1]
217 if a == b:
218 continue
219 # XXX: we could use a single variable transformation,
220 # but this is not good in practice. We get better accuracy
221 # by having 0 as an endpoint.
222 if (a, b) == (ctx.ninf, ctx.inf):
223 _f = f
224 f = lambda x: _f(-x) + _f(x)
225 a, b = (ctx.zero, ctx.inf)
226 results = []
227 err = ctx.zero
228 for degree in xrange(1, max_degree+1):
229 nodes = self.get_nodes(a, b, degree, prec, verbose)
230 if verbose:
231 print("Integrating from %s to %s (degree %s of %s)" % \
232 (ctx.nstr(a), ctx.nstr(b), degree, max_degree))
233 result = self.sum_next(f, nodes, degree, prec, results, verbose)
234 results.append(result)
235 if degree > 1:
236 err = self.estimate_error(results, prec, epsilon)
237 if verbose:
238 print("Estimated error:", ctx.nstr(err), " epsilon:", ctx.nstr(epsilon), " result: ", ctx.nstr(result))
239 if err <= epsilon:
240 break
241 I += results[-1]
242 total_err += err
243 if total_err > epsilon:
244 if verbose:
245 print("Failed to reach full accuracy. Estimated error:", ctx.nstr(total_err))
246 return I, total_err
248 def sum_next(self, f, nodes, degree, prec, previous, verbose=False):
249 r"""
250 Evaluates the step sum `\sum w_k f(x_k)` where the *nodes* list
251 contains the `(w_k, x_k)` pairs.
253 :func:`~mpmath.summation` will supply the list *results* of
254 values computed by :func:`~mpmath.sum_next` at previous degrees, in
255 case the quadrature rule is able to reuse them.
256 """
257 return self.ctx.fdot((w, f(x)) for (x,w) in nodes)
260class TanhSinh(QuadratureRule):
261 r"""
262 This class implements "tanh-sinh" or "doubly exponential"
263 quadrature. This quadrature rule is based on the Euler-Maclaurin
264 integral formula. By performing a change of variables involving
265 nested exponentials / hyperbolic functions (hence the name), the
266 derivatives at the endpoints vanish rapidly. Since the error term
267 in the Euler-Maclaurin formula depends on the derivatives at the
268 endpoints, a simple step sum becomes extremely accurate. In
269 practice, this means that doubling the number of evaluation
270 points roughly doubles the number of accurate digits.
272 Comparison to Gauss-Legendre:
273 * Initial computation of nodes is usually faster
274 * Handles endpoint singularities better
275 * Handles infinite integration intervals better
276 * Is slower for smooth integrands once nodes have been computed
278 The implementation of the tanh-sinh algorithm is based on the
279 description given in Borwein, Bailey & Girgensohn, "Experimentation
280 in Mathematics - Computational Paths to Discovery", A K Peters,
281 2003, pages 312-313. In the present implementation, a few
282 improvements have been made:
284 * A more efficient scheme is used to compute nodes (exploiting
285 recurrence for the exponential function)
286 * The nodes are computed successively instead of all at once
288 Various documents describing the algorithm are available online, e.g.:
290 * http://crd.lbl.gov/~dhbailey/dhbpapers/dhb-tanh-sinh.pdf
291 * http://users.cs.dal.ca/~jborwein/tanh-sinh.pdf
292 """
294 def sum_next(self, f, nodes, degree, prec, previous, verbose=False):
295 """
296 Step sum for tanh-sinh quadrature of degree `m`. We exploit the
297 fact that half of the abscissas at degree `m` are precisely the
298 abscissas from degree `m-1`. Thus reusing the result from
299 the previous level allows a 2x speedup.
300 """
301 h = self.ctx.mpf(2)**(-degree)
302 # Abscissas overlap, so reusing saves half of the time
303 if previous:
304 S = previous[-1]/(h*2)
305 else:
306 S = self.ctx.zero
307 S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
308 return h*S
310 def calc_nodes(self, degree, prec, verbose=False):
311 r"""
312 The abscissas and weights for tanh-sinh quadrature of degree
313 `m` are given by
315 .. math::
317 x_k = \tanh(\pi/2 \sinh(t_k))
319 w_k = \pi/2 \cosh(t_k) / \cosh(\pi/2 \sinh(t_k))^2
321 where `t_k = t_0 + hk` for a step length `h \sim 2^{-m}`. The
322 list of nodes is actually infinite, but the weights die off so
323 rapidly that only a few are needed.
324 """
325 ctx = self.ctx
326 nodes = []
328 extra = 20
329 ctx.prec += extra
330 tol = ctx.ldexp(1, -prec-10)
331 pi4 = ctx.pi/4
333 # For simplicity, we work in steps h = 1/2^n, with the first point
334 # offset so that we can reuse the sum from the previous degree
336 # We define degree 1 to include the "degree 0" steps, including
337 # the point x = 0. (It doesn't work well otherwise; not sure why.)
338 t0 = ctx.ldexp(1, -degree)
339 if degree == 1:
340 #nodes.append((mpf(0), pi4))
341 #nodes.append((-mpf(0), pi4))
342 nodes.append((ctx.zero, ctx.pi/2))
343 h = t0
344 else:
345 h = t0*2
347 # Since h is fixed, we can compute the next exponential
348 # by simply multiplying by exp(h)
349 expt0 = ctx.exp(t0)
350 a = pi4 * expt0
351 b = pi4 / expt0
352 udelta = ctx.exp(h)
353 urdelta = 1/udelta
355 for k in xrange(0, 20*2**degree+1):
356 # Reference implementation:
357 # t = t0 + k*h
358 # x = tanh(pi/2 * sinh(t))
359 # w = pi/2 * cosh(t) / cosh(pi/2 * sinh(t))**2
361 # Fast implementation. Note that c = exp(pi/2 * sinh(t))
362 c = ctx.exp(a-b)
363 d = 1/c
364 co = (c+d)/2
365 si = (c-d)/2
366 x = si / co
367 w = (a+b) / co**2
368 diff = abs(x-1)
369 if diff <= tol:
370 break
372 nodes.append((x, w))
373 nodes.append((-x, w))
375 a *= udelta
376 b *= urdelta
378 if verbose and k % 300 == 150:
379 # Note: the number displayed is rather arbitrary. Should
380 # figure out how to print something that looks more like a
381 # percentage
382 print("Calculating nodes:", ctx.nstr(-ctx.log(diff, 10) / prec))
384 ctx.prec -= extra
385 return nodes
388class GaussLegendre(QuadratureRule):
389 r"""
390 This class implements Gauss-Legendre quadrature, which is
391 exceptionally efficient for polynomials and polynomial-like (i.e.
392 very smooth) integrands.
394 The abscissas and weights are given by roots and values of
395 Legendre polynomials, which are the orthogonal polynomials
396 on `[-1, 1]` with respect to the unit weight
397 (see :func:`~mpmath.legendre`).
399 In this implementation, we take the "degree" `m` of the quadrature
400 to denote a Gauss-Legendre rule of degree `3 \cdot 2^m` (following
401 Borwein, Bailey & Girgensohn). This way we get quadratic, rather
402 than linear, convergence as the degree is incremented.
404 Comparison to tanh-sinh quadrature:
405 * Is faster for smooth integrands once nodes have been computed
406 * Initial computation of nodes is usually slower
407 * Handles endpoint singularities worse
408 * Handles infinite integration intervals worse
410 """
412 def calc_nodes(self, degree, prec, verbose=False):
413 r"""
414 Calculates the abscissas and weights for Gauss-Legendre
415 quadrature of degree of given degree (actually `3 \cdot 2^m`).
416 """
417 ctx = self.ctx
418 # It is important that the epsilon is set lower than the
419 # "real" epsilon
420 epsilon = ctx.ldexp(1, -prec-8)
421 # Fairly high precision might be required for accurate
422 # evaluation of the roots
423 orig = ctx.prec
424 ctx.prec = int(prec*1.5)
425 if degree == 1:
426 x = ctx.sqrt(ctx.mpf(3)/5)
427 w = ctx.mpf(5)/9
428 nodes = [(-x,w),(ctx.zero,ctx.mpf(8)/9),(x,w)]
429 ctx.prec = orig
430 return nodes
431 nodes = []
432 n = 3*2**(degree-1)
433 upto = n//2 + 1
434 for j in xrange(1, upto):
435 # Asymptotic formula for the roots
436 r = ctx.mpf(math.cos(math.pi*(j-0.25)/(n+0.5)))
437 # Newton iteration
438 while 1:
439 t1, t2 = 1, 0
440 # Evaluates the Legendre polynomial using its defining
441 # recurrence relation
442 for j1 in xrange(1,n+1):
443 t3, t2, t1 = t2, t1, ((2*j1-1)*r*t1 - (j1-1)*t2)/j1
444 t4 = n*(r*t1-t2)/(r**2-1)
445 a = t1/t4
446 r = r - a
447 if abs(a) < epsilon:
448 break
449 x = r
450 w = 2/((1-r**2)*t4**2)
451 if verbose and j % 30 == 15:
452 print("Computing nodes (%i of %i)" % (j, upto))
453 nodes.append((x, w))
454 nodes.append((-x, w))
455 ctx.prec = orig
456 return nodes
458class QuadratureMethods(object):
460 def __init__(ctx, *args, **kwargs):
461 ctx._gauss_legendre = GaussLegendre(ctx)
462 ctx._tanh_sinh = TanhSinh(ctx)
464 def quad(ctx, f, *points, **kwargs):
465 r"""
466 Computes a single, double or triple integral over a given
467 1D interval, 2D rectangle, or 3D cuboid. A basic example::
469 >>> from mpmath import *
470 >>> mp.dps = 15; mp.pretty = True
471 >>> quad(sin, [0, pi])
472 2.0
474 A basic 2D integral::
476 >>> f = lambda x, y: cos(x+y/2)
477 >>> quad(f, [-pi/2, pi/2], [0, pi])
478 4.0
480 **Interval format**
482 The integration range for each dimension may be specified
483 using a list or tuple. Arguments are interpreted as follows:
485 ``quad(f, [x1, x2])`` -- calculates
486 `\int_{x_1}^{x_2} f(x) \, dx`
488 ``quad(f, [x1, x2], [y1, y2])`` -- calculates
489 `\int_{x_1}^{x_2} \int_{y_1}^{y_2} f(x,y) \, dy \, dx`
491 ``quad(f, [x1, x2], [y1, y2], [z1, z2])`` -- calculates
492 `\int_{x_1}^{x_2} \int_{y_1}^{y_2} \int_{z_1}^{z_2} f(x,y,z)
493 \, dz \, dy \, dx`
495 Endpoints may be finite or infinite. An interval descriptor
496 may also contain more than two points. In this
497 case, the integration is split into subintervals, between
498 each pair of consecutive points. This is useful for
499 dealing with mid-interval discontinuities, or integrating
500 over large intervals where the function is irregular or
501 oscillates.
503 **Options**
505 :func:`~mpmath.quad` recognizes the following keyword arguments:
507 *method*
508 Chooses integration algorithm (described below).
509 *error*
510 If set to true, :func:`~mpmath.quad` returns `(v, e)` where `v` is the
511 integral and `e` is the estimated error.
512 *maxdegree*
513 Maximum degree of the quadrature rule to try before
514 quitting.
515 *verbose*
516 Print details about progress.
518 **Algorithms**
520 Mpmath presently implements two integration algorithms: tanh-sinh
521 quadrature and Gauss-Legendre quadrature. These can be selected
522 using *method='tanh-sinh'* or *method='gauss-legendre'* or by
523 passing the classes *method=TanhSinh*, *method=GaussLegendre*.
524 The functions :func:`~mpmath.quadts` and :func:`~mpmath.quadgl` are also available
525 as shortcuts.
527 Both algorithms have the property that doubling the number of
528 evaluation points roughly doubles the accuracy, so both are ideal
529 for high precision quadrature (hundreds or thousands of digits).
531 At high precision, computing the nodes and weights for the
532 integration can be expensive (more expensive than computing the
533 function values). To make repeated integrations fast, nodes
534 are automatically cached.
536 The advantages of the tanh-sinh algorithm are that it tends to
537 handle endpoint singularities well, and that the nodes are cheap
538 to compute on the first run. For these reasons, it is used by
539 :func:`~mpmath.quad` as the default algorithm.
541 Gauss-Legendre quadrature often requires fewer function
542 evaluations, and is therefore often faster for repeated use, but
543 the algorithm does not handle endpoint singularities as well and
544 the nodes are more expensive to compute. Gauss-Legendre quadrature
545 can be a better choice if the integrand is smooth and repeated
546 integrations are required (e.g. for multiple integrals).
548 See the documentation for :class:`TanhSinh` and
549 :class:`GaussLegendre` for additional details.
551 **Examples of 1D integrals**
553 Intervals may be infinite or half-infinite. The following two
554 examples evaluate the limits of the inverse tangent function
555 (`\int 1/(1+x^2) = \tan^{-1} x`), and the Gaussian integral
556 `\int_{\infty}^{\infty} \exp(-x^2)\,dx = \sqrt{\pi}`::
558 >>> mp.dps = 15
559 >>> quad(lambda x: 2/(x**2+1), [0, inf])
560 3.14159265358979
561 >>> quad(lambda x: exp(-x**2), [-inf, inf])**2
562 3.14159265358979
564 Integrals can typically be resolved to high precision.
565 The following computes 50 digits of `\pi` by integrating the
566 area of the half-circle defined by `x^2 + y^2 \le 1`,
567 `-1 \le x \le 1`, `y \ge 0`::
569 >>> mp.dps = 50
570 >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1])
571 3.1415926535897932384626433832795028841971693993751
573 One can just as well compute 1000 digits (output truncated)::
575 >>> mp.dps = 1000
576 >>> 2*quad(lambda x: sqrt(1-x**2), [-1, 1]) #doctest:+ELLIPSIS
577 3.141592653589793238462643383279502884...216420199
579 Complex integrals are supported. The following computes
580 a residue at `z = 0` by integrating counterclockwise along the
581 diamond-shaped path from `1` to `+i` to `-1` to `-i` to `1`::
583 >>> mp.dps = 15
584 >>> chop(quad(lambda z: 1/z, [1,j,-1,-j,1]))
585 (0.0 + 6.28318530717959j)
587 **Examples of 2D and 3D integrals**
589 Here are several nice examples of analytically solvable
590 2D integrals (taken from MathWorld [1]) that can be evaluated
591 to high precision fairly rapidly by :func:`~mpmath.quad`::
593 >>> mp.dps = 30
594 >>> f = lambda x, y: (x-1)/((1-x*y)*log(x*y))
595 >>> quad(f, [0, 1], [0, 1])
596 0.577215664901532860606512090082
597 >>> +euler
598 0.577215664901532860606512090082
600 >>> f = lambda x, y: 1/sqrt(1+x**2+y**2)
601 >>> quad(f, [-1, 1], [-1, 1])
602 3.17343648530607134219175646705
603 >>> 4*log(2+sqrt(3))-2*pi/3
604 3.17343648530607134219175646705
606 >>> f = lambda x, y: 1/(1-x**2 * y**2)
607 >>> quad(f, [0, 1], [0, 1])
608 1.23370055013616982735431137498
609 >>> pi**2 / 8
610 1.23370055013616982735431137498
612 >>> quad(lambda x, y: 1/(1-x*y), [0, 1], [0, 1])
613 1.64493406684822643647241516665
614 >>> pi**2 / 6
615 1.64493406684822643647241516665
617 Multiple integrals may be done over infinite ranges::
619 >>> mp.dps = 15
620 >>> print(quad(lambda x,y: exp(-x-y), [0, inf], [1, inf]))
621 0.367879441171442
622 >>> print(1/e)
623 0.367879441171442
625 For nonrectangular areas, one can call :func:`~mpmath.quad` recursively.
626 For example, we can replicate the earlier example of calculating
627 `\pi` by integrating over the unit-circle, and actually use double
628 quadrature to actually measure the area circle::
630 >>> f = lambda x: quad(lambda y: 1, [-sqrt(1-x**2), sqrt(1-x**2)])
631 >>> quad(f, [-1, 1])
632 3.14159265358979
634 Here is a simple triple integral::
636 >>> mp.dps = 15
637 >>> f = lambda x,y,z: x*y/(1+z)
638 >>> quad(f, [0,1], [0,1], [1,2], method='gauss-legendre')
639 0.101366277027041
640 >>> (log(3)-log(2))/4
641 0.101366277027041
643 **Singularities**
645 Both tanh-sinh and Gauss-Legendre quadrature are designed to
646 integrate smooth (infinitely differentiable) functions. Neither
647 algorithm copes well with mid-interval singularities (such as
648 mid-interval discontinuities in `f(x)` or `f'(x)`).
649 The best solution is to split the integral into parts::
651 >>> mp.dps = 15
652 >>> quad(lambda x: abs(sin(x)), [0, 2*pi]) # Bad
653 3.99900894176779
654 >>> quad(lambda x: abs(sin(x)), [0, pi, 2*pi]) # Good
655 4.0
657 The tanh-sinh rule often works well for integrands having a
658 singularity at one or both endpoints::
660 >>> mp.dps = 15
661 >>> quad(log, [0, 1], method='tanh-sinh') # Good
662 -1.0
663 >>> quad(log, [0, 1], method='gauss-legendre') # Bad
664 -0.999932197413801
666 However, the result may still be inaccurate for some functions::
668 >>> quad(lambda x: 1/sqrt(x), [0, 1], method='tanh-sinh')
669 1.99999999946942
671 This problem is not due to the quadrature rule per se, but to
672 numerical amplification of errors in the nodes. The problem can be
673 circumvented by temporarily increasing the precision::
675 >>> mp.dps = 30
676 >>> a = quad(lambda x: 1/sqrt(x), [0, 1], method='tanh-sinh')
677 >>> mp.dps = 15
678 >>> +a
679 2.0
681 **Highly variable functions**
683 For functions that are smooth (in the sense of being infinitely
684 differentiable) but contain sharp mid-interval peaks or many
685 "bumps", :func:`~mpmath.quad` may fail to provide full accuracy. For
686 example, with default settings, :func:`~mpmath.quad` is able to integrate
687 `\sin(x)` accurately over an interval of length 100 but not over
688 length 1000::
690 >>> quad(sin, [0, 100]); 1-cos(100) # Good
691 0.137681127712316
692 0.137681127712316
693 >>> quad(sin, [0, 1000]); 1-cos(1000) # Bad
694 -37.8587612408485
695 0.437620923709297
697 One solution is to break the integration into 10 intervals of
698 length 100::
700 >>> quad(sin, linspace(0, 1000, 10)) # Good
701 0.437620923709297
703 Another is to increase the degree of the quadrature::
705 >>> quad(sin, [0, 1000], maxdegree=10) # Also good
706 0.437620923709297
708 Whether splitting the interval or increasing the degree is
709 more efficient differs from case to case. Another example is the
710 function `1/(1+x^2)`, which has a sharp peak centered around
711 `x = 0`::
713 >>> f = lambda x: 1/(1+x**2)
714 >>> quad(f, [-100, 100]) # Bad
715 3.64804647105268
716 >>> quad(f, [-100, 100], maxdegree=10) # Good
717 3.12159332021646
718 >>> quad(f, [-100, 0, 100]) # Also good
719 3.12159332021646
721 **References**
723 1. http://mathworld.wolfram.com/DoubleIntegral.html
725 """
726 rule = kwargs.get('method', 'tanh-sinh')
727 if type(rule) is str:
728 if rule == 'tanh-sinh':
729 rule = ctx._tanh_sinh
730 elif rule == 'gauss-legendre':
731 rule = ctx._gauss_legendre
732 else:
733 raise ValueError("unknown quadrature rule: %s" % rule)
734 else:
735 rule = rule(ctx)
736 verbose = kwargs.get('verbose')
737 dim = len(points)
738 orig = prec = ctx.prec
739 epsilon = ctx.eps/8
740 m = kwargs.get('maxdegree') or rule.guess_degree(prec)
741 points = [ctx._as_points(p) for p in points]
742 try:
743 ctx.prec += 20
744 if dim == 1:
745 v, err = rule.summation(f, points[0], prec, epsilon, m, verbose)
746 elif dim == 2:
747 v, err = rule.summation(lambda x: \
748 rule.summation(lambda y: f(x,y), \
749 points[1], prec, epsilon, m)[0],
750 points[0], prec, epsilon, m, verbose)
751 elif dim == 3:
752 v, err = rule.summation(lambda x: \
753 rule.summation(lambda y: \
754 rule.summation(lambda z: f(x,y,z), \
755 points[2], prec, epsilon, m)[0],
756 points[1], prec, epsilon, m)[0],
757 points[0], prec, epsilon, m, verbose)
758 else:
759 raise NotImplementedError("quadrature must have dim 1, 2 or 3")
760 finally:
761 ctx.prec = orig
762 if kwargs.get("error"):
763 return +v, err
764 return +v
766 def quadts(ctx, *args, **kwargs):
767 """
768 Performs tanh-sinh quadrature. The call
770 quadts(func, *points, ...)
772 is simply a shortcut for:
774 quad(func, *points, ..., method=TanhSinh)
776 For example, a single integral and a double integral:
778 quadts(lambda x: exp(cos(x)), [0, 1])
779 quadts(lambda x, y: exp(cos(x+y)), [0, 1], [0, 1])
781 See the documentation for quad for information about how points
782 arguments and keyword arguments are parsed.
784 See documentation for TanhSinh for algorithmic information about
785 tanh-sinh quadrature.
786 """
787 kwargs['method'] = 'tanh-sinh'
788 return ctx.quad(*args, **kwargs)
790 def quadgl(ctx, *args, **kwargs):
791 """
792 Performs Gauss-Legendre quadrature. The call
794 quadgl(func, *points, ...)
796 is simply a shortcut for:
798 quad(func, *points, ..., method=GaussLegendre)
800 For example, a single integral and a double integral:
802 quadgl(lambda x: exp(cos(x)), [0, 1])
803 quadgl(lambda x, y: exp(cos(x+y)), [0, 1], [0, 1])
805 See the documentation for quad for information about how points
806 arguments and keyword arguments are parsed.
808 See documentation for TanhSinh for algorithmic information about
809 tanh-sinh quadrature.
810 """
811 kwargs['method'] = 'gauss-legendre'
812 return ctx.quad(*args, **kwargs)
814 def quadosc(ctx, f, interval, omega=None, period=None, zeros=None):
815 r"""
816 Calculates
818 .. math ::
820 I = \int_a^b f(x) dx
822 where at least one of `a` and `b` is infinite and where
823 `f(x) = g(x) \cos(\omega x + \phi)` for some slowly
824 decreasing function `g(x)`. With proper input, :func:`~mpmath.quadosc`
825 can also handle oscillatory integrals where the oscillation
826 rate is different from a pure sine or cosine wave.
828 In the standard case when `|a| < \infty, b = \infty`,
829 :func:`~mpmath.quadosc` works by evaluating the infinite series
831 .. math ::
833 I = \int_a^{x_1} f(x) dx +
834 \sum_{k=1}^{\infty} \int_{x_k}^{x_{k+1}} f(x) dx
836 where `x_k` are consecutive zeros (alternatively
837 some other periodic reference point) of `f(x)`.
838 Accordingly, :func:`~mpmath.quadosc` requires information about the
839 zeros of `f(x)`. For a periodic function, you can specify
840 the zeros by either providing the angular frequency `\omega`
841 (*omega*) or the *period* `2 \pi/\omega`. In general, you can
842 specify the `n`-th zero by providing the *zeros* arguments.
843 Below is an example of each::
845 >>> from mpmath import *
846 >>> mp.dps = 15; mp.pretty = True
847 >>> f = lambda x: sin(3*x)/(x**2+1)
848 >>> quadosc(f, [0,inf], omega=3)
849 0.37833007080198
850 >>> quadosc(f, [0,inf], period=2*pi/3)
851 0.37833007080198
852 >>> quadosc(f, [0,inf], zeros=lambda n: pi*n/3)
853 0.37833007080198
854 >>> (ei(3)*exp(-3)-exp(3)*ei(-3))/2 # Computed by Mathematica
855 0.37833007080198
857 Note that *zeros* was specified to multiply `n` by the
858 *half-period*, not the full period. In theory, it does not matter
859 whether each partial integral is done over a half period or a full
860 period. However, if done over half-periods, the infinite series
861 passed to :func:`~mpmath.nsum` becomes an *alternating series* and this
862 typically makes the extrapolation much more efficient.
864 Here is an example of an integration over the entire real line,
865 and a half-infinite integration starting at `-\infty`::
867 >>> quadosc(lambda x: cos(x)/(1+x**2), [-inf, inf], omega=1)
868 1.15572734979092
869 >>> pi/e
870 1.15572734979092
871 >>> quadosc(lambda x: cos(x)/x**2, [-inf, -1], period=2*pi)
872 -0.0844109505595739
873 >>> cos(1)+si(1)-pi/2
874 -0.0844109505595738
876 Of course, the integrand may contain a complex exponential just as
877 well as a real sine or cosine::
879 >>> quadosc(lambda x: exp(3*j*x)/(1+x**2), [-inf,inf], omega=3)
880 (0.156410688228254 + 0.0j)
881 >>> pi/e**3
882 0.156410688228254
883 >>> quadosc(lambda x: exp(3*j*x)/(2+x+x**2), [-inf,inf], omega=3)
884 (0.00317486988463794 - 0.0447701735209082j)
885 >>> 2*pi/sqrt(7)/exp(3*(j+sqrt(7))/2)
886 (0.00317486988463794 - 0.0447701735209082j)
888 **Non-periodic functions**
890 If `f(x) = g(x) h(x)` for some function `h(x)` that is not
891 strictly periodic, *omega* or *period* might not work, and it might
892 be necessary to use *zeros*.
894 A notable exception can be made for Bessel functions which, though not
895 periodic, are "asymptotically periodic" in a sufficiently strong sense
896 that the sum extrapolation will work out::
898 >>> quadosc(j0, [0, inf], period=2*pi)
899 1.0
900 >>> quadosc(j1, [0, inf], period=2*pi)
901 1.0
903 More properly, one should provide the exact Bessel function zeros::
905 >>> j0zero = lambda n: findroot(j0, pi*(n-0.25))
906 >>> quadosc(j0, [0, inf], zeros=j0zero)
907 1.0
909 For an example where *zeros* becomes necessary, consider the
910 complete Fresnel integrals
912 .. math ::
914 \int_0^{\infty} \cos x^2\,dx = \int_0^{\infty} \sin x^2\,dx
915 = \sqrt{\frac{\pi}{8}}.
917 Although the integrands do not decrease in magnitude as
918 `x \to \infty`, the integrals are convergent since the oscillation
919 rate increases (causing consecutive periods to asymptotically
920 cancel out). These integrals are virtually impossible to calculate
921 to any kind of accuracy using standard quadrature rules. However,
922 if one provides the correct asymptotic distribution of zeros
923 (`x_n \sim \sqrt{n}`), :func:`~mpmath.quadosc` works::
925 >>> mp.dps = 30
926 >>> f = lambda x: cos(x**2)
927 >>> quadosc(f, [0,inf], zeros=lambda n:sqrt(pi*n))
928 0.626657068657750125603941321203
929 >>> f = lambda x: sin(x**2)
930 >>> quadosc(f, [0,inf], zeros=lambda n:sqrt(pi*n))
931 0.626657068657750125603941321203
932 >>> sqrt(pi/8)
933 0.626657068657750125603941321203
935 (Interestingly, these integrals can still be evaluated if one
936 places some other constant than `\pi` in the square root sign.)
938 In general, if `f(x) \sim g(x) \cos(h(x))`, the zeros follow
939 the inverse-function distribution `h^{-1}(x)`::
941 >>> mp.dps = 15
942 >>> f = lambda x: sin(exp(x))
943 >>> quadosc(f, [1,inf], zeros=lambda n: log(n))
944 -0.25024394235267
945 >>> pi/2-si(e)
946 -0.250243942352671
948 **Non-alternating functions**
950 If the integrand oscillates around a positive value, without
951 alternating signs, the extrapolation might fail. A simple trick
952 that sometimes works is to multiply or divide the frequency by 2::
954 >>> f = lambda x: 1/x**2+sin(x)/x**4
955 >>> quadosc(f, [1,inf], omega=1) # Bad
956 1.28642190869861
957 >>> quadosc(f, [1,inf], omega=0.5) # Perfect
958 1.28652953559617
959 >>> 1+(cos(1)+ci(1)+sin(1))/6
960 1.28652953559617
962 **Fast decay**
964 :func:`~mpmath.quadosc` is primarily useful for slowly decaying
965 integrands. If the integrand decreases exponentially or faster,
966 :func:`~mpmath.quad` will likely handle it without trouble (and generally be
967 much faster than :func:`~mpmath.quadosc`)::
969 >>> quadosc(lambda x: cos(x)/exp(x), [0, inf], omega=1)
970 0.5
971 >>> quad(lambda x: cos(x)/exp(x), [0, inf])
972 0.5
974 """
975 a, b = ctx._as_points(interval)
976 a = ctx.convert(a)
977 b = ctx.convert(b)
978 if [omega, period, zeros].count(None) != 2:
979 raise ValueError( \
980 "must specify exactly one of omega, period, zeros")
981 if a == ctx.ninf and b == ctx.inf:
982 s1 = ctx.quadosc(f, [a, 0], omega=omega, zeros=zeros, period=period)
983 s2 = ctx.quadosc(f, [0, b], omega=omega, zeros=zeros, period=period)
984 return s1 + s2
985 if a == ctx.ninf:
986 if zeros:
987 return ctx.quadosc(lambda x:f(-x), [-b,-a], lambda n: zeros(-n))
988 else:
989 return ctx.quadosc(lambda x:f(-x), [-b,-a], omega=omega, period=period)
990 if b != ctx.inf:
991 raise ValueError("quadosc requires an infinite integration interval")
992 if not zeros:
993 if omega:
994 period = 2*ctx.pi/omega
995 zeros = lambda n: n*period/2
996 #for n in range(1,10):
997 # p = zeros(n)
998 # if p > a:
999 # break
1000 #if n >= 9:
1001 # raise ValueError("zeros do not appear to be correctly indexed")
1002 n = 1
1003 s = ctx.quadgl(f, [a, zeros(n)])
1004 def term(k):
1005 return ctx.quadgl(f, [zeros(k), zeros(k+1)])
1006 s += ctx.nsum(term, [n, ctx.inf])
1007 return s
1009if __name__ == '__main__':
1010 import doctest
1011 doctest.testmod()