Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-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 

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

137 

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

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

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

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

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

143 The entries are: 

144 

145 'neval' 

146 The number of function evaluations. 

147 'last' 

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

149 'alist' 

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

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

152 integration range. 

153 'blist' 

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

155 right end points of the subintervals. 

156 'rlist' 

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

158 integral approximations on the subintervals. 

159 'elist' 

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

161 moduli of the absolute error estimates on the subintervals. 

162 'iord' 

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

164 which are pointers to the error estimates over the subintervals 

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

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

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

168 decreasing sequence. 

169 

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

171 the following additional outputs are placed in the output 

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

173 

174 'pts' 

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

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

177 This is an array giving the subintervals over which integration 

178 will occur. 

179 'level' 

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

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

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

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

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

185 'ndin' 

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

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

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

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

190 corresponding to the subintervals for which this happens. 

191 

192 **Weighting the integrand** 

193 

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

195 integrand by a select list of functions. Different integration 

196 methods are used to compute the integral with these weighting 

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

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

199 

200 ========== =================================== ===================== 

201 ``weight`` Weight function used ``wvar`` 

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

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

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

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

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

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

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

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

210 ========== =================================== ===================== 

211 

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

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

214 

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

216 available. 

217 

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

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

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

221 

222 'momcom' 

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

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

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

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

227 'nnlog' 

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

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

230 array is equal to l if the corresponding subinterval is 

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

232 'chebmo' 

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

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

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

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

237 the first element. 

238 

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

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

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

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

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

244 output information dictionary contains the following entries instead of 

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

246 

247 'lst' 

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

249 'rslst' 

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

251 contain the integral contribution over the interval 

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

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

254 'erlst' 

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

256 corresponding to the interval in the same position in 

257 ``infodict['rslist']``. 

258 'ierlst' 

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

260 corresponding to the interval in the same position in 

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

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

263 

264 

265 **Details of QUADPACK level routines** 

266 

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

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

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

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

271 

272 ================ ============== ========== ===================== 

273 QUADPACK routine `weight` `points` infinite bounds 

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

275 qagse None No No 

276 qagie None No Yes 

277 qagpe None Yes No 

278 qawoe 'sin', 'cos' No No 

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

280 qawse 'alg*' No No 

281 qawce 'cauchy' No No 

282 ================ ============== ========== ===================== 

283 

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

285 routine. 

286 

287 qagse 

288 is an integrator based on globally adaptive interval 

289 subdivision in connection with extrapolation, which will 

290 eliminate the effects of integrand singularities of 

291 several types. 

292 qagie 

293 handles integration over infinite intervals. The infinite range is 

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

295 in ``QAGS`` is applied. 

296 qagpe 

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

298 user to provide explicit information about the location 

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

300 singularities, discontinuities and other difficulties of 

301 the integrand function. 

302 qawoe 

303 is an integrator for the evaluation of 

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

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

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

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

308 on the modified Clenshaw-Curtis technique 

309 

310 An adaptive subdivision scheme is used in connection 

311 with an extrapolation procedure, which is a modification 

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

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

314 qawfe 

315 calculates the Fourier transform 

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

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

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

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

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

321 to the series of integral approximations. 

322 qawse 

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

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

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

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

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

328 

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

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

331 applied, with modified Clenshaw-Curtis integration on those 

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

333 qawce 

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

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

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

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

338 point :math:`x = c`. 

339 

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

341 

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

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

344 written as 

345 

346 .. math:: 

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

348 

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

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

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

352 separately. 

353 

354 

355 References 

356 ---------- 

357 

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

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

360 QUADPACK: A subroutine package for automatic integration. 

361 Springer-Verlag. 

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

363 

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

365 Foundations of Analysis in the Complex Plane. 

366 Holt Rinehart Winston. 

367 ISBN 0-03-086370-8 

368 

369 Examples 

370 -------- 

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

372 

373 >>> from scipy import integrate 

374 >>> import numpy as np 

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

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

377 (21.333333333333332, 2.3684757858670003e-13) 

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

379 21.3333333333 

380 

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

382 

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

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

385 (1.0, 5.842605999138044e-11) 

386 

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

388 

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

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

391 >>> y 

392 0.5 

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

394 >>> y 

395 1.5 

396 

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

