Coverage for /usr/lib/python3/dist-packages/scipy/integrate/_quadpack_py.py: 12%

1# Author: Travis Oliphant 2001 

2# Author: Nathan Woods 2013 (nquad &c) 

3import sys 

4import warnings 

5from functools import partial 

6 

7from . import _quadpack 

8import numpy as np 

9from numpy import Inf 

10 

11__all__ = ["quad", "dblquad", "tplquad", "nquad", "IntegrationWarning"] 

12 

13 

14error = _quadpack.error 

15 

16class IntegrationWarning(UserWarning): 

17 """ 

18 Warning on issues during integration. 

19 """ 

20 pass 

21 

22 

23def quad(func, a, b, args=(), full_output=0, epsabs=1.49e-8, epsrel=1.49e-8, 

24 limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, 

25 limlst=50, complex_func=False): 

26 """ 

27 Compute a definite integral. 

28 

29 Integrate func from `a` to `b` (possibly infinite interval) using a 

30 technique from the Fortran library QUADPACK. 

31 

32 Parameters 

33 ---------- 

34 func : {function, scipy.LowLevelCallable} 

35 A Python function or method to integrate. If `func` takes many 

36 arguments, it is integrated along the axis corresponding to the 

37 first argument. 

38 

39 If the user desires improved integration performance, then `f` may 

40 be a `scipy.LowLevelCallable` with one of the signatures:: 

41 

42 double func(double x) 

43 double func(double x, void *user_data) 

44 double func(int n, double *xx) 

45 double func(int n, double *xx, void *user_data) 

46 

47 The ``user_data`` is the data contained in the `scipy.LowLevelCallable`. 

48 In the call forms with ``xx``, ``n`` is the length of the ``xx`` 

49 array which contains ``xx[0] == x`` and the rest of the items are 

50 numbers contained in the ``args`` argument of quad. 

51 

52 In addition, certain ctypes call signatures are supported for 

53 backward compatibility, but those should not be used in new code. 

54 a : float 

55 Lower limit of integration (use -numpy.inf for -infinity). 

56 b : float 

57 Upper limit of integration (use numpy.inf for +infinity). 

58 args : tuple, optional 

59 Extra arguments to pass to `func`. 

60 full_output : int, optional 

61 Non-zero to return a dictionary of integration information. 

62 If non-zero, warning messages are also suppressed and the 

63 message is appended to the output tuple. 

64 complex_func : bool, optional 

65 Indicate if the function's (`func`) return type is real 

66 (``complex_func=False``: default) or complex (``complex_func=True``). 

67 In both cases, the function's argument is real. 

68 If full_output is also non-zero, the `infodict`, `message`, and 

69 `explain` for the real and complex components are returned in 

70 a dictionary with keys "real output" and "imag output". 

71 

72 Returns 

73 ------- 

74 y : float 

75 The integral of func from `a` to `b`. 

76 abserr : float 

77 An estimate of the absolute error in the result. 

78 infodict : dict 

79 A dictionary containing additional information. 

80 message 

81 A convergence message. 

82 explain 

83 Appended only with 'cos' or 'sin' weighting and infinite 

84 integration limits, it contains an explanation of the codes in 

85 infodict['ierlst'] 

86 

87 Other Parameters 

88 ---------------- 

89 epsabs : float or int, optional 

90 Absolute error tolerance. Default is 1.49e-8. `quad` tries to obtain 

91 an accuracy of ``abs(i-result) <= max(epsabs, epsrel*abs(i))`` 

92 where ``i`` = integral of `func` from `a` to `b`, and ``result`` is the 

93 numerical approximation. See `epsrel` below. 

94 epsrel : float or int, optional 

95 Relative error tolerance. Default is 1.49e-8. 

96 If ``epsabs <= 0``, `epsrel` must be greater than both 5e-29 

97 and ``50 * (machine epsilon)``. See `epsabs` above. 

98 limit : float or int, optional 

99 An upper bound on the number of subintervals used in the adaptive 

100 algorithm. 

101 points : (sequence of floats,ints), optional 

102 A sequence of break points in the bounded integration interval 

103 where local difficulties of the integrand may occur (e.g., 

104 singularities, discontinuities). The sequence does not have 

105 to be sorted. Note that this option cannot be used in conjunction 

106 with ``weight``. 

107 weight : float or int, optional 

108 String indicating weighting function. Full explanation for this 

109 and the remaining arguments can be found below. 

110 wvar : optional 

111 Variables for use with weighting functions. 

112 wopts : optional 

113 Optional input for reusing Chebyshev moments. 

114 maxp1 : float or int, optional 

115 An upper bound on the number of Chebyshev moments. 

116 limlst : int, optional 

117 Upper bound on the number of cycles (>=3) for use with a sinusoidal 

118 weighting and an infinite end-point. 

119 

120 See Also 

121 -------- 

122 dblquad : double integral 

123 tplquad : triple integral 

124 nquad : n-dimensional integrals (uses `quad` recursively) 

125 fixed_quad : fixed-order Gaussian quadrature 

126 quadrature : adaptive Gaussian quadrature 

127 odeint : ODE integrator 

128 ode : ODE integrator 

129 simpson : integrator for sampled data 

130 romb : integrator for sampled data 

131 scipy.special : for coefficients and roots of orthogonal polynomials 

132 

133 Notes 

134 ----- 

135 For valid results, the integral must converge; behavior for divergent 

136 integrals is not guaranteed. 

137 

138 **Extra information for quad() inputs and outputs** 

139 

140 If full_output is non-zero, then the third output argument 

141 (infodict) is a dictionary with entries as tabulated below. For 

142 infinite limits, the range is transformed to (0,1) and the 

143 optional outputs are given with respect to this transformed range. 

144 Let M be the input argument limit and let K be infodict['last']. 

145 The entries are: 

146 

147 'neval' 

148 The number of function evaluations. 

149 'last' 

150 The number, K, of subintervals produced in the subdivision process. 

151 'alist' 

152 A rank-1 array of length M, the first K elements of which are the 

153 left end points of the subintervals in the partition of the 

154 integration range. 

155 'blist' 

156 A rank-1 array of length M, the first K elements of which are the 

157 right end points of the subintervals. 

158 'rlist' 

159 A rank-1 array of length M, the first K elements of which are the 

160 integral approximations on the subintervals. 

161 'elist' 

162 A rank-1 array of length M, the first K elements of which are the 

163 moduli of the absolute error estimates on the subintervals. 

164 'iord' 

165 A rank-1 integer array of length M, the first L elements of 

166 which are pointers to the error estimates over the subintervals 

167 with ``L=K`` if ``K<=M/2+2`` or ``L=M+1-K`` otherwise. Let I be the 

168 sequence ``infodict['iord']`` and let E be the sequence 

169 ``infodict['elist']``. Then ``E[I[1]], ..., E[I[L]]`` forms a 

170 decreasing sequence. 

171 

172 If the input argument points is provided (i.e., it is not None), 

173 the following additional outputs are placed in the output 

174 dictionary. Assume the points sequence is of length P. 

175 

176 'pts' 

177 A rank-1 array of length P+2 containing the integration limits 

178 and the break points of the intervals in ascending order. 

179 This is an array giving the subintervals over which integration 

180 will occur. 

181 'level' 

182 A rank-1 integer array of length M (=limit), containing the 

183 subdivision levels of the subintervals, i.e., if (aa,bb) is a 

184 subinterval of ``(pts[1], pts[2])`` where ``pts[0]`` and ``pts[2]`` 

185 are adjacent elements of ``infodict['pts']``, then (aa,bb) has level l 

186 if ``|bb-aa| = |pts[2]-pts[1]| * 2**(-l)``. 

187 'ndin' 

188 A rank-1 integer array of length P+2. After the first integration 

189 over the intervals (pts[1], pts[2]), the error estimates over some 

190 of the intervals may have been increased artificially in order to 

191 put their subdivision forward. This array has ones in slots 

192 corresponding to the subintervals for which this happens. 

193 

194 **Weighting the integrand** 

195 

196 The input variables, *weight* and *wvar*, are used to weight the 

197 integrand by a select list of functions. Different integration 

198 methods are used to compute the integral with these weighting 

199 functions, and these do not support specifying break points. The 

200 possible values of weight and the corresponding weighting functions are. 

201 

202 ========== =================================== ===================== 

203 ``weight`` Weight function used ``wvar`` 

204 ========== =================================== ===================== 

205 'cos' cos(w*x) wvar = w 

206 'sin' sin(w*x) wvar = w 

207 'alg' g(x) = ((x-a)**alpha)*((b-x)**beta) wvar = (alpha, beta) 

208 'alg-loga' g(x)*log(x-a) wvar = (alpha, beta) 

209 'alg-logb' g(x)*log(b-x) wvar = (alpha, beta) 

210 'alg-log' g(x)*log(x-a)*log(b-x) wvar = (alpha, beta) 

211 'cauchy' 1/(x-c) wvar = c 

212 ========== =================================== ===================== 

213 

214 wvar holds the parameter w, (alpha, beta), or c depending on the weight 

215 selected. In these expressions, a and b are the integration limits. 

216 

217 For the 'cos' and 'sin' weighting, additional inputs and outputs are 

218 available. 

219 

220 For finite integration limits, the integration is performed using a 

221 Clenshaw-Curtis method which uses Chebyshev moments. For repeated 

222 calculations, these moments are saved in the output dictionary: 

223 

224 'momcom' 

225 The maximum level of Chebyshev moments that have been computed, 

226 i.e., if ``M_c`` is ``infodict['momcom']`` then the moments have been 

227 computed for intervals of length ``|b-a| * 2**(-l)``, 

228 ``l=0,1,...,M_c``. 

229 'nnlog' 

230 A rank-1 integer array of length M(=limit), containing the 

231 subdivision levels of the subintervals, i.e., an element of this 

232 array is equal to l if the corresponding subinterval is 

233 ``|b-a|* 2**(-l)``. 

234 'chebmo' 

235 A rank-2 array of shape (25, maxp1) containing the computed 

236 Chebyshev moments. These can be passed on to an integration 

237 over the same interval by passing this array as the second 

238 element of the sequence wopts and passing infodict['momcom'] as 

239 the first element. 

240 

241 If one of the integration limits is infinite, then a Fourier integral is 

242 computed (assuming w neq 0). If full_output is 1 and a numerical error 

243 is encountered, besides the error message attached to the output tuple, 

244 a dictionary is also appended to the output tuple which translates the 

245 error codes in the array ``info['ierlst']`` to English messages. The 

246 output information dictionary contains the following entries instead of 

247 'last', 'alist', 'blist', 'rlist', and 'elist': 

248 

249 'lst' 

250 The number of subintervals needed for the integration (call it ``K_f``). 

251 'rslst' 

252 A rank-1 array of length M_f=limlst, whose first ``K_f`` elements 

253 contain the integral contribution over the interval 

254 ``(a+(k-1)c, a+kc)`` where ``c = (2*floor(|w|) + 1) * pi / |w|`` 

255 and ``k=1,2,...,K_f``. 

256 'erlst' 

257 A rank-1 array of length ``M_f`` containing the error estimate 

258 corresponding to the interval in the same position in 

259 ``infodict['rslist']``. 

260 'ierlst' 

261 A rank-1 integer array of length ``M_f`` containing an error flag 

262 corresponding to the interval in the same position in 

263 ``infodict['rslist']``. See the explanation dictionary (last entry 

264 in the output tuple) for the meaning of the codes. 

265 

266 

267 **Details of QUADPACK level routines** 

268 

269 `quad` calls routines from the FORTRAN library QUADPACK. This section 

270 provides details on the conditions for each routine to be called and a 

271 short description of each routine. The routine called depends on 

272 `weight`, `points` and the integration limits `a` and `b`. 

273 

274 ================ ============== ========== ===================== 

275 QUADPACK routine `weight` `points` infinite bounds 

276 ================ ============== ========== ===================== 

277 qagse None No No 

278 qagie None No Yes 

279 qagpe None Yes No 

280 qawoe 'sin', 'cos' No No 

281 qawfe 'sin', 'cos' No either `a` or `b` 

282 qawse 'alg*' No No 

283 qawce 'cauchy' No No 

284 ================ ============== ========== ===================== 

285 

286 The following provides a short desciption from [1]_ for each 

287 routine. 

288 

289 qagse 

290 is an integrator based on globally adaptive interval 

291 subdivision in connection with extrapolation, which will 

292 eliminate the effects of integrand singularities of 

293 several types. 

294 qagie 

295 handles integration over infinite intervals. The infinite range is 

296 mapped onto a finite interval and subsequently the same strategy as 

297 in ``QAGS`` is applied. 

298 qagpe 

299 serves the same purposes as QAGS, but also allows the 

300 user to provide explicit information about the location 

301 and type of trouble-spots i.e. the abscissae of internal 

302 singularities, discontinuities and other difficulties of 

303 the integrand function. 

304 qawoe 

305 is an integrator for the evaluation of 

306 :math:`\\int^b_a \\cos(\\omega x)f(x)dx` or 

307 :math:`\\int^b_a \\sin(\\omega x)f(x)dx` 

308 over a finite interval [a,b], where :math:`\\omega` and :math:`f` 

309 are specified by the user. The rule evaluation component is based 

310 on the modified Clenshaw-Curtis technique 

311 

312 An adaptive subdivision scheme is used in connection 

313 with an extrapolation procedure, which is a modification 

314 of that in ``QAGS`` and allows the algorithm to deal with 

315 singularities in :math:`f(x)`. 

316 qawfe 

317 calculates the Fourier transform 

318 :math:`\\int^\\infty_a \\cos(\\omega x)f(x)dx` or 

319 :math:`\\int^\\infty_a \\sin(\\omega x)f(x)dx` 

320 for user-provided :math:`\\omega` and :math:`f`. The procedure of 

321 ``QAWO`` is applied on successive finite intervals, and convergence 

322 acceleration by means of the :math:`\\varepsilon`-algorithm is applied 

323 to the series of integral approximations. 

324 qawse 

325 approximate :math:`\\int^b_a w(x)f(x)dx`, with :math:`a < b` where 

326 :math:`w(x) = (x-a)^{\\alpha}(b-x)^{\\beta}v(x)` with 

327 :math:`\\alpha,\\beta > -1`, where :math:`v(x)` may be one of the 

328 following functions: :math:`1`, :math:`\\log(x-a)`, :math:`\\log(b-x)`, 

329 :math:`\\log(x-a)\\log(b-x)`. 

330 

331 The user specifies :math:`\\alpha`, :math:`\\beta` and the type of the 

332 function :math:`v`. A globally adaptive subdivision strategy is 

333 applied, with modified Clenshaw-Curtis integration on those 

334 subintervals which contain `a` or `b`. 

335 qawce 

336 compute :math:`\\int^b_a f(x) / (x-c)dx` where the integral must be 

337 interpreted as a Cauchy principal value integral, for user specified 

338 :math:`c` and :math:`f`. The strategy is globally adaptive. Modified 

339 Clenshaw-Curtis integration is used on those intervals containing the 

340 point :math:`x = c`. 

341 

342 **Integration of Complex Function of a Real Variable** 

343 

344 A complex valued function, :math:`f`, of a real variable can be written as 

345 :math:`f = g + ih`. Similarly, the integral of :math:`f` can be 

346 written as 

347 

348 .. math:: 

349 \\int_a^b f(x) dx = \\int_a^b g(x) dx + i\\int_a^b h(x) dx 

350 

351 assuming that the integrals of :math:`g` and :math:`h` exist 

352 over the inteval :math:`[a,b]` [2]_. Therefore, ``quad`` integrates 

353 complex-valued functions by integrating the real and imaginary components 

354 separately. 

355 

356 

357 References 

358 ---------- 

359 

360 .. [1] Piessens, Robert; de Doncker-Kapenga, Elise; 

361 Überhuber, Christoph W.; Kahaner, David (1983). 

362 QUADPACK: A subroutine package for automatic integration. 

363 Springer-Verlag. 

364 ISBN 978-3-540-12553-2. 

365 

366 .. [2] McCullough, Thomas; Phillips, Keith (1973). 

367 Foundations of Analysis in the Complex Plane. 

368 Holt Rinehart Winston. 

369 ISBN 0-03-086370-8 

370 

371 Examples 

372 -------- 

373 Calculate :math:`\\int^4_0 x^2 dx` and compare with an analytic result 

374 

375 >>> from scipy import integrate 

376 >>> import numpy as np 

377 >>> x2 = lambda x: x**2 

378 >>> integrate.quad(x2, 0, 4) 

379 (21.333333333333332, 2.3684757858670003e-13) 

380 >>> print(4**3 / 3.) # analytical result 

381 21.3333333333 

382 

383 Calculate :math:`\\int^\\infty_0 e^{-x} dx` 

384 

385 >>> invexp = lambda x: np.exp(-x) 

386 >>> integrate.quad(invexp, 0, np.inf) 

387 (1.0, 5.842605999138044e-11) 

388 

389 Calculate :math:`\\int^1_0 a x \\,dx` for :math:`a = 1, 3` 

390 

391 >>> f = lambda x, a: a*x 

392 >>> y, err = integrate.quad(f, 0, 1, args=(1,)) 

393 >>> y 

394 0.5 

395 >>> y, err = integrate.quad(f, 0, 1, args=(3,)) 

396 >>> y 

397 1.5 

398 

399 Calculate :math:`\\int^1_0 x^2 + y^2 dx` with ctypes, holding 

400 y parameter as 1:: 

401 

402 testlib.c => 

403 double func(int n, double args[n]){ 

404 return args[0]*args[0] + args[1]*args[1];} 

405 compile to library testlib.* 

406 

407 :: 

408 

409 from scipy import integrate 

410 import ctypes 

411 lib = ctypes.CDLL('/home/.../testlib.*') #use absolute path 

412 lib.func.restype = ctypes.c_double 

413 lib.func.argtypes = (ctypes.c_int,ctypes.c_double) 

414 integrate.quad(lib.func,0,1,(1)) 

415 #(1.3333333333333333, 1.4802973661668752e-14) 

416 print((1.0**3/3.0 + 1.0) - (0.0**3/3.0 + 0.0)) #Analytic result 

417 # 1.3333333333333333 

418 

419 Be aware that pulse shapes and other sharp features as compared to the 

420 size of the integration interval may not be integrated correctly using 

421 this method. A simplified example of this limitation is integrating a 

422 y-axis reflected step function with many zero values within the integrals 

423 bounds. 

424 

425 >>> y = lambda x: 1 if x<=0 else 0 

426 >>> integrate.quad(y, -1, 1) 

427 (1.0, 1.1102230246251565e-14) 

428 >>> integrate.quad(y, -1, 100) 

429 (1.0000000002199108, 1.0189464580163188e-08) 

430 >>> integrate.quad(y, -1, 10000) 

431 (0.0, 0.0) 

432 

433 """ 

