Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# Author: Travis Oliphant 2001 

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

3import sys 

4import warnings 

5from functools import partial 

6 

7from . import _quadpack 

8import numpy 

9from numpy import Inf 

10 

11__all__ = ['quad', 'dblquad', 'tplquad', 'nquad', 'quad_explain', 

12 'IntegrationWarning'] 

13 

14 

15error = _quadpack.error 

16 

17class IntegrationWarning(UserWarning): 

18 """ 

19 Warning on issues during integration. 

20 """ 

21 pass 

22 

23 

24def quad_explain(output=sys.stdout): 

25 """ 

26 Print extra information about integrate.quad() parameters and returns. 

27 

28 Parameters 

29 ---------- 

30 output : instance with "write" method, optional 

31 Information about `quad` is passed to ``output.write()``. 

32 Default is ``sys.stdout``. 

33 

34 Returns 

35 ------- 

36 None 

37 

38 Examples 

39 -------- 

40 We can show detailed information of the `integrate.quad` function in stdout: 

41 

42 >>> from scipy.integrate import quad_explain 

43 >>> quad_explain() 

44 

45 """ 

46 output.write(quad.__doc__) 

47 

48 

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

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

51 limlst=50): 

52 """ 

53 Compute a definite integral. 

54 

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

56 technique from the Fortran library QUADPACK. 

57 

58 Parameters 

59 ---------- 

60 func : {function, scipy.LowLevelCallable} 

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

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

63 first argument. 

64 

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

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

67 

68 double func(double x) 

69 double func(double x, void *user_data) 

70 double func(int n, double *xx) 

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

72 

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

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

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

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

77 

78 In addition, certain ctypes call signatures are supported for 

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

80 a : float 

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

82 b : float 

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

84 args : tuple, optional 

85 Extra arguments to pass to `func`. 

86 full_output : int, optional 

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

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

89 message is appended to the output tuple. 

90 

91 Returns 

92 ------- 

93 y : float 

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

95 abserr : float 

96 An estimate of the absolute error in the result. 

97 infodict : dict 

98 A dictionary containing additional information. 

99 Run scipy.integrate.quad_explain() for more information. 

100 message 

101 A convergence message. 

102 explain 

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

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

105 infodict['ierlst'] 

106 

107 Other Parameters 

108 ---------------- 

109 epsabs : float or int, optional 

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

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

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

113 numerical approximation. See `epsrel` below. 

114 epsrel : float or int, optional 

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

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

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

118 limit : float or int, optional 

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

120 algorithm. 

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

122 A sequence of break points in the bounded integration interval 

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

124 singularities, discontinuities). The sequence does not have 

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

126 with ``weight``. 

127 weight : float or int, optional 

128 String indicating weighting function. Full explanation for this 

129 and the remaining arguments can be found below. 

130 wvar : optional 

131 Variables for use with weighting functions. 

132 wopts : optional 

133 Optional input for reusing Chebyshev moments. 

134 maxp1 : float or int, optional 

135 An upper bound on the number of Chebyshev moments. 

136 limlst : int, optional 

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

138 weighting and an infinite end-point. 

139 

140 See Also 

141 -------- 

142 dblquad : double integral 

143 tplquad : triple integral 

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

145 fixed_quad : fixed-order Gaussian quadrature 

146 quadrature : adaptive Gaussian quadrature 

147 odeint : ODE integrator 

148 ode : ODE integrator 

149 simps : integrator for sampled data 

150 romb : integrator for sampled data 

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

152 

153 Notes 

154 ----- 

155 

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

157 

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

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

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

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

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

163 The entries are: 

164 

165 'neval' 

166 The number of function evaluations. 

167 'last' 

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

169 'alist' 

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

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

172 integration range. 

173 'blist' 

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

175 right end points of the subintervals. 

176 'rlist' 

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

178 integral approximations on the subintervals. 

179 'elist' 

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

181 moduli of the absolute error estimates on the subintervals. 

182 'iord' 

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

184 which are pointers to the error estimates over the subintervals 

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

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

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

188 decreasing sequence. 

189 

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

191 the following additional outputs are placed in the output 

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

193 

194 'pts' 

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

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

197 This is an array giving the subintervals over which integration 

198 will occur. 

199 'level' 

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

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

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

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

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

205 'ndin' 

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

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

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

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

210 corresponding to the subintervals for which this happens. 

211 

212 **Weighting the integrand** 

213 

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

215 integrand by a select list of functions. Different integration 

216 methods are used to compute the integral with these weighting 

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

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

219 

220 ========== =================================== ===================== 

221 ``weight`` Weight function used ``wvar`` 

222 ========== =================================== ===================== 

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

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

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

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

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

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

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

230 ========== =================================== ===================== 

231 

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

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

234 

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

236 available. 

237 

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

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

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

241 

242 'momcom' 

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

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

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

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

247 'nnlog' 

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

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

250 array is equal to l if the corresponding subinterval is 

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

252 'chebmo' 

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

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

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

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

257 the first element. 

258 

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

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

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

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

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

264 output information dictionary contains the following entries instead of 

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

266 

267 'lst' 

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

269 'rslst' 

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

271 contain the integral contribution over the interval 

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

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

274 'erlst' 

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

276 corresponding to the interval in the same position in 

277 ``infodict['rslist']``. 

278 'ierlst' 

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

280 corresponding to the interval in the same position in 

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

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

283 

284 Examples 

285 -------- 

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

287 

288 >>> from scipy import integrate 

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

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

291 (21.333333333333332, 2.3684757858670003e-13) 

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

293 21.3333333333 

294 

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

296 

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

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

299 (1.0, 5.842605999138044e-11) 

300 

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

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

303 >>> y 

304 0.5 

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

306 >>> y 

307 1.5 

308 

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

310 y parameter as 1:: 

311 

312 testlib.c => 

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

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

315 compile to library testlib.* 

316 

317 :: 

318 

319 from scipy import integrate 

320 import ctypes 

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

322 lib.func.restype = ctypes.c_double 

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

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

325 #(1.3333333333333333, 1.4802973661668752e-14) 

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

327 # 1.3333333333333333 

328 

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

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

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

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

333 bounds. 

334 

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

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

337 (1.0, 1.1102230246251565e-14) 

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

339 (1.0000000002199108, 1.0189464580163188e-08) 

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

341 (0.0, 0.0) 

342 

343 """ 

