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

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
1# Author: Travis Oliphant 2001
2# Author: Nathan Woods 2013 (nquad &c)
3import sys
4import warnings
5from functools import partial
7from . import _quadpack
8import numpy
9from numpy import Inf
11__all__ = ['quad', 'dblquad', 'tplquad', 'nquad', 'quad_explain',
12 'IntegrationWarning']
15error = _quadpack.error
17class IntegrationWarning(UserWarning):
18 """
19 Warning on issues during integration.
20 """
21 pass
24def quad_explain(output=sys.stdout):
25 """
26 Print extra information about integrate.quad() parameters and returns.
28 Parameters
29 ----------
30 output : instance with "write" method, optional
31 Information about `quad` is passed to ``output.write()``.
32 Default is ``sys.stdout``.
34 Returns
35 -------
36 None
38 Examples
39 --------
40 We can show detailed information of the `integrate.quad` function in stdout:
42 >>> from scipy.integrate import quad_explain
43 >>> quad_explain()
45 """
46 output.write(quad.__doc__)
49def quad(func, a, b, args=(), full_output=0, epsabs=1.49e-8, epsrel=1.49e-8,
50 limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50,
51 limlst=50):
52 """
53 Compute a definite integral.
55 Integrate func from `a` to `b` (possibly infinite interval) using a
56 technique from the Fortran library QUADPACK.
58 Parameters
59 ----------
60 func : {function, scipy.LowLevelCallable}
61 A Python function or method to integrate. If `func` takes many
62 arguments, it is integrated along the axis corresponding to the
63 first argument.
65 If the user desires improved integration performance, then `f` may
66 be a `scipy.LowLevelCallable` with one of the signatures::
68 double func(double x)
69 double func(double x, void *user_data)
70 double func(int n, double *xx)
71 double func(int n, double *xx, void *user_data)
73 The ``user_data`` is the data contained in the `scipy.LowLevelCallable`.
74 In the call forms with ``xx``, ``n`` is the length of the ``xx``
75 array which contains ``xx[0] == x`` and the rest of the items are
76 numbers contained in the ``args`` argument of quad.
78 In addition, certain ctypes call signatures are supported for
79 backward compatibility, but those should not be used in new code.
80 a : float
81 Lower limit of integration (use -numpy.inf for -infinity).
82 b : float
83 Upper limit of integration (use numpy.inf for +infinity).
84 args : tuple, optional
85 Extra arguments to pass to `func`.
86 full_output : int, optional
87 Non-zero to return a dictionary of integration information.
88 If non-zero, warning messages are also suppressed and the
89 message is appended to the output tuple.
91 Returns
92 -------
93 y : float
94 The integral of func from `a` to `b`.
95 abserr : float
96 An estimate of the absolute error in the result.
97 infodict : dict
98 A dictionary containing additional information.
99 Run scipy.integrate.quad_explain() for more information.
100 message
101 A convergence message.
102 explain
103 Appended only with 'cos' or 'sin' weighting and infinite
104 integration limits, it contains an explanation of the codes in
105 infodict['ierlst']
107 Other Parameters
108 ----------------
109 epsabs : float or int, optional
110 Absolute error tolerance. Default is 1.49e-8. `quad` tries to obtain
111 an accuracy of ``abs(i-result) <= max(epsabs, epsrel*abs(i))``
112 where ``i`` = integral of `func` from `a` to `b`, and ``result`` is the
113 numerical approximation. See `epsrel` below.
114 epsrel : float or int, optional
115 Relative error tolerance. Default is 1.49e-8.
116 If ``epsabs <= 0``, `epsrel` must be greater than both 5e-29
117 and ``50 * (machine epsilon)``. See `epsabs` above.
118 limit : float or int, optional
119 An upper bound on the number of subintervals used in the adaptive
120 algorithm.
121 points : (sequence of floats,ints), optional
122 A sequence of break points in the bounded integration interval
123 where local difficulties of the integrand may occur (e.g.,
124 singularities, discontinuities). The sequence does not have
125 to be sorted. Note that this option cannot be used in conjunction
126 with ``weight``.
127 weight : float or int, optional
128 String indicating weighting function. Full explanation for this
129 and the remaining arguments can be found below.
130 wvar : optional
131 Variables for use with weighting functions.
132 wopts : optional
133 Optional input for reusing Chebyshev moments.
134 maxp1 : float or int, optional
135 An upper bound on the number of Chebyshev moments.
136 limlst : int, optional
137 Upper bound on the number of cycles (>=3) for use with a sinusoidal
138 weighting and an infinite end-point.
140 See Also
141 --------
142 dblquad : double integral
143 tplquad : triple integral
144 nquad : n-dimensional integrals (uses `quad` recursively)
145 fixed_quad : fixed-order Gaussian quadrature
146 quadrature : adaptive Gaussian quadrature
147 odeint : ODE integrator
148 ode : ODE integrator
149 simps : integrator for sampled data
150 romb : integrator for sampled data
151 scipy.special : for coefficients and roots of orthogonal polynomials
153 Notes
154 -----
156 **Extra information for quad() inputs and outputs**
158 If full_output is non-zero, then the third output argument
159 (infodict) is a dictionary with entries as tabulated below. For
160 infinite limits, the range is transformed to (0,1) and the
161 optional outputs are given with respect to this transformed range.
162 Let M be the input argument limit and let K be infodict['last'].
163 The entries are:
165 'neval'
166 The number of function evaluations.
167 'last'
168 The number, K, of subintervals produced in the subdivision process.
169 'alist'
170 A rank-1 array of length M, the first K elements of which are the
171 left end points of the subintervals in the partition of the
172 integration range.
173 'blist'
174 A rank-1 array of length M, the first K elements of which are the
175 right end points of the subintervals.
176 'rlist'
177 A rank-1 array of length M, the first K elements of which are the
178 integral approximations on the subintervals.
179 'elist'
180 A rank-1 array of length M, the first K elements of which are the
181 moduli of the absolute error estimates on the subintervals.
182 'iord'
183 A rank-1 integer array of length M, the first L elements of
184 which are pointers to the error estimates over the subintervals
185 with ``L=K`` if ``K<=M/2+2`` or ``L=M+1-K`` otherwise. Let I be the
186 sequence ``infodict['iord']`` and let E be the sequence
187 ``infodict['elist']``. Then ``E[I[1]], ..., E[I[L]]`` forms a
188 decreasing sequence.
190 If the input argument points is provided (i.e., it is not None),
191 the following additional outputs are placed in the output
192 dictionary. Assume the points sequence is of length P.
194 'pts'
195 A rank-1 array of length P+2 containing the integration limits
196 and the break points of the intervals in ascending order.
197 This is an array giving the subintervals over which integration
198 will occur.
199 'level'
200 A rank-1 integer array of length M (=limit), containing the
201 subdivision levels of the subintervals, i.e., if (aa,bb) is a
202 subinterval of ``(pts[1], pts[2])`` where ``pts[0]`` and ``pts[2]``
203 are adjacent elements of ``infodict['pts']``, then (aa,bb) has level l
204 if ``|bb-aa| = |pts[2]-pts[1]| * 2**(-l)``.
205 'ndin'
206 A rank-1 integer array of length P+2. After the first integration
207 over the intervals (pts[1], pts[2]), the error estimates over some
208 of the intervals may have been increased artificially in order to
209 put their subdivision forward. This array has ones in slots
210 corresponding to the subintervals for which this happens.
212 **Weighting the integrand**
214 The input variables, *weight* and *wvar*, are used to weight the
215 integrand by a select list of functions. Different integration
216 methods are used to compute the integral with these weighting
217 functions, and these do not support specifying break points. The
218 possible values of weight and the corresponding weighting functions are.
220 ========== =================================== =====================
221 ``weight`` Weight function used ``wvar``
222 ========== =================================== =====================
223 'cos' cos(w*x) wvar = w
224 'sin' sin(w*x) wvar = w
225 'alg' g(x) = ((x-a)**alpha)*((b-x)**beta) wvar = (alpha, beta)
226 'alg-loga' g(x)*log(x-a) wvar = (alpha, beta)
227 'alg-logb' g(x)*log(b-x) wvar = (alpha, beta)
228 'alg-log' g(x)*log(x-a)*log(b-x) wvar = (alpha, beta)
229 'cauchy' 1/(x-c) wvar = c
230 ========== =================================== =====================
232 wvar holds the parameter w, (alpha, beta), or c depending on the weight
233 selected. In these expressions, a and b are the integration limits.
235 For the 'cos' and 'sin' weighting, additional inputs and outputs are
236 available.
238 For finite integration limits, the integration is performed using a
239 Clenshaw-Curtis method which uses Chebyshev moments. For repeated
240 calculations, these moments are saved in the output dictionary:
242 'momcom'
243 The maximum level of Chebyshev moments that have been computed,
244 i.e., if ``M_c`` is ``infodict['momcom']`` then the moments have been
245 computed for intervals of length ``|b-a| * 2**(-l)``,
246 ``l=0,1,...,M_c``.
247 'nnlog'
248 A rank-1 integer array of length M(=limit), containing the
249 subdivision levels of the subintervals, i.e., an element of this
250 array is equal to l if the corresponding subinterval is
251 ``|b-a|* 2**(-l)``.
252 'chebmo'
253 A rank-2 array of shape (25, maxp1) containing the computed
254 Chebyshev moments. These can be passed on to an integration
255 over the same interval by passing this array as the second
256 element of the sequence wopts and passing infodict['momcom'] as
257 the first element.
259 If one of the integration limits is infinite, then a Fourier integral is
260 computed (assuming w neq 0). If full_output is 1 and a numerical error
261 is encountered, besides the error message attached to the output tuple,
262 a dictionary is also appended to the output tuple which translates the
263 error codes in the array ``info['ierlst']`` to English messages. The
264 output information dictionary contains the following entries instead of
265 'last', 'alist', 'blist', 'rlist', and 'elist':
267 'lst'
268 The number of subintervals needed for the integration (call it ``K_f``).
269 'rslst'
270 A rank-1 array of length M_f=limlst, whose first ``K_f`` elements
271 contain the integral contribution over the interval
272 ``(a+(k-1)c, a+kc)`` where ``c = (2*floor(|w|) + 1) * pi / |w|``
273 and ``k=1,2,...,K_f``.
274 'erlst'
275 A rank-1 array of length ``M_f`` containing the error estimate
276 corresponding to the interval in the same position in
277 ``infodict['rslist']``.
278 'ierlst'
279 A rank-1 integer array of length ``M_f`` containing an error flag
280 corresponding to the interval in the same position in
281 ``infodict['rslist']``. See the explanation dictionary (last entry
282 in the output tuple) for the meaning of the codes.
284 Examples
285 --------
286 Calculate :math:`\\int^4_0 x^2 dx` and compare with an analytic result
288 >>> from scipy import integrate
289 >>> x2 = lambda x: x**2
290 >>> integrate.quad(x2, 0, 4)
291 (21.333333333333332, 2.3684757858670003e-13)
292 >>> print(4**3 / 3.) # analytical result
293 21.3333333333
295 Calculate :math:`\\int^\\infty_0 e^{-x} dx`
297 >>> invexp = lambda x: np.exp(-x)
298 >>> integrate.quad(invexp, 0, np.inf)
299 (1.0, 5.842605999138044e-11)
301 >>> f = lambda x,a : a*x
302 >>> y, err = integrate.quad(f, 0, 1, args=(1,))
303 >>> y
304 0.5
305 >>> y, err = integrate.quad(f, 0, 1, args=(3,))
306 >>> y
307 1.5
309 Calculate :math:`\\int^1_0 x^2 + y^2 dx` with ctypes, holding
310 y parameter as 1::
312 testlib.c =>
313 double func(int n, double args[n]){
314 return args[0]*args[0] + args[1]*args[1];}
315 compile to library testlib.*
317 ::
319 from scipy import integrate
320 import ctypes
321 lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path
322 lib.func.restype = ctypes.c_double
323 lib.func.argtypes = (ctypes.c_int,ctypes.c_double)
324 integrate.quad(lib.func,0,1,(1))
325 #(1.3333333333333333, 1.4802973661668752e-14)
326 print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result
327 # 1.3333333333333333
329 Be aware that pulse shapes and other sharp features as compared to the
330 size of the integration interval may not be integrated correctly using
331 this method. A simplified example of this limitation is integrating a
332 y-axis reflected step function with many zero values within the integrals
333 bounds.
335 >>> y = lambda x: 1 if x<=0 else 0
336 >>> integrate.quad(y, -1, 1)
337 (1.0, 1.1102230246251565e-14)
338 >>> integrate.quad(y, -1, 100)
339 (1.0000000002199108, 1.0189464580163188e-08)
340 >>> integrate.quad(y, -1, 10000)
341 (0.0, 0.0)
343 """
344 if not isinstance(args, tuple):
345 args = (args,)
347 # check the limits of integration: \int_a^b, expect a < b
348 flip, a, b = b < a, min(a, b), max(a, b)
350 if weight is None:
351 retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
352 points)
353 else:
354 if points is not None:
355 msg = ("Break points cannot be specified when using weighted integrand.\n"
356 "Continuing, ignoring specified points.")
357 warnings.warn(msg, IntegrationWarning, stacklevel=2)
358 retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,
359 limlst, limit, maxp1, weight, wvar, wopts)
361 if flip:
362 retval = (-retval[0],) + retval[1:]
364 ier = retval[-1]
365 if ier == 0:
366 return retval[:-1]
368 msgs = {80: "A Python error occurred possibly while calling the function.",
369 1: "The maximum number of subdivisions (%d) has been achieved.\n If increasing the limit yields no improvement it is advised to analyze \n the integrand in order to determine the difficulties. If the position of a \n local difficulty can be determined (singularity, discontinuity) one will \n probably gain from splitting up the interval and calling the integrator \n on the subranges. Perhaps a special-purpose integrator should be used." % limit,
370 2: "The occurrence of roundoff error is detected, which prevents \n the requested tolerance from being achieved. The error may be \n underestimated.",
371 3: "Extremely bad integrand behavior occurs at some points of the\n integration interval.",
372 4: "The algorithm does not converge. Roundoff error is detected\n in the extrapolation table. It is assumed that the requested tolerance\n cannot be achieved, and that the returned result (if full_output = 1) is \n the best which can be obtained.",
373 5: "The integral is probably divergent, or slowly convergent.",
374 6: "The input is invalid.",
375 7: "Abnormal termination of the routine. The estimates for result\n and error are less reliable. It is assumed that the requested accuracy\n has not been achieved.",
376 'unknown': "Unknown error."}
378 if weight in ['cos','sin'] and (b == Inf or a == -Inf):
379 msgs[1] = "The maximum number of cycles allowed has been achieved., e.e.\n of subintervals (a+(k-1)c, a+kc) where c = (2*int(abs(omega)+1))\n *pi/abs(omega), for k = 1, 2, ..., lst. One can allow more cycles by increasing the value of limlst. Look at info['ierlst'] with full_output=1."
380 msgs[4] = "The extrapolation table constructed for convergence acceleration\n of the series formed by the integral contributions over the cycles, \n does not converge to within the requested accuracy. Look at \n info['ierlst'] with full_output=1."
381 msgs[7] = "Bad integrand behavior occurs within one or more of the cycles.\n Location and type of the difficulty involved can be determined from \n the vector info['ierlist'] obtained with full_output=1."
382 explain = {1: "The maximum number of subdivisions (= limit) has been \n achieved on this cycle.",
383 2: "The occurrence of roundoff error is detected and prevents\n the tolerance imposed on this cycle from being achieved.",
384 3: "Extremely bad integrand behavior occurs at some points of\n this cycle.",
385 4: "The integral over this cycle does not converge (to within the required accuracy) due to roundoff in the extrapolation procedure invoked on this cycle. It is assumed that the result on this interval is the best which can be obtained.",
386 5: "The integral over this cycle is probably divergent or slowly convergent."}
388 try:
389 msg = msgs[ier]
390 except KeyError:
391 msg = msgs['unknown']
393 if ier in [1,2,3,4,5,7]:
394 if full_output:
395 if weight in ['cos', 'sin'] and (b == Inf or a == -Inf):
396 return retval[:-1] + (msg, explain)
397 else:
398 return retval[:-1] + (msg,)
399 else:
400 warnings.warn(msg, IntegrationWarning, stacklevel=2)
401 return retval[:-1]
403 elif ier == 6: # Forensic decision tree when QUADPACK throws ier=6
404 if epsabs <= 0: # Small error tolerance - applies to all methods
405 if epsrel < max(50 * sys.float_info.epsilon, 5e-29):
406 msg = ("If 'epsabs'<=0, 'epsrel' must be greater than both"
407 " 5e-29 and 50*(machine epsilon).")
408 elif weight in ['sin', 'cos'] and (abs(a) + abs(b) == Inf):
409 msg = ("Sine or cosine weighted intergals with infinite domain"
410 " must have 'epsabs'>0.")
412 elif weight is None:
413 if points is None: # QAGSE/QAGIE
414 msg = ("Invalid 'limit' argument. There must be"
415 " at least one subinterval")
416 else: # QAGPE
417 if not (min(a, b) <= min(points) <= max(points) <= max(a, b)):
418 msg = ("All break points in 'points' must lie within the"
419 " integration limits.")
420 elif len(points) >= limit:
421 msg = ("Number of break points ({:d})"
422 " must be less than subinterval"
423 " limit ({:d})").format(len(points), limit)
425 else:
426 if maxp1 < 1:
427 msg = "Chebyshev moment limit maxp1 must be >=1."
429 elif weight in ('cos', 'sin') and abs(a+b) == Inf: # QAWFE
430 msg = "Cycle limit limlst must be >=3."
432 elif weight.startswith('alg'): # QAWSE
433 if min(wvar) < -1:
434 msg = "wvar parameters (alpha, beta) must both be >= -1."
435 if b < a:
436 msg = "Integration limits a, b must satistfy a<b."
438 elif weight == 'cauchy' and wvar in (a, b):
439 msg = ("Parameter 'wvar' must not equal"
440 " integration limits 'a' or 'b'.")
442 raise ValueError(msg)
445def _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points):
446 infbounds = 0
447 if (b != Inf and a != -Inf):
448 pass # standard integration
449 elif (b == Inf and a != -Inf):
450 infbounds = 1
451 bound = a
452 elif (b == Inf and a == -Inf):
453 infbounds = 2
454 bound = 0 # ignored
455 elif (b != Inf and a == -Inf):
456 infbounds = -1
457 bound = b
458 else:
459 raise RuntimeError("Infinity comparisons don't work for you.")
461 if points is None:
462 if infbounds == 0:
463 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
464 else:
465 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
466 else:
467 if infbounds != 0:
468 raise ValueError("Infinity inputs cannot be used with break points.")
469 else:
470 #Duplicates force function evaluation at singular points
471 the_points = numpy.unique(points)
472 the_points = the_points[a < the_points]
473 the_points = the_points[the_points < b]
474 the_points = numpy.concatenate((the_points, (0., 0.)))
475 return _quadpack._qagpe(func,a,b,the_points,args,full_output,epsabs,epsrel,limit)
478def _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts):
479 if weight not in ['cos','sin','alg','alg-loga','alg-logb','alg-log','cauchy']:
480 raise ValueError("%s not a recognized weighting function." % weight)
482 strdict = {'cos':1,'sin':2,'alg':1,'alg-loga':2,'alg-logb':3,'alg-log':4}
484 if weight in ['cos','sin']:
485 integr = strdict[weight]
486 if (b != Inf and a != -Inf): # finite limits
487 if wopts is None: # no precomputed Chebyshev moments
488 return _quadpack._qawoe(func, a, b, wvar, integr, args, full_output,
489 epsabs, epsrel, limit, maxp1,1)
490 else: # precomputed Chebyshev moments
491 momcom = wopts[0]
492 chebcom = wopts[1]
493 return _quadpack._qawoe(func, a, b, wvar, integr, args, full_output,
494 epsabs, epsrel, limit, maxp1, 2, momcom, chebcom)
496 elif (b == Inf and a != -Inf):
497 return _quadpack._qawfe(func, a, wvar, integr, args, full_output,
498 epsabs,limlst,limit,maxp1)
499 elif (b != Inf and a == -Inf): # remap function and interval
500 if weight == 'cos':
501 def thefunc(x,*myargs):
502 y = -x
503 func = myargs[0]
504 myargs = (y,) + myargs[1:]
505 return func(*myargs)
506 else:
507 def thefunc(x,*myargs):
508 y = -x
509 func = myargs[0]
510 myargs = (y,) + myargs[1:]
511 return -func(*myargs)
512 args = (func,) + args
513 return _quadpack._qawfe(thefunc, -b, wvar, integr, args,
514 full_output, epsabs, limlst, limit, maxp1)
515 else:
516 raise ValueError("Cannot integrate with this weight from -Inf to +Inf.")
517 else:
518 if a in [-Inf,Inf] or b in [-Inf,Inf]:
519 raise ValueError("Cannot integrate with this weight over an infinite interval.")
521 if weight.startswith('alg'):
522 integr = strdict[weight]
523 return _quadpack._qawse(func, a, b, wvar, integr, args,
524 full_output, epsabs, epsrel, limit)
525 else: # weight == 'cauchy'
526 return _quadpack._qawce(func, a, b, wvar, args, full_output,
527 epsabs, epsrel, limit)
530def dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-8, epsrel=1.49e-8):
531 """
532 Compute a double integral.
534 Return the double (definite) integral of ``func(y, x)`` from ``x = a..b``
535 and ``y = gfun(x)..hfun(x)``.
537 Parameters
538 ----------
539 func : callable
540 A Python function or method of at least two variables: y must be the
541 first argument and x the second argument.
542 a, b : float
543 The limits of integration in x: `a` < `b`
544 gfun : callable or float
545 The lower boundary curve in y which is a function taking a single
546 floating point argument (x) and returning a floating point result
547 or a float indicating a constant boundary curve.
548 hfun : callable or float
549 The upper boundary curve in y (same requirements as `gfun`).
550 args : sequence, optional
551 Extra arguments to pass to `func`.
552 epsabs : float, optional
553 Absolute tolerance passed directly to the inner 1-D quadrature
554 integration. Default is 1.49e-8. `dblquad`` tries to obtain
555 an accuracy of ``abs(i-result) <= max(epsabs, epsrel*abs(i))``
556 where ``i`` = inner integral of ``func(y, x)`` from ``gfun(x)``
557 to ``hfun(x)``, and ``result`` is the numerical approximation.
558 See `epsrel` below.
559 epsrel : float, optional
560 Relative tolerance of the inner 1-D integrals. Default is 1.49e-8.
561 If ``epsabs <= 0``, `epsrel` must be greater than both 5e-29
562 and ``50 * (machine epsilon)``. See `epsabs` above.
564 Returns
565 -------
566 y : float
567 The resultant integral.
568 abserr : float
569 An estimate of the error.
571 See also
572 --------
573 quad : single integral
574 tplquad : triple integral
575 nquad : N-dimensional integrals
576 fixed_quad : fixed-order Gaussian quadrature
577 quadrature : adaptive Gaussian quadrature
578 odeint : ODE integrator
579 ode : ODE integrator
580 simps : integrator for sampled data
581 romb : integrator for sampled data
582 scipy.special : for coefficients and roots of orthogonal polynomials
584 Examples
585 --------
587 Compute the double integral of ``x * y**2`` over the box
588 ``x`` ranging from 0 to 2 and ``y`` ranging from 0 to 1.
590 >>> from scipy import integrate
591 >>> f = lambda y, x: x*y**2
592 >>> integrate.dblquad(f, 0, 2, lambda x: 0, lambda x: 1)
593 (0.6666666666666667, 7.401486830834377e-15)
595 """
597 def temp_ranges(*args):
598 return [gfun(args[0]) if callable(gfun) else gfun,
599 hfun(args[0]) if callable(hfun) else hfun]
601 return nquad(func, [temp_ranges, [a, b]], args=args,
602 opts={"epsabs": epsabs, "epsrel": epsrel})
605def tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-8,
606 epsrel=1.49e-8):
607 """
608 Compute a triple (definite) integral.
610 Return the triple integral of ``func(z, y, x)`` from ``x = a..b``,
611 ``y = gfun(x)..hfun(x)``, and ``z = qfun(x,y)..rfun(x,y)``.
613 Parameters
614 ----------
615 func : function
616 A Python function or method of at least three variables in the
617 order (z, y, x).
618 a, b : float
619 The limits of integration in x: `a` < `b`
620 gfun : function or float
621 The lower boundary curve in y which is a function taking a single
622 floating point argument (x) and returning a floating point result
623 or a float indicating a constant boundary curve.
624 hfun : function or float
625 The upper boundary curve in y (same requirements as `gfun`).
626 qfun : function or float
627 The lower boundary surface in z. It must be a function that takes
628 two floats in the order (x, y) and returns a float or a float
629 indicating a constant boundary surface.
630 rfun : function or float
631 The upper boundary surface in z. (Same requirements as `qfun`.)
632 args : tuple, optional
633 Extra arguments to pass to `func`.
634 epsabs : float, optional
635 Absolute tolerance passed directly to the innermost 1-D quadrature
636 integration. Default is 1.49e-8.
637 epsrel : float, optional
638 Relative tolerance of the innermost 1-D integrals. Default is 1.49e-8.
640 Returns
641 -------
642 y : float
643 The resultant integral.
644 abserr : float
645 An estimate of the error.
647 See Also
648 --------
649 quad: Adaptive quadrature using QUADPACK
650 quadrature: Adaptive Gaussian quadrature
651 fixed_quad: Fixed-order Gaussian quadrature
652 dblquad: Double integrals
653 nquad : N-dimensional integrals
654 romb: Integrators for sampled data
655 simps: Integrators for sampled data
656 ode: ODE integrators
657 odeint: ODE integrators
658 scipy.special: For coefficients and roots of orthogonal polynomials
660 Examples
661 --------
663 Compute the triple integral of ``x * y * z``, over ``x`` ranging
664 from 1 to 2, ``y`` ranging from 2 to 3, ``z`` ranging from 0 to 1.
666 >>> from scipy import integrate
667 >>> f = lambda z, y, x: x*y*z
668 >>> integrate.tplquad(f, 1, 2, lambda x: 2, lambda x: 3,
669 ... lambda x, y: 0, lambda x, y: 1)
670 (1.8750000000000002, 3.324644794257407e-14)
673 """
674 # f(z, y, x)
675 # qfun/rfun (x, y)
676 # gfun/hfun(x)
677 # nquad will hand (y, x, t0, ...) to ranges0
678 # nquad will hand (x, t0, ...) to ranges1
679 # Stupid different API...
681 def ranges0(*args):
682 return [qfun(args[1], args[0]) if callable(qfun) else qfun,
683 rfun(args[1], args[0]) if callable(rfun) else rfun]
685 def ranges1(*args):
686 return [gfun(args[0]) if callable(gfun) else gfun,
687 hfun(args[0]) if callable(hfun) else hfun]
689 ranges = [ranges0, ranges1, [a, b]]
690 return nquad(func, ranges, args=args,
691 opts={"epsabs": epsabs, "epsrel": epsrel})
694def nquad(func, ranges, args=None, opts=None, full_output=False):
695 """
696 Integration over multiple variables.
698 Wraps `quad` to enable integration over multiple variables.
699 Various options allow improved integration of discontinuous functions, as
700 well as the use of weighted integration, and generally finer control of the
701 integration process.
703 Parameters
704 ----------
705 func : {callable, scipy.LowLevelCallable}
706 The function to be integrated. Has arguments of ``x0, ... xn``,
707 ``t0, tm``, where integration is carried out over ``x0, ... xn``, which
708 must be floats. Function signature should be
709 ``func(x0, x1, ..., xn, t0, t1, ..., tm)``. Integration is carried out
710 in order. That is, integration over ``x0`` is the innermost integral,
711 and ``xn`` is the outermost.
713 If the user desires improved integration performance, then `f` may
714 be a `scipy.LowLevelCallable` with one of the signatures::
716 double func(int n, double *xx)
717 double func(int n, double *xx, void *user_data)
719 where ``n`` is the number of extra parameters and args is an array
720 of doubles of the additional parameters, the ``xx`` array contains the
721 coordinates. The ``user_data`` is the data contained in the
722 `scipy.LowLevelCallable`.
723 ranges : iterable object
724 Each element of ranges may be either a sequence of 2 numbers, or else
725 a callable that returns such a sequence. ``ranges[0]`` corresponds to
726 integration over x0, and so on. If an element of ranges is a callable,
727 then it will be called with all of the integration arguments available,
728 as well as any parametric arguments. e.g., if
729 ``func = f(x0, x1, x2, t0, t1)``, then ``ranges[0]`` may be defined as
730 either ``(a, b)`` or else as ``(a, b) = range0(x1, x2, t0, t1)``.
731 args : iterable object, optional
732 Additional arguments ``t0, ..., tn``, required by `func`, `ranges`, and
733 ``opts``.
734 opts : iterable object or dict, optional
735 Options to be passed to `quad`. May be empty, a dict, or
736 a sequence of dicts or functions that return a dict. If empty, the
737 default options from scipy.integrate.quad are used. If a dict, the same
738 options are used for all levels of integraion. If a sequence, then each
739 element of the sequence corresponds to a particular integration. e.g.,
740 opts[0] corresponds to integration over x0, and so on. If a callable,
741 the signature must be the same as for ``ranges``. The available
742 options together with their default values are:
744 - epsabs = 1.49e-08
745 - epsrel = 1.49e-08
746 - limit = 50
747 - points = None
748 - weight = None
749 - wvar = None
750 - wopts = None
752 For more information on these options, see `quad` and `quad_explain`.
754 full_output : bool, optional
755 Partial implementation of ``full_output`` from scipy.integrate.quad.
756 The number of integrand function evaluations ``neval`` can be obtained
757 by setting ``full_output=True`` when calling nquad.
759 Returns
760 -------
761 result : float
762 The result of the integration.
763 abserr : float
764 The maximum of the estimates of the absolute error in the various
765 integration results.
766 out_dict : dict, optional
767 A dict containing additional information on the integration.
769 See Also
770 --------
771 quad : 1-D numerical integration
772 dblquad, tplquad : double and triple integrals
773 fixed_quad : fixed-order Gaussian quadrature
774 quadrature : adaptive Gaussian quadrature
776 Examples
777 --------
778 >>> from scipy import integrate
779 >>> func = lambda x0,x1,x2,x3 : x0**2 + x1*x2 - x3**3 + np.sin(x0) + (
780 ... 1 if (x0-.2*x3-.5-.25*x1>0) else 0)
781 >>> points = [[lambda x1,x2,x3 : 0.2*x3 + 0.5 + 0.25*x1], [], [], []]
782 >>> def opts0(*args, **kwargs):
783 ... return {'points':[0.2*args[2] + 0.5 + 0.25*args[0]]}
784 >>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]],
785 ... opts=[opts0,{},{},{}], full_output=True)
786 (1.5267454070738633, 2.9437360001402324e-14, {'neval': 388962})
788 >>> scale = .1
789 >>> def func2(x0, x1, x2, x3, t0, t1):
790 ... return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0+t1*x1-t0>0 else 0)
791 >>> def lim0(x1, x2, x3, t0, t1):
792 ... return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
793 ... scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]
794 >>> def lim1(x2, x3, t0, t1):
795 ... return [scale * (t0*x2 + t1*x3) - 1,
796 ... scale * (t0*x2 + t1*x3) + 1]
797 >>> def lim2(x3, t0, t1):
798 ... return [scale * (x3 + t0**2*t1**3) - 1,
799 ... scale * (x3 + t0**2*t1**3) + 1]
800 >>> def lim3(t0, t1):
801 ... return [scale * (t0+t1) - 1, scale * (t0+t1) + 1]
802 >>> def opts0(x1, x2, x3, t0, t1):
803 ... return {'points' : [t0 - t1*x1]}
804 >>> def opts1(x2, x3, t0, t1):
805 ... return {}
806 >>> def opts2(x3, t0, t1):
807 ... return {}
808 >>> def opts3(t0, t1):
809 ... return {}
810 >>> integrate.nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0),
811 ... opts=[opts0, opts1, opts2, opts3])
812 (25.066666666666666, 2.7829590483937256e-13)
814 """
815 depth = len(ranges)
816 ranges = [rng if callable(rng) else _RangeFunc(rng) for rng in ranges]
817 if args is None:
818 args = ()
819 if opts is None:
820 opts = [dict([])] * depth
822 if isinstance(opts, dict):
823 opts = [_OptFunc(opts)] * depth
824 else:
825 opts = [opt if callable(opt) else _OptFunc(opt) for opt in opts]
826 return _NQuad(func, ranges, opts, full_output).integrate(*args)
829class _RangeFunc(object):
830 def __init__(self, range_):
831 self.range_ = range_
833 def __call__(self, *args):
834 """Return stored value.
836 *args needed because range_ can be float or func, and is called with
837 variable number of parameters.
838 """
839 return self.range_
842class _OptFunc(object):
843 def __init__(self, opt):
844 self.opt = opt
846 def __call__(self, *args):
847 """Return stored dict."""
848 return self.opt
851class _NQuad(object):
852 def __init__(self, func, ranges, opts, full_output):
853 self.abserr = 0
854 self.func = func
855 self.ranges = ranges
856 self.opts = opts
857 self.maxdepth = len(ranges)
858 self.full_output = full_output
859 if self.full_output:
860 self.out_dict = {'neval': 0}
862 def integrate(self, *args, **kwargs):
863 depth = kwargs.pop('depth', 0)
864 if kwargs:
865 raise ValueError('unexpected kwargs')
867 # Get the integration range and options for this depth.
868 ind = -(depth + 1)
869 fn_range = self.ranges[ind]
870 low, high = fn_range(*args)
871 fn_opt = self.opts[ind]
872 opt = dict(fn_opt(*args))
874 if 'points' in opt:
875 opt['points'] = [x for x in opt['points'] if low <= x <= high]
876 if depth + 1 == self.maxdepth:
877 f = self.func
878 else:
879 f = partial(self.integrate, depth=depth+1)
880 quad_r = quad(f, low, high, args=args, full_output=self.full_output,
881 **opt)
882 value = quad_r[0]
883 abserr = quad_r[1]
884 if self.full_output:
885 infodict = quad_r[2]
886 # The 'neval' parameter in full_output returns the total
887 # number of times the integrand function was evaluated.
888 # Therefore, only the innermost integration loop counts.
889 if depth + 1 == self.maxdepth:
890 self.out_dict['neval'] += infodict['neval']
891 self.abserr = max(self.abserr, abserr)
892 if depth > 0:
893 return value
894 else:
895 # Final result of N-D integration with error
896 if self.full_output:
897 return value, self.abserr, self.out_dict
898 else:
899 return value, self.abserr