434 if not isinstance(args, tuple): 

435 args = (args,) 

436 

437 # check the limits of integration: \int_a^b, expect a < b 

438 flip, a, b = b < a, min(a, b), max(a, b) 

439 

440 if complex_func: 

441 def imfunc(x, *args): 

442 return np.imag(func(x, *args)) 

443 

444 def refunc(x, *args): 

445 return np.real(func(x, *args)) 

446 

447 re_retval = quad(refunc, a, b, args, full_output, epsabs, 

448 epsrel, limit, points, weight, wvar, wopts, 

449 maxp1, limlst, complex_func=False) 

450 im_retval = quad(imfunc, a, b, args, full_output, epsabs, 

451 epsrel, limit, points, weight, wvar, wopts, 

452 maxp1, limlst, complex_func=False) 

453 integral = re_retval[0] + 1j*im_retval[0] 

454 error_estimate = re_retval[1] + 1j*im_retval[1] 

455 retval = integral, error_estimate 

456 if full_output: 

457 msgexp = {} 

458 msgexp["real"] = re_retval[2:] 

459 msgexp["imag"] = im_retval[2:] 

460 retval = retval + (msgexp,) 

461 

462 return retval 

463 

464 if weight is None: 

465 retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit, 

466 points) 

467 else: 