344 if not isinstance(args, tuple): 

345 args = (args,) 

346 

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

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

349 

350 if weight is None: 

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

352 points) 

353 else: 

354 if points is not None: 

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

356 "Continuing, ignoring specified points.") 

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

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

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

360 

361 if flip: 

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

363 

364 ier = retval[-1] 

365 if ier == 0: 

366 return retval[:-1] 

367 

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

369 1: "The maximum number of subdivisions (%d) has been achieved.\n If increasing the limit yields no improvement it is advised to analyze \n the integrand in order to determine the difficulties. If the position of a \n local difficulty can be determined (singularity, discontinuity) one will \n probably gain from splitting up the interval and calling the integrator \n on the subranges. Perhaps a special-purpose integrator should be used." % limit, 

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

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

372 4: "The algorithm does not converge. Roundoff error is detected\n in the extrapolation table. It is assumed that the requested tolerance\n cannot be achieved, and that the returned result (if full_output = 1) is \n the best which can be obtained.", 

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

374 6: "The input is invalid.", 

375 7: "Abnormal termination of the routine. The estimates for result\n and error are less reliable. It is assumed that the requested accuracy\n has not been achieved.", 

376 'unknown': "Unknown error."} 

377 

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

379 msgs[1] = "The maximum number of cycles allowed has been achieved., e.e.\n of subintervals (a+(k-1)c, a+kc) where c = (2*int(abs(omega)+1))\n *pi/abs(omega), for k = 1, 2, ..., lst. One can allow more cycles by increasing the value of limlst. Look at info['ierlst'] with full_output=1." 

380 msgs[4] = "The extrapolation table constructed for convergence acceleration\n of the series formed by the integral contributions over the cycles, \n does not converge to within the requested accuracy. Look at \n info['ierlst'] with full_output=1." 

381 msgs[7] = "Bad integrand behavior occurs within one or more of the cycles.\n Location and type of the difficulty involved can be determined from \n the vector info['ierlist'] obtained with full_output=1." 

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

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

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

385 4: "The integral over this cycle does not converge (to within the required accuracy) due to roundoff in the extrapolation procedure invoked on this cycle. It is assumed that the result on this interval is the best which can be obtained.", 

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

387 

388 try: 

389 msg = msgs[ier] 

390 except KeyError: 

391 msg = msgs['unknown'] 