398 y parameter as 1:: 

399 

400 testlib.c => 

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

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

403 compile to library testlib.* 

404 

405 :: 

406 

407 from scipy import integrate 

408 import ctypes 

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

410 lib.func.restype = ctypes.c_double 

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

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

413 #(1.3333333333333333, 1.4802973661668752e-14) 

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

415 # 1.3333333333333333 

416 

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

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

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

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

421 bounds. 

422 

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

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

425 (1.0, 1.1102230246251565e-14) 

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

427 (1.0000000002199108, 1.0189464580163188e-08) 

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

429 (0.0, 0.0) 

430 

431 """ 

432 if not isinstance(args, tuple): 

433 args = (args,) 

434 

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

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

437 

438 if complex_func: 

439 def imfunc(x, *args): 

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

441 

442 def refunc(x, *args): 

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

444 

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

446 epsrel, limit, points, weight, wvar, wopts, 

447 maxp1, limlst, complex_func=False) 

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

449 epsrel, limit, points, weight, wvar, wopts, 

450 maxp1, limlst, complex_func=False) 

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

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

453 retval = integral, error_estimate 

454 if full_output: 

455 msgexp = {} 

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

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

458 retval = retval + (msgexp,) 

459 

460 return retval 

461 

462 if weight is None: 

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

464 points) 

465 else: 

466 if points is not None: 

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

468 "Continuing, ignoring specified points.") 

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

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

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

472 

473 if flip: 

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

475 

476 ier = retval[-1] 

477 if ier == 0: 

478 return retval[:-1] 

479 

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

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

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

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

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

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

486 6: "The input is invalid.", 

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

488 'unknown': "Unknown error."} 

489 

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

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

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

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

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

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

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

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

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

499 

500 try: 

501 msg = msgs[ier] 

502 except KeyError: 

503 msg = msgs['unknown'] 

504 

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

506 if full_output: 

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

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

509 else: 

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

511 else: 

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

513 return retval[:-1] 

514 

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

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

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

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

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

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

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

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

523 

524 elif weight is None: 

525 if points is None: # QAGSE/QAGIE 

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

527 " at least one subinterval") 

528 else: # QAGPE 

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

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

531 " integration limits.") 

532 elif len(points) >= limit: 

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

534 " must be less than subinterval" 

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

536 

537 else: 

538 if maxp1 < 1: 

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

540 

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

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

543 

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

545 if min(wvar) < -1: 

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

547 if b < a: 

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

549 

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

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

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

553 

554 raise ValueError(msg) 

555 

556 

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

558 infbounds = 0 

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

560 pass # standard integration 

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

562 infbounds = 1 

563 bound = a 

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

565 infbounds = 2 

566 bound = 0 # ignored 

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

568 infbounds = -1 

569 bound = b 

570 else: 

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

572 

573 if points is None: 

574 if infbounds == 0: 

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

576 else: 

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

578 else: 

579 if infbounds != 0: 

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

581 else: 

582 #Duplicates force function evaluation at singular points 

583 the_points = np.unique(points) 

584 the_points = the_points[a < the_points] 

585 the_points = the_points[the_points < b] 

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

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

588 

589 

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

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

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

593 

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

595 

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

597 integr = strdict[weight] 

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

599 if wopts is None: # no precomputed Chebyshev moments 

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

601 epsabs, epsrel, limit, maxp1,1) 

602 else: # precomputed Chebyshev moments 

603 momcom = wopts[0] 

604 chebcom = wopts[1] 

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

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

607 

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

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

610 epsabs,limlst,limit,maxp1) 

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

612 if weight == 'cos': 

613 def thefunc(x,*myargs): 

614 y = -x 

615 func = myargs[0] 

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

617 return func(*myargs) 

618 else: 

619 def thefunc(x,*myargs): 

620 y = -x 

621 func = myargs[0] 

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

623 return -func(*myargs) 

624 args = (func,) + args 

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

626 full_output, epsabs, limlst, limit, maxp1) 

627 else: 

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

629 else: 

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

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

632 

633 if weight.startswith('alg'): 

634 integr = strdict[weight] 

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

636 full_output, epsabs, epsrel, limit) 

637 else: # weight == 'cauchy' 

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

639 epsabs, epsrel, limit) 

640 

641 

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

643 """ 