468 if points is not None: 

469 msg = ("Break points cannot be specified when using weighted integrand.\n" 

470 "Continuing, ignoring specified points.") 

471 warnings.warn(msg, IntegrationWarning, stacklevel=2) 

472 retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel, 

473 limlst, limit, maxp1, weight, wvar, wopts) 

474 

475 if flip: 

476 retval = (-retval[0],) + retval[1:] 

477 

478 ier = retval[-1] 

479 if ier == 0: 

480 return retval[:-1] 

481 

482 msgs = {80: "A Python error occurred possibly while calling the function.", 

483 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, 

484 2: "The occurrence of roundoff error is detected, which prevents \n the requested tolerance from being achieved. The error may be \n underestimated.", 

485 3: "Extremely bad integrand behavior occurs at some points of the\n integration interval.", 

486 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.", 

487 5: "The integral is probably divergent, or slowly convergent.", 

488 6: "The input is invalid.", 

489 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.", 

490 'unknown': "Unknown error."} 

491 

492 if weight in ['cos','sin'] and (b == Inf or a == -Inf): 

493 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." 

494 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." 

495 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." 

496 explain = {1: "The maximum number of subdivisions (= limit) has been \n achieved on this cycle.", 

497 2: "The occurrence of roundoff error is detected and prevents\n the tolerance imposed on this cycle from being achieved.", 

498 3: "Extremely bad integrand behavior occurs at some points of\n this cycle.", 

499 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.", 

500 5: "The integral over this cycle is probably divergent or slowly convergent."} 