392 

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

394 if full_output: 

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

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

397 else: 

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

399 else: 

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

401 return retval[:-1] 

402 

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

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

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

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

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

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

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

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

411 

412 elif weight is None: 

413 if points is None: # QAGSE/QAGIE 

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

415 " at least one subinterval") 

416 else: # QAGPE 

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

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

419 " integration limits.") 

420 elif len(points) >= limit: 

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

422 " must be less than subinterval" 

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

424 

425 else: 

426 if maxp1 < 1: 

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

428 

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

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

431 

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

433 if min(wvar) < -1: 

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

435 if b < a: 

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

437 

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

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

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

441 

442 raise ValueError(msg) 

443 

444 

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

446 infbounds = 0 

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

448 pass # standard integration 

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

450 infbounds = 1 

451 bound = a 

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

453 infbounds = 2 

454 bound = 0 # ignored 

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

456 infbounds = -1 

457 bound = b 

458 else: 

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

460 

461 if points is None: 

462 if infbounds == 0: 

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

464 else: 

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

466 else: 

467 if infbounds != 0: 

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

469 else: 

470 #Duplicates force function evaluation at singular points 

471 the_points = numpy.unique(points) 

472 the_points = the_points[a < the_points] 

473 the_points = the_points[the_points < b] 

474 the_points = numpy.concatenate((the_points, (0., 0.))) 

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

476 

477 

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

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

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

481 

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

483 

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

485 integr = strdict[weight] 

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

487 if wopts is None: # no precomputed Chebyshev moments 

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

489 epsabs, epsrel, limit, maxp1,1) 

490 else: # precomputed Chebyshev moments 

491 momcom = wopts[0] 

492 chebcom = wopts[1] 

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

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

495 

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

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

498 epsabs,limlst,limit,maxp1) 

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

500 if weight == 'cos': 

501 def thefunc(x,*myargs): 

502 y = -x 

503 func = myargs[0] 

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

505 return func(*myargs) 

506 else: 

507 def thefunc(x,*myargs): 

508 y = -x 

509 func = myargs[0] 

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

511 return -func(*myargs) 

512 args = (func,) + args 

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

514 full_output, epsabs, limlst, limit, maxp1) 

515 else: 

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

517 else: 

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

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

520 

521 if weight.startswith('alg'): 

522 integr = strdict[weight] 

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

524 full_output, epsabs, epsrel, limit) 

525 else: # weight == 'cauchy' 

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

527 epsabs, epsrel, limit) 

528 

529 

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

531 """ 

532 Compute a double integral. 

533 

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

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

536 

537 Parameters 

538 ---------- 

539 func : callable 

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

541 first argument and x the second argument. 

542 a, b : float 

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

544 gfun : callable or float 

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

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

547 or a float indicating a constant boundary curve. 

548 hfun : callable or float 

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

550 args : sequence, optional 

551 Extra arguments to pass to `func`. 

552 epsabs : float, optional 

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

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

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

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

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

558 See `epsrel` below. 

559 epsrel : float, optional 

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

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

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

563 

564 Returns 

565 ------- 

566 y : float 

567 The resultant integral. 

568 abserr : float 

569 An estimate of the error. 

570 

571 See also 

572 -------- 

573 quad : single integral 

574 tplquad : triple integral 

575 nquad : N-dimensional integrals 

576 fixed_quad : fixed-order Gaussian quadrature 

577 quadrature : adaptive Gaussian quadrature 

578 odeint : ODE integrator 

579 ode : ODE integrator 

580 simps : integrator for sampled data 

581 romb : integrator for sampled data 

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

583 

584 Examples 

585 -------- 

586 

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

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

589 

590 >>> from scipy import integrate 

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

592 >>> integrate.dblquad(f, 0, 2, lambda x: 0, lambda x: 1) 

593 (0.6666666666666667, 7.401486830834377e-15) 

594 

595 """ 

596 

597 def temp_ranges(*args): 

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

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

600 

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

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

603 

604 

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

606 epsrel=1.49e-8): 

