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

2A collection of functions to find the weights and abscissas for 

3Gaussian Quadrature. 

4 

5These calculations are done by finding the eigenvalues of a 

6tridiagonal matrix whose entries are dependent on the coefficients 

7in the recursion formula for the orthogonal polynomials with the 

8corresponding weighting function over the interval. 

9 

10Many recursion relations for orthogonal polynomials are given: 

11 

12.. math:: 

13 

14 a1n f_{n+1} (x) = (a2n + a3n x ) f_n (x) - a4n f_{n-1} (x) 

15 

16The recursion relation of interest is 

17 

18.. math:: 

19 

20 P_{n+1} (x) = (x - A_n) P_n (x) - B_n P_{n-1} (x) 

21 

22where :math:`P` has a different normalization than :math:`f`. 

23 

24The coefficients can be found as: 

25 

26.. math:: 

27 

28 A_n = -a2n / a3n 

29 \\qquad 

30 B_n = ( a4n / a3n \\sqrt{h_n-1 / h_n})^2 

31 

32where 

33 

34.. math:: 

35 

36 h_n = \\int_a^b w(x) f_n(x)^2 

37 

38assume: 

39 

40.. math:: 

41 

42 P_0 (x) = 1 

43 \\qquad 

44 P_{-1} (x) == 0 

45 

46For the mathematical background, see [golub.welsch-1969-mathcomp]_ and 

47[abramowitz.stegun-1965]_. 

48 

49References 

50---------- 

51.. [golub.welsch-1969-mathcomp] 

52 Golub, Gene H, and John H Welsch. 1969. Calculation of Gauss 

53 Quadrature Rules. *Mathematics of Computation* 23, 221-230+s1--s10. 

54 

55.. [abramowitz.stegun-1965] 

56 Abramowitz, Milton, and Irene A Stegun. (1965) *Handbook of 

57 Mathematical Functions: with Formulas, Graphs, and Mathematical 

58 Tables*. Gaithersburg, MD: National Bureau of Standards. 

59 http://www.math.sfu.ca/~cbm/aands/ 

60 

61.. [townsend.trogdon.olver-2014] 

62 Townsend, A. and Trogdon, T. and Olver, S. (2014) 

63 *Fast computation of Gauss quadrature nodes and 

64 weights on the whole real line*. :arXiv:`1410.5286`. 

65 

66.. [townsend.trogdon.olver-2015] 

67 Townsend, A. and Trogdon, T. and Olver, S. (2015) 

68 *Fast computation of Gauss quadrature nodes and 

69 weights on the whole real line*. 

70 IMA Journal of Numerical Analysis 

71 :doi:`10.1093/imanum/drv002`. 

72""" 

73# 

74# Author: Travis Oliphant 2000 

75# Updated Sep. 2003 (fixed bugs --- tested to be accurate) 

76 

77# SciPy imports. 

78import numpy as np 

79from numpy import (exp, inf, pi, sqrt, floor, sin, cos, around, int, 

80 hstack, arccos, arange) 

81from scipy import linalg 

82from scipy.special import airy 

83 

84# Local imports. 

85from . import _ufuncs 

86from . import _ufuncs as cephes 

87_gam = cephes.gamma 

88from . import specfun 

89 

90_polyfuns = ['legendre', 'chebyt', 'chebyu', 'chebyc', 'chebys', 

91 'jacobi', 'laguerre', 'genlaguerre', 'hermite', 

92 'hermitenorm', 'gegenbauer', 'sh_legendre', 'sh_chebyt', 

93 'sh_chebyu', 'sh_jacobi'] 

94 

95# Correspondence between new and old names of root functions 

96_rootfuns_map = {'roots_legendre': 'p_roots', 

97 'roots_chebyt': 't_roots', 

98 'roots_chebyu': 'u_roots', 

99 'roots_chebyc': 'c_roots', 

100 'roots_chebys': 's_roots', 

101 'roots_jacobi': 'j_roots', 

102 'roots_laguerre': 'l_roots', 

103 'roots_genlaguerre': 'la_roots', 

104 'roots_hermite': 'h_roots', 

105 'roots_hermitenorm': 'he_roots', 

106 'roots_gegenbauer': 'cg_roots', 

107 'roots_sh_legendre': 'ps_roots', 

108 'roots_sh_chebyt': 'ts_roots', 

109 'roots_sh_chebyu': 'us_roots', 

110 'roots_sh_jacobi': 'js_roots'} 

111 

112_evalfuns = ['eval_legendre', 'eval_chebyt', 'eval_chebyu', 

113 'eval_chebyc', 'eval_chebys', 'eval_jacobi', 

114 'eval_laguerre', 'eval_genlaguerre', 'eval_hermite', 

115 'eval_hermitenorm', 'eval_gegenbauer', 

116 'eval_sh_legendre', 'eval_sh_chebyt', 'eval_sh_chebyu', 

117 'eval_sh_jacobi'] 

118 

119__all__ = _polyfuns + list(_rootfuns_map.keys()) + _evalfuns + ['poch', 'binom'] 

120 

121 

122class orthopoly1d(np.poly1d): 

123 

124 def __init__(self, roots, weights=None, hn=1.0, kn=1.0, wfunc=None, 

125 limits=None, monic=False, eval_func=None): 

126 equiv_weights = [weights[k] / wfunc(roots[k]) for 

127 k in range(len(roots))] 

128 mu = sqrt(hn) 

129 if monic: 

130 evf = eval_func 

131 if evf: 

132 knn = kn 

133 eval_func = lambda x: evf(x) / knn 

134 mu = mu / abs(kn) 

135 kn = 1.0 

136 

137 # compute coefficients from roots, then scale 

138 poly = np.poly1d(roots, r=True) 

139 np.poly1d.__init__(self, poly.coeffs * float(kn)) 

140 

141 self.weights = np.array(list(zip(roots, weights, equiv_weights))) 

142 self.weight_func = wfunc 

143 self.limits = limits 

144 self.normcoef = mu 

145 

146 # Note: eval_func will be discarded on arithmetic 

147 self._eval_func = eval_func 

148 

149 def __call__(self, v): 

150 if self._eval_func and not isinstance(v, np.poly1d): 

151 return self._eval_func(v) 

152 else: 

153 return np.poly1d.__call__(self, v) 

154 

155 def _scale(self, p): 

156 if p == 1.0: 

157 return 

158 self._coeffs *= p 

159 

160 evf = self._eval_func 

161 if evf: 

162 self._eval_func = lambda x: evf(x) * p 

163 self.normcoef *= p 

164 

165 

166def _gen_roots_and_weights(n, mu0, an_func, bn_func, f, df, symmetrize, mu): 

167 """[x,w] = gen_roots_and_weights(n,an_func,sqrt_bn_func,mu) 

168 

169 Returns the roots (x) of an nth order orthogonal polynomial, 

170 and weights (w) to use in appropriate Gaussian quadrature with that 

171 orthogonal polynomial. 

172 

173 The polynomials have the recurrence relation 

174 P_n+1(x) = (x - A_n) P_n(x) - B_n P_n-1(x) 

175 

176 an_func(n) should return A_n 

177 sqrt_bn_func(n) should return sqrt(B_n) 

178 mu ( = h_0 ) is the integral of the weight over the orthogonal 

179 interval 

180 """ 

181 k = np.arange(n, dtype='d') 

182 c = np.zeros((2, n)) 

183 c[0,1:] = bn_func(k[1:]) 

184 c[1,:] = an_func(k) 

185 x = linalg.eigvals_banded(c, overwrite_a_band=True) 

186 

187 # improve roots by one application of Newton's method 

188 y = f(n, x) 

189 dy = df(n, x) 

190 x -= y/dy 

191 

192 fm = f(n-1, x) 

193 fm /= np.abs(fm).max() 

194 dy /= np.abs(dy).max() 

195 w = 1.0 / (fm * dy) 

196 

197 if symmetrize: 

198 w = (w + w[::-1]) / 2 

199 x = (x - x[::-1]) / 2 

200 

201 w *= mu0 / w.sum() 

202 

203 if mu: 

204 return x, w, mu0 

205 else: 

206 return x, w 

207 

208# Jacobi Polynomials 1 P^(alpha,beta)_n(x) 

209 

210 

211def roots_jacobi(n, alpha, beta, mu=False): 

212 r"""Gauss-Jacobi quadrature. 

213 

214 Compute the sample points and weights for Gauss-Jacobi 

215 quadrature. The sample points are the roots of the nth degree 

216 Jacobi polynomial, :math:`P^{\alpha, \beta}_n(x)`. These sample 

217 points and weights correctly integrate polynomials of degree 

218 :math:`2n - 1` or less over the interval :math:`[-1, 1]` with 

219 weight function :math:`w(x) = (1 - x)^{\alpha} (1 + 

220 x)^{\beta}`. See 22.2.1 in [AS]_ for details. 

221 

222 Parameters 

223 ---------- 

224 n : int 

225 quadrature order 

226 alpha : float 

227 alpha must be > -1 

228 beta : float 

229 beta must be > -1 

230 mu : bool, optional 

231 If True, return the sum of the weights, optional. 

232 

233 Returns 

234 ------- 

235 x : ndarray 

236 Sample points 

237 w : ndarray 

238 Weights 

239 mu : float 

240 Sum of the weights 

241 

242 See Also 

243 -------- 

244 scipy.integrate.quadrature 

245 scipy.integrate.fixed_quad 

246 

247 References 

248 ---------- 

249 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

250 Handbook of Mathematical Functions with Formulas, 

251 Graphs, and Mathematical Tables. New York: Dover, 1972. 

252 

253 """ 

254 m = int(n) 

255 if n < 1 or n != m: 

256 raise ValueError("n must be a positive integer.") 

257 if alpha <= -1 or beta <= -1: 

258 raise ValueError("alpha and beta must be greater than -1.") 

259 

260 if alpha == 0.0 and beta == 0.0: 

261 return roots_legendre(m, mu) 

262 if alpha == beta: 

263 return roots_gegenbauer(m, alpha+0.5, mu) 

264 

265 mu0 = 2.0**(alpha+beta+1)*cephes.beta(alpha+1, beta+1) 

266 a = alpha 

267 b = beta 

268 if a + b == 0.0: 

269 an_func = lambda k: np.where(k == 0, (b-a)/(2+a+b), 0.0) 

270 else: 

271 an_func = lambda k: np.where(k == 0, (b-a)/(2+a+b), 

272 (b*b - a*a) / ((2.0*k+a+b)*(2.0*k+a+b+2))) 

273 

274 bn_func = lambda k: 2.0 / (2.0*k+a+b)*np.sqrt((k+a)*(k+b) / (2*k+a+b+1)) \ 

275 * np.where(k == 1, 1.0, np.sqrt(k*(k+a+b) / (2.0*k+a+b-1))) 

276 

277 f = lambda n, x: cephes.eval_jacobi(n, a, b, x) 

278 df = lambda n, x: 0.5 * (n + a + b + 1) \ 

279 * cephes.eval_jacobi(n-1, a+1, b+1, x) 

280 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, False, mu) 

281 

282 

283def jacobi(n, alpha, beta, monic=False): 

284 r"""Jacobi polynomial. 