501 

502 try: 

503 msg = msgs[ier] 

504 except KeyError: 

505 msg = msgs['unknown'] 

506 

507 if ier in [1,2,3,4,5,7]: 

508 if full_output: 

509 if weight in ['cos', 'sin'] and (b == Inf or a == -Inf): 

510 return retval[:-1] + (msg, explain) 

511 else: 

512 return retval[:-1] + (msg,) 

513 else: 

514 warnings.warn(msg, IntegrationWarning, stacklevel=2) 

515 return retval[:-1] 

516 

517 elif ier == 6: # Forensic decision tree when QUADPACK throws ier=6 

518 if epsabs <= 0: # Small error tolerance - applies to all methods 

519 if epsrel < max(50 * sys.float_info.epsilon, 5e-29): 

520 msg = ("If 'epsabs'<=0, 'epsrel' must be greater than both" 

521 " 5e-29 and 50*(machine epsilon).") 

522 elif weight in ['sin', 'cos'] and (abs(a) + abs(b) == Inf): 

523 msg = ("Sine or cosine weighted intergals with infinite domain" 

524 " must have 'epsabs'>0.") 

525 

526 elif weight is None: 

527 if points is None: # QAGSE/QAGIE 

528 msg = ("Invalid 'limit' argument. There must be" 

529 " at least one subinterval") 

530 else: # QAGPE 

531 if not (min(a, b) <= min(points) <= max(points) <= max(a, b)): 

532 msg = ("All break points in 'points' must lie within the" 

533 " integration limits.") 

534 elif len(points) >= limit: 

535 msg = ("Number of break points ({:d})" 

536 " must be less than subinterval" 

537 " limit ({:d})").format(len(points), limit) 

538 

539 else: 

540 if maxp1 < 1: 

541 msg = "Chebyshev moment limit maxp1 must be >=1." 

542 

543 elif weight in ('cos', 'sin') and abs(a+b) == Inf: # QAWFE 

544 msg = "Cycle limit limlst must be >=3." 

545 

546 elif weight.startswith('alg'): # QAWSE 

547 if min(wvar) < -1: 

548 msg = "wvar parameters (alpha, beta) must both be >= -1." 

549 if b < a: 

550 msg = "Integration limits a, b must satistfy a<b." 

551 

552 elif weight == 'cauchy' and wvar in (a, b): 