607 """ 

608 Compute a triple (definite) integral. 

609 

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

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

612 

613 Parameters 

614 ---------- 

615 func : function 

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

617 order (z, y, x). 

618 a, b : float 

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

620 gfun : function or float 

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

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

623 or a float indicating a constant boundary curve. 

624 hfun : function or float 

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

626 qfun : function or float 

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

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

629 indicating a constant boundary surface. 

630 rfun : function or float 

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

632 args : tuple, optional 

633 Extra arguments to pass to `func`. 

634 epsabs : float, optional 

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

636 integration. Default is 1.49e-8. 

637 epsrel : float, optional 

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

639 

640 Returns 

641 ------- 

642 y : float 

643 The resultant integral. 

644 abserr : float 

645 An estimate of the error. 

646 

647 See Also 

648 -------- 

649 quad: Adaptive quadrature using QUADPACK 

650 quadrature: Adaptive Gaussian quadrature 

651 fixed_quad: Fixed-order Gaussian quadrature 

652 dblquad: Double integrals 

653 nquad : N-dimensional integrals 

654 romb: Integrators for sampled data 

655 simps: Integrators for sampled data 

656 ode: ODE integrators 

657 odeint: ODE integrators 

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

659 

660 Examples 

661 -------- 

662 

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

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

665 

666 >>> from scipy import integrate 

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

668 >>> integrate.tplquad(f, 1, 2, lambda x: 2, lambda x: 3, 

669 ... lambda x, y: 0, lambda x, y: 1) 

670 (1.8750000000000002, 3.324644794257407e-14) 

671 

672 

673 """ 

674 # f(z, y, x) 

675 # qfun/rfun (x, y) 

676 # gfun/hfun(x) 

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

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

679 # Stupid different API... 

680 

681 def ranges0(*args): 

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

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

684 

685 def ranges1(*args): 

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

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

688 

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

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

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

692 

693 

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

695 """ 

696 Integration over multiple variables. 

697 

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

699 Various options allow improved integration of discontinuous functions, as 

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

701 integration process. 

702 

703 Parameters 

704 ---------- 

705 func : {callable, scipy.LowLevelCallable} 

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

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

708 must be floats. Function signature should be 

709 ``func(x0, x1, ..., xn, t0, t1, ..., tm)``. Integration is carried out 

710 in order. That is, integration over ``x0`` is the innermost integral, 

711 and ``xn`` is the outermost. 

712 

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

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

715 

716 double func(int n, double *xx) 

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

718 

719 where ``n`` is the number of extra parameters and args is an array 

720 of doubles of the additional parameters, the ``xx`` array contains the 

721 coordinates. The ``user_data`` is the data contained in the 

722 `scipy.LowLevelCallable`. 

723 ranges : iterable object 

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

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

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

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

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

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

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

731 args : iterable object, optional 

732 Additional arguments ``t0, ..., tn``, required by `func`, `ranges`, and 

733 ``opts``. 

734 opts : iterable object or dict, optional 

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

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

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

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

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

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

741 the signature must be the same as for ``ranges``. The available 

742 options together with their default values are: 

743 

744 - epsabs = 1.49e-08 

745 - epsrel = 1.49e-08 

746 - limit = 50 

747 - points = None 

748 - weight = None 

749 - wvar = None 

750 - wopts = None 

751 

752 For more information on these options, see `quad` and `quad_explain`. 

753 

754 full_output : bool, optional 

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

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

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

758 

759 Returns 

760 ------- 

761 result : float 

762 The result of the integration. 

763 abserr : float 

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

765 integration results. 

766 out_dict : dict, optional 

767 A dict containing additional information on the integration. 

768 

769 See Also 

770 -------- 

771 quad : 1-D numerical integration 

772 dblquad, tplquad : double and triple integrals 

773 fixed_quad : fixed-order Gaussian quadrature 

774 quadrature : adaptive Gaussian quadrature 

775 

776 Examples 

777 -------- 

778 >>> from scipy import integrate 

779 >>> func = lambda x0,x1,x2,x3 : x0**2 + x1*x2 - x3**3 + np.sin(x0) + ( 

780 ... 1 if (x0-.2*x3-.5-.25*x1>0) else 0) 

781 >>> points = [[lambda x1,x2,x3 : 0.2*x3 + 0.5 + 0.25*x1], [], [], []] 

782 >>> def opts0(*args, **kwargs): 

783 ... return {'points':[0.2*args[2] + 0.5 + 0.25*args[0]]} 

784 >>> integrate.nquad(func, [[0,1], [-1,1], [.13,.8], [-.15,1]], 

785 ... opts=[opts0,{},{},{}], full_output=True) 

786 (1.5267454070738633, 2.9437360001402324e-14, {'neval': 388962}) 

787 

788 >>> scale = .1 

789 >>> def func2(x0, x1, x2, x3, t0, t1): 

790 ... return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0+t1*x1-t0>0 else 0) 

791 >>> def lim0(x1, x2, x3, t0, t1): 

792 ... return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1, 

793 ... scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1] 

794 >>> def lim1(x2, x3, t0, t1): 

795 ... return [scale * (t0*x2 + t1*x3) - 1, 

796 ... scale * (t0*x2 + t1*x3) + 1] 

797 >>> def lim2(x3, t0, t1): 

798 ... return [scale * (x3 + t0**2*t1**3) - 1, 

799 ... scale * (x3 + t0**2*t1**3) + 1] 

800 >>> def lim3(t0, t1): 

801 ... return [scale * (t0+t1) - 1, scale * (t0+t1) + 1] 

802 >>> def opts0(x1, x2, x3, t0, t1): 

803 ... return {'points' : [t0 - t1*x1]} 

804 >>> def opts1(x2, x3, t0, t1): 

805 ... return {} 

806 >>> def opts2(x3, t0, t1): 

807 ... return {} 

808 >>> def opts3(t0, t1): 

809 ... return {} 

810 >>> integrate.nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0), 

811 ... opts=[opts0, opts1, opts2, opts3]) 

812 (25.066666666666666, 2.7829590483937256e-13) 

813 

814 """ 