285 

286 Defined to be the solution of 

287 

288 .. math:: 

289 (1 - x^2)\frac{d^2}{dx^2}P_n^{(\alpha, \beta)} 

290 + (\beta - \alpha - (\alpha + \beta + 2)x) 

291 \frac{d}{dx}P_n^{(\alpha, \beta)} 

292 + n(n + \alpha + \beta + 1)P_n^{(\alpha, \beta)} = 0 

293 

294 for :math:`\alpha, \beta > -1`; :math:`P_n^{(\alpha, \beta)}` is a 

295 polynomial of degree :math:`n`. 

296 

297 Parameters 

298 ---------- 

299 n : int 

300 Degree of the polynomial. 

301 alpha : float 

302 Parameter, must be greater than -1. 

303 beta : float 

304 Parameter, must be greater than -1. 

305 monic : bool, optional 

306 If `True`, scale the leading coefficient to be 1. Default is 

307 `False`. 

308 

309 Returns 

310 ------- 

311 P : orthopoly1d 

312 Jacobi polynomial. 

313 

314 Notes 

315 ----- 

316 For fixed :math:`\alpha, \beta`, the polynomials 

317 :math:`P_n^{(\alpha, \beta)}` are orthogonal over :math:`[-1, 1]` 

318 with weight function :math:`(1 - x)^\alpha(1 + x)^\beta`. 

319 

320 """ 

321 if n < 0: 

322 raise ValueError("n must be nonnegative.") 

323 

324 wfunc = lambda x: (1 - x)**alpha * (1 + x)**beta 

325 if n == 0: 

326 return orthopoly1d([], [], 1.0, 1.0, wfunc, (-1, 1), monic, 

327 eval_func=np.ones_like) 

328 x, w, mu = roots_jacobi(n, alpha, beta, mu=True) 

329 ab1 = alpha + beta + 1.0 

330 hn = 2**ab1 / (2 * n + ab1) * _gam(n + alpha + 1) 

331 hn *= _gam(n + beta + 1.0) / _gam(n + 1) / _gam(n + ab1) 

332 kn = _gam(2 * n + ab1) / 2.0**n / _gam(n + 1) / _gam(n + ab1) 

333 # here kn = coefficient on x^n term 

334 p = orthopoly1d(x, w, hn, kn, wfunc, (-1, 1), monic, 

335 lambda x: eval_jacobi(n, alpha, beta, x)) 

336 return p 

337 

338# Jacobi Polynomials shifted G_n(p,q,x) 

339 

340 

341def roots_sh_jacobi(n, p1, q1, mu=False): 

342 """Gauss-Jacobi (shifted) quadrature. 

343 

344 Compute the sample points and weights for Gauss-Jacobi (shifted) 

345 quadrature. The sample points are the roots of the nth degree 

346 shifted Jacobi polynomial, :math:`G^{p,q}_n(x)`. These sample 

347 points and weights correctly integrate polynomials of degree 

348 :math:`2n - 1` or less over the interval :math:`[0, 1]` with 

349 weight function :math:`w(x) = (1 - x)^{p-q} x^{q-1}`. See 22.2.2 

350 in [AS]_ for details. 

351 

352 Parameters 

353 ---------- 

354 n : int 

355 quadrature order 

356 p1 : float 

357 (p1 - q1) must be > -1 

358 q1 : float 

359 q1 must be > 0 

360 mu : bool, optional 

361 If True, return the sum of the weights, optional. 

362 

363 Returns 

364 ------- 

365 x : ndarray 

366 Sample points 

367 w : ndarray 

368 Weights 

369 mu : float 

370 Sum of the weights 

371 

372 See Also 

373 -------- 

374 scipy.integrate.quadrature 

375 scipy.integrate.fixed_quad 

376 

377 References 

378 ---------- 

379 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

380 Handbook of Mathematical Functions with Formulas, 

381 Graphs, and Mathematical Tables. New York: Dover, 1972. 

382 

383 """ 

384 if (p1-q1) <= -1 or q1 <= 0: 

385 raise ValueError("(p - q) must be greater than -1, and q must be greater than 0.") 

386 x, w, m = roots_jacobi(n, p1-q1, q1-1, True) 

387 x = (x + 1) / 2 

388 scale = 2.0**p1 

389 w /= scale 

390 m /= scale 

391 if mu: 

392 return x, w, m 

393 else: 

394 return x, w 

395 

396def sh_jacobi(n, p, q, monic=False): 

397 r"""Shifted Jacobi polynomial. 

398 

399 Defined by 

400 

401 .. math:: 

402 

403 G_n^{(p, q)}(x) 

404 = \binom{2n + p - 1}{n}^{-1}P_n^{(p - q, q - 1)}(2x - 1), 

405 

406 where :math:`P_n^{(\cdot, \cdot)}` is the nth Jacobi polynomial. 

407 

408 Parameters 

409 ---------- 

410 n : int 

411 Degree of the polynomial. 

412 p : float 

413 Parameter, must have :math:`p > q - 1`. 

414 q : float 

415 Parameter, must be greater than 0. 

416 monic : bool, optional 

417 If `True`, scale the leading coefficient to be 1. Default is 

418 `False`. 

419 

420 Returns 

421 ------- 

422 G : orthopoly1d 

423 Shifted Jacobi polynomial. 

424 

425 Notes 

426 ----- 

427 For fixed :math:`p, q`, the polynomials :math:`G_n^{(p, q)}` are 

428 orthogonal over :math:`[0, 1]` with weight function :math:`(1 - 

429 x)^{p - q}x^{q - 1}`. 

430 

431 """ 

432 if n < 0: 

433 raise ValueError("n must be nonnegative.") 

434 

435 wfunc = lambda x: (1.0 - x)**(p - q) * (x)**(q - 1.) 

436 if n == 0: 

437 return orthopoly1d([], [], 1.0, 1.0, wfunc, (-1, 1), monic, 

438 eval_func=np.ones_like) 

439 n1 = n 

440 x, w, mu0 = roots_sh_jacobi(n1, p, q, mu=True) 

441 hn = _gam(n + 1) * _gam(n + q) * _gam(n + p) * _gam(n + p - q + 1) 

442 hn /= (2 * n + p) * (_gam(2 * n + p)**2) 

443 # kn = 1.0 in standard form so monic is redundant. Kept for compatibility. 

444 kn = 1.0 

445 pp = orthopoly1d(x, w, hn, kn, wfunc=wfunc, limits=(0, 1), monic=monic, 

446 eval_func=lambda x: eval_sh_jacobi(n, p, q, x)) 

447 return pp 

448 

449# Generalized Laguerre L^(alpha)_n(x) 

450 

451 

452def roots_genlaguerre(n, alpha, mu=False): 

453 r"""Gauss-generalized Laguerre quadrature. 

454 

455 Compute the sample points and weights for Gauss-generalized 

456 Laguerre quadrature. The sample points are the roots of the nth 

457 degree generalized Laguerre polynomial, :math:`L^{\alpha}_n(x)`. 

458 These sample points and weights correctly integrate polynomials of 

459 degree :math:`2n - 1` or less over the interval :math:`[0, 

460 \infty]` with weight function :math:`w(x) = x^{\alpha} 

461 e^{-x}`. See 22.3.9 in [AS]_ for details. 

462 

463 Parameters 

464 ---------- 

465 n : int 

466 quadrature order 

467 alpha : float 

468 alpha must be > -1 

469 mu : bool, optional 

470 If True, return the sum of the weights, optional. 

471 

472 Returns 

473 ------- 

474 x : ndarray 

475 Sample points 

476 w : ndarray 

477 Weights 

478 mu : float 

479 Sum of the weights 

480 

481 See Also 

482 -------- 

483 scipy.integrate.quadrature 

484 scipy.integrate.fixed_quad 

485 

486 References 

487 ---------- 

488 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

489 Handbook of Mathematical Functions with Formulas, 

490 Graphs, and Mathematical Tables. New York: Dover, 1972. 

491 

492 """ 

493 m = int(n) 

494 if n < 1 or n != m: 

495 raise ValueError("n must be a positive integer.") 

496 if alpha < -1: 

497 raise ValueError("alpha must be greater than -1.") 

498 

499 mu0 = cephes.gamma(alpha + 1) 

500 

501 if m == 1: 

502 x = np.array([alpha+1.0], 'd') 

503 w = np.array([mu0], 'd') 

504 if mu: 

505 return x, w, mu0 

506 else: 

507 return x, w 

508 

509 an_func = lambda k: 2 * k + alpha + 1 

510 bn_func = lambda k: -np.sqrt(k * (k + alpha)) 

511 f = lambda n, x: cephes.eval_genlaguerre(n, alpha, x) 

512 df = lambda n, x: (n*cephes.eval_genlaguerre(n, alpha, x) 

513 - (n + alpha)*cephes.eval_genlaguerre(n-1, alpha, x))/x 

514 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, False, mu) 

515 

516 

517def genlaguerre(n, alpha, monic=False): 

518 r"""Generalized (associated) Laguerre polynomial. 

519 

520 Defined to be the solution of 

521 

522 .. math:: 

523 x\frac{d^2}{dx^2}L_n^{(\alpha)} 

524 + (\alpha + 1 - x)\frac{d}{dx}L_n^{(\alpha)} 

525 + nL_n^{(\alpha)} = 0, 

526 

527 where :math:`\alpha > -1`; :math:`L_n^{(\alpha)}` is a polynomial 

528 of degree :math:`n`. 

529 

530 Parameters 

531 ---------- 

532 n : int 

533 Degree of the polynomial. 

534 alpha : float 

535 Parameter, must be greater than -1. 

536 monic : bool, optional 

537 If `True`, scale the leading coefficient to be 1. Default is 

538 `False`. 

539 

540 Returns 

541 ------- 

542 L : orthopoly1d 

543 Generalized Laguerre polynomial. 

544 

545 Notes 

546 ----- 

547 For fixed :math:`\alpha`, the polynomials :math:`L_n^{(\alpha)}` 

548 are orthogonal over :math:`[0, \infty)` with weight function 

549 :math:`e^{-x}x^\alpha`. 