644 Compute a double integral. 

645 

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

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

648 

649 Parameters 

650 ---------- 

651 func : callable 

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

653 first argument and x the second argument. 

654 a, b : float 

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

656 gfun : callable or float 

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

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

659 or a float indicating a constant boundary curve. 

660 hfun : callable or float 

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

662 args : sequence, optional 

663 Extra arguments to pass to `func`. 

664 epsabs : float, optional 

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

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

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

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

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

670 See `epsrel` below. 

671 epsrel : float, optional 

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

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

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

675 

676 Returns 

677 ------- 

678 y : float 

679 The resultant integral. 

680 abserr : float 

681 An estimate of the error. 

682 

683 See Also 

684 -------- 

685 quad : single integral 

686 tplquad : triple integral 

687 nquad : N-dimensional integrals 

688 fixed_quad : fixed-order Gaussian quadrature 

689 quadrature : adaptive Gaussian quadrature 

690 odeint : ODE integrator 

691 ode : ODE integrator 

692 simpson : integrator for sampled data 

693 romb : integrator for sampled data 

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

695 

696 

697 Notes 

698 ----- 

699 

700 **Details of QUADPACK level routines** 

701 

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

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

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

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

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

707 routine. 

708 

709 qagse 

710 is an integrator based on globally adaptive interval 

711 subdivision in connection with extrapolation, which will 

712 eliminate the effects of integrand singularities of 

713 several types. 

714 qagie 

715 handles integration over infinite intervals. The infinite range is 

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

717 in ``QAGS`` is applied. 

718 

719 References 

720 ---------- 

721 

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

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

724 QUADPACK: A subroutine package for automatic integration. 

725 Springer-Verlag. 

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

727 

728 Examples 

729 -------- 

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

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

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

733 

734 >>> import numpy as np 

735 >>> from scipy import integrate 

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

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

738 (0.6666666666666667, 7.401486830834377e-15) 

739 

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

741 \\,dy \\,dx`. 

742 

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

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

745 (0.41421356237309503, 1.1083280054755938e-14) 

746 

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

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

749 

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

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

752 (0.33333333333333337, 5.551115123125783e-15) 

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

754 (0.9999999999999999, 1.6653345369377348e-14) 

755 

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

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

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

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

760 

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

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

763 (3.141592653589777, 2.5173086737433208e-08) 

764 

765 """ 

766 

767 def temp_ranges(*args): 

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

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

770 

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

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

773 

774 

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

776 epsrel=1.49e-8): 

777 """ 

778 Compute a triple (definite) integral. 

779 

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

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

782 

783 Parameters 

784 ---------- 

785 func : function 

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

787 order (z, y, x). 

788 a, b : float 

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

790 gfun : function or float 

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

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

793 or a float indicating a constant boundary curve. 

794 hfun : function or float 

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

796 qfun : function or float 

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

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

799 indicating a constant boundary surface. 

800 rfun : function or float 

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

802 args : tuple, optional 

803 Extra arguments to pass to `func`. 

804 epsabs : float, optional 

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

806 integration. Default is 1.49e-8. 

807 epsrel : float, optional 

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

809 

810 Returns 

811 ------- 

812 y : float 

813 The resultant integral. 

814 abserr : float 

815 An estimate of the error. 

816 

817 See Also 

818 -------- 

819 quad : Adaptive quadrature using QUADPACK 

820 quadrature : Adaptive Gaussian quadrature 

821 fixed_quad : Fixed-order Gaussian quadrature 

822 dblquad : Double integrals 

823 nquad : N-dimensional integrals 

824 romb : Integrators for sampled data 

825 simpson : Integrators for sampled data 

826 ode : ODE integrators 

827 odeint : ODE integrators 

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

829 

830 Notes 

831 ----- 

832 

833 **Details of QUADPACK level routines** 

834 

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

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

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

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

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

840 routine. 

841 

842 qagse 

843 is an integrator based on globally adaptive interval 

844 subdivision in connection with extrapolation, which will 

845 eliminate the effects of integrand singularities of 

846 several types. 

847 qagie 

848 handles integration over infinite intervals. The infinite range is 

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

850 in ``QAGS`` is applied. 

851 

852 References 

853 ---------- 

854 

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

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

857 QUADPACK: A subroutine package for automatic integration. 

858 Springer-Verlag. 

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

860 

861 Examples 

862 -------- 

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

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

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

866 \\,dz \\,dy \\,dx`. 

