Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/special/_orthogonal.py: 11%
527 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""
2A collection of functions to find the weights and abscissas for
3Gaussian Quadrature.
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.
10Many recursion relations for orthogonal polynomials are given:
12.. math::
14 a1n f_{n+1} (x) = (a2n + a3n x ) f_n (x) - a4n f_{n-1} (x)
16The recursion relation of interest is
18.. math::
20 P_{n+1} (x) = (x - A_n) P_n (x) - B_n P_{n-1} (x)
22where :math:`P` has a different normalization than :math:`f`.
24The coefficients can be found as:
26.. math::
28 A_n = -a2n / a3n
29 \\qquad
30 B_n = ( a4n / a3n \\sqrt{h_n-1 / h_n})^2
32where
34.. math::
36 h_n = \\int_a^b w(x) f_n(x)^2
38assume:
40.. math::
42 P_0 (x) = 1
43 \\qquad
44 P_{-1} (x) == 0
46For the mathematical background, see [golub.welsch-1969-mathcomp]_ and
47[abramowitz.stegun-1965]_.
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.
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/
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`.
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)
77# SciPy imports.
78import numpy as np
79from numpy import (exp, inf, pi, sqrt, floor, sin, cos, around,
80 hstack, arccos, arange)
81from scipy import linalg
82from scipy.special import airy
84# Local imports.
85from . import _ufuncs
86_gam = _ufuncs.gamma
87# There is no .pyi file for _specfun
88from . import _specfun # type: ignore
90_polyfuns = ['legendre', 'chebyt', 'chebyu', 'chebyc', 'chebys',
91 'jacobi', 'laguerre', 'genlaguerre', 'hermite',
92 'hermitenorm', 'gegenbauer', 'sh_legendre', 'sh_chebyt',
93 'sh_chebyu', 'sh_jacobi']
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'}
112__all__ = _polyfuns + list(_rootfuns_map.keys())
115class orthopoly1d(np.poly1d):
117 def __init__(self, roots, weights=None, hn=1.0, kn=1.0, wfunc=None,
118 limits=None, monic=False, eval_func=None):
119 equiv_weights = [weights[k] / wfunc(roots[k]) for
120 k in range(len(roots))]
121 mu = sqrt(hn)
122 if monic:
123 evf = eval_func
124 if evf:
125 knn = kn
126 eval_func = lambda x: evf(x) / knn
127 mu = mu / abs(kn)
128 kn = 1.0
130 # compute coefficients from roots, then scale
131 poly = np.poly1d(roots, r=True)
132 np.poly1d.__init__(self, poly.coeffs * float(kn))
134 self.weights = np.array(list(zip(roots, weights, equiv_weights)))
135 self.weight_func = wfunc
136 self.limits = limits
137 self.normcoef = mu
139 # Note: eval_func will be discarded on arithmetic
140 self._eval_func = eval_func
142 def __call__(self, v):
143 if self._eval_func and not isinstance(v, np.poly1d):
144 return self._eval_func(v)
145 else:
146 return np.poly1d.__call__(self, v)
148 def _scale(self, p):
149 if p == 1.0:
150 return
151 self._coeffs *= p
153 evf = self._eval_func
154 if evf:
155 self._eval_func = lambda x: evf(x) * p
156 self.normcoef *= p
159def _gen_roots_and_weights(n, mu0, an_func, bn_func, f, df, symmetrize, mu):
160 """[x,w] = gen_roots_and_weights(n,an_func,sqrt_bn_func,mu)
162 Returns the roots (x) of an nth order orthogonal polynomial,
163 and weights (w) to use in appropriate Gaussian quadrature with that
164 orthogonal polynomial.
166 The polynomials have the recurrence relation
167 P_n+1(x) = (x - A_n) P_n(x) - B_n P_n-1(x)
169 an_func(n) should return A_n
170 sqrt_bn_func(n) should return sqrt(B_n)
171 mu ( = h_0 ) is the integral of the weight over the orthogonal
172 interval
173 """
174 k = np.arange(n, dtype='d')
175 c = np.zeros((2, n))
176 c[0,1:] = bn_func(k[1:])
177 c[1,:] = an_func(k)
178 x = linalg.eigvals_banded(c, overwrite_a_band=True)
180 # improve roots by one application of Newton's method
181 y = f(n, x)
182 dy = df(n, x)
183 x -= y/dy
185 # fm and dy may contain very large/small values, so we
186 # log-normalize them to maintain precision in the product fm*dy
187 fm = f(n-1, x)
188 log_fm = np.log(np.abs(fm))
189 log_dy = np.log(np.abs(dy))
190 fm /= np.exp((log_fm.max() + log_fm.min()) / 2.)
191 dy /= np.exp((log_dy.max() + log_dy.min()) / 2.)
192 w = 1.0 / (fm * dy)
194 if symmetrize:
195 w = (w + w[::-1]) / 2
196 x = (x - x[::-1]) / 2
198 w *= mu0 / w.sum()
200 if mu:
201 return x, w, mu0
202 else:
203 return x, w
205# Jacobi Polynomials 1 P^(alpha,beta)_n(x)
208def roots_jacobi(n, alpha, beta, mu=False):
209 r"""Gauss-Jacobi quadrature.
211 Compute the sample points and weights for Gauss-Jacobi
212 quadrature. The sample points are the roots of the nth degree
213 Jacobi polynomial, :math:`P^{\alpha, \beta}_n(x)`. These sample
214 points and weights correctly integrate polynomials of degree
215 :math:`2n - 1` or less over the interval :math:`[-1, 1]` with
216 weight function :math:`w(x) = (1 - x)^{\alpha} (1 +
217 x)^{\beta}`. See 22.2.1 in [AS]_ for details.
219 Parameters
220 ----------
221 n : int
222 quadrature order
223 alpha : float
224 alpha must be > -1
225 beta : float
226 beta must be > -1
227 mu : bool, optional
228 If True, return the sum of the weights, optional.
230 Returns
231 -------
232 x : ndarray
233 Sample points
234 w : ndarray
235 Weights
236 mu : float
237 Sum of the weights
239 See Also
240 --------
241 scipy.integrate.quadrature
242 scipy.integrate.fixed_quad
244 References
245 ----------
246 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
247 Handbook of Mathematical Functions with Formulas,
248 Graphs, and Mathematical Tables. New York: Dover, 1972.
250 """
251 m = int(n)
252 if n < 1 or n != m:
253 raise ValueError("n must be a positive integer.")
254 if alpha <= -1 or beta <= -1:
255 raise ValueError("alpha and beta must be greater than -1.")
257 if alpha == 0.0 and beta == 0.0:
258 return roots_legendre(m, mu)
259 if alpha == beta:
260 return roots_gegenbauer(m, alpha+0.5, mu)
262 if (alpha + beta) <= 1000:
263 mu0 = 2.0**(alpha+beta+1) * _ufuncs.beta(alpha+1, beta+1)
264 else:
265 # Avoid overflows in pow and beta for very large parameters
266 mu0 = np.exp((alpha + beta + 1) * np.log(2.0)
267 + _ufuncs.betaln(alpha+1, beta+1))
268 a = alpha
269 b = beta
270 if a + b == 0.0:
271 an_func = lambda k: np.where(k == 0, (b-a)/(2+a+b), 0.0)
272 else:
273 an_func = lambda k: np.where(k == 0, (b-a)/(2+a+b),
274 (b*b - a*a) / ((2.0*k+a+b)*(2.0*k+a+b+2)))
276 bn_func = lambda k: 2.0 / (2.0*k+a+b)*np.sqrt((k+a)*(k+b) / (2*k+a+b+1)) \
277 * np.where(k == 1, 1.0, np.sqrt(k*(k+a+b) / (2.0*k+a+b-1)))
279 f = lambda n, x: _ufuncs.eval_jacobi(n, a, b, x)
280 df = lambda n, x: (0.5 * (n + a + b + 1)
281 * _ufuncs.eval_jacobi(n-1, a+1, b+1, x))
282 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, False, mu)
285def jacobi(n, alpha, beta, monic=False):
286 r"""Jacobi polynomial.
288 Defined to be the solution of
290 .. math::
291 (1 - x^2)\frac{d^2}{dx^2}P_n^{(\alpha, \beta)}
292 + (\beta - \alpha - (\alpha + \beta + 2)x)
293 \frac{d}{dx}P_n^{(\alpha, \beta)}
294 + n(n + \alpha + \beta + 1)P_n^{(\alpha, \beta)} = 0
296 for :math:`\alpha, \beta > -1`; :math:`P_n^{(\alpha, \beta)}` is a
297 polynomial of degree :math:`n`.
299 Parameters
300 ----------
301 n : int
302 Degree of the polynomial.
303 alpha : float
304 Parameter, must be greater than -1.
305 beta : float
306 Parameter, must be greater than -1.
307 monic : bool, optional
308 If `True`, scale the leading coefficient to be 1. Default is
309 `False`.
311 Returns
312 -------
313 P : orthopoly1d
314 Jacobi polynomial.
316 Notes
317 -----
318 For fixed :math:`\alpha, \beta`, the polynomials
319 :math:`P_n^{(\alpha, \beta)}` are orthogonal over :math:`[-1, 1]`
320 with weight function :math:`(1 - x)^\alpha(1 + x)^\beta`.
322 References
323 ----------
324 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
325 Handbook of Mathematical Functions with Formulas,
326 Graphs, and Mathematical Tables. New York: Dover, 1972.
328 Examples
329 --------
330 The Jacobi polynomials satisfy the recurrence relation:
332 .. math::
333 P_n^{(\alpha, \beta-1)}(x) - P_n^{(\alpha-1, \beta)}(x)
334 = P_{n-1}^{(\alpha, \beta)}(x)
336 This can be verified, for example, for :math:`\alpha = \beta = 2`
337 and :math:`n = 1` over the interval :math:`[-1, 1]`:
339 >>> import numpy as np
340 >>> from scipy.special import jacobi
341 >>> x = np.arange(-1.0, 1.0, 0.01)
342 >>> np.allclose(jacobi(0, 2, 2)(x),
343 ... jacobi(1, 2, 1)(x) - jacobi(1, 1, 2)(x))
344 True
346 Plot of the Jacobi polynomial :math:`P_5^{(\alpha, -0.5)}` for
347 different values of :math:`\alpha`:
349 >>> import matplotlib.pyplot as plt
350 >>> x = np.arange(-1.0, 1.0, 0.01)
351 >>> fig, ax = plt.subplots()
352 >>> ax.set_ylim(-2.0, 2.0)
353 >>> ax.set_title(r'Jacobi polynomials $P_5^{(\alpha, -0.5)}$')
354 >>> for alpha in np.arange(0, 4, 1):
355 ... ax.plot(x, jacobi(5, alpha, -0.5)(x), label=rf'$\alpha={alpha}$')
356 >>> plt.legend(loc='best')
357 >>> plt.show()
359 """
360 if n < 0:
361 raise ValueError("n must be nonnegative.")
363 wfunc = lambda x: (1 - x)**alpha * (1 + x)**beta
364 if n == 0:
365 return orthopoly1d([], [], 1.0, 1.0, wfunc, (-1, 1), monic,
366 eval_func=np.ones_like)
367 x, w, mu = roots_jacobi(n, alpha, beta, mu=True)
368 ab1 = alpha + beta + 1.0
369 hn = 2**ab1 / (2 * n + ab1) * _gam(n + alpha + 1)
370 hn *= _gam(n + beta + 1.0) / _gam(n + 1) / _gam(n + ab1)
371 kn = _gam(2 * n + ab1) / 2.0**n / _gam(n + 1) / _gam(n + ab1)
372 # here kn = coefficient on x^n term
373 p = orthopoly1d(x, w, hn, kn, wfunc, (-1, 1), monic,
374 lambda x: _ufuncs.eval_jacobi(n, alpha, beta, x))
375 return p
377# Jacobi Polynomials shifted G_n(p,q,x)
380def roots_sh_jacobi(n, p1, q1, mu=False):
381 """Gauss-Jacobi (shifted) quadrature.
383 Compute the sample points and weights for Gauss-Jacobi (shifted)
384 quadrature. The sample points are the roots of the nth degree
385 shifted Jacobi polynomial, :math:`G^{p,q}_n(x)`. These sample
386 points and weights correctly integrate polynomials of degree
387 :math:`2n - 1` or less over the interval :math:`[0, 1]` with
388 weight function :math:`w(x) = (1 - x)^{p-q} x^{q-1}`. See 22.2.2
389 in [AS]_ for details.
391 Parameters
392 ----------
393 n : int
394 quadrature order
395 p1 : float
396 (p1 - q1) must be > -1
397 q1 : float
398 q1 must be > 0
399 mu : bool, optional
400 If True, return the sum of the weights, optional.
402 Returns
403 -------
404 x : ndarray
405 Sample points
406 w : ndarray
407 Weights
408 mu : float
409 Sum of the weights
411 See Also
412 --------
413 scipy.integrate.quadrature
414 scipy.integrate.fixed_quad
416 References
417 ----------
418 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
419 Handbook of Mathematical Functions with Formulas,
420 Graphs, and Mathematical Tables. New York: Dover, 1972.
422 """
423 if (p1-q1) <= -1 or q1 <= 0:
424 raise ValueError("(p - q) must be greater than -1, and q must be greater than 0.")
425 x, w, m = roots_jacobi(n, p1-q1, q1-1, True)
426 x = (x + 1) / 2
427 scale = 2.0**p1
428 w /= scale
429 m /= scale
430 if mu:
431 return x, w, m
432 else:
433 return x, w
436def sh_jacobi(n, p, q, monic=False):
437 r"""Shifted Jacobi polynomial.
439 Defined by
441 .. math::
443 G_n^{(p, q)}(x)
444 = \binom{2n + p - 1}{n}^{-1}P_n^{(p - q, q - 1)}(2x - 1),
446 where :math:`P_n^{(\cdot, \cdot)}` is the nth Jacobi polynomial.
448 Parameters
449 ----------
450 n : int
451 Degree of the polynomial.
452 p : float
453 Parameter, must have :math:`p > q - 1`.
454 q : float
455 Parameter, must be greater than 0.
456 monic : bool, optional
457 If `True`, scale the leading coefficient to be 1. Default is
458 `False`.
460 Returns
461 -------
462 G : orthopoly1d
463 Shifted Jacobi polynomial.
465 Notes
466 -----
467 For fixed :math:`p, q`, the polynomials :math:`G_n^{(p, q)}` are
468 orthogonal over :math:`[0, 1]` with weight function :math:`(1 -
469 x)^{p - q}x^{q - 1}`.
471 """
472 if n < 0:
473 raise ValueError("n must be nonnegative.")
475 wfunc = lambda x: (1.0 - x)**(p - q) * (x)**(q - 1.)
476 if n == 0:
477 return orthopoly1d([], [], 1.0, 1.0, wfunc, (-1, 1), monic,
478 eval_func=np.ones_like)
479 n1 = n
480 x, w = roots_sh_jacobi(n1, p, q)
481 hn = _gam(n + 1) * _gam(n + q) * _gam(n + p) * _gam(n + p - q + 1)
482 hn /= (2 * n + p) * (_gam(2 * n + p)**2)
483 # kn = 1.0 in standard form so monic is redundant. Kept for compatibility.
484 kn = 1.0
485 pp = orthopoly1d(x, w, hn, kn, wfunc=wfunc, limits=(0, 1), monic=monic,
486 eval_func=lambda x: _ufuncs.eval_sh_jacobi(n, p, q, x))
487 return pp
489# Generalized Laguerre L^(alpha)_n(x)
492def roots_genlaguerre(n, alpha, mu=False):
493 r"""Gauss-generalized Laguerre quadrature.
495 Compute the sample points and weights for Gauss-generalized
496 Laguerre quadrature. The sample points are the roots of the nth
497 degree generalized Laguerre polynomial, :math:`L^{\alpha}_n(x)`.
498 These sample points and weights correctly integrate polynomials of
499 degree :math:`2n - 1` or less over the interval :math:`[0,
500 \infty]` with weight function :math:`w(x) = x^{\alpha}
501 e^{-x}`. See 22.3.9 in [AS]_ for details.
503 Parameters
504 ----------
505 n : int
506 quadrature order
507 alpha : float
508 alpha must be > -1
509 mu : bool, optional
510 If True, return the sum of the weights, optional.
512 Returns
513 -------
514 x : ndarray
515 Sample points
516 w : ndarray
517 Weights
518 mu : float
519 Sum of the weights
521 See Also
522 --------
523 scipy.integrate.quadrature
524 scipy.integrate.fixed_quad
526 References
527 ----------
528 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
529 Handbook of Mathematical Functions with Formulas,
530 Graphs, and Mathematical Tables. New York: Dover, 1972.
532 """
533 m = int(n)
534 if n < 1 or n != m:
535 raise ValueError("n must be a positive integer.")
536 if alpha < -1:
537 raise ValueError("alpha must be greater than -1.")
539 mu0 = _ufuncs.gamma(alpha + 1)
541 if m == 1:
542 x = np.array([alpha+1.0], 'd')
543 w = np.array([mu0], 'd')
544 if mu:
545 return x, w, mu0
546 else:
547 return x, w
549 an_func = lambda k: 2 * k + alpha + 1
550 bn_func = lambda k: -np.sqrt(k * (k + alpha))
551 f = lambda n, x: _ufuncs.eval_genlaguerre(n, alpha, x)
552 df = lambda n, x: (n*_ufuncs.eval_genlaguerre(n, alpha, x)
553 - (n + alpha)*_ufuncs.eval_genlaguerre(n-1, alpha, x))/x
554 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, False, mu)
557def genlaguerre(n, alpha, monic=False):
558 r"""Generalized (associated) Laguerre polynomial.
560 Defined to be the solution of
562 .. math::
563 x\frac{d^2}{dx^2}L_n^{(\alpha)}
564 + (\alpha + 1 - x)\frac{d}{dx}L_n^{(\alpha)}
565 + nL_n^{(\alpha)} = 0,
567 where :math:`\alpha > -1`; :math:`L_n^{(\alpha)}` is a polynomial
568 of degree :math:`n`.
570 Parameters
571 ----------
572 n : int
573 Degree of the polynomial.
574 alpha : float
575 Parameter, must be greater than -1.
576 monic : bool, optional
577 If `True`, scale the leading coefficient to be 1. Default is
578 `False`.
580 Returns
581 -------
582 L : orthopoly1d
583 Generalized Laguerre polynomial.
585 Notes
586 -----
587 For fixed :math:`\alpha`, the polynomials :math:`L_n^{(\alpha)}`
588 are orthogonal over :math:`[0, \infty)` with weight function
589 :math:`e^{-x}x^\alpha`.
591 The Laguerre polynomials are the special case where :math:`\alpha
592 = 0`.
594 See Also
595 --------
596 laguerre : Laguerre polynomial.
597 hyp1f1 : confluent hypergeometric function
599 References
600 ----------
601 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
602 Handbook of Mathematical Functions with Formulas,
603 Graphs, and Mathematical Tables. New York: Dover, 1972.
605 Examples
606 --------
607 The generalized Laguerre polynomials are closely related to the confluent
608 hypergeometric function :math:`{}_1F_1`:
610 .. math::
611 L_n^{(\alpha)} = \binom{n + \alpha}{n} {}_1F_1(-n, \alpha +1, x)
613 This can be verified, for example, for :math:`n = \alpha = 3` over the
614 interval :math:`[-1, 1]`:
616 >>> import numpy as np
617 >>> from scipy.special import binom
618 >>> from scipy.special import genlaguerre
619 >>> from scipy.special import hyp1f1
620 >>> x = np.arange(-1.0, 1.0, 0.01)
621 >>> np.allclose(genlaguerre(3, 3)(x), binom(6, 3) * hyp1f1(-3, 4, x))
622 True
624 This is the plot of the generalized Laguerre polynomials
625 :math:`L_3^{(\alpha)}` for some values of :math:`\alpha`:
627 >>> import matplotlib.pyplot as plt
628 >>> x = np.arange(-4.0, 12.0, 0.01)
629 >>> fig, ax = plt.subplots()
630 >>> ax.set_ylim(-5.0, 10.0)
631 >>> ax.set_title(r'Generalized Laguerre polynomials $L_3^{\alpha}$')
632 >>> for alpha in np.arange(0, 5):
633 ... ax.plot(x, genlaguerre(3, alpha)(x), label=rf'$L_3^{(alpha)}$')
634 >>> plt.legend(loc='best')
635 >>> plt.show()
637 """
638 if alpha <= -1:
639 raise ValueError("alpha must be > -1")
640 if n < 0:
641 raise ValueError("n must be nonnegative.")
643 if n == 0:
644 n1 = n + 1
645 else:
646 n1 = n
647 x, w = roots_genlaguerre(n1, alpha)
648 wfunc = lambda x: exp(-x) * x**alpha
649 if n == 0:
650 x, w = [], []
651 hn = _gam(n + alpha + 1) / _gam(n + 1)
652 kn = (-1)**n / _gam(n + 1)
653 p = orthopoly1d(x, w, hn, kn, wfunc, (0, inf), monic,
654 lambda x: _ufuncs.eval_genlaguerre(n, alpha, x))
655 return p
657# Laguerre L_n(x)
660def roots_laguerre(n, mu=False):
661 r"""Gauss-Laguerre quadrature.
663 Compute the sample points and weights for Gauss-Laguerre
664 quadrature. The sample points are the roots of the nth degree
665 Laguerre polynomial, :math:`L_n(x)`. These sample points and
666 weights correctly integrate polynomials of degree :math:`2n - 1`
667 or less over the interval :math:`[0, \infty]` with weight function
668 :math:`w(x) = e^{-x}`. See 22.2.13 in [AS]_ for details.
670 Parameters
671 ----------
672 n : int
673 quadrature order
674 mu : bool, optional
675 If True, return the sum of the weights, optional.
677 Returns
678 -------
679 x : ndarray
680 Sample points
681 w : ndarray
682 Weights
683 mu : float
684 Sum of the weights
686 See Also
687 --------
688 scipy.integrate.quadrature
689 scipy.integrate.fixed_quad
690 numpy.polynomial.laguerre.laggauss
692 References
693 ----------
694 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
695 Handbook of Mathematical Functions with Formulas,
696 Graphs, and Mathematical Tables. New York: Dover, 1972.
698 """
699 return roots_genlaguerre(n, 0.0, mu=mu)
702def laguerre(n, monic=False):
703 r"""Laguerre polynomial.
705 Defined to be the solution of
707 .. math::
708 x\frac{d^2}{dx^2}L_n + (1 - x)\frac{d}{dx}L_n + nL_n = 0;
710 :math:`L_n` is a polynomial of degree :math:`n`.
712 Parameters
713 ----------
714 n : int
715 Degree of the polynomial.
716 monic : bool, optional
717 If `True`, scale the leading coefficient to be 1. Default is
718 `False`.
720 Returns
721 -------
722 L : orthopoly1d
723 Laguerre Polynomial.
725 Notes
726 -----
727 The polynomials :math:`L_n` are orthogonal over :math:`[0,
728 \infty)` with weight function :math:`e^{-x}`.
730 See Also
731 --------
732 genlaguerre : Generalized (associated) Laguerre polynomial.
734 References
735 ----------
736 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
737 Handbook of Mathematical Functions with Formulas,
738 Graphs, and Mathematical Tables. New York: Dover, 1972.
740 Examples
741 --------
742 The Laguerre polynomials :math:`L_n` are the special case
743 :math:`\alpha = 0` of the generalized Laguerre polynomials
744 :math:`L_n^{(\alpha)}`.
745 Let's verify it on the interval :math:`[-1, 1]`:
747 >>> import numpy as np
748 >>> from scipy.special import genlaguerre
749 >>> from scipy.special import laguerre
750 >>> x = np.arange(-1.0, 1.0, 0.01)
751 >>> np.allclose(genlaguerre(3, 0)(x), laguerre(3)(x))
752 True
754 The polynomials :math:`L_n` also satisfy the recurrence relation:
756 .. math::
757 (n + 1)L_{n+1}(x) = (2n +1 -x)L_n(x) - nL_{n-1}(x)
759 This can be easily checked on :math:`[0, 1]` for :math:`n = 3`:
761 >>> x = np.arange(0.0, 1.0, 0.01)
762 >>> np.allclose(4 * laguerre(4)(x),
763 ... (7 - x) * laguerre(3)(x) - 3 * laguerre(2)(x))
764 True
766 This is the plot of the first few Laguerre polynomials :math:`L_n`:
768 >>> import matplotlib.pyplot as plt
769 >>> x = np.arange(-1.0, 5.0, 0.01)
770 >>> fig, ax = plt.subplots()
771 >>> ax.set_ylim(-5.0, 5.0)
772 >>> ax.set_title(r'Laguerre polynomials $L_n$')
773 >>> for n in np.arange(0, 5):
774 ... ax.plot(x, laguerre(n)(x), label=rf'$L_{n}$')
775 >>> plt.legend(loc='best')
776 >>> plt.show()
778 """
779 if n < 0:
780 raise ValueError("n must be nonnegative.")
782 if n == 0:
783 n1 = n + 1
784 else:
785 n1 = n
786 x, w = roots_laguerre(n1)
787 if n == 0:
788 x, w = [], []
789 hn = 1.0
790 kn = (-1)**n / _gam(n + 1)
791 p = orthopoly1d(x, w, hn, kn, lambda x: exp(-x), (0, inf), monic,
792 lambda x: _ufuncs.eval_laguerre(n, x))
793 return p
795# Hermite 1 H_n(x)
798def roots_hermite(n, mu=False):
799 r"""Gauss-Hermite (physicist's) quadrature.
801 Compute the sample points and weights for Gauss-Hermite
802 quadrature. The sample points are the roots of the nth degree
803 Hermite polynomial, :math:`H_n(x)`. These sample points and
804 weights correctly integrate polynomials of degree :math:`2n - 1`
805 or less over the interval :math:`[-\infty, \infty]` with weight
806 function :math:`w(x) = e^{-x^2}`. See 22.2.14 in [AS]_ for
807 details.
809 Parameters
810 ----------
811 n : int
812 quadrature order
813 mu : bool, optional
814 If True, return the sum of the weights, optional.
816 Returns
817 -------
818 x : ndarray
819 Sample points
820 w : ndarray
821 Weights
822 mu : float
823 Sum of the weights
825 Notes
826 -----
827 For small n up to 150 a modified version of the Golub-Welsch
828 algorithm is used. Nodes are computed from the eigenvalue
829 problem and improved by one step of a Newton iteration.
830 The weights are computed from the well-known analytical formula.
832 For n larger than 150 an optimal asymptotic algorithm is applied
833 which computes nodes and weights in a numerically stable manner.
834 The algorithm has linear runtime making computation for very
835 large n (several thousand or more) feasible.
837 See Also
838 --------
839 scipy.integrate.quadrature
840 scipy.integrate.fixed_quad
841 numpy.polynomial.hermite.hermgauss
842 roots_hermitenorm
844 References
845 ----------
846 .. [townsend.trogdon.olver-2014]
847 Townsend, A. and Trogdon, T. and Olver, S. (2014)
848 *Fast computation of Gauss quadrature nodes and
849 weights on the whole real line*. :arXiv:`1410.5286`.
850 .. [townsend.trogdon.olver-2015]
851 Townsend, A. and Trogdon, T. and Olver, S. (2015)
852 *Fast computation of Gauss quadrature nodes and
853 weights on the whole real line*.
854 IMA Journal of Numerical Analysis
855 :doi:`10.1093/imanum/drv002`.
856 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
857 Handbook of Mathematical Functions with Formulas,
858 Graphs, and Mathematical Tables. New York: Dover, 1972.
860 """
861 m = int(n)
862 if n < 1 or n != m:
863 raise ValueError("n must be a positive integer.")
865 mu0 = np.sqrt(np.pi)
866 if n <= 150:
867 an_func = lambda k: 0.0*k
868 bn_func = lambda k: np.sqrt(k/2.0)
869 f = _ufuncs.eval_hermite
870 df = lambda n, x: 2.0 * n * _ufuncs.eval_hermite(n-1, x)
871 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu)
872 else:
873 nodes, weights = _roots_hermite_asy(m)
874 if mu:
875 return nodes, weights, mu0
876 else:
877 return nodes, weights
880def _compute_tauk(n, k, maxit=5):
881 """Helper function for Tricomi initial guesses
883 For details, see formula 3.1 in lemma 3.1 in the
884 original paper.
886 Parameters
887 ----------
888 n : int
889 Quadrature order
890 k : ndarray of type int
891 Index of roots :math:`\tau_k` to compute
892 maxit : int
893 Number of Newton maxit performed, the default
894 value of 5 is sufficient.
896 Returns
897 -------
898 tauk : ndarray
899 Roots of equation 3.1
901 See Also
902 --------
903 initial_nodes_a
904 roots_hermite_asy
905 """
906 a = n % 2 - 0.5
907 c = (4.0*floor(n/2.0) - 4.0*k + 3.0)*pi / (4.0*floor(n/2.0) + 2.0*a + 2.0)
908 f = lambda x: x - sin(x) - c
909 df = lambda x: 1.0 - cos(x)
910 xi = 0.5*pi
911 for i in range(maxit):
912 xi = xi - f(xi)/df(xi)
913 return xi
916def _initial_nodes_a(n, k):
917 r"""Tricomi initial guesses
919 Computes an initial approximation to the square of the `k`-th
920 (positive) root :math:`x_k` of the Hermite polynomial :math:`H_n`
921 of order :math:`n`. The formula is the one from lemma 3.1 in the
922 original paper. The guesses are accurate except in the region
923 near :math:`\sqrt{2n + 1}`.
925 Parameters
926 ----------
927 n : int
928 Quadrature order
929 k : ndarray of type int
930 Index of roots to compute
932 Returns
933 -------
934 xksq : ndarray
935 Square of the approximate roots
937 See Also
938 --------
939 initial_nodes
940 roots_hermite_asy
941 """
942 tauk = _compute_tauk(n, k)
943 sigk = cos(0.5*tauk)**2
944 a = n % 2 - 0.5
945 nu = 4.0*floor(n/2.0) + 2.0*a + 2.0
946 # Initial approximation of Hermite roots (square)
947 xksq = nu*sigk - 1.0/(3.0*nu) * (5.0/(4.0*(1.0-sigk)**2) - 1.0/(1.0-sigk) - 0.25)
948 return xksq
951def _initial_nodes_b(n, k):
952 r"""Gatteschi initial guesses
954 Computes an initial approximation to the square of the kth
955 (positive) root :math:`x_k` of the Hermite polynomial :math:`H_n`
956 of order :math:`n`. The formula is the one from lemma 3.2 in the
957 original paper. The guesses are accurate in the region just
958 below :math:`\sqrt{2n + 1}`.
960 Parameters
961 ----------
962 n : int
963 Quadrature order
964 k : ndarray of type int
965 Index of roots to compute
967 Returns
968 -------
969 xksq : ndarray
970 Square of the approximate root
972 See Also
973 --------
974 initial_nodes
975 roots_hermite_asy
976 """
977 a = n % 2 - 0.5
978 nu = 4.0*floor(n/2.0) + 2.0*a + 2.0
979 # Airy roots by approximation
980 ak = _specfun.airyzo(k.max(), 1)[0][::-1]
981 # Initial approximation of Hermite roots (square)
982 xksq = (nu +
983 2.0**(2.0/3.0) * ak * nu**(1.0/3.0) +
984 1.0/5.0 * 2.0**(4.0/3.0) * ak**2 * nu**(-1.0/3.0) +
985 (9.0/140.0 - 12.0/175.0 * ak**3) * nu**(-1.0) +
986 (16.0/1575.0 * ak + 92.0/7875.0 * ak**4) * 2.0**(2.0/3.0) * nu**(-5.0/3.0) -
987 (15152.0/3031875.0 * ak**5 + 1088.0/121275.0 * ak**2) * 2.0**(1.0/3.0) * nu**(-7.0/3.0))
988 return xksq
991def _initial_nodes(n):
992 """Initial guesses for the Hermite roots
994 Computes an initial approximation to the non-negative
995 roots :math:`x_k` of the Hermite polynomial :math:`H_n`
996 of order :math:`n`. The Tricomi and Gatteschi initial
997 guesses are used in the region where they are accurate.
999 Parameters
1000 ----------
1001 n : int
1002 Quadrature order
1004 Returns
1005 -------
1006 xk : ndarray
1007 Approximate roots
1009 See Also
1010 --------
1011 roots_hermite_asy
1012 """
1013 # Turnover point
1014 # linear polynomial fit to error of 10, 25, 40, ..., 1000 point rules
1015 fit = 0.49082003*n - 4.37859653
1016 turnover = around(fit).astype(int)
1017 # Compute all approximations
1018 ia = arange(1, int(floor(n*0.5)+1))
1019 ib = ia[::-1]
1020 xasq = _initial_nodes_a(n, ia[:turnover+1])
1021 xbsq = _initial_nodes_b(n, ib[turnover+1:])
1022 # Combine
1023 iv = sqrt(hstack([xasq, xbsq]))
1024 # Central node is always zero
1025 if n % 2 == 1:
1026 iv = hstack([0.0, iv])
1027 return iv
1030def _pbcf(n, theta):
1031 r"""Asymptotic series expansion of parabolic cylinder function
1033 The implementation is based on sections 3.2 and 3.3 from the
1034 original paper. Compared to the published version this code
1035 adds one more term to the asymptotic series. The detailed
1036 formulas can be found at [parabolic-asymptotics]_. The evaluation
1037 is done in a transformed variable :math:`\theta := \arccos(t)`
1038 where :math:`t := x / \mu` and :math:`\mu := \sqrt{2n + 1}`.
1040 Parameters
1041 ----------
1042 n : int
1043 Quadrature order
1044 theta : ndarray
1045 Transformed position variable
1047 Returns
1048 -------
1049 U : ndarray
1050 Value of the parabolic cylinder function :math:`U(a, \theta)`.
1051 Ud : ndarray
1052 Value of the derivative :math:`U^{\prime}(a, \theta)` of
1053 the parabolic cylinder function.
1055 See Also
1056 --------
1057 roots_hermite_asy
1059 References
1060 ----------
1061 .. [parabolic-asymptotics]
1062 https://dlmf.nist.gov/12.10#vii
1063 """
1064 st = sin(theta)
1065 ct = cos(theta)
1066 # https://dlmf.nist.gov/12.10#vii
1067 mu = 2.0*n + 1.0
1068 # https://dlmf.nist.gov/12.10#E23
1069 eta = 0.5*theta - 0.5*st*ct
1070 # https://dlmf.nist.gov/12.10#E39
1071 zeta = -(3.0*eta/2.0) ** (2.0/3.0)
1072 # https://dlmf.nist.gov/12.10#E40
1073 phi = (-zeta / st**2) ** (0.25)
1074 # Coefficients
1075 # https://dlmf.nist.gov/12.10#E43
1076 a0 = 1.0
1077 a1 = 0.10416666666666666667
1078 a2 = 0.08355034722222222222
1079 a3 = 0.12822657455632716049
1080 a4 = 0.29184902646414046425
1081 a5 = 0.88162726744375765242
1082 b0 = 1.0
1083 b1 = -0.14583333333333333333
1084 b2 = -0.09874131944444444444
1085 b3 = -0.14331205391589506173
1086 b4 = -0.31722720267841354810
1087 b5 = -0.94242914795712024914
1088 # Polynomials
1089 # https://dlmf.nist.gov/12.10#E9
1090 # https://dlmf.nist.gov/12.10#E10
1091 ctp = ct ** arange(16).reshape((-1,1))
1092 u0 = 1.0
1093 u1 = (1.0*ctp[3,:] - 6.0*ct) / 24.0
1094 u2 = (-9.0*ctp[4,:] + 249.0*ctp[2,:] + 145.0) / 1152.0
1095 u3 = (-4042.0*ctp[9,:] + 18189.0*ctp[7,:] - 28287.0*ctp[5,:] - 151995.0*ctp[3,:] - 259290.0*ct) / 414720.0
1096 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
1097 u5 = (82393456.0*ctp[15,:] - 617950920.0*ctp[13,:] + 1994971575.0*ctp[11,:] - 3630137104.0*ctp[9,:] + 4433574213.0*ctp[7,:]
1098 - 37370295816.0*ctp[5,:] - 119582875013.0*ctp[3,:] - 34009066266.0*ct) / 6688604160.0
1099 v0 = 1.0
1100 v1 = (1.0*ctp[3,:] + 6.0*ct) / 24.0
1101 v2 = (15.0*ctp[4,:] - 327.0*ctp[2,:] - 143.0) / 1152.0
1102 v3 = (-4042.0*ctp[9,:] + 18189.0*ctp[7,:] - 36387.0*ctp[5,:] + 238425.0*ctp[3,:] + 259290.0*ct) / 414720.0
1103 v4 = (-121260.0*ctp[10,:] + 551733.0*ctp[8,:] - 151958.0*ctp[6,:] - 57484425.0*ctp[4,:] - 132752238.0*ctp[2,:] - 12118727) / 39813120.0
1104 v5 = (82393456.0*ctp[15,:] - 617950920.0*ctp[13,:] + 2025529095.0*ctp[11,:] - 3750839308.0*ctp[9,:] + 3832454253.0*ctp[7,:]
1105 + 35213253348.0*ctp[5,:] + 130919230435.0*ctp[3,:] + 34009066266*ct) / 6688604160.0
1106 # Airy Evaluation (Bi and Bip unused)
1107 Ai, Aip, Bi, Bip = airy(mu**(4.0/6.0) * zeta)
1108 # Prefactor for U
1109 P = 2.0*sqrt(pi) * mu**(1.0/6.0) * phi
1110 # Terms for U
1111 # https://dlmf.nist.gov/12.10#E42
1112 phip = phi ** arange(6, 31, 6).reshape((-1,1))
1113 A0 = b0*u0
1114 A1 = (b2*u0 + phip[0,:]*b1*u1 + phip[1,:]*b0*u2) / zeta**3
1115 A2 = (b4*u0 + phip[0,:]*b3*u1 + phip[1,:]*b2*u2 + phip[2,:]*b1*u3 + phip[3,:]*b0*u4) / zeta**6
1116 B0 = -(a1*u0 + phip[0,:]*a0*u1) / zeta**2
1117 B1 = -(a3*u0 + phip[0,:]*a2*u1 + phip[1,:]*a1*u2 + phip[2,:]*a0*u3) / zeta**5
1118 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
1119 # U
1120 # https://dlmf.nist.gov/12.10#E35
1121 U = P * (Ai * (A0 + A1/mu**2.0 + A2/mu**4.0) +
1122 Aip * (B0 + B1/mu**2.0 + B2/mu**4.0) / mu**(8.0/6.0))
1123 # Prefactor for derivative of U
1124 Pd = sqrt(2.0*pi) * mu**(2.0/6.0) / phi
1125 # Terms for derivative of U
1126 # https://dlmf.nist.gov/12.10#E46
1127 C0 = -(b1*v0 + phip[0,:]*b0*v1) / zeta
1128 C1 = -(b3*v0 + phip[0,:]*b2*v1 + phip[1,:]*b1*v2 + phip[2,:]*b0*v3) / zeta**4
1129 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
1130 D0 = a0*v0
1131 D1 = (a2*v0 + phip[0,:]*a1*v1 + phip[1,:]*a0*v2) / zeta**3
1132 D2 = (a4*v0 + phip[0,:]*a3*v1 + phip[1,:]*a2*v2 + phip[2,:]*a1*v3 + phip[3,:]*a0*v4) / zeta**6
1133 # Derivative of U
1134 # https://dlmf.nist.gov/12.10#E36
1135 Ud = Pd * (Ai * (C0 + C1/mu**2.0 + C2/mu**4.0) / mu**(4.0/6.0) +
1136 Aip * (D0 + D1/mu**2.0 + D2/mu**4.0))
1137 return U, Ud
1140def _newton(n, x_initial, maxit=5):
1141 """Newton iteration for polishing the asymptotic approximation
1142 to the zeros of the Hermite polynomials.
1144 Parameters
1145 ----------
1146 n : int
1147 Quadrature order
1148 x_initial : ndarray
1149 Initial guesses for the roots
1150 maxit : int
1151 Maximal number of Newton iterations.
1152 The default 5 is sufficient, usually
1153 only one or two steps are needed.
1155 Returns
1156 -------
1157 nodes : ndarray
1158 Quadrature nodes
1159 weights : ndarray
1160 Quadrature weights
1162 See Also
1163 --------
1164 roots_hermite_asy
1165 """
1166 # Variable transformation
1167 mu = sqrt(2.0*n + 1.0)
1168 t = x_initial / mu
1169 theta = arccos(t)
1170 # Newton iteration
1171 for i in range(maxit):
1172 u, ud = _pbcf(n, theta)
1173 dtheta = u / (sqrt(2.0) * mu * sin(theta) * ud)
1174 theta = theta + dtheta
1175 if max(abs(dtheta)) < 1e-14:
1176 break
1177 # Undo variable transformation
1178 x = mu * cos(theta)
1179 # Central node is always zero
1180 if n % 2 == 1:
1181 x[0] = 0.0
1182 # Compute weights
1183 w = exp(-x**2) / (2.0*ud**2)
1184 return x, w
1187def _roots_hermite_asy(n):
1188 r"""Gauss-Hermite (physicist's) quadrature for large n.
1190 Computes the sample points and weights for Gauss-Hermite quadrature.
1191 The sample points are the roots of the nth degree Hermite polynomial,
1192 :math:`H_n(x)`. These sample points and weights correctly integrate
1193 polynomials of degree :math:`2n - 1` or less over the interval
1194 :math:`[-\infty, \infty]` with weight function :math:`f(x) = e^{-x^2}`.
1196 This method relies on asymptotic expansions which work best for n > 150.
1197 The algorithm has linear runtime making computation for very large n
1198 feasible.
1200 Parameters
1201 ----------
1202 n : int
1203 quadrature order
1205 Returns
1206 -------
1207 nodes : ndarray
1208 Quadrature nodes
1209 weights : ndarray
1210 Quadrature weights
1212 See Also
1213 --------
1214 roots_hermite
1216 References
1217 ----------
1218 .. [townsend.trogdon.olver-2014]
1219 Townsend, A. and Trogdon, T. and Olver, S. (2014)
1220 *Fast computation of Gauss quadrature nodes and
1221 weights on the whole real line*. :arXiv:`1410.5286`.
1223 .. [townsend.trogdon.olver-2015]
1224 Townsend, A. and Trogdon, T. and Olver, S. (2015)
1225 *Fast computation of Gauss quadrature nodes and
1226 weights on the whole real line*.
1227 IMA Journal of Numerical Analysis
1228 :doi:`10.1093/imanum/drv002`.
1229 """
1230 iv = _initial_nodes(n)
1231 nodes, weights = _newton(n, iv)
1232 # Combine with negative parts
1233 if n % 2 == 0:
1234 nodes = hstack([-nodes[::-1], nodes])
1235 weights = hstack([weights[::-1], weights])
1236 else:
1237 nodes = hstack([-nodes[-1:0:-1], nodes])
1238 weights = hstack([weights[-1:0:-1], weights])
1239 # Scale weights
1240 weights *= sqrt(pi) / sum(weights)
1241 return nodes, weights
1244def hermite(n, monic=False):
1245 r"""Physicist's Hermite polynomial.
1247 Defined by
1249 .. math::
1251 H_n(x) = (-1)^ne^{x^2}\frac{d^n}{dx^n}e^{-x^2};
1253 :math:`H_n` is a polynomial of degree :math:`n`.
1255 Parameters
1256 ----------
1257 n : int
1258 Degree of the polynomial.
1259 monic : bool, optional
1260 If `True`, scale the leading coefficient to be 1. Default is
1261 `False`.
1263 Returns
1264 -------
1265 H : orthopoly1d
1266 Hermite polynomial.
1268 Notes
1269 -----
1270 The polynomials :math:`H_n` are orthogonal over :math:`(-\infty,
1271 \infty)` with weight function :math:`e^{-x^2}`.
1273 Examples
1274 --------
1275 >>> from scipy import special
1276 >>> import matplotlib.pyplot as plt
1277 >>> import numpy as np
1279 >>> p_monic = special.hermite(3, monic=True)
1280 >>> p_monic
1281 poly1d([ 1. , 0. , -1.5, 0. ])
1282 >>> p_monic(1)
1283 -0.49999999999999983
1284 >>> x = np.linspace(-3, 3, 400)
1285 >>> y = p_monic(x)
1286 >>> plt.plot(x, y)
1287 >>> plt.title("Monic Hermite polynomial of degree 3")
1288 >>> plt.xlabel("x")
1289 >>> plt.ylabel("H_3(x)")
1290 >>> plt.show()
1292 """
1293 if n < 0:
1294 raise ValueError("n must be nonnegative.")
1296 if n == 0:
1297 n1 = n + 1
1298 else:
1299 n1 = n
1300 x, w = roots_hermite(n1)
1301 wfunc = lambda x: exp(-x * x)
1302 if n == 0:
1303 x, w = [], []
1304 hn = 2**n * _gam(n + 1) * sqrt(pi)
1305 kn = 2**n
1306 p = orthopoly1d(x, w, hn, kn, wfunc, (-inf, inf), monic,
1307 lambda x: _ufuncs.eval_hermite(n, x))
1308 return p
1310# Hermite 2 He_n(x)
1313def roots_hermitenorm(n, mu=False):
1314 r"""Gauss-Hermite (statistician's) quadrature.
1316 Compute the sample points and weights for Gauss-Hermite
1317 quadrature. The sample points are the roots of the nth degree
1318 Hermite polynomial, :math:`He_n(x)`. These sample points and
1319 weights correctly integrate polynomials of degree :math:`2n - 1`
1320 or less over the interval :math:`[-\infty, \infty]` with weight
1321 function :math:`w(x) = e^{-x^2/2}`. See 22.2.15 in [AS]_ for more
1322 details.
1324 Parameters
1325 ----------
1326 n : int
1327 quadrature order
1328 mu : bool, optional
1329 If True, return the sum of the weights, optional.
1331 Returns
1332 -------
1333 x : ndarray
1334 Sample points
1335 w : ndarray
1336 Weights
1337 mu : float
1338 Sum of the weights
1340 Notes
1341 -----
1342 For small n up to 150 a modified version of the Golub-Welsch
1343 algorithm is used. Nodes are computed from the eigenvalue
1344 problem and improved by one step of a Newton iteration.
1345 The weights are computed from the well-known analytical formula.
1347 For n larger than 150 an optimal asymptotic algorithm is used
1348 which computes nodes and weights in a numerical stable manner.
1349 The algorithm has linear runtime making computation for very
1350 large n (several thousand or more) feasible.
1352 See Also
1353 --------
1354 scipy.integrate.quadrature
1355 scipy.integrate.fixed_quad
1356 numpy.polynomial.hermite_e.hermegauss
1358 References
1359 ----------
1360 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1361 Handbook of Mathematical Functions with Formulas,
1362 Graphs, and Mathematical Tables. New York: Dover, 1972.
1364 """
1365 m = int(n)
1366 if n < 1 or n != m:
1367 raise ValueError("n must be a positive integer.")
1369 mu0 = np.sqrt(2.0*np.pi)
1370 if n <= 150:
1371 an_func = lambda k: 0.0*k
1372 bn_func = lambda k: np.sqrt(k)
1373 f = _ufuncs.eval_hermitenorm
1374 df = lambda n, x: n * _ufuncs.eval_hermitenorm(n-1, x)
1375 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu)
1376 else:
1377 nodes, weights = _roots_hermite_asy(m)
1378 # Transform
1379 nodes *= sqrt(2)
1380 weights *= sqrt(2)
1381 if mu:
1382 return nodes, weights, mu0
1383 else:
1384 return nodes, weights
1387def hermitenorm(n, monic=False):
1388 r"""Normalized (probabilist's) Hermite polynomial.
1390 Defined by
1392 .. math::
1394 He_n(x) = (-1)^ne^{x^2/2}\frac{d^n}{dx^n}e^{-x^2/2};
1396 :math:`He_n` is a polynomial of degree :math:`n`.
1398 Parameters
1399 ----------
1400 n : int
1401 Degree of the polynomial.
1402 monic : bool, optional
1403 If `True`, scale the leading coefficient to be 1. Default is
1404 `False`.
1406 Returns
1407 -------
1408 He : orthopoly1d
1409 Hermite polynomial.
1411 Notes
1412 -----
1414 The polynomials :math:`He_n` are orthogonal over :math:`(-\infty,
1415 \infty)` with weight function :math:`e^{-x^2/2}`.
1417 """
1418 if n < 0:
1419 raise ValueError("n must be nonnegative.")
1421 if n == 0:
1422 n1 = n + 1
1423 else:
1424 n1 = n
1425 x, w = roots_hermitenorm(n1)
1426 wfunc = lambda x: exp(-x * x / 2.0)
1427 if n == 0:
1428 x, w = [], []
1429 hn = sqrt(2 * pi) * _gam(n + 1)
1430 kn = 1.0
1431 p = orthopoly1d(x, w, hn, kn, wfunc=wfunc, limits=(-inf, inf), monic=monic,
1432 eval_func=lambda x: _ufuncs.eval_hermitenorm(n, x))
1433 return p
1435# The remainder of the polynomials can be derived from the ones above.
1437# Ultraspherical (Gegenbauer) C^(alpha)_n(x)
1440def roots_gegenbauer(n, alpha, mu=False):
1441 r"""Gauss-Gegenbauer quadrature.
1443 Compute the sample points and weights for Gauss-Gegenbauer
1444 quadrature. The sample points are the roots of the nth degree
1445 Gegenbauer polynomial, :math:`C^{\alpha}_n(x)`. These sample
1446 points and weights correctly integrate polynomials of degree
1447 :math:`2n - 1` or less over the interval :math:`[-1, 1]` with
1448 weight function :math:`w(x) = (1 - x^2)^{\alpha - 1/2}`. See
1449 22.2.3 in [AS]_ for more details.
1451 Parameters
1452 ----------
1453 n : int
1454 quadrature order
1455 alpha : float
1456 alpha must be > -0.5
1457 mu : bool, optional
1458 If True, return the sum of the weights, optional.
1460 Returns
1461 -------
1462 x : ndarray
1463 Sample points
1464 w : ndarray
1465 Weights
1466 mu : float
1467 Sum of the weights
1469 See Also
1470 --------
1471 scipy.integrate.quadrature
1472 scipy.integrate.fixed_quad
1474 References
1475 ----------
1476 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1477 Handbook of Mathematical Functions with Formulas,
1478 Graphs, and Mathematical Tables. New York: Dover, 1972.
1480 """
1481 m = int(n)
1482 if n < 1 or n != m:
1483 raise ValueError("n must be a positive integer.")
1484 if alpha < -0.5:
1485 raise ValueError("alpha must be greater than -0.5.")
1486 elif alpha == 0.0:
1487 # C(n,0,x) == 0 uniformly, however, as alpha->0, C(n,alpha,x)->T(n,x)
1488 # strictly, we should just error out here, since the roots are not
1489 # really defined, but we used to return something useful, so let's
1490 # keep doing so.
1491 return roots_chebyt(n, mu)
1493 if alpha <= 170:
1494 mu0 = (np.sqrt(np.pi) * _ufuncs.gamma(alpha + 0.5)) \
1495 / _ufuncs.gamma(alpha + 1)
1496 else:
1497 # For large alpha we use a Taylor series expansion around inf,
1498 # expressed as a 6th order polynomial of a^-1 and using Horner's
1499 # method to minimize computation and maximize precision
1500 inv_alpha = 1. / alpha
1501 coeffs = np.array([0.000207186, -0.00152206, -0.000640869,
1502 0.00488281, 0.0078125, -0.125, 1.])
1503 mu0 = coeffs[0]
1504 for term in range(1, len(coeffs)):
1505 mu0 = mu0 * inv_alpha + coeffs[term]
1506 mu0 = mu0 * np.sqrt(np.pi / alpha)
1507 an_func = lambda k: 0.0 * k
1508 bn_func = lambda k: np.sqrt(k * (k + 2 * alpha - 1)
1509 / (4 * (k + alpha) * (k + alpha - 1)))
1510 f = lambda n, x: _ufuncs.eval_gegenbauer(n, alpha, x)
1511 df = lambda n, x: ((-n*x*_ufuncs.eval_gegenbauer(n, alpha, x)
1512 + ((n + 2*alpha - 1)
1513 * _ufuncs.eval_gegenbauer(n - 1, alpha, x)))
1514 / (1 - x**2))
1515 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu)
1518def gegenbauer(n, alpha, monic=False):
1519 r"""Gegenbauer (ultraspherical) polynomial.
1521 Defined to be the solution of
1523 .. math::
1524 (1 - x^2)\frac{d^2}{dx^2}C_n^{(\alpha)}
1525 - (2\alpha + 1)x\frac{d}{dx}C_n^{(\alpha)}
1526 + n(n + 2\alpha)C_n^{(\alpha)} = 0
1528 for :math:`\alpha > -1/2`; :math:`C_n^{(\alpha)}` is a polynomial
1529 of degree :math:`n`.
1531 Parameters
1532 ----------
1533 n : int
1534 Degree of the polynomial.
1535 alpha : float
1536 Parameter, must be greater than -0.5.
1537 monic : bool, optional
1538 If `True`, scale the leading coefficient to be 1. Default is
1539 `False`.
1541 Returns
1542 -------
1543 C : orthopoly1d
1544 Gegenbauer polynomial.
1546 Notes
1547 -----
1548 The polynomials :math:`C_n^{(\alpha)}` are orthogonal over
1549 :math:`[-1,1]` with weight function :math:`(1 - x^2)^{(\alpha -
1550 1/2)}`.
1552 Examples
1553 --------
1554 >>> import numpy as np
1555 >>> from scipy import special
1556 >>> import matplotlib.pyplot as plt
1558 We can initialize a variable ``p`` as a Gegenbauer polynomial using the
1559 `gegenbauer` function and evaluate at a point ``x = 1``.
1561 >>> p = special.gegenbauer(3, 0.5, monic=False)
1562 >>> p
1563 poly1d([ 2.5, 0. , -1.5, 0. ])
1564 >>> p(1)
1565 1.0
1567 To evaluate ``p`` at various points ``x`` in the interval ``(-3, 3)``,
1568 simply pass an array ``x`` to ``p`` as follows:
1570 >>> x = np.linspace(-3, 3, 400)
1571 >>> y = p(x)
1573 We can then visualize ``x, y`` using `matplotlib.pyplot`.
1575 >>> fig, ax = plt.subplots()
1576 >>> ax.plot(x, y)
1577 >>> ax.set_title("Gegenbauer (ultraspherical) polynomial of degree 3")
1578 >>> ax.set_xlabel("x")
1579 >>> ax.set_ylabel("G_3(x)")
1580 >>> plt.show()
1582 """
1583 base = jacobi(n, alpha - 0.5, alpha - 0.5, monic=monic)
1584 if monic:
1585 return base
1586 # Abrahmowitz and Stegan 22.5.20
1587 factor = (_gam(2*alpha + n) * _gam(alpha + 0.5) /
1588 _gam(2*alpha) / _gam(alpha + 0.5 + n))
1589 base._scale(factor)
1590 base.__dict__['_eval_func'] = lambda x: _ufuncs.eval_gegenbauer(float(n),
1591 alpha, x)
1592 return base
1594# Chebyshev of the first kind: T_n(x) =
1595# n! sqrt(pi) / _gam(n+1./2)* P^(-1/2,-1/2)_n(x)
1596# Computed anew.
1599def roots_chebyt(n, mu=False):
1600 r"""Gauss-Chebyshev (first kind) quadrature.
1602 Computes the sample points and weights for Gauss-Chebyshev
1603 quadrature. The sample points are the roots of the nth degree
1604 Chebyshev polynomial of the first kind, :math:`T_n(x)`. These
1605 sample points and weights correctly integrate polynomials of
1606 degree :math:`2n - 1` or less over the interval :math:`[-1, 1]`
1607 with weight function :math:`w(x) = 1/\sqrt{1 - x^2}`. See 22.2.4
1608 in [AS]_ for more details.
1610 Parameters
1611 ----------
1612 n : int
1613 quadrature order
1614 mu : bool, optional
1615 If True, return the sum of the weights, optional.
1617 Returns
1618 -------
1619 x : ndarray
1620 Sample points
1621 w : ndarray
1622 Weights
1623 mu : float
1624 Sum of the weights
1626 See Also
1627 --------
1628 scipy.integrate.quadrature
1629 scipy.integrate.fixed_quad
1630 numpy.polynomial.chebyshev.chebgauss
1632 References
1633 ----------
1634 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1635 Handbook of Mathematical Functions with Formulas,
1636 Graphs, and Mathematical Tables. New York: Dover, 1972.
1638 """
1639 m = int(n)
1640 if n < 1 or n != m:
1641 raise ValueError('n must be a positive integer.')
1642 x = _ufuncs._sinpi(np.arange(-m + 1, m, 2) / (2*m))
1643 w = np.full_like(x, pi/m)
1644 if mu:
1645 return x, w, pi
1646 else:
1647 return x, w
1650def chebyt(n, monic=False):
1651 r"""Chebyshev polynomial of the first kind.
1653 Defined to be the solution of
1655 .. math::
1656 (1 - x^2)\frac{d^2}{dx^2}T_n - x\frac{d}{dx}T_n + n^2T_n = 0;
1658 :math:`T_n` is a polynomial of degree :math:`n`.
1660 Parameters
1661 ----------
1662 n : int
1663 Degree of the polynomial.
1664 monic : bool, optional
1665 If `True`, scale the leading coefficient to be 1. Default is
1666 `False`.
1668 Returns
1669 -------
1670 T : orthopoly1d
1671 Chebyshev polynomial of the first kind.
1673 Notes
1674 -----
1675 The polynomials :math:`T_n` are orthogonal over :math:`[-1, 1]`
1676 with weight function :math:`(1 - x^2)^{-1/2}`.
1678 See Also
1679 --------
1680 chebyu : Chebyshev polynomial of the second kind.
1682 References
1683 ----------
1684 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1685 Handbook of Mathematical Functions with Formulas,
1686 Graphs, and Mathematical Tables. New York: Dover, 1972.
1688 Examples
1689 --------
1690 Chebyshev polynomials of the first kind of order :math:`n` can
1691 be obtained as the determinant of specific :math:`n \times n`
1692 matrices. As an example we can check how the points obtained from
1693 the determinant of the following :math:`3 \times 3` matrix
1694 lay exacty on :math:`T_3`:
1696 >>> import numpy as np
1697 >>> import matplotlib.pyplot as plt
1698 >>> from scipy.linalg import det
1699 >>> from scipy.special import chebyt
1700 >>> x = np.arange(-1.0, 1.0, 0.01)
1701 >>> fig, ax = plt.subplots()
1702 >>> ax.set_ylim(-2.0, 2.0)
1703 >>> ax.set_title(r'Chebyshev polynomial $T_3$')
1704 >>> ax.plot(x, chebyt(3)(x), label=rf'$T_3$')
1705 >>> for p in np.arange(-1.0, 1.0, 0.1):
1706 ... ax.plot(p,
1707 ... det(np.array([[p, 1, 0], [1, 2*p, 1], [0, 1, 2*p]])),
1708 ... 'rx')
1709 >>> plt.legend(loc='best')
1710 >>> plt.show()
1712 They are also related to the Jacobi Polynomials
1713 :math:`P_n^{(-0.5, -0.5)}` through the relation:
1715 .. math::
1716 P_n^{(-0.5, -0.5)}(x) = \frac{1}{4^n} \binom{2n}{n} T_n(x)
1718 Let's verify it for :math:`n = 3`:
1720 >>> from scipy.special import binom
1721 >>> from scipy.special import jacobi
1722 >>> x = np.arange(-1.0, 1.0, 0.01)
1723 >>> np.allclose(jacobi(3, -0.5, -0.5)(x),
1724 ... 1/64 * binom(6, 3) * chebyt(3)(x))
1725 True
1727 We can plot the Chebyshev polynomials :math:`T_n` for some values
1728 of :math:`n`:
1730 >>> x = np.arange(-1.5, 1.5, 0.01)
1731 >>> fig, ax = plt.subplots()
1732 >>> ax.set_ylim(-4.0, 4.0)
1733 >>> ax.set_title(r'Chebyshev polynomials $T_n$')
1734 >>> for n in np.arange(2,5):
1735 ... ax.plot(x, chebyt(n)(x), label=rf'$T_n={n}$')
1736 >>> plt.legend(loc='best')
1737 >>> plt.show()
1739 """
1740 if n < 0:
1741 raise ValueError("n must be nonnegative.")
1743 wfunc = lambda x: 1.0 / sqrt(1 - x * x)
1744 if n == 0:
1745 return orthopoly1d([], [], pi, 1.0, wfunc, (-1, 1), monic,
1746 lambda x: _ufuncs.eval_chebyt(n, x))
1747 n1 = n
1748 x, w, mu = roots_chebyt(n1, mu=True)
1749 hn = pi / 2
1750 kn = 2**(n - 1)
1751 p = orthopoly1d(x, w, hn, kn, wfunc, (-1, 1), monic,
1752 lambda x: _ufuncs.eval_chebyt(n, x))
1753 return p
1755# Chebyshev of the second kind
1756# U_n(x) = (n+1)! sqrt(pi) / (2*_gam(n+3./2)) * P^(1/2,1/2)_n(x)
1759def roots_chebyu(n, mu=False):
1760 r"""Gauss-Chebyshev (second kind) quadrature.
1762 Computes the sample points and weights for Gauss-Chebyshev
1763 quadrature. The sample points are the roots of the nth degree
1764 Chebyshev polynomial of the second kind, :math:`U_n(x)`. These
1765 sample points and weights correctly integrate polynomials of
1766 degree :math:`2n - 1` or less over the interval :math:`[-1, 1]`
1767 with weight function :math:`w(x) = \sqrt{1 - x^2}`. See 22.2.5 in
1768 [AS]_ for details.
1770 Parameters
1771 ----------
1772 n : int
1773 quadrature order
1774 mu : bool, optional
1775 If True, return the sum of the weights, optional.
1777 Returns
1778 -------
1779 x : ndarray
1780 Sample points
1781 w : ndarray
1782 Weights
1783 mu : float
1784 Sum of the weights
1786 See Also
1787 --------
1788 scipy.integrate.quadrature
1789 scipy.integrate.fixed_quad
1791 References
1792 ----------
1793 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1794 Handbook of Mathematical Functions with Formulas,
1795 Graphs, and Mathematical Tables. New York: Dover, 1972.
1797 """
1798 m = int(n)
1799 if n < 1 or n != m:
1800 raise ValueError('n must be a positive integer.')
1801 t = np.arange(m, 0, -1) * pi / (m + 1)
1802 x = np.cos(t)
1803 w = pi * np.sin(t)**2 / (m + 1)
1804 if mu:
1805 return x, w, pi / 2
1806 else:
1807 return x, w
1810def chebyu(n, monic=False):
1811 r"""Chebyshev polynomial of the second kind.
1813 Defined to be the solution of
1815 .. math::
1816 (1 - x^2)\frac{d^2}{dx^2}U_n - 3x\frac{d}{dx}U_n
1817 + n(n + 2)U_n = 0;
1819 :math:`U_n` is a polynomial of degree :math:`n`.
1821 Parameters
1822 ----------
1823 n : int
1824 Degree of the polynomial.
1825 monic : bool, optional
1826 If `True`, scale the leading coefficient to be 1. Default is
1827 `False`.
1829 Returns
1830 -------
1831 U : orthopoly1d
1832 Chebyshev polynomial of the second kind.
1834 Notes
1835 -----
1836 The polynomials :math:`U_n` are orthogonal over :math:`[-1, 1]`
1837 with weight function :math:`(1 - x^2)^{1/2}`.
1839 See Also
1840 --------
1841 chebyt : Chebyshev polynomial of the first kind.
1843 References
1844 ----------
1845 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1846 Handbook of Mathematical Functions with Formulas,
1847 Graphs, and Mathematical Tables. New York: Dover, 1972.
1849 Examples
1850 --------
1851 Chebyshev polynomials of the second kind of order :math:`n` can
1852 be obtained as the determinant of specific :math:`n \times n`
1853 matrices. As an example we can check how the points obtained from
1854 the determinant of the following :math:`3 \times 3` matrix
1855 lay exacty on :math:`U_3`:
1857 >>> import numpy as np
1858 >>> import matplotlib.pyplot as plt
1859 >>> from scipy.linalg import det
1860 >>> from scipy.special import chebyu
1861 >>> x = np.arange(-1.0, 1.0, 0.01)
1862 >>> fig, ax = plt.subplots()
1863 >>> ax.set_ylim(-2.0, 2.0)
1864 >>> ax.set_title(r'Chebyshev polynomial $U_3$')
1865 >>> ax.plot(x, chebyu(3)(x), label=rf'$U_3$')
1866 >>> for p in np.arange(-1.0, 1.0, 0.1):
1867 ... ax.plot(p,
1868 ... det(np.array([[2*p, 1, 0], [1, 2*p, 1], [0, 1, 2*p]])),
1869 ... 'rx')
1870 >>> plt.legend(loc='best')
1871 >>> plt.show()
1873 They satisfy the recurrence relation:
1875 .. math::
1876 U_{2n-1}(x) = 2 T_n(x)U_{n-1}(x)
1878 where the :math:`T_n` are the Chebyshev polynomial of the first kind.
1879 Let's verify it for :math:`n = 2`:
1881 >>> from scipy.special import chebyt
1882 >>> x = np.arange(-1.0, 1.0, 0.01)
1883 >>> np.allclose(chebyu(3)(x), 2 * chebyt(2)(x) * chebyu(1)(x))
1884 True
1886 We can plot the Chebyshev polynomials :math:`U_n` for some values
1887 of :math:`n`:
1889 >>> x = np.arange(-1.0, 1.0, 0.01)
1890 >>> fig, ax = plt.subplots()
1891 >>> ax.set_ylim(-1.5, 1.5)
1892 >>> ax.set_title(r'Chebyshev polynomials $U_n$')
1893 >>> for n in np.arange(1,5):
1894 ... ax.plot(x, chebyu(n)(x), label=rf'$U_n={n}$')
1895 >>> plt.legend(loc='best')
1896 >>> plt.show()
1898 """
1899 base = jacobi(n, 0.5, 0.5, monic=monic)
1900 if monic:
1901 return base
1902 factor = sqrt(pi) / 2.0 * _gam(n + 2) / _gam(n + 1.5)
1903 base._scale(factor)
1904 return base
1906# Chebyshev of the first kind C_n(x)
1909def roots_chebyc(n, mu=False):
1910 r"""Gauss-Chebyshev (first kind) quadrature.
1912 Compute the sample points and weights for Gauss-Chebyshev
1913 quadrature. The sample points are the roots of the nth degree
1914 Chebyshev polynomial of the first kind, :math:`C_n(x)`. These
1915 sample points and weights correctly integrate polynomials of
1916 degree :math:`2n - 1` or less over the interval :math:`[-2, 2]`
1917 with weight function :math:`w(x) = 1 / \sqrt{1 - (x/2)^2}`. See
1918 22.2.6 in [AS]_ for more details.
1920 Parameters
1921 ----------
1922 n : int
1923 quadrature order
1924 mu : bool, optional
1925 If True, return the sum of the weights, optional.
1927 Returns
1928 -------
1929 x : ndarray
1930 Sample points
1931 w : ndarray
1932 Weights
1933 mu : float
1934 Sum of the weights
1936 See Also
1937 --------
1938 scipy.integrate.quadrature
1939 scipy.integrate.fixed_quad
1941 References
1942 ----------
1943 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1944 Handbook of Mathematical Functions with Formulas,
1945 Graphs, and Mathematical Tables. New York: Dover, 1972.
1947 """
1948 x, w, m = roots_chebyt(n, True)
1949 x *= 2
1950 w *= 2
1951 m *= 2
1952 if mu:
1953 return x, w, m
1954 else:
1955 return x, w
1958def chebyc(n, monic=False):
1959 r"""Chebyshev polynomial of the first kind on :math:`[-2, 2]`.
1961 Defined as :math:`C_n(x) = 2T_n(x/2)`, where :math:`T_n` is the
1962 nth Chebychev polynomial of the first kind.
1964 Parameters
1965 ----------
1966 n : int
1967 Degree of the polynomial.
1968 monic : bool, optional
1969 If `True`, scale the leading coefficient to be 1. Default is
1970 `False`.
1972 Returns
1973 -------
1974 C : orthopoly1d
1975 Chebyshev polynomial of the first kind on :math:`[-2, 2]`.
1977 Notes
1978 -----
1979 The polynomials :math:`C_n(x)` are orthogonal over :math:`[-2, 2]`
1980 with weight function :math:`1/\sqrt{1 - (x/2)^2}`.
1982 See Also
1983 --------
1984 chebyt : Chebyshev polynomial of the first kind.
1986 References
1987 ----------
1988 .. [1] Abramowitz and Stegun, "Handbook of Mathematical Functions"
1989 Section 22. National Bureau of Standards, 1972.
1991 """
1992 if n < 0:
1993 raise ValueError("n must be nonnegative.")
1995 if n == 0:
1996 n1 = n + 1
1997 else:
1998 n1 = n
1999 x, w = roots_chebyc(n1)
2000 if n == 0:
2001 x, w = [], []
2002 hn = 4 * pi * ((n == 0) + 1)
2003 kn = 1.0
2004 p = orthopoly1d(x, w, hn, kn,
2005 wfunc=lambda x: 1.0 / sqrt(1 - x * x / 4.0),
2006 limits=(-2, 2), monic=monic)
2007 if not monic:
2008 p._scale(2.0 / p(2))
2009 p.__dict__['_eval_func'] = lambda x: _ufuncs.eval_chebyc(n, x)
2010 return p
2012# Chebyshev of the second kind S_n(x)
2015def roots_chebys(n, mu=False):
2016 r"""Gauss-Chebyshev (second kind) quadrature.
2018 Compute the sample points and weights for Gauss-Chebyshev
2019 quadrature. The sample points are the roots of the nth degree
2020 Chebyshev polynomial of the second kind, :math:`S_n(x)`. These
2021 sample points and weights correctly integrate polynomials of
2022 degree :math:`2n - 1` or less over the interval :math:`[-2, 2]`
2023 with weight function :math:`w(x) = \sqrt{1 - (x/2)^2}`. See 22.2.7
2024 in [AS]_ for more details.
2026 Parameters
2027 ----------
2028 n : int
2029 quadrature order
2030 mu : bool, optional
2031 If True, return the sum of the weights, optional.
2033 Returns
2034 -------
2035 x : ndarray
2036 Sample points
2037 w : ndarray
2038 Weights
2039 mu : float
2040 Sum of the weights
2042 See Also
2043 --------
2044 scipy.integrate.quadrature
2045 scipy.integrate.fixed_quad
2047 References
2048 ----------
2049 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2050 Handbook of Mathematical Functions with Formulas,
2051 Graphs, and Mathematical Tables. New York: Dover, 1972.
2053 """
2054 x, w, m = roots_chebyu(n, True)
2055 x *= 2
2056 w *= 2
2057 m *= 2
2058 if mu:
2059 return x, w, m
2060 else:
2061 return x, w
2064def chebys(n, monic=False):
2065 r"""Chebyshev polynomial of the second kind on :math:`[-2, 2]`.
2067 Defined as :math:`S_n(x) = U_n(x/2)` where :math:`U_n` is the
2068 nth Chebychev polynomial of the second kind.
2070 Parameters
2071 ----------
2072 n : int
2073 Degree of the polynomial.
2074 monic : bool, optional
2075 If `True`, scale the leading coefficient to be 1. Default is
2076 `False`.
2078 Returns
2079 -------
2080 S : orthopoly1d
2081 Chebyshev polynomial of the second kind on :math:`[-2, 2]`.
2083 Notes
2084 -----
2085 The polynomials :math:`S_n(x)` are orthogonal over :math:`[-2, 2]`
2086 with weight function :math:`\sqrt{1 - (x/2)}^2`.
2088 See Also
2089 --------
2090 chebyu : Chebyshev polynomial of the second kind
2092 References
2093 ----------
2094 .. [1] Abramowitz and Stegun, "Handbook of Mathematical Functions"
2095 Section 22. National Bureau of Standards, 1972.
2097 """
2098 if n < 0:
2099 raise ValueError("n must be nonnegative.")
2101 if n == 0:
2102 n1 = n + 1
2103 else:
2104 n1 = n
2105 x, w = roots_chebys(n1)
2106 if n == 0:
2107 x, w = [], []
2108 hn = pi
2109 kn = 1.0
2110 p = orthopoly1d(x, w, hn, kn,
2111 wfunc=lambda x: sqrt(1 - x * x / 4.0),
2112 limits=(-2, 2), monic=monic)
2113 if not monic:
2114 factor = (n + 1.0) / p(2)
2115 p._scale(factor)
2116 p.__dict__['_eval_func'] = lambda x: _ufuncs.eval_chebys(n, x)
2117 return p
2119# Shifted Chebyshev of the first kind T^*_n(x)
2122def roots_sh_chebyt(n, mu=False):
2123 r"""Gauss-Chebyshev (first kind, shifted) quadrature.
2125 Compute the sample points and weights for Gauss-Chebyshev
2126 quadrature. The sample points are the roots of the nth degree
2127 shifted Chebyshev polynomial of the first kind, :math:`T_n(x)`.
2128 These sample points and weights correctly integrate polynomials of
2129 degree :math:`2n - 1` or less over the interval :math:`[0, 1]`
2130 with weight function :math:`w(x) = 1/\sqrt{x - x^2}`. See 22.2.8
2131 in [AS]_ for more details.
2133 Parameters
2134 ----------
2135 n : int
2136 quadrature order
2137 mu : bool, optional
2138 If True, return the sum of the weights, optional.
2140 Returns
2141 -------
2142 x : ndarray
2143 Sample points
2144 w : ndarray
2145 Weights
2146 mu : float
2147 Sum of the weights
2149 See Also
2150 --------
2151 scipy.integrate.quadrature
2152 scipy.integrate.fixed_quad
2154 References
2155 ----------
2156 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2157 Handbook of Mathematical Functions with Formulas,
2158 Graphs, and Mathematical Tables. New York: Dover, 1972.
2160 """
2161 xw = roots_chebyt(n, mu)
2162 return ((xw[0] + 1) / 2,) + xw[1:]
2165def sh_chebyt(n, monic=False):
2166 r"""Shifted Chebyshev polynomial of the first kind.
2168 Defined as :math:`T^*_n(x) = T_n(2x - 1)` for :math:`T_n` the nth
2169 Chebyshev polynomial of the first kind.
2171 Parameters
2172 ----------
2173 n : int
2174 Degree of the polynomial.
2175 monic : bool, optional
2176 If `True`, scale the leading coefficient to be 1. Default is
2177 `False`.
2179 Returns
2180 -------
2181 T : orthopoly1d
2182 Shifted Chebyshev polynomial of the first kind.
2184 Notes
2185 -----
2186 The polynomials :math:`T^*_n` are orthogonal over :math:`[0, 1]`
2187 with weight function :math:`(x - x^2)^{-1/2}`.
2189 """
2190 base = sh_jacobi(n, 0.0, 0.5, monic=monic)
2191 if monic:
2192 return base
2193 if n > 0:
2194 factor = 4**n / 2.0
2195 else:
2196 factor = 1.0
2197 base._scale(factor)
2198 return base
2201# Shifted Chebyshev of the second kind U^*_n(x)
2202def roots_sh_chebyu(n, mu=False):
2203 r"""Gauss-Chebyshev (second kind, shifted) quadrature.
2205 Computes the sample points and weights for Gauss-Chebyshev
2206 quadrature. The sample points are the roots of the nth degree
2207 shifted Chebyshev polynomial of the second kind, :math:`U_n(x)`.
2208 These sample points and weights correctly integrate polynomials of
2209 degree :math:`2n - 1` or less over the interval :math:`[0, 1]`
2210 with weight function :math:`w(x) = \sqrt{x - x^2}`. See 22.2.9 in
2211 [AS]_ for more details.
2213 Parameters
2214 ----------
2215 n : int
2216 quadrature order
2217 mu : bool, optional
2218 If True, return the sum of the weights, optional.
2220 Returns
2221 -------
2222 x : ndarray
2223 Sample points
2224 w : ndarray
2225 Weights
2226 mu : float
2227 Sum of the weights
2229 See Also
2230 --------
2231 scipy.integrate.quadrature
2232 scipy.integrate.fixed_quad
2234 References
2235 ----------
2236 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2237 Handbook of Mathematical Functions with Formulas,
2238 Graphs, and Mathematical Tables. New York: Dover, 1972.
2240 """
2241 x, w, m = roots_chebyu(n, True)
2242 x = (x + 1) / 2
2243 m_us = _ufuncs.beta(1.5, 1.5)
2244 w *= m_us / m
2245 if mu:
2246 return x, w, m_us
2247 else:
2248 return x, w
2251def sh_chebyu(n, monic=False):
2252 r"""Shifted Chebyshev polynomial of the second kind.
2254 Defined as :math:`U^*_n(x) = U_n(2x - 1)` for :math:`U_n` the nth
2255 Chebyshev polynomial of the second kind.
2257 Parameters
2258 ----------
2259 n : int
2260 Degree of the polynomial.
2261 monic : bool, optional
2262 If `True`, scale the leading coefficient to be 1. Default is
2263 `False`.
2265 Returns
2266 -------
2267 U : orthopoly1d
2268 Shifted Chebyshev polynomial of the second kind.
2270 Notes
2271 -----
2272 The polynomials :math:`U^*_n` are orthogonal over :math:`[0, 1]`
2273 with weight function :math:`(x - x^2)^{1/2}`.
2275 """
2276 base = sh_jacobi(n, 2.0, 1.5, monic=monic)
2277 if monic:
2278 return base
2279 factor = 4**n
2280 base._scale(factor)
2281 return base
2283# Legendre
2286def roots_legendre(n, mu=False):
2287 r"""Gauss-Legendre quadrature.
2289 Compute the sample points and weights for Gauss-Legendre
2290 quadrature [GL]_. The sample points are the roots of the nth degree
2291 Legendre polynomial :math:`P_n(x)`. These sample points and
2292 weights correctly integrate polynomials of degree :math:`2n - 1`
2293 or less over the interval :math:`[-1, 1]` with weight function
2294 :math:`w(x) = 1`. See 2.2.10 in [AS]_ for more details.
2296 Parameters
2297 ----------
2298 n : int
2299 quadrature order
2300 mu : bool, optional
2301 If True, return the sum of the weights, optional.
2303 Returns
2304 -------
2305 x : ndarray
2306 Sample points
2307 w : ndarray
2308 Weights
2309 mu : float
2310 Sum of the weights
2312 See Also
2313 --------
2314 scipy.integrate.quadrature
2315 scipy.integrate.fixed_quad
2316 numpy.polynomial.legendre.leggauss
2318 References
2319 ----------
2320 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2321 Handbook of Mathematical Functions with Formulas,
2322 Graphs, and Mathematical Tables. New York: Dover, 1972.
2323 .. [GL] Gauss-Legendre quadrature, Wikipedia,
2324 https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature
2326 Examples
2327 --------
2328 >>> import numpy as np
2329 >>> from scipy.special import roots_legendre, eval_legendre
2330 >>> roots, weights = roots_legendre(9)
2332 ``roots`` holds the roots, and ``weights`` holds the weights for
2333 Gauss-Legendre quadrature.
2335 >>> roots
2336 array([-0.96816024, -0.83603111, -0.61337143, -0.32425342, 0. ,
2337 0.32425342, 0.61337143, 0.83603111, 0.96816024])
2338 >>> weights
2339 array([0.08127439, 0.18064816, 0.2606107 , 0.31234708, 0.33023936,
2340 0.31234708, 0.2606107 , 0.18064816, 0.08127439])
2342 Verify that we have the roots by evaluating the degree 9 Legendre
2343 polynomial at ``roots``. All the values are approximately zero:
2345 >>> eval_legendre(9, roots)
2346 array([-8.88178420e-16, -2.22044605e-16, 1.11022302e-16, 1.11022302e-16,
2347 0.00000000e+00, -5.55111512e-17, -1.94289029e-16, 1.38777878e-16,
2348 -8.32667268e-17])
2350 Here we'll show how the above values can be used to estimate the
2351 integral from 1 to 2 of f(t) = t + 1/t with Gauss-Legendre
2352 quadrature [GL]_. First define the function and the integration
2353 limits.
2355 >>> def f(t):
2356 ... return t + 1/t
2357 ...
2358 >>> a = 1
2359 >>> b = 2
2361 We'll use ``integral(f(t), t=a, t=b)`` to denote the definite integral
2362 of f from t=a to t=b. The sample points in ``roots`` are from the
2363 interval [-1, 1], so we'll rewrite the integral with the simple change
2364 of variable::
2366 x = 2/(b - a) * t - (a + b)/(b - a)
2368 with inverse::
2370 t = (b - a)/2 * x + (a + 2)/2
2372 Then::
2374 integral(f(t), a, b) =
2375 (b - a)/2 * integral(f((b-a)/2*x + (a+b)/2), x=-1, x=1)
2377 We can approximate the latter integral with the values returned
2378 by `roots_legendre`.
2380 Map the roots computed above from [-1, 1] to [a, b].
2382 >>> t = (b - a)/2 * roots + (a + b)/2
2384 Approximate the integral as the weighted sum of the function values.
2386 >>> (b - a)/2 * f(t).dot(weights)
2387 2.1931471805599276
2389 Compare that to the exact result, which is 3/2 + log(2):
2391 >>> 1.5 + np.log(2)
2392 2.1931471805599454
2394 """
2395 m = int(n)
2396 if n < 1 or n != m:
2397 raise ValueError("n must be a positive integer.")
2399 mu0 = 2.0
2400 an_func = lambda k: 0.0 * k
2401 bn_func = lambda k: k * np.sqrt(1.0 / (4 * k * k - 1))
2402 f = _ufuncs.eval_legendre
2403 df = lambda n, x: (-n*x*_ufuncs.eval_legendre(n, x)
2404 + n*_ufuncs.eval_legendre(n-1, x))/(1-x**2)
2405 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu)
2408def legendre(n, monic=False):
2409 r"""Legendre polynomial.
2411 Defined to be the solution of
2413 .. math::
2414 \frac{d}{dx}\left[(1 - x^2)\frac{d}{dx}P_n(x)\right]
2415 + n(n + 1)P_n(x) = 0;
2417 :math:`P_n(x)` is a polynomial of degree :math:`n`.
2419 Parameters
2420 ----------
2421 n : int
2422 Degree of the polynomial.
2423 monic : bool, optional
2424 If `True`, scale the leading coefficient to be 1. Default is
2425 `False`.
2427 Returns
2428 -------
2429 P : orthopoly1d
2430 Legendre polynomial.
2432 Notes
2433 -----
2434 The polynomials :math:`P_n` are orthogonal over :math:`[-1, 1]`
2435 with weight function 1.
2437 Examples
2438 --------
2439 Generate the 3rd-order Legendre polynomial 1/2*(5x^3 + 0x^2 - 3x + 0):
2441 >>> from scipy.special import legendre
2442 >>> legendre(3)
2443 poly1d([ 2.5, 0. , -1.5, 0. ])
2445 """
2446 if n < 0:
2447 raise ValueError("n must be nonnegative.")
2449 if n == 0:
2450 n1 = n + 1
2451 else:
2452 n1 = n
2453 x, w = roots_legendre(n1)
2454 if n == 0:
2455 x, w = [], []
2456 hn = 2.0 / (2 * n + 1)
2457 kn = _gam(2 * n + 1) / _gam(n + 1)**2 / 2.0**n
2458 p = orthopoly1d(x, w, hn, kn, wfunc=lambda x: 1.0, limits=(-1, 1),
2459 monic=monic,
2460 eval_func=lambda x: _ufuncs.eval_legendre(n, x))
2461 return p
2463# Shifted Legendre P^*_n(x)
2466def roots_sh_legendre(n, mu=False):
2467 r"""Gauss-Legendre (shifted) quadrature.
2469 Compute the sample points and weights for Gauss-Legendre
2470 quadrature. The sample points are the roots of the nth degree
2471 shifted Legendre polynomial :math:`P^*_n(x)`. These sample points
2472 and weights correctly integrate polynomials of degree :math:`2n -
2473 1` or less over the interval :math:`[0, 1]` with weight function
2474 :math:`w(x) = 1.0`. See 2.2.11 in [AS]_ for details.
2476 Parameters
2477 ----------
2478 n : int
2479 quadrature order
2480 mu : bool, optional
2481 If True, return the sum of the weights, optional.
2483 Returns
2484 -------
2485 x : ndarray
2486 Sample points
2487 w : ndarray
2488 Weights
2489 mu : float
2490 Sum of the weights
2492 See Also
2493 --------
2494 scipy.integrate.quadrature
2495 scipy.integrate.fixed_quad
2497 References
2498 ----------
2499 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2500 Handbook of Mathematical Functions with Formulas,
2501 Graphs, and Mathematical Tables. New York: Dover, 1972.
2503 """
2504 x, w = roots_legendre(n)
2505 x = (x + 1) / 2
2506 w /= 2
2507 if mu:
2508 return x, w, 1.0
2509 else:
2510 return x, w
2513def sh_legendre(n, monic=False):
2514 r"""Shifted Legendre polynomial.
2516 Defined as :math:`P^*_n(x) = P_n(2x - 1)` for :math:`P_n` the nth
2517 Legendre polynomial.
2519 Parameters
2520 ----------
2521 n : int
2522 Degree of the polynomial.
2523 monic : bool, optional
2524 If `True`, scale the leading coefficient to be 1. Default is
2525 `False`.
2527 Returns
2528 -------
2529 P : orthopoly1d
2530 Shifted Legendre polynomial.
2532 Notes
2533 -----
2534 The polynomials :math:`P^*_n` are orthogonal over :math:`[0, 1]`
2535 with weight function 1.
2537 """
2538 if n < 0:
2539 raise ValueError("n must be nonnegative.")
2541 wfunc = lambda x: 0.0 * x + 1.0
2542 if n == 0:
2543 return orthopoly1d([], [], 1.0, 1.0, wfunc, (0, 1), monic,
2544 lambda x: _ufuncs.eval_sh_legendre(n, x))
2545 x, w = roots_sh_legendre(n)
2546 hn = 1.0 / (2 * n + 1.0)
2547 kn = _gam(2 * n + 1) / _gam(n + 1)**2
2548 p = orthopoly1d(x, w, hn, kn, wfunc, limits=(0, 1), monic=monic,
2549 eval_func=lambda x: _ufuncs.eval_sh_legendre(n, x))
2550 return p
2553# Make the old root function names an alias for the new ones
2554_modattrs = globals()
2555for newfun, oldfun in _rootfuns_map.items():
2556 _modattrs[oldfun] = _modattrs[newfun]
2557 __all__.append(oldfun)