550 

551 The Laguerre polynomials are the special case where :math:`\alpha 

552 = 0`. 

553 

554 See Also 

555 -------- 

556 laguerre : Laguerre polynomial. 

557 

558 """ 

559 if alpha <= -1: 

560 raise ValueError("alpha must be > -1") 

561 if n < 0: 

562 raise ValueError("n must be nonnegative.") 

563 

564 if n == 0: 

565 n1 = n + 1 

566 else: 

567 n1 = n 

568 x, w, mu0 = roots_genlaguerre(n1, alpha, mu=True) 

569 wfunc = lambda x: exp(-x) * x**alpha 

570 if n == 0: 

571 x, w = [], [] 

572 hn = _gam(n + alpha + 1) / _gam(n + 1) 

573 kn = (-1)**n / _gam(n + 1) 

574 p = orthopoly1d(x, w, hn, kn, wfunc, (0, inf), monic, 

575 lambda x: eval_genlaguerre(n, alpha, x)) 

576 return p 

577 

578# Laguerre L_n(x) 

579 

580 

581def roots_laguerre(n, mu=False): 

582 r"""Gauss-Laguerre quadrature. 

583 

584 Compute the sample points and weights for Gauss-Laguerre 

585 quadrature. The sample points are the roots of the nth degree 

586 Laguerre polynomial, :math:`L_n(x)`. These sample points and 

587 weights correctly integrate polynomials of degree :math:`2n - 1` 

588 or less over the interval :math:`[0, \infty]` with weight function 

589 :math:`w(x) = e^{-x}`. See 22.2.13 in [AS]_ for details. 

590 

591 Parameters 

592 ---------- 

593 n : int 

594 quadrature order 

595 mu : bool, optional 

596 If True, return the sum of the weights, optional. 

597 

598 Returns 

599 ------- 

600 x : ndarray 

601 Sample points 

602 w : ndarray 

603 Weights 

604 mu : float 

605 Sum of the weights 

606 

607 See Also 

608 -------- 

609 scipy.integrate.quadrature 

610 scipy.integrate.fixed_quad 

611 numpy.polynomial.laguerre.laggauss 

612 

613 References 

614 ---------- 

615 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

616 Handbook of Mathematical Functions with Formulas, 

617 Graphs, and Mathematical Tables. New York: Dover, 1972. 

618 

619 """ 

620 return roots_genlaguerre(n, 0.0, mu=mu) 

621 

622 

623def laguerre(n, monic=False): 

624 r"""Laguerre polynomial. 

625 

626 Defined to be the solution of 

627 

628 .. math:: 

629 x\frac{d^2}{dx^2}L_n + (1 - x)\frac{d}{dx}L_n + nL_n = 0; 

630 

631 :math:`L_n` is a polynomial of degree :math:`n`. 

632 

633 Parameters 

634 ---------- 

635 n : int 

636 Degree of the polynomial. 

637 monic : bool, optional 

638 If `True`, scale the leading coefficient to be 1. Default is 

639 `False`. 

640 

641 Returns 

642 ------- 

643 L : orthopoly1d 

644 Laguerre Polynomial. 

645 

646 Notes 

647 ----- 

648 The polynomials :math:`L_n` are orthogonal over :math:`[0, 

649 \infty)` with weight function :math:`e^{-x}`. 

650 