553 msg = ("Parameter 'wvar' must not equal" 

554 " integration limits 'a' or 'b'.") 

555 

556 raise ValueError(msg) 

557 

558 

559def _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points): 

560 infbounds = 0 

561 if (b != Inf and a != -Inf): 

562 pass # standard integration 

563 elif (b == Inf and a != -Inf): 

564 infbounds = 1 

565 bound = a 

566 elif (b == Inf and a == -Inf): 

567 infbounds = 2 

568 bound = 0 # ignored 

569 elif (b != Inf and a == -Inf): 

570 infbounds = -1 

571 bound = b 

572 else: 

573 raise RuntimeError("Infinity comparisons don't work for you.") 

574 

575 if points is None: 

576 if infbounds == 0: 

577 return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit) 

578 else: 

579 return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit) 

580 else: 

581 if infbounds != 0: 

582 raise ValueError("Infinity inputs cannot be used with break points.") 

583 else: 

584 #Duplicates force function evaluation at singular points 

585 the_points = np.unique(points) 

586 the_points = the_points[a < the_points] 

587 the_points = the_points[the_points < b] 

588 the_points = np.concatenate((the_points, (0., 0.))) 

589 return _quadpack._qagpe(func,a,b,the_points,args,full_output,epsabs,epsrel,limit) 

590 

591 

592def _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts): 

593 if weight not in ['cos','sin','alg','alg-loga','alg-logb','alg-log','cauchy']: 

594 raise ValueError("%s not a recognized weighting function." % weight) 

595 

596 strdict = {'cos':1,'sin':2,'alg':1,'alg-loga':2,'alg-logb':3,'alg-log':4} 

597 

598 if weight in ['cos','sin']: 

599 integr = strdict[weight] 

600 if (b != Inf and a != -Inf): # finite limits 

601 if wopts is None: # no precomputed Chebyshev moments 

602 return _quadpack._qawoe(func, a, b, wvar, integr, args, full_output, 

603 epsabs, epsrel, limit, maxp1,1) 

604 else: # precomputed Chebyshev moments 

605 momcom = wopts[0] 

606 chebcom = wopts[1] 

607 return _quadpack._qawoe(func, a, b, wvar, integr, args, full_output, 

608 epsabs, epsrel, limit, maxp1, 2, momcom, chebcom) 

609 

610 elif (b == Inf and a != -Inf): 

611 return _quadpack._qawfe(func, a, wvar, integr, args, full_output, 

612 epsabs,limlst,limit,maxp1) 

613 elif (b != Inf and a == -Inf): # remap function and interval 

614 if weight == 'cos': 

615 def thefunc(x,*myargs): 

616 y = -x 

617 func = myargs[0] 

618 myargs = (y,) + myargs[1:] 

619 return func(*myargs) 

620 else: 

621 def thefunc(x,*myargs): 

622 y = -x 

623 func = myargs[0] 

624 myargs = (y,) + myargs[1:] 

625 return -func(*myargs) 

626 args = (func,) + args 

627 return _quadpack._qawfe(thefunc, -b, wvar, integr, args, 

628 full_output, epsabs, limlst, limit, maxp1) 

629 else: 

630 raise ValueError("Cannot integrate with this weight from -Inf to +Inf.") 

631 else: 

632 if a in [-Inf,Inf] or b in [-Inf,Inf]: 

633 raise ValueError("Cannot integrate with this weight over an infinite interval.") 

634 

635 if weight.startswith('alg'): 

636 integr = strdict[weight] 

637 return _quadpack._qawse(func, a, b, wvar, integr, args, 

638 full_output, epsabs, epsrel, limit) 

639 else: # weight == 'cauchy' 

640 return _quadpack._qawce(func, a, b, wvar, args, full_output, 

641 epsabs, epsrel, limit) 

642 

643 

644def dblquad(func, a, b, gfun, hfun, args=(), epsabs=1.49e-8, epsrel=1.49e-8): 