815 depth = len(ranges) 

816 ranges = [rng if callable(rng) else _RangeFunc(rng) for rng in ranges] 

817 if args is None: 

818 args = () 

819 if opts is None: 

820 opts = [dict([])] * depth 

821 

822 if isinstance(opts, dict): 

823 opts = [_OptFunc(opts)] * depth 

824 else: 

825 opts = [opt if callable(opt) else _OptFunc(opt) for opt in opts] 

826 return _NQuad(func, ranges, opts, full_output).integrate(*args) 

827 

828 

829class _RangeFunc(object): 

830 def __init__(self, range_): 

831 self.range_ = range_ 

832 

833 def __call__(self, *args): 

834 """Return stored value. 

835 

836 *args needed because range_ can be float or func, and is called with 

837 variable number of parameters. 

838 """ 

839 return self.range_ 

840 

841 

842class _OptFunc(object): 

843 def __init__(self, opt): 

844 self.opt = opt 

845 

846 def __call__(self, *args): 

847 """Return stored dict.""" 

848 return self.opt 

849 

850 

851class _NQuad(object): 

852 def __init__(self, func, ranges, opts, full_output): 

853 self.abserr = 0 

854 self.func = func 

855 self.ranges = ranges 

856 self.opts = opts 

857 self.maxdepth = len(ranges) 

858 self.full_output = full_output 

859 if self.full_output: 

860 self.out_dict = {'neval': 0} 

861 

862 def integrate(self, *args, **kwargs): 

863 depth = kwargs.pop('depth', 0) 

864 if kwargs: 

865 raise ValueError('unexpected kwargs') 

866 

867 # Get the integration range and options for this depth. 

868 ind = -(depth + 1) 

869 fn_range = self.ranges[ind] 

870 low, high = fn_range(*args) 

871 fn_opt = self.opts[ind] 

872 opt = dict(fn_opt(*args)) 

873 

874 if 'points' in opt: 

875 opt['points'] = [x for x in opt['points'] if low <= x <= high] 

876 if depth + 1 == self.maxdepth: 

877 f = self.func 

878 else: 

879 f = partial(self.integrate, depth=depth+1) 

880 quad_r = quad(f, low, high, args=args, full_output=self.full_output, 

881 **opt) 

882 value = quad_r[0] 

883 abserr = quad_r[1] 

884 if self.full_output: 

885 infodict = quad_r[2] 

886 # The 'neval' parameter in full_output returns the total 

887 # number of times the integrand function was evaluated. 

888 # Therefore, only the innermost integration loop counts. 

889 if depth + 1 == self.maxdepth: 

890 self.out_dict['neval'] += infodict['neval'] 

891 self.abserr = max(self.abserr, abserr) 

892 if depth > 0: 

893 return value 

894 else: 

895 # Final result of N-D integration with error 

896 if self.full_output: 

897 return value, self.abserr, self.out_dict 

898 else: 

899 return value, self.abserr