Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/integrate/_quadrature.py: 11%
397 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1from __future__ import annotations
2from typing import TYPE_CHECKING, Callable, Dict, Tuple, Any, cast
3import functools
4import numpy as np
5import math
6import types
7import warnings
8from collections import namedtuple
10from scipy.special import roots_legendre
11from scipy.special import gammaln, logsumexp
12from scipy._lib._util import _rng_spawn
15__all__ = ['fixed_quad', 'quadrature', 'romberg', 'romb',
16 'trapezoid', 'trapz', 'simps', 'simpson',
17 'cumulative_trapezoid', 'cumtrapz', 'newton_cotes',
18 'AccuracyWarning']
21def trapezoid(y, x=None, dx=1.0, axis=-1):
22 r"""
23 Integrate along the given axis using the composite trapezoidal rule.
25 If `x` is provided, the integration happens in sequence along its
26 elements - they are not sorted.
28 Integrate `y` (`x`) along each 1d slice on the given axis, compute
29 :math:`\int y(x) dx`.
30 When `x` is specified, this integrates along the parametric curve,
31 computing :math:`\int_t y(t) dt =
32 \int_t y(t) \left.\frac{dx}{dt}\right|_{x=x(t)} dt`.
34 Parameters
35 ----------
36 y : array_like
37 Input array to integrate.
38 x : array_like, optional
39 The sample points corresponding to the `y` values. If `x` is None,
40 the sample points are assumed to be evenly spaced `dx` apart. The
41 default is None.
42 dx : scalar, optional
43 The spacing between sample points when `x` is None. The default is 1.
44 axis : int, optional
45 The axis along which to integrate.
47 Returns
48 -------
49 trapezoid : float or ndarray
50 Definite integral of `y` = n-dimensional array as approximated along
51 a single axis by the trapezoidal rule. If `y` is a 1-dimensional array,
52 then the result is a float. If `n` is greater than 1, then the result
53 is an `n`-1 dimensional array.
55 See Also
56 --------
57 cumulative_trapezoid, simpson, romb
59 Notes
60 -----
61 Image [2]_ illustrates trapezoidal rule -- y-axis locations of points
62 will be taken from `y` array, by default x-axis distances between
63 points will be 1.0, alternatively they can be provided with `x` array
64 or with `dx` scalar. Return value will be equal to combined area under
65 the red lines.
67 References
68 ----------
69 .. [1] Wikipedia page: https://en.wikipedia.org/wiki/Trapezoidal_rule
71 .. [2] Illustration image:
72 https://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
74 Examples
75 --------
76 Use the trapezoidal rule on evenly spaced points:
78 >>> import numpy as np
79 >>> from scipy import integrate
80 >>> integrate.trapezoid([1, 2, 3])
81 4.0
83 The spacing between sample points can be selected by either the
84 ``x`` or ``dx`` arguments:
86 >>> integrate.trapezoid([1, 2, 3], x=[4, 6, 8])
87 8.0
88 >>> integrate.trapezoid([1, 2, 3], dx=2)
89 8.0
91 Using a decreasing ``x`` corresponds to integrating in reverse:
93 >>> integrate.trapezoid([1, 2, 3], x=[8, 6, 4])
94 -8.0
96 More generally ``x`` is used to integrate along a parametric curve. We can
97 estimate the integral :math:`\int_0^1 x^2 = 1/3` using:
99 >>> x = np.linspace(0, 1, num=50)
100 >>> y = x**2
101 >>> integrate.trapezoid(y, x)
102 0.33340274885464394
104 Or estimate the area of a circle, noting we repeat the sample which closes
105 the curve:
107 >>> theta = np.linspace(0, 2 * np.pi, num=1000, endpoint=True)
108 >>> integrate.trapezoid(np.cos(theta), x=np.sin(theta))
109 3.141571941375841
111 ``trapezoid`` can be applied along a specified axis to do multiple
112 computations in one call:
114 >>> a = np.arange(6).reshape(2, 3)
115 >>> a
116 array([[0, 1, 2],
117 [3, 4, 5]])
118 >>> integrate.trapezoid(a, axis=0)
119 array([1.5, 2.5, 3.5])
120 >>> integrate.trapezoid(a, axis=1)
121 array([2., 8.])
122 """
123 # Future-proofing, in case NumPy moves from trapz to trapezoid for the same
124 # reasons as SciPy
125 if hasattr(np, 'trapezoid'):
126 return np.trapezoid(y, x=x, dx=dx, axis=axis)
127 else:
128 return np.trapz(y, x=x, dx=dx, axis=axis)
131# Note: alias kept for backwards compatibility. Rename was done
132# because trapz is a slur in colloquial English (see gh-12924).
133def trapz(y, x=None, dx=1.0, axis=-1):
134 """An alias of `trapezoid`.
136 `trapz` is kept for backwards compatibility. For new code, prefer
137 `trapezoid` instead.
138 """
139 return trapezoid(y, x=x, dx=dx, axis=axis)
142class AccuracyWarning(Warning):
143 pass
146if TYPE_CHECKING:
147 # workaround for mypy function attributes see:
148 # https://github.com/python/mypy/issues/2087#issuecomment-462726600
149 from typing import Protocol
151 class CacheAttributes(Protocol):
152 cache: Dict[int, Tuple[Any, Any]]
153else:
154 CacheAttributes = Callable
157def cache_decorator(func: Callable) -> CacheAttributes:
158 return cast(CacheAttributes, func)
161@cache_decorator
162def _cached_roots_legendre(n):
163 """
164 Cache roots_legendre results to speed up calls of the fixed_quad
165 function.
166 """
167 if n in _cached_roots_legendre.cache:
168 return _cached_roots_legendre.cache[n]
170 _cached_roots_legendre.cache[n] = roots_legendre(n)
171 return _cached_roots_legendre.cache[n]
174_cached_roots_legendre.cache = dict()
177def fixed_quad(func, a, b, args=(), n=5):
178 """
179 Compute a definite integral using fixed-order Gaussian quadrature.
181 Integrate `func` from `a` to `b` using Gaussian quadrature of
182 order `n`.
184 Parameters
185 ----------
186 func : callable
187 A Python function or method to integrate (must accept vector inputs).
188 If integrating a vector-valued function, the returned array must have
189 shape ``(..., len(x))``.
190 a : float
191 Lower limit of integration.
192 b : float
193 Upper limit of integration.
194 args : tuple, optional
195 Extra arguments to pass to function, if any.
196 n : int, optional
197 Order of quadrature integration. Default is 5.
199 Returns
200 -------
201 val : float
202 Gaussian quadrature approximation to the integral
203 none : None
204 Statically returned value of None
206 See Also
207 --------
208 quad : adaptive quadrature using QUADPACK
209 dblquad : double integrals
210 tplquad : triple integrals
211 romberg : adaptive Romberg quadrature
212 quadrature : adaptive Gaussian quadrature
213 romb : integrators for sampled data
214 simpson : integrators for sampled data
215 cumulative_trapezoid : cumulative integration for sampled data
216 ode : ODE integrator
217 odeint : ODE integrator
219 Examples
220 --------
221 >>> from scipy import integrate
222 >>> import numpy as np
223 >>> f = lambda x: x**8
224 >>> integrate.fixed_quad(f, 0.0, 1.0, n=4)
225 (0.1110884353741496, None)
226 >>> integrate.fixed_quad(f, 0.0, 1.0, n=5)
227 (0.11111111111111102, None)
228 >>> print(1/9.0) # analytical result
229 0.1111111111111111
231 >>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=4)
232 (0.9999999771971152, None)
233 >>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=5)
234 (1.000000000039565, None)
235 >>> np.sin(np.pi/2)-np.sin(0) # analytical result
236 1.0
238 """
239 x, w = _cached_roots_legendre(n)
240 x = np.real(x)
241 if np.isinf(a) or np.isinf(b):
242 raise ValueError("Gaussian quadrature is only available for "
243 "finite limits.")
244 y = (b-a)*(x+1)/2.0 + a
245 return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None
248def vectorize1(func, args=(), vec_func=False):
249 """Vectorize the call to a function.
251 This is an internal utility function used by `romberg` and
252 `quadrature` to create a vectorized version of a function.
254 If `vec_func` is True, the function `func` is assumed to take vector
255 arguments.
257 Parameters
258 ----------
259 func : callable
260 User defined function.
261 args : tuple, optional
262 Extra arguments for the function.
263 vec_func : bool, optional
264 True if the function func takes vector arguments.
266 Returns
267 -------
268 vfunc : callable
269 A function that will take a vector argument and return the
270 result.
272 """
273 if vec_func:
274 def vfunc(x):
275 return func(x, *args)
276 else:
277 def vfunc(x):
278 if np.isscalar(x):
279 return func(x, *args)
280 x = np.asarray(x)
281 # call with first point to get output type
282 y0 = func(x[0], *args)
283 n = len(x)
284 dtype = getattr(y0, 'dtype', type(y0))
285 output = np.empty((n,), dtype=dtype)
286 output[0] = y0
287 for i in range(1, n):
288 output[i] = func(x[i], *args)
289 return output
290 return vfunc
293def quadrature(func, a, b, args=(), tol=1.49e-8, rtol=1.49e-8, maxiter=50,
294 vec_func=True, miniter=1):
295 """
296 Compute a definite integral using fixed-tolerance Gaussian quadrature.
298 Integrate `func` from `a` to `b` using Gaussian quadrature
299 with absolute tolerance `tol`.
301 Parameters
302 ----------
303 func : function
304 A Python function or method to integrate.
305 a : float
306 Lower limit of integration.
307 b : float
308 Upper limit of integration.
309 args : tuple, optional
310 Extra arguments to pass to function.
311 tol, rtol : float, optional
312 Iteration stops when error between last two iterates is less than
313 `tol` OR the relative change is less than `rtol`.
314 maxiter : int, optional
315 Maximum order of Gaussian quadrature.
316 vec_func : bool, optional
317 True or False if func handles arrays as arguments (is
318 a "vector" function). Default is True.
319 miniter : int, optional
320 Minimum order of Gaussian quadrature.
322 Returns
323 -------
324 val : float
325 Gaussian quadrature approximation (within tolerance) to integral.
326 err : float
327 Difference between last two estimates of the integral.
329 See Also
330 --------
331 romberg : adaptive Romberg quadrature
332 fixed_quad : fixed-order Gaussian quadrature
333 quad : adaptive quadrature using QUADPACK
334 dblquad : double integrals
335 tplquad : triple integrals
336 romb : integrator for sampled data
337 simpson : integrator for sampled data
338 cumulative_trapezoid : cumulative integration for sampled data
339 ode : ODE integrator
340 odeint : ODE integrator
342 Examples
343 --------
344 >>> from scipy import integrate
345 >>> import numpy as np
346 >>> f = lambda x: x**8
347 >>> integrate.quadrature(f, 0.0, 1.0)
348 (0.11111111111111106, 4.163336342344337e-17)
349 >>> print(1/9.0) # analytical result
350 0.1111111111111111
352 >>> integrate.quadrature(np.cos, 0.0, np.pi/2)
353 (0.9999999999999536, 3.9611425250996035e-11)
354 >>> np.sin(np.pi/2)-np.sin(0) # analytical result
355 1.0
357 """
358 if not isinstance(args, tuple):
359 args = (args,)
360 vfunc = vectorize1(func, args, vec_func=vec_func)
361 val = np.inf
362 err = np.inf
363 maxiter = max(miniter+1, maxiter)
364 for n in range(miniter, maxiter+1):
365 newval = fixed_quad(vfunc, a, b, (), n)[0]
366 err = abs(newval-val)
367 val = newval
369 if err < tol or err < rtol*abs(val):
370 break
371 else:
372 warnings.warn(
373 "maxiter (%d) exceeded. Latest difference = %e" % (maxiter, err),
374 AccuracyWarning)
375 return val, err
378def tupleset(t, i, value):
379 l = list(t)
380 l[i] = value
381 return tuple(l)
384# Note: alias kept for backwards compatibility. Rename was done
385# because cumtrapz is a slur in colloquial English (see gh-12924).
386def cumtrapz(y, x=None, dx=1.0, axis=-1, initial=None):
387 """An alias of `cumulative_trapezoid`.
389 `cumtrapz` is kept for backwards compatibility. For new code, prefer
390 `cumulative_trapezoid` instead.
391 """
392 return cumulative_trapezoid(y, x=x, dx=dx, axis=axis, initial=initial)
395def cumulative_trapezoid(y, x=None, dx=1.0, axis=-1, initial=None):
396 """
397 Cumulatively integrate y(x) using the composite trapezoidal rule.
399 Parameters
400 ----------
401 y : array_like
402 Values to integrate.
403 x : array_like, optional
404 The coordinate to integrate along. If None (default), use spacing `dx`
405 between consecutive elements in `y`.
406 dx : float, optional
407 Spacing between elements of `y`. Only used if `x` is None.
408 axis : int, optional
409 Specifies the axis to cumulate. Default is -1 (last axis).
410 initial : scalar, optional
411 If given, insert this value at the beginning of the returned result.
412 Typically this value should be 0. Default is None, which means no
413 value at ``x[0]`` is returned and `res` has one element less than `y`
414 along the axis of integration.
416 Returns
417 -------
418 res : ndarray
419 The result of cumulative integration of `y` along `axis`.
420 If `initial` is None, the shape is such that the axis of integration
421 has one less value than `y`. If `initial` is given, the shape is equal
422 to that of `y`.
424 See Also
425 --------
426 numpy.cumsum, numpy.cumprod
427 quad : adaptive quadrature using QUADPACK
428 romberg : adaptive Romberg quadrature
429 quadrature : adaptive Gaussian quadrature
430 fixed_quad : fixed-order Gaussian quadrature
431 dblquad : double integrals
432 tplquad : triple integrals
433 romb : integrators for sampled data
434 ode : ODE integrators
435 odeint : ODE integrators
437 Examples
438 --------
439 >>> from scipy import integrate
440 >>> import numpy as np
441 >>> import matplotlib.pyplot as plt
443 >>> x = np.linspace(-2, 2, num=20)
444 >>> y = x
445 >>> y_int = integrate.cumulative_trapezoid(y, x, initial=0)
446 >>> plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-')
447 >>> plt.show()
449 """
450 y = np.asarray(y)
451 if x is None:
452 d = dx
453 else:
454 x = np.asarray(x)
455 if x.ndim == 1:
456 d = np.diff(x)
457 # reshape to correct shape
458 shape = [1] * y.ndim
459 shape[axis] = -1
460 d = d.reshape(shape)
461 elif len(x.shape) != len(y.shape):
462 raise ValueError("If given, shape of x must be 1-D or the "
463 "same as y.")
464 else:
465 d = np.diff(x, axis=axis)
467 if d.shape[axis] != y.shape[axis] - 1:
468 raise ValueError("If given, length of x along axis must be the "
469 "same as y.")
471 nd = len(y.shape)
472 slice1 = tupleset((slice(None),)*nd, axis, slice(1, None))
473 slice2 = tupleset((slice(None),)*nd, axis, slice(None, -1))
474 res = np.cumsum(d * (y[slice1] + y[slice2]) / 2.0, axis=axis)
476 if initial is not None:
477 if not np.isscalar(initial):
478 raise ValueError("`initial` parameter should be a scalar.")
480 shape = list(res.shape)
481 shape[axis] = 1
482 res = np.concatenate([np.full(shape, initial, dtype=res.dtype), res],
483 axis=axis)
485 return res
488def _basic_simpson(y, start, stop, x, dx, axis):
489 nd = len(y.shape)
490 if start is None:
491 start = 0
492 step = 2
493 slice_all = (slice(None),)*nd
494 slice0 = tupleset(slice_all, axis, slice(start, stop, step))
495 slice1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
496 slice2 = tupleset(slice_all, axis, slice(start+2, stop+2, step))
498 if x is None: # Even-spaced Simpson's rule.
499 result = np.sum(y[slice0] + 4.0*y[slice1] + y[slice2], axis=axis)
500 result *= dx / 3.0
501 else:
502 # Account for possibly different spacings.
503 # Simpson's rule changes a bit.
504 h = np.diff(x, axis=axis)
505 sl0 = tupleset(slice_all, axis, slice(start, stop, step))
506 sl1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
507 h0 = np.float64(h[sl0])
508 h1 = np.float64(h[sl1])
509 hsum = h0 + h1
510 hprod = h0 * h1
511 h0divh1 = np.true_divide(h0, h1, out=np.zeros_like(h0), where=h1 != 0)
512 tmp = hsum/6.0 * (y[slice0] *
513 (2.0 - np.true_divide(1.0, h0divh1,
514 out=np.zeros_like(h0divh1),
515 where=h0divh1 != 0)) +
516 y[slice1] * (hsum *
517 np.true_divide(hsum, hprod,
518 out=np.zeros_like(hsum),
519 where=hprod != 0)) +
520 y[slice2] * (2.0 - h0divh1))
521 result = np.sum(tmp, axis=axis)
522 return result
525# Note: alias kept for backwards compatibility. simps was renamed to simpson
526# because the former is a slur in colloquial English (see gh-12924).
527def simps(y, x=None, dx=1.0, axis=-1, even='avg'):
528 """An alias of `simpson`.
530 `simps` is kept for backwards compatibility. For new code, prefer
531 `simpson` instead.
532 """
533 return simpson(y, x=x, dx=dx, axis=axis, even=even)
536def simpson(y, x=None, dx=1.0, axis=-1, even='avg'):
537 """
538 Integrate y(x) using samples along the given axis and the composite
539 Simpson's rule. If x is None, spacing of dx is assumed.
541 If there are an even number of samples, N, then there are an odd
542 number of intervals (N-1), but Simpson's rule requires an even number
543 of intervals. The parameter 'even' controls how this is handled.
545 Parameters
546 ----------
547 y : array_like
548 Array to be integrated.
549 x : array_like, optional
550 If given, the points at which `y` is sampled.
551 dx : float, optional
552 Spacing of integration points along axis of `x`. Only used when
553 `x` is None. Default is 1.
554 axis : int, optional
555 Axis along which to integrate. Default is the last axis.
556 even : str {'avg', 'first', 'last'}, optional
557 'avg' : Average two results:1) use the first N-2 intervals with
558 a trapezoidal rule on the last interval and 2) use the last
559 N-2 intervals with a trapezoidal rule on the first interval.
561 'first' : Use Simpson's rule for the first N-2 intervals with
562 a trapezoidal rule on the last interval.
564 'last' : Use Simpson's rule for the last N-2 intervals with a
565 trapezoidal rule on the first interval.
567 Returns
568 -------
569 float
570 The estimated integral computed with the composite Simpson's rule.
572 See Also
573 --------
574 quad : adaptive quadrature using QUADPACK
575 romberg : adaptive Romberg quadrature
576 quadrature : adaptive Gaussian quadrature
577 fixed_quad : fixed-order Gaussian quadrature
578 dblquad : double integrals
579 tplquad : triple integrals
580 romb : integrators for sampled data
581 cumulative_trapezoid : cumulative integration for sampled data
582 ode : ODE integrators
583 odeint : ODE integrators
585 Notes
586 -----
587 For an odd number of samples that are equally spaced the result is
588 exact if the function is a polynomial of order 3 or less. If
589 the samples are not equally spaced, then the result is exact only
590 if the function is a polynomial of order 2 or less.
592 Examples
593 --------
594 >>> from scipy import integrate
595 >>> import numpy as np
596 >>> x = np.arange(0, 10)
597 >>> y = np.arange(0, 10)
599 >>> integrate.simpson(y, x)
600 40.5
602 >>> y = np.power(x, 3)
603 >>> integrate.simpson(y, x)
604 1642.5
605 >>> integrate.quad(lambda x: x**3, 0, 9)[0]
606 1640.25
608 >>> integrate.simpson(y, x, even='first')
609 1644.5
611 """
612 y = np.asarray(y)
613 nd = len(y.shape)
614 N = y.shape[axis]
615 last_dx = dx
616 first_dx = dx
617 returnshape = 0
618 if x is not None:
619 x = np.asarray(x)
620 if len(x.shape) == 1:
621 shapex = [1] * nd
622 shapex[axis] = x.shape[0]
623 saveshape = x.shape
624 returnshape = 1
625 x = x.reshape(tuple(shapex))
626 elif len(x.shape) != len(y.shape):
627 raise ValueError("If given, shape of x must be 1-D or the "
628 "same as y.")
629 if x.shape[axis] != N:
630 raise ValueError("If given, length of x along axis must be the "
631 "same as y.")
632 if N % 2 == 0:
633 val = 0.0
634 result = 0.0
635 slice1 = (slice(None),)*nd
636 slice2 = (slice(None),)*nd
637 if even not in ['avg', 'last', 'first']:
638 raise ValueError("Parameter 'even' must be "
639 "'avg', 'last', or 'first'.")
640 # Compute using Simpson's rule on first intervals
641 if even in ['avg', 'first']:
642 slice1 = tupleset(slice1, axis, -1)
643 slice2 = tupleset(slice2, axis, -2)
644 if x is not None:
645 last_dx = x[slice1] - x[slice2]
646 val += 0.5*last_dx*(y[slice1]+y[slice2])
647 result = _basic_simpson(y, 0, N-3, x, dx, axis)
648 # Compute using Simpson's rule on last set of intervals
649 if even in ['avg', 'last']:
650 slice1 = tupleset(slice1, axis, 0)
651 slice2 = tupleset(slice2, axis, 1)
652 if x is not None:
653 first_dx = x[tuple(slice2)] - x[tuple(slice1)]
654 val += 0.5*first_dx*(y[slice2]+y[slice1])
655 result += _basic_simpson(y, 1, N-2, x, dx, axis)
656 if even == 'avg':
657 val /= 2.0
658 result /= 2.0
659 result = result + val
660 else:
661 result = _basic_simpson(y, 0, N-2, x, dx, axis)
662 if returnshape:
663 x = x.reshape(saveshape)
664 return result
667def romb(y, dx=1.0, axis=-1, show=False):
668 """
669 Romberg integration using samples of a function.
671 Parameters
672 ----------
673 y : array_like
674 A vector of ``2**k + 1`` equally-spaced samples of a function.
675 dx : float, optional
676 The sample spacing. Default is 1.
677 axis : int, optional
678 The axis along which to integrate. Default is -1 (last axis).
679 show : bool, optional
680 When `y` is a single 1-D array, then if this argument is True
681 print the table showing Richardson extrapolation from the
682 samples. Default is False.
684 Returns
685 -------
686 romb : ndarray
687 The integrated result for `axis`.
689 See Also
690 --------
691 quad : adaptive quadrature using QUADPACK
692 romberg : adaptive Romberg quadrature
693 quadrature : adaptive Gaussian quadrature
694 fixed_quad : fixed-order Gaussian quadrature
695 dblquad : double integrals
696 tplquad : triple integrals
697 simpson : integrators for sampled data
698 cumulative_trapezoid : cumulative integration for sampled data
699 ode : ODE integrators
700 odeint : ODE integrators
702 Examples
703 --------
704 >>> from scipy import integrate
705 >>> import numpy as np
706 >>> x = np.arange(10, 14.25, 0.25)
707 >>> y = np.arange(3, 12)
709 >>> integrate.romb(y)
710 56.0
712 >>> y = np.sin(np.power(x, 2.5))
713 >>> integrate.romb(y)
714 -0.742561336672229
716 >>> integrate.romb(y, show=True)
717 Richardson Extrapolation Table for Romberg Integration
718 ======================================================
719 -0.81576
720 4.63862 6.45674
721 -1.10581 -3.02062 -3.65245
722 -2.57379 -3.06311 -3.06595 -3.05664
723 -1.34093 -0.92997 -0.78776 -0.75160 -0.74256
724 ======================================================
725 -0.742561336672229 # may vary
727 """
728 y = np.asarray(y)
729 nd = len(y.shape)
730 Nsamps = y.shape[axis]
731 Ninterv = Nsamps-1
732 n = 1
733 k = 0
734 while n < Ninterv:
735 n <<= 1
736 k += 1
737 if n != Ninterv:
738 raise ValueError("Number of samples must be one plus a "
739 "non-negative power of 2.")
741 R = {}
742 slice_all = (slice(None),) * nd
743 slice0 = tupleset(slice_all, axis, 0)
744 slicem1 = tupleset(slice_all, axis, -1)
745 h = Ninterv * np.asarray(dx, dtype=float)
746 R[(0, 0)] = (y[slice0] + y[slicem1])/2.0*h
747 slice_R = slice_all
748 start = stop = step = Ninterv
749 for i in range(1, k+1):
750 start >>= 1
751 slice_R = tupleset(slice_R, axis, slice(start, stop, step))
752 step >>= 1
753 R[(i, 0)] = 0.5*(R[(i-1, 0)] + h*y[slice_R].sum(axis=axis))
754 for j in range(1, i+1):
755 prev = R[(i, j-1)]
756 R[(i, j)] = prev + (prev-R[(i-1, j-1)]) / ((1 << (2*j))-1)
757 h /= 2.0
759 if show:
760 if not np.isscalar(R[(0, 0)]):
761 print("*** Printing table only supported for integrals" +
762 " of a single data set.")
763 else:
764 try:
765 precis = show[0]
766 except (TypeError, IndexError):
767 precis = 5
768 try:
769 width = show[1]
770 except (TypeError, IndexError):
771 width = 8
772 formstr = "%%%d.%df" % (width, precis)
774 title = "Richardson Extrapolation Table for Romberg Integration"
775 print(title, "=" * len(title), sep="\n", end="\n")
776 for i in range(k+1):
777 for j in range(i+1):
778 print(formstr % R[(i, j)], end=" ")
779 print()
780 print("=" * len(title))
782 return R[(k, k)]
784# Romberg quadratures for numeric integration.
785#
786# Written by Scott M. Ransom <ransom@cfa.harvard.edu>
787# last revision: 14 Nov 98
788#
789# Cosmetic changes by Konrad Hinsen <hinsen@cnrs-orleans.fr>
790# last revision: 1999-7-21
791#
792# Adapted to SciPy by Travis Oliphant <oliphant.travis@ieee.org>
793# last revision: Dec 2001
796def _difftrap(function, interval, numtraps):
797 """
798 Perform part of the trapezoidal rule to integrate a function.
799 Assume that we had called difftrap with all lower powers-of-2
800 starting with 1. Calling difftrap only returns the summation
801 of the new ordinates. It does _not_ multiply by the width
802 of the trapezoids. This must be performed by the caller.
803 'function' is the function to evaluate (must accept vector arguments).
804 'interval' is a sequence with lower and upper limits
805 of integration.
806 'numtraps' is the number of trapezoids to use (must be a
807 power-of-2).
808 """
809 if numtraps <= 0:
810 raise ValueError("numtraps must be > 0 in difftrap().")
811 elif numtraps == 1:
812 return 0.5*(function(interval[0])+function(interval[1]))
813 else:
814 numtosum = numtraps/2
815 h = float(interval[1]-interval[0])/numtosum
816 lox = interval[0] + 0.5 * h
817 points = lox + h * np.arange(numtosum)
818 s = np.sum(function(points), axis=0)
819 return s
822def _romberg_diff(b, c, k):
823 """
824 Compute the differences for the Romberg quadrature corrections.
825 See Forman Acton's "Real Computing Made Real," p 143.
826 """
827 tmp = 4.0**k
828 return (tmp * c - b)/(tmp - 1.0)
831def _printresmat(function, interval, resmat):
832 # Print the Romberg result matrix.
833 i = j = 0
834 print('Romberg integration of', repr(function), end=' ')
835 print('from', interval)
836 print('')
837 print('%6s %9s %9s' % ('Steps', 'StepSize', 'Results'))
838 for i in range(len(resmat)):
839 print('%6d %9f' % (2**i, (interval[1]-interval[0])/(2.**i)), end=' ')
840 for j in range(i+1):
841 print('%9f' % (resmat[i][j]), end=' ')
842 print('')
843 print('')
844 print('The final result is', resmat[i][j], end=' ')
845 print('after', 2**(len(resmat)-1)+1, 'function evaluations.')
848def romberg(function, a, b, args=(), tol=1.48e-8, rtol=1.48e-8, show=False,
849 divmax=10, vec_func=False):
850 """
851 Romberg integration of a callable function or method.
853 Returns the integral of `function` (a function of one variable)
854 over the interval (`a`, `b`).
856 If `show` is 1, the triangular array of the intermediate results
857 will be printed. If `vec_func` is True (default is False), then
858 `function` is assumed to support vector arguments.
860 Parameters
861 ----------
862 function : callable
863 Function to be integrated.
864 a : float
865 Lower limit of integration.
866 b : float
867 Upper limit of integration.
869 Returns
870 -------
871 results : float
872 Result of the integration.
874 Other Parameters
875 ----------------
876 args : tuple, optional
877 Extra arguments to pass to function. Each element of `args` will
878 be passed as a single argument to `func`. Default is to pass no
879 extra arguments.
880 tol, rtol : float, optional
881 The desired absolute and relative tolerances. Defaults are 1.48e-8.
882 show : bool, optional
883 Whether to print the results. Default is False.
884 divmax : int, optional
885 Maximum order of extrapolation. Default is 10.
886 vec_func : bool, optional
887 Whether `func` handles arrays as arguments (i.e., whether it is a
888 "vector" function). Default is False.
890 See Also
891 --------
892 fixed_quad : Fixed-order Gaussian quadrature.
893 quad : Adaptive quadrature using QUADPACK.
894 dblquad : Double integrals.
895 tplquad : Triple integrals.
896 romb : Integrators for sampled data.
897 simpson : Integrators for sampled data.
898 cumulative_trapezoid : Cumulative integration for sampled data.
899 ode : ODE integrator.
900 odeint : ODE integrator.
902 References
903 ----------
904 .. [1] 'Romberg's method' https://en.wikipedia.org/wiki/Romberg%27s_method
906 Examples
907 --------
908 Integrate a gaussian from 0 to 1 and compare to the error function.
910 >>> from scipy import integrate
911 >>> from scipy.special import erf
912 >>> import numpy as np
913 >>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
914 >>> result = integrate.romberg(gaussian, 0, 1, show=True)
915 Romberg integration of <function vfunc at ...> from [0, 1]
917 ::
919 Steps StepSize Results
920 1 1.000000 0.385872
921 2 0.500000 0.412631 0.421551
922 4 0.250000 0.419184 0.421368 0.421356
923 8 0.125000 0.420810 0.421352 0.421350 0.421350
924 16 0.062500 0.421215 0.421350 0.421350 0.421350 0.421350
925 32 0.031250 0.421317 0.421350 0.421350 0.421350 0.421350 0.421350
927 The final result is 0.421350396475 after 33 function evaluations.
929 >>> print("%g %g" % (2*result, erf(1)))
930 0.842701 0.842701
932 """
933 if np.isinf(a) or np.isinf(b):
934 raise ValueError("Romberg integration only available "
935 "for finite limits.")
936 vfunc = vectorize1(function, args, vec_func=vec_func)
937 n = 1
938 interval = [a, b]
939 intrange = b - a
940 ordsum = _difftrap(vfunc, interval, n)
941 result = intrange * ordsum
942 resmat = [[result]]
943 err = np.inf
944 last_row = resmat[0]
945 for i in range(1, divmax+1):
946 n *= 2
947 ordsum += _difftrap(vfunc, interval, n)
948 row = [intrange * ordsum / n]
949 for k in range(i):
950 row.append(_romberg_diff(last_row[k], row[k], k+1))
951 result = row[i]
952 lastresult = last_row[i-1]
953 if show:
954 resmat.append(row)
955 err = abs(result - lastresult)
956 if err < tol or err < rtol * abs(result):
957 break
958 last_row = row
959 else:
960 warnings.warn(
961 "divmax (%d) exceeded. Latest difference = %e" % (divmax, err),
962 AccuracyWarning)
964 if show:
965 _printresmat(vfunc, interval, resmat)
966 return result
969# Coefficients for Newton-Cotes quadrature
970#
971# These are the points being used
972# to construct the local interpolating polynomial
973# a are the weights for Newton-Cotes integration
974# B is the error coefficient.
975# error in these coefficients grows as N gets larger.
976# or as samples are closer and closer together
978# You can use maxima to find these rational coefficients
979# for equally spaced data using the commands
980# a(i,N) := integrate(product(r-j,j,0,i-1) * product(r-j,j,i+1,N),r,0,N) / ((N-i)! * i!) * (-1)^(N-i);
981# Be(N) := N^(N+2)/(N+2)! * (N/(N+3) - sum((i/N)^(N+2)*a(i,N),i,0,N));
982# Bo(N) := N^(N+1)/(N+1)! * (N/(N+2) - sum((i/N)^(N+1)*a(i,N),i,0,N));
983# B(N) := (if (mod(N,2)=0) then Be(N) else Bo(N));
984#
985# pre-computed for equally-spaced weights
986#
987# num_a, den_a, int_a, num_B, den_B = _builtincoeffs[N]
988#
989# a = num_a*array(int_a)/den_a
990# B = num_B*1.0 / den_B
991#
992# integrate(f(x),x,x_0,x_N) = dx*sum(a*f(x_i)) + B*(dx)^(2k+3) f^(2k+2)(x*)
993# where k = N // 2
994#
995_builtincoeffs = {
996 1: (1,2,[1,1],-1,12),
997 2: (1,3,[1,4,1],-1,90),
998 3: (3,8,[1,3,3,1],-3,80),
999 4: (2,45,[7,32,12,32,7],-8,945),
1000 5: (5,288,[19,75,50,50,75,19],-275,12096),
1001 6: (1,140,[41,216,27,272,27,216,41],-9,1400),
1002 7: (7,17280,[751,3577,1323,2989,2989,1323,3577,751],-8183,518400),
1003 8: (4,14175,[989,5888,-928,10496,-4540,10496,-928,5888,989],
1004 -2368,467775),
1005 9: (9,89600,[2857,15741,1080,19344,5778,5778,19344,1080,
1006 15741,2857], -4671, 394240),
1007 10: (5,299376,[16067,106300,-48525,272400,-260550,427368,
1008 -260550,272400,-48525,106300,16067],
1009 -673175, 163459296),
1010 11: (11,87091200,[2171465,13486539,-3237113, 25226685,-9595542,
1011 15493566,15493566,-9595542,25226685,-3237113,
1012 13486539,2171465], -2224234463, 237758976000),
1013 12: (1, 5255250, [1364651,9903168,-7587864,35725120,-51491295,
1014 87516288,-87797136,87516288,-51491295,35725120,
1015 -7587864,9903168,1364651], -3012, 875875),
1016 13: (13, 402361344000,[8181904909, 56280729661, -31268252574,
1017 156074417954,-151659573325,206683437987,
1018 -43111992612,-43111992612,206683437987,
1019 -151659573325,156074417954,-31268252574,
1020 56280729661,8181904909], -2639651053,
1021 344881152000),
1022 14: (7, 2501928000, [90241897,710986864,-770720657,3501442784,
1023 -6625093363,12630121616,-16802270373,19534438464,
1024 -16802270373,12630121616,-6625093363,3501442784,
1025 -770720657,710986864,90241897], -3740727473,
1026 1275983280000)
1027 }
1030def newton_cotes(rn, equal=0):
1031 r"""
1032 Return weights and error coefficient for Newton-Cotes integration.
1034 Suppose we have (N+1) samples of f at the positions
1035 x_0, x_1, ..., x_N. Then an N-point Newton-Cotes formula for the
1036 integral between x_0 and x_N is:
1038 :math:`\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i)
1039 + B_N (\Delta x)^{N+2} f^{N+1} (\xi)`
1041 where :math:`\xi \in [x_0,x_N]`
1042 and :math:`\Delta x = \frac{x_N-x_0}{N}` is the average samples spacing.
1044 If the samples are equally-spaced and N is even, then the error
1045 term is :math:`B_N (\Delta x)^{N+3} f^{N+2}(\xi)`.
1047 Parameters
1048 ----------
1049 rn : int
1050 The integer order for equally-spaced data or the relative positions of
1051 the samples with the first sample at 0 and the last at N, where N+1 is
1052 the length of `rn`. N is the order of the Newton-Cotes integration.
1053 equal : int, optional
1054 Set to 1 to enforce equally spaced data.
1056 Returns
1057 -------
1058 an : ndarray
1059 1-D array of weights to apply to the function at the provided sample
1060 positions.
1061 B : float
1062 Error coefficient.
1064 Notes
1065 -----
1066 Normally, the Newton-Cotes rules are used on smaller integration
1067 regions and a composite rule is used to return the total integral.
1069 Examples
1070 --------
1071 Compute the integral of sin(x) in [0, :math:`\pi`]:
1073 >>> from scipy.integrate import newton_cotes
1074 >>> import numpy as np
1075 >>> def f(x):
1076 ... return np.sin(x)
1077 >>> a = 0
1078 >>> b = np.pi
1079 >>> exact = 2
1080 >>> for N in [2, 4, 6, 8, 10]:
1081 ... x = np.linspace(a, b, N + 1)
1082 ... an, B = newton_cotes(N, 1)
1083 ... dx = (b - a) / N
1084 ... quad = dx * np.sum(an * f(x))
1085 ... error = abs(quad - exact)
1086 ... print('{:2d} {:10.9f} {:.5e}'.format(N, quad, error))
1087 ...
1088 2 2.094395102 9.43951e-02
1089 4 1.998570732 1.42927e-03
1090 6 2.000017814 1.78136e-05
1091 8 1.999999835 1.64725e-07
1092 10 2.000000001 1.14677e-09
1094 """
1095 try:
1096 N = len(rn)-1
1097 if equal:
1098 rn = np.arange(N+1)
1099 elif np.all(np.diff(rn) == 1):
1100 equal = 1
1101 except Exception:
1102 N = rn
1103 rn = np.arange(N+1)
1104 equal = 1
1106 if equal and N in _builtincoeffs:
1107 na, da, vi, nb, db = _builtincoeffs[N]
1108 an = na * np.array(vi, dtype=float) / da
1109 return an, float(nb)/db
1111 if (rn[0] != 0) or (rn[-1] != N):
1112 raise ValueError("The sample positions must start at 0"
1113 " and end at N")
1114 yi = rn / float(N)
1115 ti = 2 * yi - 1
1116 nvec = np.arange(N+1)
1117 C = ti ** nvec[:, np.newaxis]
1118 Cinv = np.linalg.inv(C)
1119 # improve precision of result
1120 for i in range(2):
1121 Cinv = 2*Cinv - Cinv.dot(C).dot(Cinv)
1122 vec = 2.0 / (nvec[::2]+1)
1123 ai = Cinv[:, ::2].dot(vec) * (N / 2.)
1125 if (N % 2 == 0) and equal:
1126 BN = N/(N+3.)
1127 power = N+2
1128 else:
1129 BN = N/(N+2.)
1130 power = N+1
1132 BN = BN - np.dot(yi**power, ai)
1133 p1 = power+1
1134 fac = power*math.log(N) - gammaln(p1)
1135 fac = math.exp(fac)
1136 return ai, BN*fac
1139def _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log):
1141 # lazy import to avoid issues with partially-initialized submodule
1142 if not hasattr(_qmc_quad, 'qmc'):
1143 from scipy import stats
1144 _qmc_quad.stats = stats
1145 else:
1146 stats = _qmc_quad.stats
1148 if not callable(func):
1149 message = "`func` must be callable."
1150 raise TypeError(message)
1152 # a, b will be modified, so copy. Oh well if it's copied twice.
1153 a = np.atleast_1d(a).copy()
1154 b = np.atleast_1d(b).copy()
1155 a, b = np.broadcast_arrays(a, b)
1156 dim = a.shape[0]
1158 try:
1159 func((a + b) / 2)
1160 except Exception as e:
1161 message = ("`func` must evaluate the integrand at points within "
1162 "the integration range; e.g. `func( (a + b) / 2)` "
1163 "must return the integrand at the centroid of the "
1164 "integration volume.")
1165 raise ValueError(message) from e
1167 try:
1168 func(np.array([a, b]))
1169 vfunc = func
1170 except Exception as e:
1171 message = ("Exception encountered when attempting vectorized call to "
1172 f"`func`: {e}. `func` should accept two-dimensional array "
1173 "with shape `(n_points, len(a))` and return an array with "
1174 "the integrand value at each of the `n_points` for better "
1175 "performance.")
1176 warnings.warn(message, stacklevel=3)
1178 def vfunc(x):
1179 return np.apply_along_axis(func, axis=-1, arr=x)
1181 n_points_int = np.int64(n_points)
1182 if n_points != n_points_int:
1183 message = "`n_points` must be an integer."
1184 raise TypeError(message)
1186 n_estimates_int = np.int64(n_estimates)
1187 if n_estimates != n_estimates_int:
1188 message = "`n_estimates` must be an integer."
1189 raise TypeError(message)
1191 if qrng is None:
1192 qrng = stats.qmc.Halton(dim)
1193 elif not isinstance(qrng, stats.qmc.QMCEngine):
1194 message = "`qrng` must be an instance of scipy.stats.qmc.QMCEngine."
1195 raise TypeError(message)
1197 if qrng.d != a.shape[0]:
1198 message = ("`qrng` must be initialized with dimensionality equal to "
1199 "the number of variables in `a`, i.e., "
1200 "`qrng.random().shape[-1]` must equal `a.shape[0]`.")
1201 raise ValueError(message)
1203 rng_seed = getattr(qrng, 'rng_seed', None)
1204 rng = stats._qmc.check_random_state(rng_seed)
1206 if log not in {True, False}:
1207 message = "`log` must be boolean (`True` or `False`)."
1208 raise TypeError(message)
1210 return (vfunc, a, b, n_points_int, n_estimates_int, qrng, rng, log, stats)
1213QMCQuadResult = namedtuple('QMCQuadResult', ['integral', 'standard_error'])
1216def _qmc_quad(func, a, b, *, n_points=1024, n_estimates=8, qrng=None,
1217 log=False, args=None):
1218 """
1219 Compute an integral in N-dimensions using Quasi-Monte Carlo quadrature.
1221 Parameters
1222 ----------
1223 func : callable
1224 The integrand. Must accept a single arguments `x`, an array which
1225 specifies the point at which to evaluate the integrand. For efficiency,
1226 the function should be vectorized to compute the integrand for each
1227 element an array of shape ``(n_points, n)``, where ``n`` is number of
1228 variables.
1229 a, b : array-like
1230 One-dimensional arrays specifying the lower and upper integration
1231 limits, respectively, of each of the ``n`` variables.
1232 n_points, n_estimates : int, optional
1233 One QMC sample of `n_points` (default: 256) points will be generated
1234 by `qrng`, and `n_estimates` (default: 8) statistically independent
1235 estimates of the integral will be produced. The total number of points
1236 at which the integrand `func` will be evaluated is
1237 ``n_points * n_estimates``. See Notes for details.
1238 qrng : `~scipy.stats.qmc.QMCEngine`, optional
1239 An instance of the QMCEngine from which to sample QMC points.
1240 The QMCEngine must be initialized to a number of dimensions
1241 corresponding with the number of variables ``x0, ..., xn`` passed to
1242 `func`.
1243 The provided QMCEngine is used to produce the first integral estimate.
1244 If `n_estimates` is greater than one, additional QMCEngines are
1245 spawned from the first (with scrambling enabled, if it is an option.)
1246 If a QMCEngine is not provided, the default `scipy.stats.qmc.Halton`
1247 will be initialized with the number of dimensions determine from
1248 `a`.
1249 log : boolean, default: False
1250 When set to True, `func` returns the log of the integrand, and
1251 the result object contains the log of the integral.
1253 Returns
1254 -------
1255 result : object
1256 A result object with attributes:
1258 integral : float
1259 The estimate of the integral.
1260 standard_error :
1261 The error estimate. See Notes for interpretation.
1263 Notes
1264 -----
1265 Values of the integrand at each of the `n_points` points of a QMC sample
1266 are used to produce an estimate of the integral. This estimate is drawn
1267 from a population of possible estimates of the integral, the value of
1268 which we obtain depends on the particular points at which the integral
1269 was evaluated. We perform this process `n_estimates` times, each time
1270 evaluating the integrand at different scrambled QMC points, effectively
1271 drawing i.i.d. random samples from the population of integral estimates.
1272 The sample mean :math:`m` of these integral estimates is an
1273 unbiased estimator of the true value of the integral, and the standard
1274 error of the mean :math:`s` of these estimates may be used to generate
1275 confidence intervals using the t distribution with ``n_estimates - 1``
1276 degrees of freedom. Perhaps counter-intuitively, increasing `n_points`
1277 while keeping the total number of function evaluation points
1278 ``n_points * n_estimates`` fixed tends to reduce the actual error, whereas
1279 increasing `n_estimates` tends to decrease the error estimate.
1281 Examples
1282 --------
1283 QMC quadrature is particularly useful for computing integrals in higher
1284 dimensions. An example integrand is the probability density function
1285 of a multivariate normal distribution.
1287 >>> import numpy as np
1288 >>> from scipy import stats
1289 >>> dim = 8
1290 >>> mean = np.zeros(dim)
1291 >>> cov = np.eye(dim)
1292 >>> def func(x):
1293 ... return stats.multivariate_normal.pdf(x, mean, cov)
1295 To compute the integral over the unit hypercube:
1297 >>> from scipy.integrate import qmc_quad
1298 >>> a = np.zeros(dim)
1299 >>> b = np.ones(dim)
1300 >>> rng = np.random.default_rng()
1301 >>> qrng = stats.qmc.Halton(d=dim, seed=rng)
1302 >>> n_estimates = 8
1303 >>> res = qmc_quad(func, a, b, n_estimates=n_estimates, qrng=qrng)
1304 >>> res.integral, res.standard_error
1305 (0.00018441088533413305, 1.1255608140911588e-07)
1307 A two-sided, 99% confidence interval for the integral may be estimated
1308 as:
1310 >>> t = stats.t(df=n_estimates-1, loc=res.integral,
1311 ... scale=res.standard_error)
1312 >>> t.interval(0.99)
1313 (0.00018401699720722663, 0.00018480477346103947)
1315 Indeed, the value reported by `scipy.stats.multivariate_normal` is
1316 within this range.
1318 >>> stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
1319 0.00018430867675187443
1321 """
1322 args = _qmc_quad_iv(func, a, b, n_points, n_estimates, qrng, log)
1323 func, a, b, n_points, n_estimates, qrng, rng, log, stats = args
1325 # The sign of the integral depends on the order of the limits. Fix this by
1326 # ensuring that lower bounds are indeed lower and setting sign of resulting
1327 # integral manually
1328 if np.any(a == b):
1329 message = ("A lower limit was equal to an upper limit, so the value "
1330 "of the integral is zero by definition.")
1331 warnings.warn(message, stacklevel=2)
1332 return QMCQuadResult(-np.inf if log else 0, 0)
1334 i_swap = b < a
1335 sign = (-1)**(i_swap.sum(axis=-1)) # odd # of swaps -> negative
1336 a[i_swap], b[i_swap] = b[i_swap], a[i_swap]
1338 A = np.prod(b - a)
1339 dA = A / n_points
1341 estimates = np.zeros(n_estimates)
1342 rngs = _rng_spawn(qrng.rng, n_estimates)
1343 for i in range(n_estimates):
1344 # Generate integral estimate
1345 sample = qrng.random(n_points)
1346 x = stats.qmc.scale(sample, a, b)
1347 integrands = func(x)
1348 if log:
1349 estimate = logsumexp(integrands) + np.log(dA)
1350 else:
1351 estimate = np.sum(integrands * dA)
1352 estimates[i] = estimate
1354 # Get a new, independently-scrambled QRNG for next time
1355 qrng = type(qrng)(seed=rngs[i], **qrng._init_quad)
1357 integral = np.mean(estimates)
1358 integral = integral + np.pi*1j if (log and sign < 0) else integral*sign
1359 standard_error = stats.sem(estimates)
1360 return QMCQuadResult(integral, standard_error)