645 """ 

646 Compute a double integral. 

647 

648 Return the double (definite) integral of ``func(y, x)`` from ``x = a..b`` 

649 and ``y = gfun(x)..hfun(x)``. 

650 

651 Parameters 

652 ---------- 

653 func : callable 

654 A Python function or method of at least two variables: y must be the 

655 first argument and x the second argument. 

656 a, b : float 

657 The limits of integration in x: `a` < `b` 

658 gfun : callable or float 

659 The lower boundary curve in y which is a function taking a single 

660 floating point argument (x) and returning a floating point result 

661 or a float indicating a constant boundary curve. 

662 hfun : callable or float 

663 The upper boundary curve in y (same requirements as `gfun`). 

664 args : sequence, optional 

665 Extra arguments to pass to `func`. 

666 epsabs : float, optional 

667 Absolute tolerance passed directly to the inner 1-D quadrature 

668 integration. Default is 1.49e-8. ``dblquad`` tries to obtain 

669 an accuracy of ``abs(i-result) <= max(epsabs, epsrel*abs(i))`` 

670 where ``i`` = inner integral of ``func(y, x)`` from ``gfun(x)`` 

671 to ``hfun(x)``, and ``result`` is the numerical approximation. 

672 See `epsrel` below. 

673 epsrel : float, optional 

674 Relative tolerance of the inner 1-D integrals. Default is 1.49e-8. 

675 If ``epsabs <= 0``, `epsrel` must be greater than both 5e-29 

676 and ``50 * (machine epsilon)``. See `epsabs` above. 

677 

678 Returns 

679 ------- 

680 y : float 

681 The resultant integral. 

682 abserr : float 

683 An estimate of the error. 

684 

685 See Also 

686 -------- 

687 quad : single integral 

688 tplquad : triple integral 

689 nquad : N-dimensional integrals 

690 fixed_quad : fixed-order Gaussian quadrature 

691 quadrature : adaptive Gaussian quadrature 

692 odeint : ODE integrator 

693 ode : ODE integrator 

694 simpson : integrator for sampled data 

695 romb : integrator for sampled data 

696 scipy.special : for coefficients and roots of orthogonal polynomials 

697 

698 

699 Notes 

700 ----- 

701 For valid results, the integral must converge; behavior for divergent 

702 integrals is not guaranteed. 

703 

704 **Details of QUADPACK level routines** 

705 

706 `quad` calls routines from the FORTRAN library QUADPACK. This section 

707 provides details on the conditions for each routine to be called and a 

708 short description of each routine. For each level of integration, ``qagse`` 

709 is used for finite limits or ``qagie`` is used if either limit (or both!) 

710 are infinite. The following provides a short description from [1]_ for each 

711 routine. 

712 

713 qagse 

714 is an integrator based on globally adaptive interval 

715 subdivision in connection with extrapolation, which will 

716 eliminate the effects of integrand singularities of 

717 several types. 

718 qagie 

719 handles integration over infinite intervals. The infinite range is 

720 mapped onto a finite interval and subsequently the same strategy as 

721 in ``QAGS`` is applied. 

722 

723 References 

724 ---------- 

725 

726 .. [1] Piessens, Robert; de Doncker-Kapenga, Elise; 

727 Überhuber, Christoph W.; Kahaner, David (1983). 

728 QUADPACK: A subroutine package for automatic integration. 

729 Springer-Verlag. 

730 ISBN 978-3-540-12553-2. 

731 

732 Examples 

733 -------- 

734 Compute the double integral of ``x * y**2`` over the box 

735 ``x`` ranging from 0 to 2 and ``y`` ranging from 0 to 1. 

736 That is, :math:`\\int^{x=2}_{x=0} \\int^{y=1}_{y=0} x y^2 \\,dy \\,dx`. 

737 

738 >>> import numpy as np 

739 >>> from scipy import integrate 

740 >>> f = lambda y, x: x*y**2 

741 >>> integrate.dblquad(f, 0, 2, 0, 1) 

742 (0.6666666666666667, 7.401486830834377e-15) 

743 

744 Calculate :math:`\\int^{x=\\pi/4}_{x=0} \\int^{y=\\cos(x)}_{y=\\sin(x)} 1 

745 \\,dy \\,dx`. 

746 

747 >>> f = lambda y, x: 1 

748 >>> integrate.dblquad(f, 0, np.pi/4, np.sin, np.cos) 

749 (0.41421356237309503, 1.1083280054755938e-14) 

750 

751 Calculate :math:`\\int^{x=1}_{x=0} \\int^{y=2-x}_{y=x} a x y \\,dy \\,dx` 

752 for :math:`a=1, 3`. 

753 

754 >>> f = lambda y, x, a: a*x*y 

755 >>> integrate.dblquad(f, 0, 1, lambda x: x, lambda x: 2-x, args=(1,)) 

756 (0.33333333333333337, 5.551115123125783e-15) 

757 >>> integrate.dblquad(f, 0, 1, lambda x: x, lambda x: 2-x, args=(3,)) 

758 (0.9999999999999999, 1.6653345369377348e-14) 

759 

760 Compute the two-dimensional Gaussian Integral, which is the integral of the 

761 Gaussian function :math:`f(x,y) = e^{-(x^{2} + y^{2})}`, over 

762 :math:`(-\\infty,+\\infty)`. That is, compute the integral 

763 :math:`\\iint^{+\\infty}_{-\\infty} e^{-(x^{2} + y^{2})} \\,dy\\,dx`. 

764 

765 >>> f = lambda x, y: np.exp(-(x ** 2 + y ** 2)) 

766 >>> integrate.dblquad(f, -np.inf, np.inf, -np.inf, np.inf) 

767 (3.141592653589777, 2.5173086737433208e-08) 

768 

769 """ 

770 

771 def temp_ranges(*args): 

772 return [gfun(args[0]) if callable(gfun) else gfun, 

773 hfun(args[0]) if callable(hfun) else hfun] 

774 

775 return nquad(func, [temp_ranges, [a, b]], args=args, 

776 opts={"epsabs": epsabs, "epsrel": epsrel}) 

777 

778 

779def tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-8, 

780 epsrel=1.49e-8): 