651 """ 

652 if n < 0: 

653 raise ValueError("n must be nonnegative.") 

654 

655 if n == 0: 

656 n1 = n + 1 

657 else: 

658 n1 = n 

659 x, w, mu0 = roots_laguerre(n1, mu=True) 

660 if n == 0: 

661 x, w = [], [] 

662 hn = 1.0 

663 kn = (-1)**n / _gam(n + 1) 

664 p = orthopoly1d(x, w, hn, kn, lambda x: exp(-x), (0, inf), monic, 

665 lambda x: eval_laguerre(n, x)) 

666 return p 

667 

668# Hermite 1 H_n(x) 

669 

670 

671def roots_hermite(n, mu=False): 

672 r"""Gauss-Hermite (physicist's) quadrature. 

673 

674 Compute the sample points and weights for Gauss-Hermite 

675 quadrature. The sample points are the roots of the nth degree 

676 Hermite polynomial, :math:`H_n(x)`. These sample points and 

677 weights correctly integrate polynomials of degree :math:`2n - 1` 

678 or less over the interval :math:`[-\infty, \infty]` with weight 

679 function :math:`w(x) = e^{-x^2}`. See 22.2.14 in [AS]_ for 

680 details. 

681 

682 Parameters 

683 ---------- 

684 n : int 

685 quadrature order 

686 mu : bool, optional 

687 If True, return the sum of the weights, optional. 

688 

689 Returns 

690 ------- 

691 x : ndarray 

692 Sample points 

693 w : ndarray 

694 Weights 

695 mu : float 

696 Sum of the weights 

697 

698 Notes 

699 ----- 

700 For small n up to 150 a modified version of the Golub-Welsch 

701 algorithm is used. Nodes are computed from the eigenvalue 

702 problem and improved by one step of a Newton iteration. 

703 The weights are computed from the well-known analytical formula. 

704 

705 For n larger than 150 an optimal asymptotic algorithm is applied 

706 which computes nodes and weights in a numerically stable manner. 

707 The algorithm has linear runtime making computation for very 

708 large n (several thousand or more) feasible. 

709 

710 See Also 

711 -------- 

712 scipy.integrate.quadrature 

713 scipy.integrate.fixed_quad 

714 numpy.polynomial.hermite.hermgauss 

715 roots_hermitenorm 

716 

717 References 

718 ---------- 

719 .. [townsend.trogdon.olver-2014] 

720 Townsend, A. and Trogdon, T. and Olver, S. (2014) 

721 *Fast computation of Gauss quadrature nodes and 

722 weights on the whole real line*. :arXiv:`1410.5286`. 

723 .. [townsend.trogdon.olver-2015] 

724 Townsend, A. and Trogdon, T. and Olver, S. (2015) 

725 *Fast computation of Gauss quadrature nodes and 

726 weights on the whole real line*. 

727 IMA Journal of Numerical Analysis 

728 :doi:`10.1093/imanum/drv002`. 

729 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

730 Handbook of Mathematical Functions with Formulas, 

731 Graphs, and Mathematical Tables. New York: Dover, 1972. 

732 

733 """ 

734 m = int(n) 

735 if n < 1 or n != m: 

736 raise ValueError("n must be a positive integer.") 

737 

738 mu0 = np.sqrt(np.pi) 

739 if n <= 150: 

740 an_func = lambda k: 0.0*k 

741 bn_func = lambda k: np.sqrt(k/2.0) 

742 f = cephes.eval_hermite 

743 df = lambda n, x: 2.0 * n * cephes.eval_hermite(n-1, x) 

744 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu) 

745 else: 

746 nodes, weights = _roots_hermite_asy(m) 

747 if mu: 

748 return nodes, weights, mu0 

749 else: 

750 return nodes, weights 

751 

752 

753def _compute_tauk(n, k, maxit=5): 

754 """Helper function for Tricomi initial guesses 

755 

756 For details, see formula 3.1 in lemma 3.1 in the 

757 original paper. 

758 

759 Parameters 

760 ---------- 

761 n : int 

762 Quadrature order 

763 k : ndarray of type int 

764 Index of roots :math:`\tau_k` to compute 

765 maxit : int 

766 Number of Newton maxit performed, the default 

767 value of 5 is sufficient. 

768 

769 Returns 

770 ------- 

771 tauk : ndarray 

772 Roots of equation 3.1 

773 

774 See Also 

775 -------- 

776 initial_nodes_a 

777 roots_hermite_asy 

778 """ 

779 a = n % 2 - 0.5 

780 c = (4.0*floor(n/2.0) - 4.0*k + 3.0)*pi / (4.0*floor(n/2.0) + 2.0*a + 2.0) 

781 f = lambda x: x - sin(x) - c 

782 df = lambda x: 1.0 - cos(x) 

783 xi = 0.5*pi 

784 for i in range(maxit): 

785 xi = xi - f(xi)/df(xi) 

786 return xi 

787 

788 

789def _initial_nodes_a(n, k): 

790 r"""Tricomi initial guesses 

791 

792 Computes an initial approximation to the square of the `k`-th 

793 (positive) root :math:`x_k` of the Hermite polynomial :math:`H_n` 

794 of order :math:`n`. The formula is the one from lemma 3.1 in the 

795 original paper. The guesses are accurate except in the region 

796 near :math:`\sqrt{2n + 1}`. 

797 

798 Parameters 

799 ---------- 

800 n : int 

801 Quadrature order 

802 k : ndarray of type int 

803 Index of roots to compute 

804 

805 Returns 

806 ------- 

807 xksq : ndarray 

808 Square of the approximate roots 

809 

810 See Also 

811 -------- 

812 initial_nodes 

813 roots_hermite_asy 

814 """ 

815 tauk = _compute_tauk(n, k) 

816 sigk = cos(0.5*tauk)**2 

817 a = n % 2 - 0.5 

818 nu = 4.0*floor(n/2.0) + 2.0*a + 2.0 

819 # Initial approximation of Hermite roots (square) 

820 xksq = nu*sigk - 1.0/(3.0*nu) * (5.0/(4.0*(1.0-sigk)**2) - 1.0/(1.0-sigk) - 0.25) 

821 return xksq 

822 

823 

824def _initial_nodes_b(n, k): 

825 r"""Gatteschi initial guesses 

826 

827 Computes an initial approximation to the square of the kth 

828 (positive) root :math:`x_k` of the Hermite polynomial :math:`H_n` 

829 of order :math:`n`. The formula is the one from lemma 3.2 in the 

830 original paper. The guesses are accurate in the region just 

831 below :math:`\sqrt{2n + 1}`. 

832 

833 Parameters 

834 ---------- 

835 n : int 

836 Quadrature order 

837 k : ndarray of type int 

838 Index of roots to compute 

839 

840 Returns 

841 ------- 

842 xksq : ndarray 

843 Square of the approximate root 

844 

845 See Also 

846 -------- 

847 initial_nodes 

848 roots_hermite_asy 

849 """ 

850 a = n % 2 - 0.5 

851 nu = 4.0*floor(n/2.0) + 2.0*a + 2.0 

852 # Airy roots by approximation 

853 ak = specfun.airyzo(k.max(), 1)[0][::-1] 

854 # Initial approximation of Hermite roots (square) 

855 xksq = (nu + 

856 2.0**(2.0/3.0) * ak * nu**(1.0/3.0) + 

857 1.0/5.0 * 2.0**(4.0/3.0) * ak**2 * nu**(-1.0/3.0) + 

858 (9.0/140.0 - 12.0/175.0 * ak**3) * nu**(-1.0) + 

859 (16.0/1575.0 * ak + 92.0/7875.0 * ak**4) * 2.0**(2.0/3.0) * nu**(-5.0/3.0) - 

860 (15152.0/3031875.0 * ak**5 + 1088.0/121275.0 * ak**2) * 2.0**(1.0/3.0) * nu**(-7.0/3.0)) 

861 return xksq 

862 

863 

864def _initial_nodes(n): 

865 """Initial guesses for the Hermite roots 

866 

867 Computes an initial approximation to the non-negative 

868 roots :math:`x_k` of the Hermite polynomial :math:`H_n` 

869 of order :math:`n`. The Tricomi and Gatteschi initial 

870 guesses are used in the region where they are accurate. 

871 

872 Parameters 

873 ---------- 

874 n : int 

875 Quadrature order 

876 

877 Returns 

878 ------- 

879 xk : ndarray 

880 Approximate roots 

881 

882 See Also 

883 -------- 

884 roots_hermite_asy 

885 """ 

886 # Turnover point 

887 # linear polynomial fit to error of 10, 25, 40, ..., 1000 point rules 

888 fit = 0.49082003*n - 4.37859653 

889 turnover = around(fit).astype(int) 

890 # Compute all approximations 

891 ia = arange(1, int(floor(n*0.5)+1)) 

892 ib = ia[::-1] 

893 xasq = _initial_nodes_a(n, ia[:turnover+1]) 

894 xbsq = _initial_nodes_b(n, ib[turnover+1:]) 

895 # Combine 

896 iv = sqrt(hstack([xasq, xbsq])) 

897 # Central node is always zero 

898 if n % 2 == 1: 

899 iv = hstack([0.0, iv]) 

900 return iv 

901 

902 

903def _pbcf(n, theta): 

904 r"""Asymptotic series expansion of parabolic cylinder function 

905 

906 The implementation is based on sections 3.2 and 3.3 from the 

907 original paper. Compared to the published version this code 

908 adds one more term to the asymptotic series. The detailed 

909 formulas can be found at [parabolic-asymptotics]_. The evaluation 

910 is done in a transformed variable :math:`\theta := \arccos(t)` 

911 where :math:`t := x / \mu` and :math:`\mu := \sqrt{2n + 1}`. 

912 

913 Parameters 

914 ---------- 

915 n : int 

916 Quadrature order 

917 theta : ndarray 

918 Transformed position variable 

919 

920 Returns 

921 ------- 

922 U : ndarray 

923 Value of the parabolic cylinder function :math:`U(a, \theta)`. 

924 Ud : ndarray 

925 Value of the derivative :math:`U^{\prime}(a, \theta)` of 

926 the parabolic cylinder function. 

927 

928 See Also 

929 -------- 

930 roots_hermite_asy 

931 

932 References 

933 ---------- 

934 .. [parabolic-asymptotics] 

935 https://dlmf.nist.gov/12.10#vii 

936 """ 

937 st = sin(theta) 

938 ct = cos(theta) 

939 # https://dlmf.nist.gov/12.10#vii 

940 mu = 2.0*n + 1.0 

941 # https://dlmf.nist.gov/12.10#E23 

942 eta = 0.5*theta - 0.5*st*ct 

943 # https://dlmf.nist.gov/12.10#E39 

944 zeta = -(3.0*eta/2.0) ** (2.0/3.0) 

945 # https://dlmf.nist.gov/12.10#E40 

946 phi = (-zeta / st**2) ** (0.25) 

947 # Coefficients 

948 # https://dlmf.nist.gov/12.10#E43 

949 a0 = 1.0 

950 a1 = 0.10416666666666666667 

951 a2 = 0.08355034722222222222 

952 a3 = 0.12822657455632716049 

953 a4 = 0.29184902646414046425 

954 a5 = 0.88162726744375765242 

955 b0 = 1.0 

956 b1 = -0.14583333333333333333 

957 b2 = -0.09874131944444444444 

958 b3 = -0.14331205391589506173 

959 b4 = -0.31722720267841354810 

960 b5 = -0.94242914795712024914 

961 # Polynomials 

962 # https://dlmf.nist.gov/12.10#E9 

963 # https://dlmf.nist.gov/12.10#E10 

964 ctp = ct ** arange(16).reshape((-1,1)) 

965 u0 = 1.0 

966 u1 = (1.0*ctp[3,:] - 6.0*ct) / 24.0 

967 u2 = (-9.0*ctp[4,:] + 249.0*ctp[2,:] + 145.0) / 1152.0 

968 u3 = (-4042.0*ctp[9,:] + 18189.0*ctp[7,:] - 28287.0*ctp[5,:] - 151995.0*ctp[3,:] - 259290.0*ct) / 414720.0 

969 u4 = (72756.0*ctp[10,:] - 321339.0*ctp[8,:] - 154982.0*ctp[6,:] + 50938215.0*ctp[4,:] + 122602962.0*ctp[2,:] + 12773113.0) / 39813120.0 

970 u5 = (82393456.0*ctp[15,:] - 617950920.0*ctp[13,:] + 1994971575.0*ctp[11,:] - 3630137104.0*ctp[9,:] + 4433574213.0*ctp[7,:] 

971 - 37370295816.0*ctp[5,:] - 119582875013.0*ctp[3,:] - 34009066266.0*ct) / 6688604160.0 

972 v0 = 1.0 

973 v1 = (1.0*ctp[3,:] + 6.0*ct) / 24.0 

974 v2 = (15.0*ctp[4,:] - 327.0*ctp[2,:] - 143.0) / 1152.0 

975 v3 = (-4042.0*ctp[9,:] + 18189.0*ctp[7,:] - 36387.0*ctp[5,:] + 238425.0*ctp[3,:] + 259290.0*ct) / 414720.0 

976 v4 = (-121260.0*ctp[10,:] + 551733.0*ctp[8,:] - 151958.0*ctp[6,:] - 57484425.0*ctp[4,:] - 132752238.0*ctp[2,:] - 12118727) / 39813120.0 

977 v5 = (82393456.0*ctp[15,:] - 617950920.0*ctp[13,:] + 2025529095.0*ctp[11,:] - 3750839308.0*ctp[9,:] + 3832454253.0*ctp[7,:] 

978 + 35213253348.0*ctp[5,:] + 130919230435.0*ctp[3,:] + 34009066266*ct) / 6688604160.0 

979 # Airy Evaluation (Bi and Bip unused) 

980 Ai, Aip, Bi, Bip = airy(mu**(4.0/6.0) * zeta) 

981 # Prefactor for U 

982 P = 2.0*sqrt(pi) * mu**(1.0/6.0) * phi 

983 # Terms for U 

984 # https://dlmf.nist.gov/12.10#E42 

985 phip = phi ** arange(6, 31, 6).reshape((-1,1)) 

986 A0 = b0*u0 

987 A1 = (b2*u0 + phip[0,:]*b1*u1 + phip[1,:]*b0*u2) / zeta**3 

988 A2 = (b4*u0 + phip[0,:]*b3*u1 + phip[1,:]*b2*u2 + phip[2,:]*b1*u3 + phip[3,:]*b0*u4) / zeta**6 

989 B0 = -(a1*u0 + phip[0,:]*a0*u1) / zeta**2 

990 B1 = -(a3*u0 + phip[0,:]*a2*u1 + phip[1,:]*a1*u2 + phip[2,:]*a0*u3) / zeta**5 

991 B2 = -(a5*u0 + phip[0,:]*a4*u1 + phip[1,:]*a3*u2 + phip[2,:]*a2*u3 + phip[3,:]*a1*u4 + phip[4,:]*a0*u5) / zeta**8 

992 # U 

993 # https://dlmf.nist.gov/12.10#E35 

994 U = P * (Ai * (A0 + A1/mu**2.0 + A2/mu**4.0) + 

995 Aip * (B0 + B1/mu**2.0 + B2/mu**4.0) / mu**(8.0/6.0)) 

996 # Prefactor for derivative of U 

997 Pd = sqrt(2.0*pi) * mu**(2.0/6.0) / phi 

998 # Terms for derivative of U 

999 # https://dlmf.nist.gov/12.10#E46 

1000 C0 = -(b1*v0 + phip[0,:]*b0*v1) / zeta 

1001 C1 = -(b3*v0 + phip[0,:]*b2*v1 + phip[1,:]*b1*v2 + phip[2,:]*b0*v3) / zeta**4 

1002 C2 = -(b5*v0 + phip[0,:]*b4*v1 + phip[1,:]*b3*v2 + phip[2,:]*b2*v3 + phip[3,:]*b1*v4 + phip[4,:]*b0*v5) / zeta**7 

1003 D0 = a0*v0 

1004 D1 = (a2*v0 + phip[0,:]*a1*v1 + phip[1,:]*a0*v2) / zeta**3 

1005 D2 = (a4*v0 + phip[0,:]*a3*v1 + phip[1,:]*a2*v2 + phip[2,:]*a1*v3 + phip[3,:]*a0*v4) / zeta**6 

1006 # Derivative of U 

1007 # https://dlmf.nist.gov/12.10#E36 

1008 Ud = Pd * (Ai * (C0 + C1/mu**2.0 + C2/mu**4.0) / mu**(4.0/6.0) + 

1009 Aip * (D0 + D1/mu**2.0 + D2/mu**4.0)) 

1010 return U, Ud 

1011 

1012 

1013def _newton(n, x_initial, maxit=5): 

1014 """Newton iteration for polishing the asymptotic approximation 

1015 to the zeros of the Hermite polynomials. 

1016 

1017 Parameters 

1018 ---------- 

1019 n : int 

1020 Quadrature order 

1021 x_initial : ndarray 

1022 Initial guesses for the roots 

1023 maxit : int 

1024 Maximal number of Newton iterations. 

1025 The default 5 is sufficient, usually 

1026 only one or two steps are needed. 

1027 

1028 Returns 

1029 ------- 

1030 nodes : ndarray 

1031 Quadrature nodes 

1032 weights : ndarray 

1033 Quadrature weights 

1034 

1035 See Also 

1036 -------- 

1037 roots_hermite_asy 

1038 """ 

1039 # Variable transformation 

1040 mu = sqrt(2.0*n + 1.0) 

1041 t = x_initial / mu 

1042 theta = arccos(t) 

1043 # Newton iteration 

1044 for i in range(maxit): 

1045 u, ud = _pbcf(n, theta) 

1046 dtheta = u / (sqrt(2.0) * mu * sin(theta) * ud) 

1047 theta = theta + dtheta 

1048 if max(abs(dtheta)) < 1e-14: 

1049 break 

1050 # Undo variable transformation 

1051 x = mu * cos(theta) 

1052 # Central node is always zero 

1053 if n % 2 == 1: 

1054 x[0] = 0.0 

1055 # Compute weights 

1056 w = exp(-x**2) / (2.0*ud**2) 

1057 return x, w 

1058 

1059 

1060def _roots_hermite_asy(n): 

1061 r"""Gauss-Hermite (physicist's) quadrature for large n. 

1062 

1063 Computes the sample points and weights for Gauss-Hermite quadrature. 

1064 The sample points are the roots of the nth degree Hermite polynomial, 

1065 :math:`H_n(x)`. These sample points and weights correctly integrate 

1066 polynomials of degree :math:`2n - 1` or less over the interval 

1067 :math:`[-\infty, \infty]` with weight function :math:`f(x) = e^{-x^2}`. 

1068 

1069 This method relies on asymptotic expansions which work best for n > 150. 

1070 The algorithm has linear runtime making computation for very large n 

1071 feasible. 

1072 

1073 Parameters 

1074 ---------- 

1075 n : int 

1076 quadrature order 

1077 

1078 Returns 

1079 ------- 

1080 nodes : ndarray 

1081 Quadrature nodes 

1082 weights : ndarray 

1083 Quadrature weights 

1084 

1085 See Also 

1086 -------- 

1087 roots_hermite 

1088 

1089 References 

1090 ---------- 

1091 .. [townsend.trogdon.olver-2014] 

1092 Townsend, A. and Trogdon, T. and Olver, S. (2014) 

1093 *Fast computation of Gauss quadrature nodes and 

1094 weights on the whole real line*. :arXiv:`1410.5286`. 

1095 

1096 .. [townsend.trogdon.olver-2015] 

1097 Townsend, A. and Trogdon, T. and Olver, S. (2015) 

1098 *Fast computation of Gauss quadrature nodes and 

1099 weights on the whole real line*. 

1100 IMA Journal of Numerical Analysis 

1101 :doi:`10.1093/imanum/drv002`. 

1102 """ 

1103 iv = _initial_nodes(n) 

1104 nodes, weights = _newton(n, iv) 

1105 # Combine with negative parts 

1106 if n % 2 == 0: 

1107 nodes = hstack([-nodes[::-1], nodes]) 

1108 weights = hstack([weights[::-1], weights]) 

1109 else: 

1110 nodes = hstack([-nodes[-1:0:-1], nodes]) 

1111 weights = hstack([weights[-1:0:-1], weights]) 

1112 # Scale weights 

1113 weights *= sqrt(pi) / sum(weights) 

1114 return nodes, weights 

1115 

1116 

1117def hermite(n, monic=False): 

1118 r"""Physicist's Hermite polynomial. 

1119 

1120 Defined by 

1121 

1122 .. math:: 

1123 

1124 H_n(x) = (-1)^ne^{x^2}\frac{d^n}{dx^n}e^{-x^2}; 

1125 

1126 :math:`H_n` is a polynomial of degree :math:`n`. 

1127 

1128 Parameters 

1129 ---------- 

1130 n : int 

1131 Degree of the polynomial. 

1132 monic : bool, optional 

1133 If `True`, scale the leading coefficient to be 1. Default is 

1134 `False`. 

1135 

1136 Returns 

1137 ------- 

1138 H : orthopoly1d 

1139 Hermite polynomial. 

1140 

1141 Notes 

1142 ----- 

1143 The polynomials :math:`H_n` are orthogonal over :math:`(-\infty, 

1144 \infty)` with weight function :math:`e^{-x^2}`. 

1145 

1146 """ 

1147 if n < 0: 

1148 raise ValueError("n must be nonnegative.") 

1149 

1150 if n == 0: 

1151 n1 = n + 1 

1152 else: 

1153 n1 = n 

1154 x, w, mu0 = roots_hermite(n1, mu=True) 

1155 wfunc = lambda x: exp(-x * x) 

1156 if n == 0: 

1157 x, w = [], [] 

1158 hn = 2**n * _gam(n + 1) * sqrt(pi) 

1159 kn = 2**n 

1160 p = orthopoly1d(x, w, hn, kn, wfunc, (-inf, inf), monic, 

1161 lambda x: eval_hermite(n, x)) 

1162 return p 

1163 

1164# Hermite 2 He_n(x) 

1165 

1166 

1167def roots_hermitenorm(n, mu=False): 

1168 r"""Gauss-Hermite (statistician's) quadrature. 

1169 

1170 Compute the sample points and weights for Gauss-Hermite 

1171 quadrature. The sample points are the roots of the nth degree 

1172 Hermite polynomial, :math:`He_n(x)`. These sample points and 

1173 weights correctly integrate polynomials of degree :math:`2n - 1` 

1174 or less over the interval :math:`[-\infty, \infty]` with weight 

1175 function :math:`w(x) = e^{-x^2/2}`. See 22.2.15 in [AS]_ for more 

1176 details. 

1177 

1178 Parameters 

1179 ---------- 

1180 n : int 

1181 quadrature order 

1182 mu : bool, optional 

1183 If True, return the sum of the weights, optional. 

1184 

1185 Returns 

1186 ------- 

1187 x : ndarray 

1188 Sample points 

1189 w : ndarray 

1190 Weights 

1191 mu : float 

1192 Sum of the weights 

1193 

1194 Notes 

1195 ----- 

1196 For small n up to 150 a modified version of the Golub-Welsch 

1197 algorithm is used. Nodes are computed from the eigenvalue 

1198 problem and improved by one step of a Newton iteration. 

1199 The weights are computed from the well-known analytical formula. 

1200 

1201 For n larger than 150 an optimal asymptotic algorithm is used 

1202 which computes nodes and weights in a numerical stable manner. 

1203 The algorithm has linear runtime making computation for very 

1204 large n (several thousand or more) feasible. 

1205 

1206 See Also 

1207 -------- 

1208 scipy.integrate.quadrature 

1209 scipy.integrate.fixed_quad 

1210 numpy.polynomial.hermite_e.hermegauss 

1211 

1212 References 

1213 ---------- 

1214 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

1215 Handbook of Mathematical Functions with Formulas, 

1216 Graphs, and Mathematical Tables. New York: Dover, 1972. 

1217 

1218 """ 

1219 m = int(n) 

1220 if n < 1 or n != m: 

1221 raise ValueError("n must be a positive integer.") 

1222 

1223 mu0 = np.sqrt(2.0*np.pi) 

1224 if n <= 150: 

1225 an_func = lambda k: 0.0*k 

1226 bn_func = lambda k: np.sqrt(k) 

1227 f = cephes.eval_hermitenorm 

1228 df = lambda n, x: n * cephes.eval_hermitenorm(n-1, x) 

1229 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu) 

1230 else: 

1231 nodes, weights = _roots_hermite_asy(m) 

1232 # Transform 

1233 nodes *= sqrt(2) 

1234 weights *= sqrt(2) 

1235 if mu: 

1236 return nodes, weights, mu0 

1237 else: 

1238 return nodes, weights 

1239 

1240 

1241def hermitenorm(n, monic=False): 

1242 r"""Normalized (probabilist's) Hermite polynomial. 

1243 

1244 Defined by 

1245 

1246 .. math:: 

1247 

1248 He_n(x) = (-1)^ne^{x^2/2}\frac{d^n}{dx^n}e^{-x^2/2}; 

1249 

1250 :math:`He_n` is a polynomial of degree :math:`n`. 

1251 

1252 Parameters 

1253 ---------- 

1254 n : int 

1255 Degree of the polynomial. 

1256 monic : bool, optional 

1257 If `True`, scale the leading coefficient to be 1. Default is 

1258 `False`. 

1259 

1260 Returns 

1261 ------- 

1262 He : orthopoly1d 

1263 Hermite polynomial. 

1264 

1265 Notes 

1266 ----- 

1267 

1268 The polynomials :math:`He_n` are orthogonal over :math:`(-\infty, 

1269 \infty)` with weight function :math:`e^{-x^2/2}`. 

1270 

1271 """ 

1272 if n < 0: 

1273 raise ValueError("n must be nonnegative.") 

1274 

1275 if n == 0: 

1276 n1 = n + 1 

1277 else: 

1278 n1 = n 

1279 x, w, mu0 = roots_hermitenorm(n1, mu=True) 

1280 wfunc = lambda x: exp(-x * x / 2.0) 

1281 if n == 0: 

1282 x, w = [], [] 

1283 hn = sqrt(2 * pi) * _gam(n + 1) 

1284 kn = 1.0 

1285 p = orthopoly1d(x, w, hn, kn, wfunc=wfunc, limits=(-inf, inf), monic=monic, 

1286 eval_func=lambda x: eval_hermitenorm(n, x)) 

1287 return p 

1288 

1289# The remainder of the polynomials can be derived from the ones above. 

1290 

1291# Ultraspherical (Gegenbauer) C^(alpha)_n(x) 

1292 

1293 

1294def roots_gegenbauer(n, alpha, mu=False): 

1295 r"""Gauss-Gegenbauer quadrature. 

1296 

1297 Compute the sample points and weights for Gauss-Gegenbauer 

1298 quadrature. The sample points are the roots of the nth degree 

1299 Gegenbauer polynomial, :math:`C^{\alpha}_n(x)`. These sample 

1300 points and weights correctly integrate polynomials of degree 

1301 :math:`2n - 1` or less over the interval :math:`[-1, 1]` with 

1302 weight function :math:`w(x) = (1 - x^2)^{\alpha - 1/2}`. See 

1303 22.2.3 in [AS]_ for more details. 

1304 

1305 Parameters 

1306 ---------- 

1307 n : int 

1308 quadrature order 

1309 alpha : float 

1310 alpha must be > -0.5 

1311 mu : bool, optional 

1312 If True, return the sum of the weights, optional. 

1313 

1314 Returns 

1315 ------- 

1316 x : ndarray 

1317 Sample points 

1318 w : ndarray 

1319 Weights 

1320 mu : float 

1321 Sum of the weights 

1322 

1323 See Also 

1324 -------- 

1325 scipy.integrate.quadrature 

1326 scipy.integrate.fixed_quad 

1327 

1328 References 

1329 ---------- 

1330 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

1331 Handbook of Mathematical Functions with Formulas, 

1332 Graphs, and Mathematical Tables. New York: Dover, 1972. 

1333 

1334 """ 

1335 m = int(n) 

1336 if n < 1 or n != m: 

1337 raise ValueError("n must be a positive integer.") 

1338 if alpha < -0.5: 

1339 raise ValueError("alpha must be greater than -0.5.") 

1340 elif alpha == 0.0: 

1341 # C(n,0,x) == 0 uniformly, however, as alpha->0, C(n,alpha,x)->T(n,x) 

1342 # strictly, we should just error out here, since the roots are not 

1343 # really defined, but we used to return something useful, so let's 

1344 # keep doing so. 

1345 return roots_chebyt(n, mu) 

1346 

1347 mu0 = np.sqrt(np.pi) * cephes.gamma(alpha + 0.5) / cephes.gamma(alpha + 1) 

1348 an_func = lambda k: 0.0 * k 

1349 bn_func = lambda k: np.sqrt(k * (k + 2 * alpha - 1) 

1350 / (4 * (k + alpha) * (k + alpha - 1))) 

1351 f = lambda n, x: cephes.eval_gegenbauer(n, alpha, x) 

1352 df = lambda n, x: (-n*x*cephes.eval_gegenbauer(n, alpha, x) 

1353 + (n + 2*alpha - 1)*cephes.eval_gegenbauer(n-1, alpha, x))/(1-x**2) 

1354 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu) 

1355 

1356 

1357def gegenbauer(n, alpha, monic=False): 

1358 r"""Gegenbauer (ultraspherical) polynomial. 

1359 

1360 Defined to be the solution of 

1361 

1362 .. math:: 

1363 (1 - x^2)\frac{d^2}{dx^2}C_n^{(\alpha)} 

1364 - (2\alpha + 1)x\frac{d}{dx}C_n^{(\alpha)} 

1365 + n(n + 2\alpha)C_n^{(\alpha)} = 0 

1366 

1367 for :math:`\alpha > -1/2`; :math:`C_n^{(\alpha)}` is a polynomial 

1368 of degree :math:`n`. 

1369 

1370 Parameters 

1371 ---------- 

1372 n : int 

1373 Degree of the polynomial. 

1374 monic : bool, optional 

1375 If `True`, scale the leading coefficient to be 1. Default is 

1376 `False`. 

1377 

1378 Returns 

1379 ------- 

1380 C : orthopoly1d 

1381 Gegenbauer polynomial. 

1382 

1383 Notes 

1384 ----- 

1385 The polynomials :math:`C_n^{(\alpha)}` are orthogonal over 

1386 :math:`[-1,1]` with weight function :math:`(1 - x^2)^{(\alpha - 

1387 1/2)}`. 

1388 

1389 """ 

1390 base = jacobi(n, alpha - 0.5, alpha - 0.5, monic=monic) 

1391 if monic: 

1392 return base 

1393 # Abrahmowitz and Stegan 22.5.20 

1394 factor = (_gam(2*alpha + n) * _gam(alpha + 0.5) / 

1395 _gam(2*alpha) / _gam(alpha + 0.5 + n)) 

1396 base._scale(factor) 

1397 base.__dict__['_eval_func'] = lambda x: eval_gegenbauer(float(n), alpha, x) 

1398 return base 

1399 

1400# Chebyshev of the first kind: T_n(x) = 

1401# n! sqrt(pi) / _gam(n+1./2)* P^(-1/2,-1/2)_n(x) 

1402# Computed anew. 

1403 

1404 

1405def roots_chebyt(n, mu=False): 

1406 r"""Gauss-Chebyshev (first kind) quadrature. 

1407 

1408 Computes the sample points and weights for Gauss-Chebyshev 

1409 quadrature. The sample points are the roots of the nth degree 

1410 Chebyshev polynomial of the first kind, :math:`T_n(x)`. These 

1411 sample points and weights correctly integrate polynomials of 

1412 degree :math:`2n - 1` or less over the interval :math:`[-1, 1]` 

1413 with weight function :math:`w(x) = 1/\sqrt{1 - x^2}`. See 22.2.4 

1414 in [AS]_ for more details. 

1415 

1416 Parameters 

1417 ---------- 

1418 n : int 

1419 quadrature order 

1420 mu : bool, optional 

1421 If True, return the sum of the weights, optional. 

1422 

1423 Returns 

1424 ------- 

1425 x : ndarray 

1426 Sample points 

1427 w : ndarray 

1428 Weights 

1429 mu : float 

1430 Sum of the weights 

1431 

1432 See Also 

1433 -------- 

1434 scipy.integrate.quadrature 

1435 scipy.integrate.fixed_quad 

1436 numpy.polynomial.chebyshev.chebgauss 

1437 

1438 References 

1439 ---------- 

1440 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

1441 Handbook of Mathematical Functions with Formulas, 

1442 Graphs, and Mathematical Tables. New York: Dover, 1972. 

1443 

1444 """ 

1445 m = int(n) 

1446 if n < 1 or n != m: 

1447 raise ValueError('n must be a positive integer.') 

1448 x = _ufuncs._sinpi(np.arange(-m + 1, m, 2) / (2*m)) 

1449 w = np.full_like(x, pi/m) 

1450 if mu: 

1451 return x, w, pi 

1452 else: 

1453 return x, w 

1454 

1455 

1456def chebyt(n, monic=False): 

1457 r"""Chebyshev polynomial of the first kind. 

1458 

1459 Defined to be the solution of 

1460 

1461 .. math:: 

1462 (1 - x^2)\frac{d^2}{dx^2}T_n - x\frac{d}{dx}T_n + n^2T_n = 0; 

1463 

1464 :math:`T_n` is a polynomial of degree :math:`n`. 

1465 

1466 Parameters 

1467 ---------- 

1468 n : int 

1469 Degree of the polynomial. 

1470 monic : bool, optional 

1471 If `True`, scale the leading coefficient to be 1. Default is 

1472 `False`. 

1473 

1474 Returns 

1475 ------- 

1476 T : orthopoly1d 

1477 Chebyshev polynomial of the first kind. 

1478 

1479 Notes 

1480 ----- 

1481 The polynomials :math:`T_n` are orthogonal over :math:`[-1, 1]` 

1482 with weight function :math:`(1 - x^2)^{-1/2}`. 

1483 

1484 See Also 

1485 -------- 

1486 chebyu : Chebyshev polynomial of the second kind. 

1487 

1488 """ 

1489 if n < 0: 

1490 raise ValueError("n must be nonnegative.") 

1491 

1492 wfunc = lambda x: 1.0 / sqrt(1 - x * x) 

1493 if n == 0: 

1494 return orthopoly1d([], [], pi, 1.0, wfunc, (-1, 1), monic, 

1495 lambda x: eval_chebyt(n, x)) 

1496 n1 = n 

1497 x, w, mu = roots_chebyt(n1, mu=True) 

1498 hn = pi / 2 

1499 kn = 2**(n - 1) 

1500 p = orthopoly1d(x, w, hn, kn, wfunc, (-1, 1), monic, 

1501 lambda x: eval_chebyt(n, x)) 

1502 return p 

1503 

1504# Chebyshev of the second kind 

1505# U_n(x) = (n+1)! sqrt(pi) / (2*_gam(n+3./2)) * P^(1/2,1/2)_n(x) 

1506 

1507 

1508def roots_chebyu(n, mu=False): 

1509 r"""Gauss-Chebyshev (second kind) quadrature. 

1510 

1511 Computes the sample points and weights for Gauss-Chebyshev 

1512 quadrature. The sample points are the roots of the nth degree 

1513 Chebyshev polynomial of the second kind, :math:`U_n(x)`. These 

1514 sample points and weights correctly integrate polynomials of 

1515 degree :math:`2n - 1` or less over the interval :math:`[-1, 1]` 

1516 with weight function :math:`w(x) = \sqrt{1 - x^2}`. See 22.2.5 in 

1517 [AS]_ for details. 

1518 

1519 Parameters 

1520 ---------- 

1521 n : int 

1522 quadrature order 

1523 mu : bool, optional 

1524 If True, return the sum of the weights, optional. 

1525 

1526 Returns 

1527 ------- 

1528 x : ndarray 

1529 Sample points 

1530 w : ndarray 

1531 Weights 

1532 mu : float 

1533 Sum of the weights 

1534 

1535 See Also 

1536 -------- 

1537 scipy.integrate.quadrature 

1538 scipy.integrate.fixed_quad 

1539 

1540 References 

1541 ---------- 

1542 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

1543 Handbook of Mathematical Functions with Formulas, 

1544 Graphs, and Mathematical Tables. New York: Dover, 1972. 

1545 

1546 """ 

1547 m = int(n) 

1548 if n < 1 or n != m: 

1549 raise ValueError('n must be a positive integer.') 

1550 t = np.arange(m, 0, -1) * pi / (m + 1) 

1551 x = np.cos(t) 

1552 w = pi * np.sin(t)**2 / (m + 1) 

1553 if mu: 

1554 return x, w, pi / 2 

1555 else: 

1556 return x, w 

1557 

1558 

1559def chebyu(n, monic=False): 

1560 r"""Chebyshev polynomial of the second kind. 

1561 

1562 Defined to be the solution of 

1563 

1564 .. math:: 

1565 (1 - x^2)\frac{d^2}{dx^2}U_n - 3x\frac{d}{dx}U_n 

1566 + n(n + 2)U_n = 0; 

1567 

1568 :math:`U_n` is a polynomial of degree :math:`n`. 

1569 

1570 Parameters 

1571 ---------- 

1572 n : int 

1573 Degree of the polynomial. 

1574 monic : bool, optional 

1575 If `True`, scale the leading coefficient to be 1. Default is 

1576 `False`. 

1577 

1578 Returns 

1579 ------- 

1580 U : orthopoly1d 

1581 Chebyshev polynomial of the second kind. 

1582 

1583 Notes 

1584 ----- 

1585 The polynomials :math:`U_n` are orthogonal over :math:`[-1, 1]` 

1586 with weight function :math:`(1 - x^2)^{1/2}`. 

1587 

1588 See Also 

1589 -------- 

1590 chebyt : Chebyshev polynomial of the first kind. 

1591 

1592 """ 

1593 base = jacobi(n, 0.5, 0.5, monic=monic) 

1594 if monic: 

1595 return base 

1596 factor = sqrt(pi) / 2.0 * _gam(n + 2) / _gam(n + 1.5) 

1597 base._scale(factor) 

1598 return base 

1599 

1600# Chebyshev of the first kind C_n(x) 

1601 

1602 

1603def roots_chebyc(n, mu=False): 

1604 r"""Gauss-Chebyshev (first kind) quadrature. 

1605 

1606 Compute the sample points and weights for Gauss-Chebyshev 

1607 quadrature. The sample points are the roots of the nth degree 

1608 Chebyshev polynomial of the first kind, :math:`C_n(x)`. These 

1609 sample points and weights correctly integrate polynomials of 

1610 degree :math:`2n - 1` or less over the interval :math:`[-2, 2]` 

1611 with weight function :math:`w(x) = 1 / \sqrt{1 - (x/2)^2}`. See 

1612 22.2.6 in [AS]_ for more details. 

1613 

1614 Parameters 

1615 ---------- 

1616 n : int 

1617 quadrature order 

1618 mu : bool, optional 

1619 If True, return the sum of the weights, optional. 

1620 

1621 Returns 

1622 ------- 

1623 x : ndarray 

1624 Sample points 

1625 w : ndarray 

1626 Weights 

1627 mu : float 

1628 Sum of the weights 

1629 

1630 See Also 

1631 -------- 

1632 scipy.integrate.quadrature 

1633 scipy.integrate.fixed_quad 

1634 

1635 References 

1636 ---------- 

1637 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

1638 Handbook of Mathematical Functions with Formulas, 

1639 Graphs, and Mathematical Tables. New York: Dover, 1972. 

1640 

1641 """ 

1642 x, w, m = roots_chebyt(n, True) 

1643 x *= 2 

1644 w *= 2 

1645 m *= 2 

1646 if mu: 

1647 return x, w, m 

1648 else: 

1649 return x, w 

1650 

1651 

1652def chebyc(n, monic=False): 

1653 r"""Chebyshev polynomial of the first kind on :math:`[-2, 2]`. 

1654 

1655 Defined as :math:`C_n(x) = 2T_n(x/2)`, where :math:`T_n` is the 

1656 nth Chebychev polynomial of the first kind. 

1657 

1658 Parameters 

1659 ---------- 

1660 n : int 

1661 Degree of the polynomial. 

1662 monic : bool, optional 

1663 If `True`, scale the leading coefficient to be 1. Default is 

1664 `False`. 

1665 

1666 Returns 

1667 ------- 

1668 C : orthopoly1d 

1669 Chebyshev polynomial of the first kind on :math:`[-2, 2]`. 

1670 

1671 Notes 

1672 ----- 

1673 The polynomials :math:`C_n(x)` are orthogonal over :math:`[-2, 2]` 

1674 with weight function :math:`1/\sqrt{1 - (x/2)^2}`. 

1675 

1676 See Also 

1677 -------- 

1678 chebyt : Chebyshev polynomial of the first kind. 

1679 

1680 References 

1681 ---------- 

1682 .. [1] Abramowitz and Stegun, "Handbook of Mathematical Functions" 

1683 Section 22. National Bureau of Standards, 1972. 

1684 

1685 """ 

1686 if n < 0: 

1687 raise ValueError("n must be nonnegative.") 

1688 

1689 if n == 0: 

1690 n1 = n + 1 

1691 else: 

1692 n1 = n 

1693 x, w, mu0 = roots_chebyc(n1, mu=True) 

1694 if n == 0: 

1695 x, w = [], [] 

1696 hn = 4 * pi * ((n == 0) + 1) 

1697 kn = 1.0 

1698 p = orthopoly1d(x, w, hn, kn, 

1699 wfunc=lambda x: 1.0 / sqrt(1 - x * x / 4.0), 

1700 limits=(-2, 2), monic=monic) 

1701 if not monic: 

1702 p._scale(2.0 / p(2)) 

1703 p.__dict__['_eval_func'] = lambda x: eval_chebyc(n, x) 

1704 return p 

1705 

1706# Chebyshev of the second kind S_n(x) 

1707 

1708 

1709def roots_chebys(n, mu=False): 

1710 r"""Gauss-Chebyshev (second kind) quadrature. 

1711 

1712 Compute the sample points and weights for Gauss-Chebyshev 

1713 quadrature. The sample points are the roots of the nth degree 

1714 Chebyshev polynomial of the second kind, :math:`S_n(x)`. These 

1715 sample points and weights correctly integrate polynomials of 

1716 degree :math:`2n - 1` or less over the interval :math:`[-2, 2]` 

1717 with weight function :math:`w(x) = \sqrt{1 - (x/2)^2}`. See 22.2.7 

1718 in [AS]_ for more details. 

1719 

1720 Parameters 

1721 ---------- 

1722 n : int 

1723 quadrature order 

1724 mu : bool, optional 

1725 If True, return the sum of the weights, optional. 

1726 

1727 Returns 

1728 ------- 

1729 x : ndarray 

1730 Sample points 

1731 w : ndarray 

1732 Weights 

1733 mu : float 

1734 Sum of the weights 

1735 

1736 See Also 

1737 -------- 

1738 scipy.integrate.quadrature 

1739 scipy.integrate.fixed_quad 

1740 

1741 References 

1742 ---------- 

1743 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

1744 Handbook of Mathematical Functions with Formulas, 

1745 Graphs, and Mathematical Tables. New York: Dover, 1972. 

1746 

1747 """ 

1748 x, w, m = roots_chebyu(n, True) 

1749 x *= 2 

1750 w *= 2 

1751 m *= 2 

1752 if mu: 

1753 return x, w, m 

1754 else: 

1755 return x, w 

1756 

1757 

1758def chebys(n, monic=False): 

1759 r"""Chebyshev polynomial of the second kind on :math:`[-2, 2]`. 

1760 

1761 Defined as :math:`S_n(x) = U_n(x/2)` where :math:`U_n` is the 

1762 nth Chebychev polynomial of the second kind. 

1763 

1764 Parameters 

1765 ---------- 

1766 n : int 

1767 Degree of the polynomial. 

1768 monic : bool, optional 

1769 If `True`, scale the leading coefficient to be 1. Default is 

1770 `False`. 

1771 

1772 Returns 

1773 ------- 

1774 S : orthopoly1d 

1775 Chebyshev polynomial of the second kind on :math:`[-2, 2]`. 

1776 

1777 Notes 

1778 ----- 

1779 The polynomials :math:`S_n(x)` are orthogonal over :math:`[-2, 2]` 

1780 with weight function :math:`\sqrt{1 - (x/2)}^2`. 

1781 

1782 See Also 

1783 -------- 

1784 chebyu : Chebyshev polynomial of the second kind 

1785 

1786 References 

1787 ---------- 

1788 .. [1] Abramowitz and Stegun, "Handbook of Mathematical Functions" 

1789 Section 22. National Bureau of Standards, 1972. 

1790 

1791 """ 

1792 if n < 0: 

1793 raise ValueError("n must be nonnegative.") 

1794 

1795 if n == 0: 

1796 n1 = n + 1 

1797 else: 

1798 n1 = n 

1799 x, w, mu0 = roots_chebys(n1, mu=True) 

1800 if n == 0: 

1801 x, w = [], [] 

1802 hn = pi 

1803 kn = 1.0 

1804 p = orthopoly1d(x, w, hn, kn, 

1805 wfunc=lambda x: sqrt(1 - x * x / 4.0), 

1806 limits=(-2, 2), monic=monic) 

1807 if not monic: 

1808 factor = (n + 1.0) / p(2) 

1809 p._scale(factor) 

1810 p.__dict__['_eval_func'] = lambda x: eval_chebys(n, x) 

1811 return p 

1812 

1813# Shifted Chebyshev of the first kind T^*_n(x) 

1814 

1815 

1816def roots_sh_chebyt(n, mu=False): 

1817 r"""Gauss-Chebyshev (first kind, shifted) quadrature. 

1818 

1819 Compute the sample points and weights for Gauss-Chebyshev 

1820 quadrature. The sample points are the roots of the nth degree 

1821 shifted Chebyshev polynomial of the first kind, :math:`T_n(x)`. 

1822 These sample points and weights correctly integrate polynomials of 

1823 degree :math:`2n - 1` or less over the interval :math:`[0, 1]` 

1824 with weight function :math:`w(x) = 1/\sqrt{x - x^2}`. See 22.2.8 

1825 in [AS]_ for more details. 

1826 

1827 Parameters 

1828 ---------- 

1829 n : int 

1830 quadrature order 

1831 mu : bool, optional 

1832 If True, return the sum of the weights, optional. 

1833 

1834 Returns 

1835 ------- 

1836 x : ndarray 

1837 Sample points 

1838 w : ndarray 

1839 Weights 

1840 mu : float 

1841 Sum of the weights 

1842 

1843 See Also 

1844 -------- 

1845 scipy.integrate.quadrature 

1846 scipy.integrate.fixed_quad 

1847 

1848 References 

1849 ---------- 

1850 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

1851 Handbook of Mathematical Functions with Formulas, 

1852 Graphs, and Mathematical Tables. New York: Dover, 1972. 

1853 

1854 """ 

1855 xw = roots_chebyt(n, mu) 

1856 return ((xw[0] + 1) / 2,) + xw[1:] 

1857 

1858 

1859def sh_chebyt(n, monic=False): 

1860 r"""Shifted Chebyshev polynomial of the first kind. 

1861 

1862 Defined as :math:`T^*_n(x) = T_n(2x - 1)` for :math:`T_n` the nth 

1863 Chebyshev polynomial of the first kind. 

1864 

1865 Parameters 

1866 ---------- 

1867 n : int 

1868 Degree of the polynomial. 

1869 monic : bool, optional 

1870 If `True`, scale the leading coefficient to be 1. Default is 

1871 `False`. 

1872 

1873 Returns 

1874 ------- 

1875 T : orthopoly1d 

1876 Shifted Chebyshev polynomial of the first kind. 

1877 

1878 Notes 

1879 ----- 

1880 The polynomials :math:`T^*_n` are orthogonal over :math:`[0, 1]` 

1881 with weight function :math:`(x - x^2)^{-1/2}`. 

1882 

1883 """ 

1884 base = sh_jacobi(n, 0.0, 0.5, monic=monic) 

1885 if monic: 

1886 return base 

1887 if n > 0: 

1888 factor = 4**n / 2.0 

1889 else: 

1890 factor = 1.0 

1891 base._scale(factor) 

1892 return base 

1893 

1894 

1895# Shifted Chebyshev of the second kind U^*_n(x) 

1896def roots_sh_chebyu(n, mu=False): 

1897 r"""Gauss-Chebyshev (second kind, shifted) quadrature. 

1898 

1899 Computes the sample points and weights for Gauss-Chebyshev 

1900 quadrature. The sample points are the roots of the nth degree 

1901 shifted Chebyshev polynomial of the second kind, :math:`U_n(x)`. 

1902 These sample points and weights correctly integrate polynomials of 

1903 degree :math:`2n - 1` or less over the interval :math:`[0, 1]` 

1904 with weight function :math:`w(x) = \sqrt{x - x^2}`. See 22.2.9 in 

1905 [AS]_ for more details. 

1906 

1907 Parameters 

1908 ---------- 

1909 n : int 

1910 quadrature order 

1911 mu : bool, optional 

1912 If True, return the sum of the weights, optional. 

1913 

1914 Returns 

1915 ------- 

1916 x : ndarray 

1917 Sample points 

1918 w : ndarray 

1919 Weights 

1920 mu : float 

1921 Sum of the weights 

1922 

1923 See Also 

1924 -------- 

1925 scipy.integrate.quadrature 

1926 scipy.integrate.fixed_quad 

1927 

1928 References 

1929 ---------- 

1930 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

1931 Handbook of Mathematical Functions with Formulas, 

1932 Graphs, and Mathematical Tables. New York: Dover, 1972. 

1933 

1934 """ 

1935 x, w, m = roots_chebyu(n, True) 

1936 x = (x + 1) / 2 

1937 m_us = cephes.beta(1.5, 1.5) 

1938 w *= m_us / m 

1939 if mu: 

1940 return x, w, m_us 

1941 else: 

1942 return x, w 

1943 

1944 

1945def sh_chebyu(n, monic=False): 

1946 r"""Shifted Chebyshev polynomial of the second kind. 

1947 

1948 Defined as :math:`U^*_n(x) = U_n(2x - 1)` for :math:`U_n` the nth 

1949 Chebyshev polynomial of the second kind. 

1950 

1951 Parameters 

1952 ---------- 

1953 n : int 

1954 Degree of the polynomial. 

1955 monic : bool, optional 

1956 If `True`, scale the leading coefficient to be 1. Default is 

1957 `False`. 

1958 

1959 Returns 

1960 ------- 

1961 U : orthopoly1d 

1962 Shifted Chebyshev polynomial of the second kind. 

1963 

1964 Notes 

1965 ----- 

1966 The polynomials :math:`U^*_n` are orthogonal over :math:`[0, 1]` 

1967 with weight function :math:`(x - x^2)^{1/2}`. 

1968 

1969 """ 

1970 base = sh_jacobi(n, 2.0, 1.5, monic=monic) 

1971 if monic: 

1972 return base 

1973 factor = 4**n 

1974 base._scale(factor) 

1975 return base 

1976 

1977# Legendre 

1978 

1979 

1980def roots_legendre(n, mu=False): 

1981 r"""Gauss-Legendre quadrature. 

1982 

1983 Compute the sample points and weights for Gauss-Legendre 

1984 quadrature. The sample points are the roots of the nth degree 

1985 Legendre polynomial :math:`P_n(x)`. These sample points and 

1986 weights correctly integrate polynomials of degree :math:`2n - 1` 

1987 or less over the interval :math:`[-1, 1]` with weight function 

1988 :math:`w(x) = 1.0`. See 2.2.10 in [AS]_ for more details. 

1989 

1990 Parameters 

1991 ---------- 

1992 n : int 

1993 quadrature order 

1994 mu : bool, optional 

1995 If True, return the sum of the weights, optional. 

1996 

1997 Returns 

1998 ------- 

1999 x : ndarray 

2000 Sample points 

2001 w : ndarray 

2002 Weights 

2003 mu : float 

2004 Sum of the weights 

2005 

2006 See Also 

2007 -------- 

2008 scipy.integrate.quadrature 

2009 scipy.integrate.fixed_quad 

2010 numpy.polynomial.legendre.leggauss 

2011 

2012 References 

2013 ---------- 

2014 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

2015 Handbook of Mathematical Functions with Formulas, 

2016 Graphs, and Mathematical Tables. New York: Dover, 1972. 

2017 

2018 """ 

2019 m = int(n) 

2020 if n < 1 or n != m: 

2021 raise ValueError("n must be a positive integer.") 

2022 

2023 mu0 = 2.0 

2024 an_func = lambda k: 0.0 * k 

2025 bn_func = lambda k: k * np.sqrt(1.0 / (4 * k * k - 1)) 

2026 f = cephes.eval_legendre 

2027 df = lambda n, x: (-n*x*cephes.eval_legendre(n, x) 

2028 + n*cephes.eval_legendre(n-1, x))/(1-x**2) 

2029 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu) 

2030 

2031 

2032def legendre(n, monic=False): 

2033 r"""Legendre polynomial. 

2034 

2035 Defined to be the solution of 

2036 

2037 .. math:: 

2038 \frac{d}{dx}\left[(1 - x^2)\frac{d}{dx}P_n(x)\right] 

2039 + n(n + 1)P_n(x) = 0; 

2040 

2041 :math:`P_n(x)` is a polynomial of degree :math:`n`. 

2042 

2043 Parameters 

2044 ---------- 

2045 n : int 

2046 Degree of the polynomial. 

2047 monic : bool, optional 

2048 If `True`, scale the leading coefficient to be 1. Default is 

2049 `False`. 

2050 

2051 Returns 

2052 ------- 

2053 P : orthopoly1d 

2054 Legendre polynomial. 

2055 

2056 Notes 

2057 ----- 

2058 The polynomials :math:`P_n` are orthogonal over :math:`[-1, 1]` 

2059 with weight function 1. 

2060 

2061 Examples 

2062 -------- 

2063 Generate the 3rd-order Legendre polynomial 1/2*(5x^3 + 0x^2 - 3x + 0): 

2064 

2065 >>> from scipy.special import legendre 

2066 >>> legendre(3) 

2067 poly1d([ 2.5, 0. , -1.5, 0. ]) 

2068 

2069 """ 

2070 if n < 0: 

2071 raise ValueError("n must be nonnegative.") 

2072 

2073 if n == 0: 

2074 n1 = n + 1 

2075 else: 

2076 n1 = n 

2077 x, w, mu0 = roots_legendre(n1, mu=True) 

2078 if n == 0: 

2079 x, w = [], [] 

2080 hn = 2.0 / (2 * n + 1) 

2081 kn = _gam(2 * n + 1) / _gam(n + 1)**2 / 2.0**n 

2082 p = orthopoly1d(x, w, hn, kn, wfunc=lambda x: 1.0, limits=(-1, 1), 

2083 monic=monic, eval_func=lambda x: eval_legendre(n, x)) 

2084 return p 

2085 

2086# Shifted Legendre P^*_n(x) 

2087 

2088 

2089def roots_sh_legendre(n, mu=False): 

2090 r"""Gauss-Legendre (shifted) quadrature. 

2091 

2092 Compute the sample points and weights for Gauss-Legendre 

2093 quadrature. The sample points are the roots of the nth degree 

2094 shifted Legendre polynomial :math:`P^*_n(x)`. These sample points 

2095 and weights correctly integrate polynomials of degree :math:`2n - 

2096 1` or less over the interval :math:`[0, 1]` with weight function 

2097 :math:`w(x) = 1.0`. See 2.2.11 in [AS]_ for details. 

2098 

2099 Parameters 

2100 ---------- 

2101 n : int 

2102 quadrature order 

2103 mu : bool, optional 

2104 If True, return the sum of the weights, optional. 

2105 

2106 Returns 

2107 ------- 

2108 x : ndarray 

2109 Sample points 

2110 w : ndarray 

2111 Weights 

2112 mu : float 

2113 Sum of the weights 

2114 

2115 See Also 

2116 -------- 

2117 scipy.integrate.quadrature 

2118 scipy.integrate.fixed_quad 

2119 

2120 References 

2121 ---------- 

2122 .. [AS] Milton Abramowitz and Irene A. Stegun, eds. 

2123 Handbook of Mathematical Functions with Formulas, 

2124 Graphs, and Mathematical Tables. New York: Dover, 1972. 

2125 

2126 """ 

2127 x, w = roots_legendre(n) 

2128 x = (x + 1) / 2 

2129 w /= 2 

2130 if mu: 

2131 return x, w, 1.0 

2132 else: 

2133 return x, w 

2134 

2135def sh_legendre(n, monic=False): 

2136 r"""Shifted Legendre polynomial. 

2137 

2138 Defined as :math:`P^*_n(x) = P_n(2x - 1)` for :math:`P_n` the nth 

2139 Legendre polynomial. 

2140 

2141 Parameters 

2142 ---------- 

2143 n : int 

2144 Degree of the polynomial. 

2145 monic : bool, optional 

2146 If `True`, scale the leading coefficient to be 1. Default is 

2147 `False`. 

2148 

2149 Returns 

2150 ------- 

2151 P : orthopoly1d 

2152 Shifted Legendre polynomial. 

2153 

2154 Notes 

2155 ----- 

2156 The polynomials :math:`P^*_n` are orthogonal over :math:`[0, 1]` 

2157 with weight function 1. 

2158 

2159 """ 

2160 if n < 0: 

2161 raise ValueError("n must be nonnegative.") 

2162 

2163 wfunc = lambda x: 0.0 * x + 1.0 

2164 if n == 0: 

2165 return orthopoly1d([], [], 1.0, 1.0, wfunc, (0, 1), monic, 

2166 lambda x: eval_sh_legendre(n, x)) 

2167 x, w, mu0 = roots_sh_legendre(n, mu=True) 

2168 hn = 1.0 / (2 * n + 1.0) 

2169 kn = _gam(2 * n + 1) / _gam(n + 1)**2 

2170 p = orthopoly1d(x, w, hn, kn, wfunc, limits=(0, 1), monic=monic, 

2171 eval_func=lambda x: eval_sh_legendre(n, x)) 

2172 return p 

2173 

2174 

2175# ----------------------------------------------------------------------------- 

2176# Code for backwards compatibility 

2177# ----------------------------------------------------------------------------- 

2178 

2179# Import functions in case someone is still calling the orthogonal 

2180# module directly. (They shouldn't be; it's not in the public API). 

2181poch = cephes.poch 

2182 

2183# eval_chebyu, eval_sh_chebyt and eval_sh_chebyu: These functions are not 

2184# used in orthogonal.py, they are not in _rootfuns_map, but their names 

2185# do appear in _evalfuns, so they must be kept. 

2186from ._ufuncs import (binom, eval_jacobi, eval_sh_jacobi, eval_gegenbauer, 

2187 eval_chebyt, eval_chebyu, eval_chebys, eval_chebyc, 

2188 eval_sh_chebyt, eval_sh_chebyu, eval_legendre, 

2189 eval_sh_legendre, eval_genlaguerre, eval_laguerre, 

2190 eval_hermite, eval_hermitenorm) 

2191 

2192# Make the old root function names an alias for the new ones 

2193_modattrs = globals() 

2194for newfun, oldfun in _rootfuns_map.items(): 

2195 _modattrs[oldfun] = _modattrs[newfun] 

2196 __all__.append(oldfun)