867 

868 >>> import numpy as np 

869 >>> from scipy import integrate 

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

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

872 (1.8749999999999998, 3.3246447942574074e-14) 

873 

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

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

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

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

878 

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

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

881 (0.05416666666666668, 2.1774196738157757e-14) 

882 

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

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

885 

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

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

888 (0.125, 5.527033708952211e-15) 

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

890 (0.375, 1.6581101126856635e-14) 

891 

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

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

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

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

896 \\,dy\\,dx`. 

897 

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

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

900 (5.568327996830833, 4.4619078828029765e-08) 

901 

902 """ 

903 # f(z, y, x) 

904 # qfun/rfun(x, y) 

905 # gfun/hfun(x) 

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

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

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

909 

910 def ranges0(*args): 

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

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

913 

914 def ranges1(*args): 

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

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

917 

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

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

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

921 

922 

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

924 r""" 

925 Integration over multiple variables. 

926 

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

928 Various options allow improved integration of discontinuous functions, as 

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

930 integration process. 

931 

932 Parameters 

933 ---------- 

934 func : {callable, scipy.LowLevelCallable} 

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

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

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

938 passed in args. 

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

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

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

942 

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

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

945 

946 double func(int n, double *xx) 

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

948 

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

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

951 contained in the `scipy.LowLevelCallable`. 

952 ranges : iterable object 

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

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

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

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

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

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

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

960 args : iterable object, optional 

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

962 and ``opts``. 

963 opts : iterable object or dict, optional 

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

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

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

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

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

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

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

971 available options together with their default values are: 

972 

973 - epsabs = 1.49e-08 

974 - epsrel = 1.49e-08 

975 - limit = 50 

976 - points = None 

977 - weight = None 

978 - wvar = None 

979 - wopts = None 

980 

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

982 

983 full_output : bool, optional 

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

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

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

987 

988 Returns 

989 ------- 

990 result : float 

991 The result of the integration. 

992 abserr : float 

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

994 integration results. 

995 out_dict : dict, optional 

996 A dict containing additional information on the integration. 

997 

998 See Also 

999 -------- 

1000 quad : 1-D numerical integration 

1001 dblquad, tplquad : double and triple integrals 

1002 fixed_quad : fixed-order Gaussian quadrature 

1003 quadrature : adaptive Gaussian quadrature 

1004 

1005 Notes 

1006 ----- 

1007 

1008 **Details of QUADPACK level routines** 

1009 

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

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

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

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

1014 

1015 ================ ============== ========== ===================== 

1016 QUADPACK routine `weight` `points` infinite bounds 

1017 ================ ============== ========== ===================== 

1018 qagse None No No 

1019 qagie None No Yes 

1020 qagpe None Yes No 

1021 qawoe 'sin', 'cos' No No 

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

1023 qawse 'alg*' No No 

1024 qawce 'cauchy' No No 

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

1026 

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

1028 routine. 

1029 

1030 qagse 

1031 is an integrator based on globally adaptive interval 

1032 subdivision in connection with extrapolation, which will 

1033 eliminate the effects of integrand singularities of 

1034 several types. 

1035 qagie 

1036 handles integration over infinite intervals. The infinite range is 

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

1038 in ``QAGS`` is applied. 

1039 qagpe 

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

1041 user to provide explicit information about the location 

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

1043 singularities, discontinuities and other difficulties of 

1044 the integrand function. 

1045 qawoe 

1046 is an integrator for the evaluation of 

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

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

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

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

1051 on the modified Clenshaw-Curtis technique 

1052 

1053 An adaptive subdivision scheme is used in connection 

1054 with an extrapolation procedure, which is a modification 

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

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

1057 qawfe 

1058 calculates the Fourier transform 

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

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

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

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

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

1064 to the series of integral approximations. 

1065 qawse 

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

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

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

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

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

1071 

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

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

1074 applied, with modified Clenshaw-Curtis integration on those 

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

1076 qawce 

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

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

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

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

1081 point :math:`x = c`. 

1082 

1083 References 

1084 ---------- 

1085 

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

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