781 """ 

782 Compute a triple (definite) integral. 

783 

784 Return the triple integral of ``func(z, y, x)`` from ``x = a..b``, 

785 ``y = gfun(x)..hfun(x)``, and ``z = qfun(x,y)..rfun(x,y)``. 

786 

787 Parameters 

788 ---------- 

789 func : function 

790 A Python function or method of at least three variables in the 

791 order (z, y, x). 

792 a, b : float 

793 The limits of integration in x: `a` < `b` 

794 gfun : function or float 

795 The lower boundary curve in y which is a function taking a single 

796 floating point argument (x) and returning a floating point result 

797 or a float indicating a constant boundary curve. 

798 hfun : function or float 

799 The upper boundary curve in y (same requirements as `gfun`). 

800 qfun : function or float 

801 The lower boundary surface in z. It must be a function that takes 

802 two floats in the order (x, y) and returns a float or a float 

803 indicating a constant boundary surface. 

804 rfun : function or float 

805 The upper boundary surface in z. (Same requirements as `qfun`.) 

806 args : tuple, optional 

807 Extra arguments to pass to `func`. 

808 epsabs : float, optional 

809 Absolute tolerance passed directly to the innermost 1-D quadrature 

810 integration. Default is 1.49e-8. 

811 epsrel : float, optional 

812 Relative tolerance of the innermost 1-D integrals. Default is 1.49e-8. 

813 

814 Returns 

815 ------- 

816 y : float 

817 The resultant integral. 

818 abserr : float 

819 An estimate of the error. 

820 

821 See Also 

822 -------- 

823 quad : Adaptive quadrature using QUADPACK 

824 quadrature : Adaptive Gaussian quadrature 

825 fixed_quad : Fixed-order Gaussian quadrature 

826 dblquad : Double integrals 

827 nquad : N-dimensional integrals 

828 romb : Integrators for sampled data 

829 simpson : Integrators for sampled data 

830 ode : ODE integrators 

831 odeint : ODE integrators 

832 scipy.special : For coefficients and roots of orthogonal polynomials 

833 

834 Notes 

835 ----- 

836 For valid results, the integral must converge; behavior for divergent 

837 integrals is not guaranteed. 

838 

839 **Details of QUADPACK level routines** 

840 

841 `quad` calls routines from the FORTRAN library QUADPACK. This section 

842 provides details on the conditions for each routine to be called and a 

843 short description of each routine. For each level of integration, ``qagse`` 

844 is used for finite limits or ``qagie`` is used, if either limit (or both!) 

845 are infinite. The following provides a short description from [1]_ for each 

846 routine. 

847 

848 qagse 

849 is an integrator based on globally adaptive interval 

850 subdivision in connection with extrapolation, which will 

851 eliminate the effects of integrand singularities of 

852 several types. 

853 qagie 

854 handles integration over infinite intervals. The infinite range is 

855 mapped onto a finite interval and subsequently the same strategy as 

856 in ``QAGS`` is applied. 

857 

858 References 

859 ---------- 

860 

861 .. [1] Piessens, Robert; de Doncker-Kapenga, Elise; 

862 Überhuber, Christoph W.; Kahaner, David (1983). 

863 QUADPACK: A subroutine package for automatic integration. 

864 Springer-Verlag. 

865 ISBN 978-3-540-12553-2. 

866 

867 Examples 

868 -------- 

869 Compute the triple integral of ``x * y * z``, over ``x`` ranging 

870 from 1 to 2, ``y`` ranging from 2 to 3, ``z`` ranging from 0 to 1. 

871 That is, :math:`\\int^{x=2}_{x=1} \\int^{y=3}_{y=2} \\int^{z=1}_{z=0} x y z 

872 \\,dz \\,dy \\,dx`. 

873 

874 >>> import numpy as np 

875 >>> from scipy import integrate 

876 >>> f = lambda z, y, x: x*y*z 

877 >>> integrate.tplquad(f, 1, 2, 2, 3, 0, 1) 

878 (1.8749999999999998, 3.3246447942574074e-14) 

879 

880 Calculate :math:`\\int^{x=1}_{x=0} \\int^{y=1-2x}_{y=0} 

881 \\int^{z=1-x-2y}_{z=0} x y z \\,dz \\,dy \\,dx`. 

882 Note: `qfun`/`rfun` takes arguments in the order (x, y), even though ``f`` 

883 takes arguments in the order (z, y, x). 

884 

885 >>> f = lambda z, y, x: x*y*z 

886 >>> integrate.tplquad(f, 0, 1, 0, lambda x: 1-2*x, 0, lambda x, y: 1-x-2*y) 

887 (0.05416666666666668, 2.1774196738157757e-14) 

888 

889 Calculate :math:`\\int^{x=1}_{x=0} \\int^{y=1}_{y=0} \\int^{z=1}_{z=0} 

890 a x y z \\,dz \\,dy \\,dx` for :math:`a=1, 3`. 

891 

892 >>> f = lambda z, y, x, a: a*x*y*z 

893 >>> integrate.tplquad(f, 0, 1, 0, 1, 0, 1, args=(1,)) 

894 (0.125, 5.527033708952211e-15) 

895 >>> integrate.tplquad(f, 0, 1, 0, 1, 0, 1, args=(3,)) 

896 (0.375, 1.6581101126856635e-14) 

897 

898 Compute the three-dimensional Gaussian Integral, which is the integral of 

899 the Gaussian function :math:`f(x,y,z) = e^{-(x^{2} + y^{2} + z^{2})}`, over 

900 :math:`(-\\infty,+\\infty)`. That is, compute the integral 

901 :math:`\\iiint^{+\\infty}_{-\\infty} e^{-(x^{2} + y^{2} + z^{2})} \\,dz 

902 \\,dy\\,dx`. 

903 

904 >>> f = lambda x, y, z: np.exp(-(x ** 2 + y ** 2 + z ** 2)) 

905 >>> integrate.tplquad(f, -np.inf, np.inf, -np.inf, np.inf, -np.inf, np.inf) 

906 (5.568327996830833, 4.4619078828029765e-08) 

907 

908 """ 

909 # f(z, y, x) 

910 # qfun/rfun(x, y) 

911 # gfun/hfun(x) 

912 # nquad will hand (y, x, t0, ...) to ranges0 

913 # nquad will hand (x, t0, ...) to ranges1 

914 # Only qfun / rfun is different API... 

915 

916 def ranges0(*args): 

917 return [qfun(args[1], args[0]) if callable(qfun) else qfun, 

918 rfun(args[1], args[0]) if callable(rfun) else rfun] 

919 

920 def ranges1(*args): 

921 return [gfun(args[0]) if callable(gfun) else gfun, 

922 hfun(args[0]) if callable(hfun) else hfun] 

923 

924 ranges = [ranges0, ranges1, [a, b]] 

925 return nquad(func, ranges, args=args, 

926 opts={"epsabs": epsabs, "epsrel": epsrel}) 

927 

928 

929def nquad(func, ranges, args=None, opts=None, full_output=False): 

