Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/integrate/_quadrature.py : 11%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import functools
2import numpy as np
3import math
4import types
5import warnings
7# trapz is a public function for scipy.integrate,
8# even though it's actually a NumPy function.
9from numpy import trapz
10from scipy.special import roots_legendre
11from scipy.special import gammaln
13__all__ = ['fixed_quad', 'quadrature', 'romberg', 'trapz', 'simps', 'romb',
14 'cumtrapz', 'newton_cotes', 'AccuracyWarning']
17# Make See Also linking for our local copy work properly
18def _copy_func(f):
19 """Based on http://stackoverflow.com/a/6528148/190597 (Glenn Maynard)"""
20 g = types.FunctionType(f.__code__, f.__globals__, name=f.__name__,
21 argdefs=f.__defaults__, closure=f.__closure__)
22 g = functools.update_wrapper(g, f)
23 g.__kwdefaults__ = f.__kwdefaults__
24 return g
27trapz = _copy_func(trapz)
28if trapz.__doc__:
29 trapz.__doc__ = trapz.__doc__.replace('sum, cumsum', 'numpy.cumsum')
32class AccuracyWarning(Warning):
33 pass
36def _cached_roots_legendre(n):
37 """
38 Cache roots_legendre results to speed up calls of the fixed_quad
39 function.
40 """
41 if n in _cached_roots_legendre.cache:
42 return _cached_roots_legendre.cache[n]
44 _cached_roots_legendre.cache[n] = roots_legendre(n)
45 return _cached_roots_legendre.cache[n]
48_cached_roots_legendre.cache = dict()
51def fixed_quad(func, a, b, args=(), n=5):
52 """
53 Compute a definite integral using fixed-order Gaussian quadrature.
55 Integrate `func` from `a` to `b` using Gaussian quadrature of
56 order `n`.
58 Parameters
59 ----------
60 func : callable
61 A Python function or method to integrate (must accept vector inputs).
62 If integrating a vector-valued function, the returned array must have
63 shape ``(..., len(x))``.
64 a : float
65 Lower limit of integration.
66 b : float
67 Upper limit of integration.
68 args : tuple, optional
69 Extra arguments to pass to function, if any.
70 n : int, optional
71 Order of quadrature integration. Default is 5.
73 Returns
74 -------
75 val : float
76 Gaussian quadrature approximation to the integral
77 none : None
78 Statically returned value of None
81 See Also
82 --------
83 quad : adaptive quadrature using QUADPACK
84 dblquad : double integrals
85 tplquad : triple integrals
86 romberg : adaptive Romberg quadrature
87 quadrature : adaptive Gaussian quadrature
88 romb : integrators for sampled data
89 simps : integrators for sampled data
90 cumtrapz : cumulative integration for sampled data
91 ode : ODE integrator
92 odeint : ODE integrator
94 Examples
95 --------
96 >>> from scipy import integrate
97 >>> f = lambda x: x**8
98 >>> integrate.fixed_quad(f, 0.0, 1.0, n=4)
99 (0.1110884353741496, None)
100 >>> integrate.fixed_quad(f, 0.0, 1.0, n=5)
101 (0.11111111111111102, None)
102 >>> print(1/9.0) # analytical result
103 0.1111111111111111
105 >>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=4)
106 (0.9999999771971152, None)
107 >>> integrate.fixed_quad(np.cos, 0.0, np.pi/2, n=5)
108 (1.000000000039565, None)
109 >>> np.sin(np.pi/2)-np.sin(0) # analytical result
110 1.0
112 """
113 x, w = _cached_roots_legendre(n)
114 x = np.real(x)
115 if np.isinf(a) or np.isinf(b):
116 raise ValueError("Gaussian quadrature is only available for "
117 "finite limits.")
118 y = (b-a)*(x+1)/2.0 + a
119 return (b-a)/2.0 * np.sum(w*func(y, *args), axis=-1), None
122def vectorize1(func, args=(), vec_func=False):
123 """Vectorize the call to a function.
125 This is an internal utility function used by `romberg` and
126 `quadrature` to create a vectorized version of a function.
128 If `vec_func` is True, the function `func` is assumed to take vector
129 arguments.
131 Parameters
132 ----------
133 func : callable
134 User defined function.
135 args : tuple, optional
136 Extra arguments for the function.
137 vec_func : bool, optional
138 True if the function func takes vector arguments.
140 Returns
141 -------
142 vfunc : callable
143 A function that will take a vector argument and return the
144 result.
146 """
147 if vec_func:
148 def vfunc(x):
149 return func(x, *args)
150 else:
151 def vfunc(x):
152 if np.isscalar(x):
153 return func(x, *args)
154 x = np.asarray(x)
155 # call with first point to get output type
156 y0 = func(x[0], *args)
157 n = len(x)
158 dtype = getattr(y0, 'dtype', type(y0))
159 output = np.empty((n,), dtype=dtype)
160 output[0] = y0
161 for i in range(1, n):
162 output[i] = func(x[i], *args)
163 return output
164 return vfunc
167def quadrature(func, a, b, args=(), tol=1.49e-8, rtol=1.49e-8, maxiter=50,
168 vec_func=True, miniter=1):
169 """
170 Compute a definite integral using fixed-tolerance Gaussian quadrature.
172 Integrate `func` from `a` to `b` using Gaussian quadrature
173 with absolute tolerance `tol`.
175 Parameters
176 ----------
177 func : function
178 A Python function or method to integrate.
179 a : float
180 Lower limit of integration.
181 b : float
182 Upper limit of integration.
183 args : tuple, optional
184 Extra arguments to pass to function.
185 tol, rtol : float, optional
186 Iteration stops when error between last two iterates is less than
187 `tol` OR the relative change is less than `rtol`.
188 maxiter : int, optional
189 Maximum order of Gaussian quadrature.
190 vec_func : bool, optional
191 True or False if func handles arrays as arguments (is
192 a "vector" function). Default is True.
193 miniter : int, optional
194 Minimum order of Gaussian quadrature.
196 Returns
197 -------
198 val : float
199 Gaussian quadrature approximation (within tolerance) to integral.
200 err : float
201 Difference between last two estimates of the integral.
203 See also
204 --------
205 romberg: adaptive Romberg quadrature
206 fixed_quad: fixed-order Gaussian quadrature
207 quad: adaptive quadrature using QUADPACK
208 dblquad: double integrals
209 tplquad: triple integrals
210 romb: integrator for sampled data
211 simps: integrator for sampled data
212 cumtrapz: cumulative integration for sampled data
213 ode: ODE integrator
214 odeint: ODE integrator
216 Examples
217 --------
218 >>> from scipy import integrate
219 >>> f = lambda x: x**8
220 >>> integrate.quadrature(f, 0.0, 1.0)
221 (0.11111111111111106, 4.163336342344337e-17)
222 >>> print(1/9.0) # analytical result
223 0.1111111111111111
225 >>> integrate.quadrature(np.cos, 0.0, np.pi/2)
226 (0.9999999999999536, 3.9611425250996035e-11)
227 >>> np.sin(np.pi/2)-np.sin(0) # analytical result
228 1.0
230 """
231 if not isinstance(args, tuple):
232 args = (args,)
233 vfunc = vectorize1(func, args, vec_func=vec_func)
234 val = np.inf
235 err = np.inf
236 maxiter = max(miniter+1, maxiter)
237 for n in range(miniter, maxiter+1):
238 newval = fixed_quad(vfunc, a, b, (), n)[0]
239 err = abs(newval-val)
240 val = newval
242 if err < tol or err < rtol*abs(val):
243 break
244 else:
245 warnings.warn(
246 "maxiter (%d) exceeded. Latest difference = %e" % (maxiter, err),
247 AccuracyWarning)
248 return val, err
251def tupleset(t, i, value):
252 l = list(t)
253 l[i] = value
254 return tuple(l)
257def cumtrapz(y, x=None, dx=1.0, axis=-1, initial=None):
258 """
259 Cumulatively integrate y(x) using the composite trapezoidal rule.
261 Parameters
262 ----------
263 y : array_like
264 Values to integrate.
265 x : array_like, optional
266 The coordinate to integrate along. If None (default), use spacing `dx`
267 between consecutive elements in `y`.
268 dx : float, optional
269 Spacing between elements of `y`. Only used if `x` is None.
270 axis : int, optional
271 Specifies the axis to cumulate. Default is -1 (last axis).
272 initial : scalar, optional
273 If given, insert this value at the beginning of the returned result.
274 Typically this value should be 0. Default is None, which means no
275 value at ``x[0]`` is returned and `res` has one element less than `y`
276 along the axis of integration.
278 Returns
279 -------
280 res : ndarray
281 The result of cumulative integration of `y` along `axis`.
282 If `initial` is None, the shape is such that the axis of integration
283 has one less value than `y`. If `initial` is given, the shape is equal
284 to that of `y`.
286 See Also
287 --------
288 numpy.cumsum, numpy.cumprod
289 quad: adaptive quadrature using QUADPACK
290 romberg: adaptive Romberg quadrature
291 quadrature: adaptive Gaussian quadrature
292 fixed_quad: fixed-order Gaussian quadrature
293 dblquad: double integrals
294 tplquad: triple integrals
295 romb: integrators for sampled data
296 ode: ODE integrators
297 odeint: ODE integrators
299 Examples
300 --------
301 >>> from scipy import integrate
302 >>> import matplotlib.pyplot as plt
304 >>> x = np.linspace(-2, 2, num=20)
305 >>> y = x
306 >>> y_int = integrate.cumtrapz(y, x, initial=0)
307 >>> plt.plot(x, y_int, 'ro', x, y[0] + 0.5 * x**2, 'b-')
308 >>> plt.show()
310 """
311 y = np.asarray(y)
312 if x is None:
313 d = dx
314 else:
315 x = np.asarray(x)
316 if x.ndim == 1:
317 d = np.diff(x)
318 # reshape to correct shape
319 shape = [1] * y.ndim
320 shape[axis] = -1
321 d = d.reshape(shape)
322 elif len(x.shape) != len(y.shape):
323 raise ValueError("If given, shape of x must be 1-D or the "
324 "same as y.")
325 else:
326 d = np.diff(x, axis=axis)
328 if d.shape[axis] != y.shape[axis] - 1:
329 raise ValueError("If given, length of x along axis must be the "
330 "same as y.")
332 nd = len(y.shape)
333 slice1 = tupleset((slice(None),)*nd, axis, slice(1, None))
334 slice2 = tupleset((slice(None),)*nd, axis, slice(None, -1))
335 res = np.cumsum(d * (y[slice1] + y[slice2]) / 2.0, axis=axis)
337 if initial is not None:
338 if not np.isscalar(initial):
339 raise ValueError("`initial` parameter should be a scalar.")
341 shape = list(res.shape)
342 shape[axis] = 1
343 res = np.concatenate([np.full(shape, initial, dtype=res.dtype), res],
344 axis=axis)
346 return res
349def _basic_simps(y, start, stop, x, dx, axis):
350 nd = len(y.shape)
351 if start is None:
352 start = 0
353 step = 2
354 slice_all = (slice(None),)*nd
355 slice0 = tupleset(slice_all, axis, slice(start, stop, step))
356 slice1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
357 slice2 = tupleset(slice_all, axis, slice(start+2, stop+2, step))
359 if x is None: # Even-spaced Simpson's rule.
360 result = np.sum(dx/3.0 * (y[slice0]+4*y[slice1]+y[slice2]),
361 axis=axis)
362 else:
363 # Account for possibly different spacings.
364 # Simpson's rule changes a bit.
365 h = np.diff(x, axis=axis)
366 sl0 = tupleset(slice_all, axis, slice(start, stop, step))
367 sl1 = tupleset(slice_all, axis, slice(start+1, stop+1, step))
368 h0 = h[sl0]
369 h1 = h[sl1]
370 hsum = h0 + h1
371 hprod = h0 * h1
372 h0divh1 = h0 / h1
373 tmp = hsum/6.0 * (y[slice0]*(2-1.0/h0divh1) +
374 y[slice1]*hsum*hsum/hprod +
375 y[slice2]*(2-h0divh1))
376 result = np.sum(tmp, axis=axis)
377 return result
380def simps(y, x=None, dx=1, axis=-1, even='avg'):
381 """
382 Integrate y(x) using samples along the given axis and the composite
383 Simpson's rule. If x is None, spacing of dx is assumed.
385 If there are an even number of samples, N, then there are an odd
386 number of intervals (N-1), but Simpson's rule requires an even number
387 of intervals. The parameter 'even' controls how this is handled.
389 Parameters
390 ----------
391 y : array_like
392 Array to be integrated.
393 x : array_like, optional
394 If given, the points at which `y` is sampled.
395 dx : int, optional
396 Spacing of integration points along axis of `x`. Only used when
397 `x` is None. Default is 1.
398 axis : int, optional
399 Axis along which to integrate. Default is the last axis.
400 even : str {'avg', 'first', 'last'}, optional
401 'avg' : Average two results:1) use the first N-2 intervals with
402 a trapezoidal rule on the last interval and 2) use the last
403 N-2 intervals with a trapezoidal rule on the first interval.
405 'first' : Use Simpson's rule for the first N-2 intervals with
406 a trapezoidal rule on the last interval.
408 'last' : Use Simpson's rule for the last N-2 intervals with a
409 trapezoidal rule on the first interval.
411 See Also
412 --------
413 quad: adaptive quadrature using QUADPACK
414 romberg: adaptive Romberg quadrature
415 quadrature: adaptive Gaussian quadrature
416 fixed_quad: fixed-order Gaussian quadrature
417 dblquad: double integrals
418 tplquad: triple integrals
419 romb: integrators for sampled data
420 cumtrapz: cumulative integration for sampled data
421 ode: ODE integrators
422 odeint: ODE integrators
424 Notes
425 -----
426 For an odd number of samples that are equally spaced the result is
427 exact if the function is a polynomial of order 3 or less. If
428 the samples are not equally spaced, then the result is exact only
429 if the function is a polynomial of order 2 or less.
431 Examples
432 --------
433 >>> from scipy import integrate
434 >>> x = np.arange(0, 10)
435 >>> y = np.arange(0, 10)
437 >>> integrate.simps(y, x)
438 40.5
440 >>> y = np.power(x, 3)
441 >>> integrate.simps(y, x)
442 1642.5
443 >>> integrate.quad(lambda x: x**3, 0, 9)[0]
444 1640.25
446 >>> integrate.simps(y, x, even='first')
447 1644.5
449 """
450 y = np.asarray(y)
451 nd = len(y.shape)
452 N = y.shape[axis]
453 last_dx = dx
454 first_dx = dx
455 returnshape = 0
456 if x is not None:
457 x = np.asarray(x)
458 if len(x.shape) == 1:
459 shapex = [1] * nd
460 shapex[axis] = x.shape[0]
461 saveshape = x.shape
462 returnshape = 1
463 x = x.reshape(tuple(shapex))
464 elif len(x.shape) != len(y.shape):
465 raise ValueError("If given, shape of x must be 1-D or the "
466 "same as y.")
467 if x.shape[axis] != N:
468 raise ValueError("If given, length of x along axis must be the "
469 "same as y.")
470 if N % 2 == 0:
471 val = 0.0
472 result = 0.0
473 slice1 = (slice(None),)*nd
474 slice2 = (slice(None),)*nd
475 if even not in ['avg', 'last', 'first']:
476 raise ValueError("Parameter 'even' must be "
477 "'avg', 'last', or 'first'.")
478 # Compute using Simpson's rule on first intervals
479 if even in ['avg', 'first']:
480 slice1 = tupleset(slice1, axis, -1)
481 slice2 = tupleset(slice2, axis, -2)
482 if x is not None:
483 last_dx = x[slice1] - x[slice2]
484 val += 0.5*last_dx*(y[slice1]+y[slice2])
485 result = _basic_simps(y, 0, N-3, x, dx, axis)
486 # Compute using Simpson's rule on last set of intervals
487 if even in ['avg', 'last']:
488 slice1 = tupleset(slice1, axis, 0)
489 slice2 = tupleset(slice2, axis, 1)
490 if x is not None:
491 first_dx = x[tuple(slice2)] - x[tuple(slice1)]
492 val += 0.5*first_dx*(y[slice2]+y[slice1])
493 result += _basic_simps(y, 1, N-2, x, dx, axis)
494 if even == 'avg':
495 val /= 2.0
496 result /= 2.0
497 result = result + val
498 else:
499 result = _basic_simps(y, 0, N-2, x, dx, axis)
500 if returnshape:
501 x = x.reshape(saveshape)
502 return result
505def romb(y, dx=1.0, axis=-1, show=False):
506 """
507 Romberg integration using samples of a function.
509 Parameters
510 ----------
511 y : array_like
512 A vector of ``2**k + 1`` equally-spaced samples of a function.
513 dx : float, optional
514 The sample spacing. Default is 1.
515 axis : int, optional
516 The axis along which to integrate. Default is -1 (last axis).
517 show : bool, optional
518 When `y` is a single 1-D array, then if this argument is True
519 print the table showing Richardson extrapolation from the
520 samples. Default is False.
522 Returns
523 -------
524 romb : ndarray
525 The integrated result for `axis`.
527 See also
528 --------
529 quad : adaptive quadrature using QUADPACK
530 romberg : adaptive Romberg quadrature
531 quadrature : adaptive Gaussian quadrature
532 fixed_quad : fixed-order Gaussian quadrature
533 dblquad : double integrals
534 tplquad : triple integrals
535 simps : integrators for sampled data
536 cumtrapz : cumulative integration for sampled data
537 ode : ODE integrators
538 odeint : ODE integrators
540 Examples
541 --------
542 >>> from scipy import integrate
543 >>> x = np.arange(10, 14.25, 0.25)
544 >>> y = np.arange(3, 12)
546 >>> integrate.romb(y)
547 56.0
549 >>> y = np.sin(np.power(x, 2.5))
550 >>> integrate.romb(y)
551 -0.742561336672229
553 >>> integrate.romb(y, show=True)
554 Richardson Extrapolation Table for Romberg Integration
555 ====================================================================
556 -0.81576
557 4.63862 6.45674
558 -1.10581 -3.02062 -3.65245
559 -2.57379 -3.06311 -3.06595 -3.05664
560 -1.34093 -0.92997 -0.78776 -0.75160 -0.74256
561 ====================================================================
562 -0.742561336672229
563 """
564 y = np.asarray(y)
565 nd = len(y.shape)
566 Nsamps = y.shape[axis]
567 Ninterv = Nsamps-1
568 n = 1
569 k = 0
570 while n < Ninterv:
571 n <<= 1
572 k += 1
573 if n != Ninterv:
574 raise ValueError("Number of samples must be one plus a "
575 "non-negative power of 2.")
577 R = {}
578 slice_all = (slice(None),) * nd
579 slice0 = tupleset(slice_all, axis, 0)
580 slicem1 = tupleset(slice_all, axis, -1)
581 h = Ninterv * np.asarray(dx, dtype=float)
582 R[(0, 0)] = (y[slice0] + y[slicem1])/2.0*h
583 slice_R = slice_all
584 start = stop = step = Ninterv
585 for i in range(1, k+1):
586 start >>= 1
587 slice_R = tupleset(slice_R, axis, slice(start, stop, step))
588 step >>= 1
589 R[(i, 0)] = 0.5*(R[(i-1, 0)] + h*y[slice_R].sum(axis=axis))
590 for j in range(1, i+1):
591 prev = R[(i, j-1)]
592 R[(i, j)] = prev + (prev-R[(i-1, j-1)]) / ((1 << (2*j))-1)
593 h /= 2.0
595 if show:
596 if not np.isscalar(R[(0, 0)]):
597 print("*** Printing table only supported for integrals" +
598 " of a single data set.")
599 else:
600 try:
601 precis = show[0]
602 except (TypeError, IndexError):
603 precis = 5
604 try:
605 width = show[1]
606 except (TypeError, IndexError):
607 width = 8
608 formstr = "%%%d.%df" % (width, precis)
610 title = "Richardson Extrapolation Table for Romberg Integration"
611 print("", title.center(68), "=" * 68, sep="\n", end="\n")
612 for i in range(k+1):
613 for j in range(i+1):
614 print(formstr % R[(i, j)], end=" ")
615 print()
616 print("=" * 68)
617 print()
619 return R[(k, k)]
621# Romberg quadratures for numeric integration.
622#
623# Written by Scott M. Ransom <ransom@cfa.harvard.edu>
624# last revision: 14 Nov 98
625#
626# Cosmetic changes by Konrad Hinsen <hinsen@cnrs-orleans.fr>
627# last revision: 1999-7-21
628#
629# Adapted to SciPy by Travis Oliphant <oliphant.travis@ieee.org>
630# last revision: Dec 2001
633def _difftrap(function, interval, numtraps):
634 """
635 Perform part of the trapezoidal rule to integrate a function.
636 Assume that we had called difftrap with all lower powers-of-2
637 starting with 1. Calling difftrap only returns the summation
638 of the new ordinates. It does _not_ multiply by the width
639 of the trapezoids. This must be performed by the caller.
640 'function' is the function to evaluate (must accept vector arguments).
641 'interval' is a sequence with lower and upper limits
642 of integration.
643 'numtraps' is the number of trapezoids to use (must be a
644 power-of-2).
645 """
646 if numtraps <= 0:
647 raise ValueError("numtraps must be > 0 in difftrap().")
648 elif numtraps == 1:
649 return 0.5*(function(interval[0])+function(interval[1]))
650 else:
651 numtosum = numtraps/2
652 h = float(interval[1]-interval[0])/numtosum
653 lox = interval[0] + 0.5 * h
654 points = lox + h * np.arange(numtosum)
655 s = np.sum(function(points), axis=0)
656 return s
659def _romberg_diff(b, c, k):
660 """
661 Compute the differences for the Romberg quadrature corrections.
662 See Forman Acton's "Real Computing Made Real," p 143.
663 """
664 tmp = 4.0**k
665 return (tmp * c - b)/(tmp - 1.0)
668def _printresmat(function, interval, resmat):
669 # Print the Romberg result matrix.
670 i = j = 0
671 print('Romberg integration of', repr(function), end=' ')
672 print('from', interval)
673 print('')
674 print('%6s %9s %9s' % ('Steps', 'StepSize', 'Results'))
675 for i in range(len(resmat)):
676 print('%6d %9f' % (2**i, (interval[1]-interval[0])/(2.**i)), end=' ')
677 for j in range(i+1):
678 print('%9f' % (resmat[i][j]), end=' ')
679 print('')
680 print('')
681 print('The final result is', resmat[i][j], end=' ')
682 print('after', 2**(len(resmat)-1)+1, 'function evaluations.')
685def romberg(function, a, b, args=(), tol=1.48e-8, rtol=1.48e-8, show=False,
686 divmax=10, vec_func=False):
687 """
688 Romberg integration of a callable function or method.
690 Returns the integral of `function` (a function of one variable)
691 over the interval (`a`, `b`).
693 If `show` is 1, the triangular array of the intermediate results
694 will be printed. If `vec_func` is True (default is False), then
695 `function` is assumed to support vector arguments.
697 Parameters
698 ----------
699 function : callable
700 Function to be integrated.
701 a : float
702 Lower limit of integration.
703 b : float
704 Upper limit of integration.
706 Returns
707 -------
708 results : float
709 Result of the integration.
711 Other Parameters
712 ----------------
713 args : tuple, optional
714 Extra arguments to pass to function. Each element of `args` will
715 be passed as a single argument to `func`. Default is to pass no
716 extra arguments.
717 tol, rtol : float, optional
718 The desired absolute and relative tolerances. Defaults are 1.48e-8.
719 show : bool, optional
720 Whether to print the results. Default is False.
721 divmax : int, optional
722 Maximum order of extrapolation. Default is 10.
723 vec_func : bool, optional
724 Whether `func` handles arrays as arguments (i.e., whether it is a
725 "vector" function). Default is False.
727 See Also
728 --------
729 fixed_quad : Fixed-order Gaussian quadrature.
730 quad : Adaptive quadrature using QUADPACK.
731 dblquad : Double integrals.
732 tplquad : Triple integrals.
733 romb : Integrators for sampled data.
734 simps : Integrators for sampled data.
735 cumtrapz : Cumulative integration for sampled data.
736 ode : ODE integrator.
737 odeint : ODE integrator.
739 References
740 ----------
741 .. [1] 'Romberg's method' https://en.wikipedia.org/wiki/Romberg%27s_method
743 Examples
744 --------
745 Integrate a gaussian from 0 to 1 and compare to the error function.
747 >>> from scipy import integrate
748 >>> from scipy.special import erf
749 >>> gaussian = lambda x: 1/np.sqrt(np.pi) * np.exp(-x**2)
750 >>> result = integrate.romberg(gaussian, 0, 1, show=True)
751 Romberg integration of <function vfunc at ...> from [0, 1]
753 ::
755 Steps StepSize Results
756 1 1.000000 0.385872
757 2 0.500000 0.412631 0.421551
758 4 0.250000 0.419184 0.421368 0.421356
759 8 0.125000 0.420810 0.421352 0.421350 0.421350
760 16 0.062500 0.421215 0.421350 0.421350 0.421350 0.421350
761 32 0.031250 0.421317 0.421350 0.421350 0.421350 0.421350 0.421350
763 The final result is 0.421350396475 after 33 function evaluations.
765 >>> print("%g %g" % (2*result, erf(1)))
766 0.842701 0.842701
768 """
769 if np.isinf(a) or np.isinf(b):
770 raise ValueError("Romberg integration only available "
771 "for finite limits.")
772 vfunc = vectorize1(function, args, vec_func=vec_func)
773 n = 1
774 interval = [a, b]
775 intrange = b - a
776 ordsum = _difftrap(vfunc, interval, n)
777 result = intrange * ordsum
778 resmat = [[result]]
779 err = np.inf
780 last_row = resmat[0]
781 for i in range(1, divmax+1):
782 n *= 2
783 ordsum += _difftrap(vfunc, interval, n)
784 row = [intrange * ordsum / n]
785 for k in range(i):
786 row.append(_romberg_diff(last_row[k], row[k], k+1))
787 result = row[i]
788 lastresult = last_row[i-1]
789 if show:
790 resmat.append(row)
791 err = abs(result - lastresult)
792 if err < tol or err < rtol * abs(result):
793 break
794 last_row = row
795 else:
796 warnings.warn(
797 "divmax (%d) exceeded. Latest difference = %e" % (divmax, err),
798 AccuracyWarning)
800 if show:
801 _printresmat(vfunc, interval, resmat)
802 return result
805# Coefficients for Newton-Cotes quadrature
806#
807# These are the points being used
808# to construct the local interpolating polynomial
809# a are the weights for Newton-Cotes integration
810# B is the error coefficient.
811# error in these coefficients grows as N gets larger.
812# or as samples are closer and closer together
814# You can use maxima to find these rational coefficients
815# for equally spaced data using the commands
816# 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);
817# Be(N) := N^(N+2)/(N+2)! * (N/(N+3) - sum((i/N)^(N+2)*a(i,N),i,0,N));
818# Bo(N) := N^(N+1)/(N+1)! * (N/(N+2) - sum((i/N)^(N+1)*a(i,N),i,0,N));
819# B(N) := (if (mod(N,2)=0) then Be(N) else Bo(N));
820#
821# pre-computed for equally-spaced weights
822#
823# num_a, den_a, int_a, num_B, den_B = _builtincoeffs[N]
824#
825# a = num_a*array(int_a)/den_a
826# B = num_B*1.0 / den_B
827#
828# integrate(f(x),x,x_0,x_N) = dx*sum(a*f(x_i)) + B*(dx)^(2k+3) f^(2k+2)(x*)
829# where k = N // 2
830#
831_builtincoeffs = {
832 1: (1,2,[1,1],-1,12),
833 2: (1,3,[1,4,1],-1,90),
834 3: (3,8,[1,3,3,1],-3,80),
835 4: (2,45,[7,32,12,32,7],-8,945),
836 5: (5,288,[19,75,50,50,75,19],-275,12096),
837 6: (1,140,[41,216,27,272,27,216,41],-9,1400),
838 7: (7,17280,[751,3577,1323,2989,2989,1323,3577,751],-8183,518400),
839 8: (4,14175,[989,5888,-928,10496,-4540,10496,-928,5888,989],
840 -2368,467775),
841 9: (9,89600,[2857,15741,1080,19344,5778,5778,19344,1080,
842 15741,2857], -4671, 394240),
843 10: (5,299376,[16067,106300,-48525,272400,-260550,427368,
844 -260550,272400,-48525,106300,16067],
845 -673175, 163459296),
846 11: (11,87091200,[2171465,13486539,-3237113, 25226685,-9595542,
847 15493566,15493566,-9595542,25226685,-3237113,
848 13486539,2171465], -2224234463, 237758976000),
849 12: (1, 5255250, [1364651,9903168,-7587864,35725120,-51491295,
850 87516288,-87797136,87516288,-51491295,35725120,
851 -7587864,9903168,1364651], -3012, 875875),
852 13: (13, 402361344000,[8181904909, 56280729661, -31268252574,
853 156074417954,-151659573325,206683437987,
854 -43111992612,-43111992612,206683437987,
855 -151659573325,156074417954,-31268252574,
856 56280729661,8181904909], -2639651053,
857 344881152000),
858 14: (7, 2501928000, [90241897,710986864,-770720657,3501442784,
859 -6625093363,12630121616,-16802270373,19534438464,
860 -16802270373,12630121616,-6625093363,3501442784,
861 -770720657,710986864,90241897], -3740727473,
862 1275983280000)
863 }
866def newton_cotes(rn, equal=0):
867 r"""
868 Return weights and error coefficient for Newton-Cotes integration.
870 Suppose we have (N+1) samples of f at the positions
871 x_0, x_1, ..., x_N. Then an N-point Newton-Cotes formula for the
872 integral between x_0 and x_N is:
874 :math:`\int_{x_0}^{x_N} f(x)dx = \Delta x \sum_{i=0}^{N} a_i f(x_i)
875 + B_N (\Delta x)^{N+2} f^{N+1} (\xi)`
877 where :math:`\xi \in [x_0,x_N]`
878 and :math:`\Delta x = \frac{x_N-x_0}{N}` is the average samples spacing.
880 If the samples are equally-spaced and N is even, then the error
881 term is :math:`B_N (\Delta x)^{N+3} f^{N+2}(\xi)`.
883 Parameters
884 ----------
885 rn : int
886 The integer order for equally-spaced data or the relative positions of
887 the samples with the first sample at 0 and the last at N, where N+1 is
888 the length of `rn`. N is the order of the Newton-Cotes integration.
889 equal : int, optional
890 Set to 1 to enforce equally spaced data.
892 Returns
893 -------
894 an : ndarray
895 1-D array of weights to apply to the function at the provided sample
896 positions.
897 B : float
898 Error coefficient.
900 Examples
901 --------
902 Compute the integral of sin(x) in [0, :math:`\pi`]:
904 >>> from scipy.integrate import newton_cotes
905 >>> def f(x):
906 ... return np.sin(x)
907 >>> a = 0
908 >>> b = np.pi
909 >>> exact = 2
910 >>> for N in [2, 4, 6, 8, 10]:
911 ... x = np.linspace(a, b, N + 1)
912 ... an, B = newton_cotes(N, 1)
913 ... dx = (b - a) / N
914 ... quad = dx * np.sum(an * f(x))
915 ... error = abs(quad - exact)
916 ... print('{:2d} {:10.9f} {:.5e}'.format(N, quad, error))
917 ...
918 2 2.094395102 9.43951e-02
919 4 1.998570732 1.42927e-03
920 6 2.000017814 1.78136e-05
921 8 1.999999835 1.64725e-07
922 10 2.000000001 1.14677e-09
924 Notes
925 -----
926 Normally, the Newton-Cotes rules are used on smaller integration
927 regions and a composite rule is used to return the total integral.
929 """
930 try:
931 N = len(rn)-1
932 if equal:
933 rn = np.arange(N+1)
934 elif np.all(np.diff(rn) == 1):
935 equal = 1
936 except Exception:
937 N = rn
938 rn = np.arange(N+1)
939 equal = 1
941 if equal and N in _builtincoeffs:
942 na, da, vi, nb, db = _builtincoeffs[N]
943 an = na * np.array(vi, dtype=float) / da
944 return an, float(nb)/db
946 if (rn[0] != 0) or (rn[-1] != N):
947 raise ValueError("The sample positions must start at 0"
948 " and end at N")
949 yi = rn / float(N)
950 ti = 2 * yi - 1
951 nvec = np.arange(N+1)
952 C = ti ** nvec[:, np.newaxis]
953 Cinv = np.linalg.inv(C)
954 # improve precision of result
955 for i in range(2):
956 Cinv = 2*Cinv - Cinv.dot(C).dot(Cinv)
957 vec = 2.0 / (nvec[::2]+1)
958 ai = Cinv[:, ::2].dot(vec) * (N / 2.)
960 if (N % 2 == 0) and equal:
961 BN = N/(N+3.)
962 power = N+2
963 else:
964 BN = N/(N+2.)
965 power = N+1
967 BN = BN - np.dot(yi**power, ai)
968 p1 = power+1
969 fac = power*math.log(N) - gammaln(p1)
970 fac = math.exp(fac)
971 return ai, BN*fac