930 r""" 

931 Integration over multiple variables. 

932 

933 Wraps `quad` to enable integration over multiple variables. 

934 Various options allow improved integration of discontinuous functions, as 

935 well as the use of weighted integration, and generally finer control of the 

936 integration process. 

937 

938 Parameters 

939 ---------- 

940 func : {callable, scipy.LowLevelCallable} 

941 The function to be integrated. Has arguments of ``x0, ... xn``, 

942 ``t0, ... tm``, where integration is carried out over ``x0, ... xn``, 

943 which must be floats. Where ``t0, ... tm`` are extra arguments 

944 passed in args. 

945 Function signature should be ``func(x0, x1, ..., xn, t0, t1, ..., tm)``. 

946 Integration is carried out in order. That is, integration over ``x0`` 

947 is the innermost integral, and ``xn`` is the outermost. 

948 

949 If the user desires improved integration performance, then `f` may 

950 be a `scipy.LowLevelCallable` with one of the signatures:: 

951 

952 double func(int n, double *xx) 

953 double func(int n, double *xx, void *user_data) 

954 

955 where ``n`` is the number of variables and args. The ``xx`` array 

956 contains the coordinates and extra arguments. ``user_data`` is the data 

957 contained in the `scipy.LowLevelCallable`. 

958 ranges : iterable object 

959 Each element of ranges may be either a sequence of 2 numbers, or else 

960 a callable that returns such a sequence. ``ranges[0]`` corresponds to 

961 integration over x0, and so on. If an element of ranges is a callable, 

962 then it will be called with all of the integration arguments available, 

963 as well as any parametric arguments. e.g., if 

964 ``func = f(x0, x1, x2, t0, t1)``, then ``ranges[0]`` may be defined as 

965 either ``(a, b)`` or else as ``(a, b) = range0(x1, x2, t0, t1)``. 

966 args : iterable object, optional 

967 Additional arguments ``t0, ... tn``, required by ``func``, ``ranges``, 

968 and ``opts``. 

969 opts : iterable object or dict, optional 

970 Options to be passed to `quad`. May be empty, a dict, or 

971 a sequence of dicts or functions that return a dict. If empty, the 

972 default options from scipy.integrate.quad are used. If a dict, the same 

973 options are used for all levels of integraion. If a sequence, then each 

974 element of the sequence corresponds to a particular integration. e.g., 

975 ``opts[0]`` corresponds to integration over ``x0``, and so on. If a 

976 callable, the signature must be the same as for ``ranges``. The 

977 available options together with their default values are: 

978 

979 - epsabs = 1.49e-08 

980 - epsrel = 1.49e-08 

981 - limit = 50 

982 - points = None 

983 - weight = None 

984 - wvar = None 

985 - wopts = None 

986 

987 For more information on these options, see `quad`. 

988 

989 full_output : bool, optional 

990 Partial implementation of ``full_output`` from scipy.integrate.quad. 

991 The number of integrand function evaluations ``neval`` can be obtained 

992 by setting ``full_output=True`` when calling nquad. 

993 

994 Returns 

995 ------- 

996 result : float 

997 The result of the integration. 

998 abserr : float 

999 The maximum of the estimates of the absolute error in the various 

1000 integration results. 

1001 out_dict : dict, optional 

1002 A dict containing additional information on the integration. 

1003 

1004 See Also 

1005 -------- 

1006 quad : 1-D numerical integration 

1007 dblquad, tplquad : double and triple integrals 

1008 fixed_quad : fixed-order Gaussian quadrature 

1009 quadrature : adaptive Gaussian quadrature 

1010 

1011 Notes 

1012 ----- 

1013 For valid results, the integral must converge; behavior for divergent 

1014 integrals is not guaranteed. 

1015 

1016 **Details of QUADPACK level routines** 

1017 

1018 `nquad` calls routines from the FORTRAN library QUADPACK. This section 

1019 provides details on the conditions for each routine to be called and a 

1020 short description of each routine. The routine called depends on 

1021 `weight`, `points` and the integration limits `a` and `b`. 

1022 

1023 ================ ============== ========== ===================== 

1024 QUADPACK routine `weight` `points` infinite bounds 

1025 ================ ============== ========== ===================== 

1026 qagse None No No 

1027 qagie None No Yes 

1028 qagpe None Yes No 

1029 qawoe 'sin', 'cos' No No 

1030 qawfe 'sin', 'cos' No either `a` or `b` 

1031 qawse 'alg*' No No 

1032 qawce 'cauchy' No No 

1033 ================ ============== ========== ===================== 

1034 

1035 The following provides a short desciption from [1]_ for each 

1036 routine. 

1037 

1038 qagse 

1039 is an integrator based on globally adaptive interval 

1040 subdivision in connection with extrapolation, which will 

1041 eliminate the effects of integrand singularities of 

1042 several types. 

1043 qagie 

1044 handles integration over infinite intervals. The infinite range is 

1045 mapped onto a finite interval and subsequently the same strategy as 

1046 in ``QAGS`` is applied. 

1047 qagpe 

1048 serves the same purposes as QAGS, but also allows the 

1049 user to provide explicit information about the location 

1050 and type of trouble-spots i.e. the abscissae of internal 

1051 singularities, discontinuities and other difficulties of 

1052 the integrand function. 

1053 qawoe 

1054 is an integrator for the evaluation of 

1055 :math:`\int^b_a \cos(\omega x)f(x)dx` or 

1056 :math:`\int^b_a \sin(\omega x)f(x)dx` 

1057 over a finite interval [a,b], where :math:`\omega` and :math:`f` 

1058 are specified by the user. The rule evaluation component is based 

1059 on the modified Clenshaw-Curtis technique 

1060 

1061 An adaptive subdivision scheme is used in connection 

1062 with an extrapolation procedure, which is a modification 

1063 of that in ``QAGS`` and allows the algorithm to deal with 

1064 singularities in :math:`f(x)`. 

1065 qawfe 

1066 calculates the Fourier transform 

1067 :math:`\int^\infty_a \cos(\omega x)f(x)dx` or 

1068 :math:`\int^\infty_a \sin(\omega x)f(x)dx` 

1069 for user-provided :math:`\omega` and :math:`f`. The procedure of 

1070 ``QAWO`` is applied on successive finite intervals, and convergence 

1071 acceleration by means of the :math:`\varepsilon`-algorithm is applied 

1072 to the series of integral approximations. 

1073 qawse 

1074 approximate :math:`\int^b_a w(x)f(x)dx`, with :math:`a < b` where 

1075 :math:`w(x) = (x-a)^{\alpha}(b-x)^{\beta}v(x)` with 

1076 :math:`\alpha,\beta > -1`, where :math:`v(x)` may be one of the 

1077 following functions: :math:`1`, :math:`\log(x-a)`, :math:`\log(b-x)`, 

1078 :math:`\log(x-a)\log(b-x)`. 

1079 

1080 The user specifies :math:`\alpha`, :math:`\beta` and the type of the 

1081 function :math:`v`. A globally adaptive subdivision strategy is 

1082 applied, with modified Clenshaw-Curtis integration on those 

1083 subintervals which contain `a` or `b`. 

1084 qawce 

1085 compute :math:`\int^b_a f(x) / (x-c)dx` where the integral must be 

1086 interpreted as a Cauchy principal value integral, for user specified