Coverage for /usr/lib/python3/dist-packages/scipy/special/_orthogonal.py: 10%
559 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
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.
85# There is no .pyi file for _specfun
86from . import _specfun # type: ignore
87from . import _ufuncs
88_gam = _ufuncs.gamma
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 def eval_func(x):
127 return evf(x) / knn
128 mu = mu / abs(kn)
129 kn = 1.0
131 # compute coefficients from roots, then scale
132 poly = np.poly1d(roots, r=True)
133 np.poly1d.__init__(self, poly.coeffs * float(kn))
135 self.weights = np.array(list(zip(roots, weights, equiv_weights)))
136 self.weight_func = wfunc
137 self.limits = limits
138 self.normcoef = mu
140 # Note: eval_func will be discarded on arithmetic
141 self._eval_func = eval_func
143 def __call__(self, v):
144 if self._eval_func and not isinstance(v, np.poly1d):
145 return self._eval_func(v)
146 else:
147 return np.poly1d.__call__(self, v)
149 def _scale(self, p):
150 if p == 1.0:
151 return
152 self._coeffs *= p
154 evf = self._eval_func
155 if evf:
156 self._eval_func = lambda x: evf(x) * p
157 self.normcoef *= p
160def _gen_roots_and_weights(n, mu0, an_func, bn_func, f, df, symmetrize, mu):
161 """[x,w] = gen_roots_and_weights(n,an_func,sqrt_bn_func,mu)
163 Returns the roots (x) of an nth order orthogonal polynomial,
164 and weights (w) to use in appropriate Gaussian quadrature with that
165 orthogonal polynomial.
167 The polynomials have the recurrence relation
168 P_n+1(x) = (x - A_n) P_n(x) - B_n P_n-1(x)
170 an_func(n) should return A_n
171 sqrt_bn_func(n) should return sqrt(B_n)
172 mu ( = h_0 ) is the integral of the weight over the orthogonal
173 interval
174 """
175 k = np.arange(n, dtype='d')
176 c = np.zeros((2, n))
177 c[0,1:] = bn_func(k[1:])
178 c[1,:] = an_func(k)
179 x = linalg.eigvals_banded(c, overwrite_a_band=True)
181 # improve roots by one application of Newton's method
182 y = f(n, x)
183 dy = df(n, x)
184 x -= y/dy
186 # fm and dy may contain very large/small values, so we
187 # log-normalize them to maintain precision in the product fm*dy
188 fm = f(n-1, x)
189 log_fm = np.log(np.abs(fm))
190 log_dy = np.log(np.abs(dy))
191 fm /= np.exp((log_fm.max() + log_fm.min()) / 2.)
192 dy /= np.exp((log_dy.max() + log_dy.min()) / 2.)
193 w = 1.0 / (fm * dy)
195 if symmetrize:
196 w = (w + w[::-1]) / 2
197 x = (x - x[::-1]) / 2
199 w *= mu0 / w.sum()
201 if mu:
202 return x, w, mu0
203 else:
204 return x, w
206# Jacobi Polynomials 1 P^(alpha,beta)_n(x)
209def roots_jacobi(n, alpha, beta, mu=False):
210 r"""Gauss-Jacobi quadrature.
212 Compute the sample points and weights for Gauss-Jacobi
213 quadrature. The sample points are the roots of the nth degree
214 Jacobi polynomial, :math:`P^{\alpha, \beta}_n(x)`. These sample
215 points and weights correctly integrate polynomials of degree
216 :math:`2n - 1` or less over the interval :math:`[-1, 1]` with
217 weight function :math:`w(x) = (1 - x)^{\alpha} (1 +
218 x)^{\beta}`. See 22.2.1 in [AS]_ for details.
220 Parameters
221 ----------
222 n : int
223 quadrature order
224 alpha : float
225 alpha must be > -1
226 beta : float
227 beta must be > -1
228 mu : bool, optional
229 If True, return the sum of the weights, optional.
231 Returns
232 -------
233 x : ndarray
234 Sample points
235 w : ndarray
236 Weights
237 mu : float
238 Sum of the weights
240 See Also
241 --------
242 scipy.integrate.quadrature
243 scipy.integrate.fixed_quad
245 References
246 ----------
247 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
248 Handbook of Mathematical Functions with Formulas,
249 Graphs, and Mathematical Tables. New York: Dover, 1972.
251 """
252 m = int(n)
253 if n < 1 or n != m:
254 raise ValueError("n must be a positive integer.")
255 if alpha <= -1 or beta <= -1:
256 raise ValueError("alpha and beta must be greater than -1.")
258 if alpha == 0.0 and beta == 0.0:
259 return roots_legendre(m, mu)
260 if alpha == beta:
261 return roots_gegenbauer(m, alpha+0.5, mu)
263 if (alpha + beta) <= 1000:
264 mu0 = 2.0**(alpha+beta+1) * _ufuncs.beta(alpha+1, beta+1)
265 else:
266 # Avoid overflows in pow and beta for very large parameters
267 mu0 = np.exp((alpha + beta + 1) * np.log(2.0)
268 + _ufuncs.betaln(alpha+1, beta+1))
269 a = alpha
270 b = beta
271 if a + b == 0.0:
272 def an_func(k):
273 return np.where(k == 0, (b - a) / (2 + a + b), 0.0)
274 else:
275 def an_func(k):
276 return np.where(k == 0, (b - a) / (2 + a + b), (b * b - a * a) / ((2.0 * k + a + b) * (2.0 * k + a + b + 2)))
278 def bn_func(k):
279 return 2.0 / (2.0 * k + a + b) * np.sqrt((k + a) * (k + b) / (2 * k + a + b + 1)) * np.where(k == 1, 1.0, np.sqrt(k * (k + a + b) / (2.0 * k + a + b - 1)))
281 def f(n, x):
282 return _ufuncs.eval_jacobi(n, a, b, x)
283 def df(n, x):
284 return 0.5 * (n + a + b + 1) * _ufuncs.eval_jacobi(n - 1, a + 1, b + 1, x)
285 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, False, mu)
288def jacobi(n, alpha, beta, monic=False):
289 r"""Jacobi polynomial.
291 Defined to be the solution of
293 .. math::
294 (1 - x^2)\frac{d^2}{dx^2}P_n^{(\alpha, \beta)}
295 + (\beta - \alpha - (\alpha + \beta + 2)x)
296 \frac{d}{dx}P_n^{(\alpha, \beta)}
297 + n(n + \alpha + \beta + 1)P_n^{(\alpha, \beta)} = 0
299 for :math:`\alpha, \beta > -1`; :math:`P_n^{(\alpha, \beta)}` is a
300 polynomial of degree :math:`n`.
302 Parameters
303 ----------
304 n : int
305 Degree of the polynomial.
306 alpha : float
307 Parameter, must be greater than -1.
308 beta : float
309 Parameter, must be greater than -1.
310 monic : bool, optional
311 If `True`, scale the leading coefficient to be 1. Default is
312 `False`.
314 Returns
315 -------
316 P : orthopoly1d
317 Jacobi polynomial.
319 Notes
320 -----
321 For fixed :math:`\alpha, \beta`, the polynomials
322 :math:`P_n^{(\alpha, \beta)}` are orthogonal over :math:`[-1, 1]`
323 with weight function :math:`(1 - x)^\alpha(1 + x)^\beta`.
325 References
326 ----------
327 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
328 Handbook of Mathematical Functions with Formulas,
329 Graphs, and Mathematical Tables. New York: Dover, 1972.
331 Examples
332 --------
333 The Jacobi polynomials satisfy the recurrence relation:
335 .. math::
336 P_n^{(\alpha, \beta-1)}(x) - P_n^{(\alpha-1, \beta)}(x)
337 = P_{n-1}^{(\alpha, \beta)}(x)
339 This can be verified, for example, for :math:`\alpha = \beta = 2`
340 and :math:`n = 1` over the interval :math:`[-1, 1]`:
342 >>> import numpy as np
343 >>> from scipy.special import jacobi
344 >>> x = np.arange(-1.0, 1.0, 0.01)
345 >>> np.allclose(jacobi(0, 2, 2)(x),
346 ... jacobi(1, 2, 1)(x) - jacobi(1, 1, 2)(x))
347 True
349 Plot of the Jacobi polynomial :math:`P_5^{(\alpha, -0.5)}` for
350 different values of :math:`\alpha`:
352 >>> import matplotlib.pyplot as plt
353 >>> x = np.arange(-1.0, 1.0, 0.01)
354 >>> fig, ax = plt.subplots()
355 >>> ax.set_ylim(-2.0, 2.0)
356 >>> ax.set_title(r'Jacobi polynomials $P_5^{(\alpha, -0.5)}$')
357 >>> for alpha in np.arange(0, 4, 1):
358 ... ax.plot(x, jacobi(5, alpha, -0.5)(x), label=rf'$\alpha={alpha}$')
359 >>> plt.legend(loc='best')
360 >>> plt.show()
362 """
363 if n < 0:
364 raise ValueError("n must be nonnegative.")
366 def wfunc(x):
367 return (1 - x) ** alpha * (1 + x) ** beta
368 if n == 0:
369 return orthopoly1d([], [], 1.0, 1.0, wfunc, (-1, 1), monic,
370 eval_func=np.ones_like)
371 x, w, mu = roots_jacobi(n, alpha, beta, mu=True)
372 ab1 = alpha + beta + 1.0
373 hn = 2**ab1 / (2 * n + ab1) * _gam(n + alpha + 1)
374 hn *= _gam(n + beta + 1.0) / _gam(n + 1) / _gam(n + ab1)
375 kn = _gam(2 * n + ab1) / 2.0**n / _gam(n + 1) / _gam(n + ab1)
376 # here kn = coefficient on x^n term
377 p = orthopoly1d(x, w, hn, kn, wfunc, (-1, 1), monic,
378 lambda x: _ufuncs.eval_jacobi(n, alpha, beta, x))
379 return p
381# Jacobi Polynomials shifted G_n(p,q,x)
384def roots_sh_jacobi(n, p1, q1, mu=False):
385 """Gauss-Jacobi (shifted) quadrature.
387 Compute the sample points and weights for Gauss-Jacobi (shifted)
388 quadrature. The sample points are the roots of the nth degree
389 shifted Jacobi polynomial, :math:`G^{p,q}_n(x)`. These sample
390 points and weights correctly integrate polynomials of degree
391 :math:`2n - 1` or less over the interval :math:`[0, 1]` with
392 weight function :math:`w(x) = (1 - x)^{p-q} x^{q-1}`. See 22.2.2
393 in [AS]_ for details.
395 Parameters
396 ----------
397 n : int
398 quadrature order
399 p1 : float
400 (p1 - q1) must be > -1
401 q1 : float
402 q1 must be > 0
403 mu : bool, optional
404 If True, return the sum of the weights, optional.
406 Returns
407 -------
408 x : ndarray
409 Sample points
410 w : ndarray
411 Weights
412 mu : float
413 Sum of the weights
415 See Also
416 --------
417 scipy.integrate.quadrature
418 scipy.integrate.fixed_quad
420 References
421 ----------
422 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
423 Handbook of Mathematical Functions with Formulas,
424 Graphs, and Mathematical Tables. New York: Dover, 1972.
426 """
427 if (p1-q1) <= -1 or q1 <= 0:
428 raise ValueError("(p - q) must be greater than -1, and q must be greater than 0.")
429 x, w, m = roots_jacobi(n, p1-q1, q1-1, True)
430 x = (x + 1) / 2
431 scale = 2.0**p1
432 w /= scale
433 m /= scale
434 if mu:
435 return x, w, m
436 else:
437 return x, w
440def sh_jacobi(n, p, q, monic=False):
441 r"""Shifted Jacobi polynomial.
443 Defined by
445 .. math::
447 G_n^{(p, q)}(x)
448 = \binom{2n + p - 1}{n}^{-1}P_n^{(p - q, q - 1)}(2x - 1),
450 where :math:`P_n^{(\cdot, \cdot)}` is the nth Jacobi polynomial.
452 Parameters
453 ----------
454 n : int
455 Degree of the polynomial.
456 p : float
457 Parameter, must have :math:`p > q - 1`.
458 q : float
459 Parameter, must be greater than 0.
460 monic : bool, optional
461 If `True`, scale the leading coefficient to be 1. Default is
462 `False`.
464 Returns
465 -------
466 G : orthopoly1d
467 Shifted Jacobi polynomial.
469 Notes
470 -----
471 For fixed :math:`p, q`, the polynomials :math:`G_n^{(p, q)}` are
472 orthogonal over :math:`[0, 1]` with weight function :math:`(1 -
473 x)^{p - q}x^{q - 1}`.
475 """
476 if n < 0:
477 raise ValueError("n must be nonnegative.")
479 def wfunc(x):
480 return (1.0 - x) ** (p - q) * x ** (q - 1.0)
481 if n == 0:
482 return orthopoly1d([], [], 1.0, 1.0, wfunc, (-1, 1), monic,
483 eval_func=np.ones_like)
484 n1 = n
485 x, w = roots_sh_jacobi(n1, p, q)
486 hn = _gam(n + 1) * _gam(n + q) * _gam(n + p) * _gam(n + p - q + 1)
487 hn /= (2 * n + p) * (_gam(2 * n + p)**2)
488 # kn = 1.0 in standard form so monic is redundant. Kept for compatibility.
489 kn = 1.0
490 pp = orthopoly1d(x, w, hn, kn, wfunc=wfunc, limits=(0, 1), monic=monic,
491 eval_func=lambda x: _ufuncs.eval_sh_jacobi(n, p, q, x))
492 return pp
494# Generalized Laguerre L^(alpha)_n(x)
497def roots_genlaguerre(n, alpha, mu=False):
498 r"""Gauss-generalized Laguerre quadrature.
500 Compute the sample points and weights for Gauss-generalized
501 Laguerre quadrature. The sample points are the roots of the nth
502 degree generalized Laguerre polynomial, :math:`L^{\alpha}_n(x)`.
503 These sample points and weights correctly integrate polynomials of
504 degree :math:`2n - 1` or less over the interval :math:`[0,
505 \infty]` with weight function :math:`w(x) = x^{\alpha}
506 e^{-x}`. See 22.3.9 in [AS]_ for details.
508 Parameters
509 ----------
510 n : int
511 quadrature order
512 alpha : float
513 alpha must be > -1
514 mu : bool, optional
515 If True, return the sum of the weights, optional.
517 Returns
518 -------
519 x : ndarray
520 Sample points
521 w : ndarray
522 Weights
523 mu : float
524 Sum of the weights
526 See Also
527 --------
528 scipy.integrate.quadrature
529 scipy.integrate.fixed_quad
531 References
532 ----------
533 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
534 Handbook of Mathematical Functions with Formulas,
535 Graphs, and Mathematical Tables. New York: Dover, 1972.
537 """
538 m = int(n)
539 if n < 1 or n != m:
540 raise ValueError("n must be a positive integer.")
541 if alpha < -1:
542 raise ValueError("alpha must be greater than -1.")
544 mu0 = _ufuncs.gamma(alpha + 1)
546 if m == 1:
547 x = np.array([alpha+1.0], 'd')
548 w = np.array([mu0], 'd')
549 if mu:
550 return x, w, mu0
551 else:
552 return x, w
554 def an_func(k):
555 return 2 * k + alpha + 1
556 def bn_func(k):
557 return -np.sqrt(k * (k + alpha))
558 def f(n, x):
559 return _ufuncs.eval_genlaguerre(n, alpha, x)
560 def df(n, x):
561 return (n * _ufuncs.eval_genlaguerre(n, alpha, x) - (n + alpha) * _ufuncs.eval_genlaguerre(n - 1, alpha, x)) / x
562 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, False, mu)
565def genlaguerre(n, alpha, monic=False):
566 r"""Generalized (associated) Laguerre polynomial.
568 Defined to be the solution of
570 .. math::
571 x\frac{d^2}{dx^2}L_n^{(\alpha)}
572 + (\alpha + 1 - x)\frac{d}{dx}L_n^{(\alpha)}
573 + nL_n^{(\alpha)} = 0,
575 where :math:`\alpha > -1`; :math:`L_n^{(\alpha)}` is a polynomial
576 of degree :math:`n`.
578 Parameters
579 ----------
580 n : int
581 Degree of the polynomial.
582 alpha : float
583 Parameter, must be greater than -1.
584 monic : bool, optional
585 If `True`, scale the leading coefficient to be 1. Default is
586 `False`.
588 Returns
589 -------
590 L : orthopoly1d
591 Generalized Laguerre polynomial.
593 Notes
594 -----
595 For fixed :math:`\alpha`, the polynomials :math:`L_n^{(\alpha)}`
596 are orthogonal over :math:`[0, \infty)` with weight function
597 :math:`e^{-x}x^\alpha`.
599 The Laguerre polynomials are the special case where :math:`\alpha
600 = 0`.
602 See Also
603 --------
604 laguerre : Laguerre polynomial.
605 hyp1f1 : confluent hypergeometric function
607 References
608 ----------
609 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
610 Handbook of Mathematical Functions with Formulas,
611 Graphs, and Mathematical Tables. New York: Dover, 1972.
613 Examples
614 --------
615 The generalized Laguerre polynomials are closely related to the confluent
616 hypergeometric function :math:`{}_1F_1`:
618 .. math::
619 L_n^{(\alpha)} = \binom{n + \alpha}{n} {}_1F_1(-n, \alpha +1, x)
621 This can be verified, for example, for :math:`n = \alpha = 3` over the
622 interval :math:`[-1, 1]`:
624 >>> import numpy as np
625 >>> from scipy.special import binom
626 >>> from scipy.special import genlaguerre
627 >>> from scipy.special import hyp1f1
628 >>> x = np.arange(-1.0, 1.0, 0.01)
629 >>> np.allclose(genlaguerre(3, 3)(x), binom(6, 3) * hyp1f1(-3, 4, x))
630 True
632 This is the plot of the generalized Laguerre polynomials
633 :math:`L_3^{(\alpha)}` for some values of :math:`\alpha`:
635 >>> import matplotlib.pyplot as plt
636 >>> x = np.arange(-4.0, 12.0, 0.01)
637 >>> fig, ax = plt.subplots()
638 >>> ax.set_ylim(-5.0, 10.0)
639 >>> ax.set_title(r'Generalized Laguerre polynomials $L_3^{\alpha}$')
640 >>> for alpha in np.arange(0, 5):
641 ... ax.plot(x, genlaguerre(3, alpha)(x), label=rf'$L_3^{(alpha)}$')
642 >>> plt.legend(loc='best')
643 >>> plt.show()
645 """
646 if alpha <= -1:
647 raise ValueError("alpha must be > -1")
648 if n < 0:
649 raise ValueError("n must be nonnegative.")
651 if n == 0:
652 n1 = n + 1
653 else:
654 n1 = n
655 x, w = roots_genlaguerre(n1, alpha)
656 def wfunc(x):
657 return exp(-x) * x ** alpha
658 if n == 0:
659 x, w = [], []
660 hn = _gam(n + alpha + 1) / _gam(n + 1)
661 kn = (-1)**n / _gam(n + 1)
662 p = orthopoly1d(x, w, hn, kn, wfunc, (0, inf), monic,
663 lambda x: _ufuncs.eval_genlaguerre(n, alpha, x))
664 return p
666# Laguerre L_n(x)
669def roots_laguerre(n, mu=False):
670 r"""Gauss-Laguerre quadrature.
672 Compute the sample points and weights for Gauss-Laguerre
673 quadrature. The sample points are the roots of the nth degree
674 Laguerre polynomial, :math:`L_n(x)`. These sample points and
675 weights correctly integrate polynomials of degree :math:`2n - 1`
676 or less over the interval :math:`[0, \infty]` with weight function
677 :math:`w(x) = e^{-x}`. See 22.2.13 in [AS]_ for details.
679 Parameters
680 ----------
681 n : int
682 quadrature order
683 mu : bool, optional
684 If True, return the sum of the weights, optional.
686 Returns
687 -------
688 x : ndarray
689 Sample points
690 w : ndarray
691 Weights
692 mu : float
693 Sum of the weights
695 See Also
696 --------
697 scipy.integrate.quadrature
698 scipy.integrate.fixed_quad
699 numpy.polynomial.laguerre.laggauss
701 References
702 ----------
703 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
704 Handbook of Mathematical Functions with Formulas,
705 Graphs, and Mathematical Tables. New York: Dover, 1972.
707 """
708 return roots_genlaguerre(n, 0.0, mu=mu)
711def laguerre(n, monic=False):
712 r"""Laguerre polynomial.
714 Defined to be the solution of
716 .. math::
717 x\frac{d^2}{dx^2}L_n + (1 - x)\frac{d}{dx}L_n + nL_n = 0;
719 :math:`L_n` is a polynomial of degree :math:`n`.
721 Parameters
722 ----------
723 n : int
724 Degree of the polynomial.
725 monic : bool, optional
726 If `True`, scale the leading coefficient to be 1. Default is
727 `False`.
729 Returns
730 -------
731 L : orthopoly1d
732 Laguerre Polynomial.
734 Notes
735 -----
736 The polynomials :math:`L_n` are orthogonal over :math:`[0,
737 \infty)` with weight function :math:`e^{-x}`.
739 See Also
740 --------
741 genlaguerre : Generalized (associated) Laguerre polynomial.
743 References
744 ----------
745 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
746 Handbook of Mathematical Functions with Formulas,
747 Graphs, and Mathematical Tables. New York: Dover, 1972.
749 Examples
750 --------
751 The Laguerre polynomials :math:`L_n` are the special case
752 :math:`\alpha = 0` of the generalized Laguerre polynomials
753 :math:`L_n^{(\alpha)}`.
754 Let's verify it on the interval :math:`[-1, 1]`:
756 >>> import numpy as np
757 >>> from scipy.special import genlaguerre
758 >>> from scipy.special import laguerre
759 >>> x = np.arange(-1.0, 1.0, 0.01)
760 >>> np.allclose(genlaguerre(3, 0)(x), laguerre(3)(x))
761 True
763 The polynomials :math:`L_n` also satisfy the recurrence relation:
765 .. math::
766 (n + 1)L_{n+1}(x) = (2n +1 -x)L_n(x) - nL_{n-1}(x)
768 This can be easily checked on :math:`[0, 1]` for :math:`n = 3`:
770 >>> x = np.arange(0.0, 1.0, 0.01)
771 >>> np.allclose(4 * laguerre(4)(x),
772 ... (7 - x) * laguerre(3)(x) - 3 * laguerre(2)(x))
773 True
775 This is the plot of the first few Laguerre polynomials :math:`L_n`:
777 >>> import matplotlib.pyplot as plt
778 >>> x = np.arange(-1.0, 5.0, 0.01)
779 >>> fig, ax = plt.subplots()
780 >>> ax.set_ylim(-5.0, 5.0)
781 >>> ax.set_title(r'Laguerre polynomials $L_n$')
782 >>> for n in np.arange(0, 5):
783 ... ax.plot(x, laguerre(n)(x), label=rf'$L_{n}$')
784 >>> plt.legend(loc='best')
785 >>> plt.show()
787 """
788 if n < 0:
789 raise ValueError("n must be nonnegative.")
791 if n == 0:
792 n1 = n + 1
793 else:
794 n1 = n
795 x, w = roots_laguerre(n1)
796 if n == 0:
797 x, w = [], []
798 hn = 1.0
799 kn = (-1)**n / _gam(n + 1)
800 p = orthopoly1d(x, w, hn, kn, lambda x: exp(-x), (0, inf), monic,
801 lambda x: _ufuncs.eval_laguerre(n, x))
802 return p
804# Hermite 1 H_n(x)
807def roots_hermite(n, mu=False):
808 r"""Gauss-Hermite (physicist's) quadrature.
810 Compute the sample points and weights for Gauss-Hermite
811 quadrature. The sample points are the roots of the nth degree
812 Hermite polynomial, :math:`H_n(x)`. These sample points and
813 weights correctly integrate polynomials of degree :math:`2n - 1`
814 or less over the interval :math:`[-\infty, \infty]` with weight
815 function :math:`w(x) = e^{-x^2}`. See 22.2.14 in [AS]_ for
816 details.
818 Parameters
819 ----------
820 n : int
821 quadrature order
822 mu : bool, optional
823 If True, return the sum of the weights, optional.
825 Returns
826 -------
827 x : ndarray
828 Sample points
829 w : ndarray
830 Weights
831 mu : float
832 Sum of the weights
834 Notes
835 -----
836 For small n up to 150 a modified version of the Golub-Welsch
837 algorithm is used. Nodes are computed from the eigenvalue
838 problem and improved by one step of a Newton iteration.
839 The weights are computed from the well-known analytical formula.
841 For n larger than 150 an optimal asymptotic algorithm is applied
842 which computes nodes and weights in a numerically stable manner.
843 The algorithm has linear runtime making computation for very
844 large n (several thousand or more) feasible.
846 See Also
847 --------
848 scipy.integrate.quadrature
849 scipy.integrate.fixed_quad
850 numpy.polynomial.hermite.hermgauss
851 roots_hermitenorm
853 References
854 ----------
855 .. [townsend.trogdon.olver-2014]
856 Townsend, A. and Trogdon, T. and Olver, S. (2014)
857 *Fast computation of Gauss quadrature nodes and
858 weights on the whole real line*. :arXiv:`1410.5286`.
859 .. [townsend.trogdon.olver-2015]
860 Townsend, A. and Trogdon, T. and Olver, S. (2015)
861 *Fast computation of Gauss quadrature nodes and
862 weights on the whole real line*.
863 IMA Journal of Numerical Analysis
864 :doi:`10.1093/imanum/drv002`.
865 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
866 Handbook of Mathematical Functions with Formulas,
867 Graphs, and Mathematical Tables. New York: Dover, 1972.
869 """
870 m = int(n)
871 if n < 1 or n != m:
872 raise ValueError("n must be a positive integer.")
874 mu0 = np.sqrt(np.pi)
875 if n <= 150:
876 def an_func(k):
877 return 0.0 * k
878 def bn_func(k):
879 return np.sqrt(k / 2.0)
880 f = _ufuncs.eval_hermite
881 def df(n, x):
882 return 2.0 * n * _ufuncs.eval_hermite(n - 1, x)
883 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu)
884 else:
885 nodes, weights = _roots_hermite_asy(m)
886 if mu:
887 return nodes, weights, mu0
888 else:
889 return nodes, weights
892def _compute_tauk(n, k, maxit=5):
893 """Helper function for Tricomi initial guesses
895 For details, see formula 3.1 in lemma 3.1 in the
896 original paper.
898 Parameters
899 ----------
900 n : int
901 Quadrature order
902 k : ndarray of type int
903 Index of roots :math:`\tau_k` to compute
904 maxit : int
905 Number of Newton maxit performed, the default
906 value of 5 is sufficient.
908 Returns
909 -------
910 tauk : ndarray
911 Roots of equation 3.1
913 See Also
914 --------
915 initial_nodes_a
916 roots_hermite_asy
917 """
918 a = n % 2 - 0.5
919 c = (4.0*floor(n/2.0) - 4.0*k + 3.0)*pi / (4.0*floor(n/2.0) + 2.0*a + 2.0)
920 def f(x):
921 return x - sin(x) - c
922 def df(x):
923 return 1.0 - cos(x)
924 xi = 0.5*pi
925 for i in range(maxit):
926 xi = xi - f(xi)/df(xi)
927 return xi
930def _initial_nodes_a(n, k):
931 r"""Tricomi initial guesses
933 Computes an initial approximation to the square of the `k`-th
934 (positive) root :math:`x_k` of the Hermite polynomial :math:`H_n`
935 of order :math:`n`. The formula is the one from lemma 3.1 in the
936 original paper. The guesses are accurate except in the region
937 near :math:`\sqrt{2n + 1}`.
939 Parameters
940 ----------
941 n : int
942 Quadrature order
943 k : ndarray of type int
944 Index of roots to compute
946 Returns
947 -------
948 xksq : ndarray
949 Square of the approximate roots
951 See Also
952 --------
953 initial_nodes
954 roots_hermite_asy
955 """
956 tauk = _compute_tauk(n, k)
957 sigk = cos(0.5*tauk)**2
958 a = n % 2 - 0.5
959 nu = 4.0*floor(n/2.0) + 2.0*a + 2.0
960 # Initial approximation of Hermite roots (square)
961 xksq = nu*sigk - 1.0/(3.0*nu) * (5.0/(4.0*(1.0-sigk)**2) - 1.0/(1.0-sigk) - 0.25)
962 return xksq
965def _initial_nodes_b(n, k):
966 r"""Gatteschi initial guesses
968 Computes an initial approximation to the square of the kth
969 (positive) root :math:`x_k` of the Hermite polynomial :math:`H_n`
970 of order :math:`n`. The formula is the one from lemma 3.2 in the
971 original paper. The guesses are accurate in the region just
972 below :math:`\sqrt{2n + 1}`.
974 Parameters
975 ----------
976 n : int
977 Quadrature order
978 k : ndarray of type int
979 Index of roots to compute
981 Returns
982 -------
983 xksq : ndarray
984 Square of the approximate root
986 See Also
987 --------
988 initial_nodes
989 roots_hermite_asy
990 """
991 a = n % 2 - 0.5
992 nu = 4.0*floor(n/2.0) + 2.0*a + 2.0
993 # Airy roots by approximation
994 ak = _specfun.airyzo(k.max(), 1)[0][::-1]
995 # Initial approximation of Hermite roots (square)
996 xksq = (nu +
997 2.0**(2.0/3.0) * ak * nu**(1.0/3.0) +
998 1.0/5.0 * 2.0**(4.0/3.0) * ak**2 * nu**(-1.0/3.0) +
999 (9.0/140.0 - 12.0/175.0 * ak**3) * nu**(-1.0) +
1000 (16.0/1575.0 * ak + 92.0/7875.0 * ak**4) * 2.0**(2.0/3.0) * nu**(-5.0/3.0) -
1001 (15152.0/3031875.0 * ak**5 + 1088.0/121275.0 * ak**2) * 2.0**(1.0/3.0) * nu**(-7.0/3.0))
1002 return xksq
1005def _initial_nodes(n):
1006 """Initial guesses for the Hermite roots
1008 Computes an initial approximation to the non-negative
1009 roots :math:`x_k` of the Hermite polynomial :math:`H_n`
1010 of order :math:`n`. The Tricomi and Gatteschi initial
1011 guesses are used in the region where they are accurate.
1013 Parameters
1014 ----------
1015 n : int
1016 Quadrature order
1018 Returns
1019 -------
1020 xk : ndarray
1021 Approximate roots
1023 See Also
1024 --------
1025 roots_hermite_asy
1026 """
1027 # Turnover point
1028 # linear polynomial fit to error of 10, 25, 40, ..., 1000 point rules
1029 fit = 0.49082003*n - 4.37859653
1030 turnover = around(fit).astype(int)
1031 # Compute all approximations
1032 ia = arange(1, int(floor(n*0.5)+1))
1033 ib = ia[::-1]
1034 xasq = _initial_nodes_a(n, ia[:turnover+1])
1035 xbsq = _initial_nodes_b(n, ib[turnover+1:])
1036 # Combine
1037 iv = sqrt(hstack([xasq, xbsq]))
1038 # Central node is always zero
1039 if n % 2 == 1:
1040 iv = hstack([0.0, iv])
1041 return iv
1044def _pbcf(n, theta):
1045 r"""Asymptotic series expansion of parabolic cylinder function
1047 The implementation is based on sections 3.2 and 3.3 from the
1048 original paper. Compared to the published version this code
1049 adds one more term to the asymptotic series. The detailed
1050 formulas can be found at [parabolic-asymptotics]_. The evaluation
1051 is done in a transformed variable :math:`\theta := \arccos(t)`
1052 where :math:`t := x / \mu` and :math:`\mu := \sqrt{2n + 1}`.
1054 Parameters
1055 ----------
1056 n : int
1057 Quadrature order
1058 theta : ndarray
1059 Transformed position variable
1061 Returns
1062 -------
1063 U : ndarray
1064 Value of the parabolic cylinder function :math:`U(a, \theta)`.
1065 Ud : ndarray
1066 Value of the derivative :math:`U^{\prime}(a, \theta)` of
1067 the parabolic cylinder function.
1069 See Also
1070 --------
1071 roots_hermite_asy
1073 References
1074 ----------
1075 .. [parabolic-asymptotics]
1076 https://dlmf.nist.gov/12.10#vii
1077 """
1078 st = sin(theta)
1079 ct = cos(theta)
1080 # https://dlmf.nist.gov/12.10#vii
1081 mu = 2.0*n + 1.0
1082 # https://dlmf.nist.gov/12.10#E23
1083 eta = 0.5*theta - 0.5*st*ct
1084 # https://dlmf.nist.gov/12.10#E39
1085 zeta = -(3.0*eta/2.0) ** (2.0/3.0)
1086 # https://dlmf.nist.gov/12.10#E40
1087 phi = (-zeta / st**2) ** (0.25)
1088 # Coefficients
1089 # https://dlmf.nist.gov/12.10#E43
1090 a0 = 1.0
1091 a1 = 0.10416666666666666667
1092 a2 = 0.08355034722222222222
1093 a3 = 0.12822657455632716049
1094 a4 = 0.29184902646414046425
1095 a5 = 0.88162726744375765242
1096 b0 = 1.0
1097 b1 = -0.14583333333333333333
1098 b2 = -0.09874131944444444444
1099 b3 = -0.14331205391589506173
1100 b4 = -0.31722720267841354810
1101 b5 = -0.94242914795712024914
1102 # Polynomials
1103 # https://dlmf.nist.gov/12.10#E9
1104 # https://dlmf.nist.gov/12.10#E10
1105 ctp = ct ** arange(16).reshape((-1,1))
1106 u0 = 1.0
1107 u1 = (1.0*ctp[3,:] - 6.0*ct) / 24.0
1108 u2 = (-9.0*ctp[4,:] + 249.0*ctp[2,:] + 145.0) / 1152.0
1109 u3 = (-4042.0*ctp[9,:] + 18189.0*ctp[7,:] - 28287.0*ctp[5,:] - 151995.0*ctp[3,:] - 259290.0*ct) / 414720.0
1110 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
1111 u5 = (82393456.0*ctp[15,:] - 617950920.0*ctp[13,:] + 1994971575.0*ctp[11,:] - 3630137104.0*ctp[9,:] + 4433574213.0*ctp[7,:]
1112 - 37370295816.0*ctp[5,:] - 119582875013.0*ctp[3,:] - 34009066266.0*ct) / 6688604160.0
1113 v0 = 1.0
1114 v1 = (1.0*ctp[3,:] + 6.0*ct) / 24.0
1115 v2 = (15.0*ctp[4,:] - 327.0*ctp[2,:] - 143.0) / 1152.0
1116 v3 = (-4042.0*ctp[9,:] + 18189.0*ctp[7,:] - 36387.0*ctp[5,:] + 238425.0*ctp[3,:] + 259290.0*ct) / 414720.0
1117 v4 = (-121260.0*ctp[10,:] + 551733.0*ctp[8,:] - 151958.0*ctp[6,:] - 57484425.0*ctp[4,:] - 132752238.0*ctp[2,:] - 12118727) / 39813120.0
1118 v5 = (82393456.0*ctp[15,:] - 617950920.0*ctp[13,:] + 2025529095.0*ctp[11,:] - 3750839308.0*ctp[9,:] + 3832454253.0*ctp[7,:]
1119 + 35213253348.0*ctp[5,:] + 130919230435.0*ctp[3,:] + 34009066266*ct) / 6688604160.0
1120 # Airy Evaluation (Bi and Bip unused)
1121 Ai, Aip, Bi, Bip = airy(mu**(4.0/6.0) * zeta)
1122 # Prefactor for U
1123 P = 2.0*sqrt(pi) * mu**(1.0/6.0) * phi
1124 # Terms for U
1125 # https://dlmf.nist.gov/12.10#E42
1126 phip = phi ** arange(6, 31, 6).reshape((-1,1))
1127 A0 = b0*u0
1128 A1 = (b2*u0 + phip[0,:]*b1*u1 + phip[1,:]*b0*u2) / zeta**3
1129 A2 = (b4*u0 + phip[0,:]*b3*u1 + phip[1,:]*b2*u2 + phip[2,:]*b1*u3 + phip[3,:]*b0*u4) / zeta**6
1130 B0 = -(a1*u0 + phip[0,:]*a0*u1) / zeta**2
1131 B1 = -(a3*u0 + phip[0,:]*a2*u1 + phip[1,:]*a1*u2 + phip[2,:]*a0*u3) / zeta**5
1132 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
1133 # U
1134 # https://dlmf.nist.gov/12.10#E35
1135 U = P * (Ai * (A0 + A1/mu**2.0 + A2/mu**4.0) +
1136 Aip * (B0 + B1/mu**2.0 + B2/mu**4.0) / mu**(8.0/6.0))
1137 # Prefactor for derivative of U
1138 Pd = sqrt(2.0*pi) * mu**(2.0/6.0) / phi
1139 # Terms for derivative of U
1140 # https://dlmf.nist.gov/12.10#E46
1141 C0 = -(b1*v0 + phip[0,:]*b0*v1) / zeta
1142 C1 = -(b3*v0 + phip[0,:]*b2*v1 + phip[1,:]*b1*v2 + phip[2,:]*b0*v3) / zeta**4
1143 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
1144 D0 = a0*v0
1145 D1 = (a2*v0 + phip[0,:]*a1*v1 + phip[1,:]*a0*v2) / zeta**3
1146 D2 = (a4*v0 + phip[0,:]*a3*v1 + phip[1,:]*a2*v2 + phip[2,:]*a1*v3 + phip[3,:]*a0*v4) / zeta**6
1147 # Derivative of U
1148 # https://dlmf.nist.gov/12.10#E36
1149 Ud = Pd * (Ai * (C0 + C1/mu**2.0 + C2/mu**4.0) / mu**(4.0/6.0) +
1150 Aip * (D0 + D1/mu**2.0 + D2/mu**4.0))
1151 return U, Ud
1154def _newton(n, x_initial, maxit=5):
1155 """Newton iteration for polishing the asymptotic approximation
1156 to the zeros of the Hermite polynomials.
1158 Parameters
1159 ----------
1160 n : int
1161 Quadrature order
1162 x_initial : ndarray
1163 Initial guesses for the roots
1164 maxit : int
1165 Maximal number of Newton iterations.
1166 The default 5 is sufficient, usually
1167 only one or two steps are needed.
1169 Returns
1170 -------
1171 nodes : ndarray
1172 Quadrature nodes
1173 weights : ndarray
1174 Quadrature weights
1176 See Also
1177 --------
1178 roots_hermite_asy
1179 """
1180 # Variable transformation
1181 mu = sqrt(2.0*n + 1.0)
1182 t = x_initial / mu
1183 theta = arccos(t)
1184 # Newton iteration
1185 for i in range(maxit):
1186 u, ud = _pbcf(n, theta)
1187 dtheta = u / (sqrt(2.0) * mu * sin(theta) * ud)
1188 theta = theta + dtheta
1189 if max(abs(dtheta)) < 1e-14:
1190 break
1191 # Undo variable transformation
1192 x = mu * cos(theta)
1193 # Central node is always zero
1194 if n % 2 == 1:
1195 x[0] = 0.0
1196 # Compute weights
1197 w = exp(-x**2) / (2.0*ud**2)
1198 return x, w
1201def _roots_hermite_asy(n):
1202 r"""Gauss-Hermite (physicist's) quadrature for large n.
1204 Computes the sample points and weights for Gauss-Hermite quadrature.
1205 The sample points are the roots of the nth degree Hermite polynomial,
1206 :math:`H_n(x)`. These sample points and weights correctly integrate
1207 polynomials of degree :math:`2n - 1` or less over the interval
1208 :math:`[-\infty, \infty]` with weight function :math:`f(x) = e^{-x^2}`.
1210 This method relies on asymptotic expansions which work best for n > 150.
1211 The algorithm has linear runtime making computation for very large n
1212 feasible.
1214 Parameters
1215 ----------
1216 n : int
1217 quadrature order
1219 Returns
1220 -------
1221 nodes : ndarray
1222 Quadrature nodes
1223 weights : ndarray
1224 Quadrature weights
1226 See Also
1227 --------
1228 roots_hermite
1230 References
1231 ----------
1232 .. [townsend.trogdon.olver-2014]
1233 Townsend, A. and Trogdon, T. and Olver, S. (2014)
1234 *Fast computation of Gauss quadrature nodes and
1235 weights on the whole real line*. :arXiv:`1410.5286`.
1237 .. [townsend.trogdon.olver-2015]
1238 Townsend, A. and Trogdon, T. and Olver, S. (2015)
1239 *Fast computation of Gauss quadrature nodes and
1240 weights on the whole real line*.
1241 IMA Journal of Numerical Analysis
1242 :doi:`10.1093/imanum/drv002`.
1243 """
1244 iv = _initial_nodes(n)
1245 nodes, weights = _newton(n, iv)
1246 # Combine with negative parts
1247 if n % 2 == 0:
1248 nodes = hstack([-nodes[::-1], nodes])
1249 weights = hstack([weights[::-1], weights])
1250 else:
1251 nodes = hstack([-nodes[-1:0:-1], nodes])
1252 weights = hstack([weights[-1:0:-1], weights])
1253 # Scale weights
1254 weights *= sqrt(pi) / sum(weights)
1255 return nodes, weights
1258def hermite(n, monic=False):
1259 r"""Physicist's Hermite polynomial.
1261 Defined by
1263 .. math::
1265 H_n(x) = (-1)^ne^{x^2}\frac{d^n}{dx^n}e^{-x^2};
1267 :math:`H_n` is a polynomial of degree :math:`n`.
1269 Parameters
1270 ----------
1271 n : int
1272 Degree of the polynomial.
1273 monic : bool, optional
1274 If `True`, scale the leading coefficient to be 1. Default is
1275 `False`.
1277 Returns
1278 -------
1279 H : orthopoly1d
1280 Hermite polynomial.
1282 Notes
1283 -----
1284 The polynomials :math:`H_n` are orthogonal over :math:`(-\infty,
1285 \infty)` with weight function :math:`e^{-x^2}`.
1287 Examples
1288 --------
1289 >>> from scipy import special
1290 >>> import matplotlib.pyplot as plt
1291 >>> import numpy as np
1293 >>> p_monic = special.hermite(3, monic=True)
1294 >>> p_monic
1295 poly1d([ 1. , 0. , -1.5, 0. ])
1296 >>> p_monic(1)
1297 -0.49999999999999983
1298 >>> x = np.linspace(-3, 3, 400)
1299 >>> y = p_monic(x)
1300 >>> plt.plot(x, y)
1301 >>> plt.title("Monic Hermite polynomial of degree 3")
1302 >>> plt.xlabel("x")
1303 >>> plt.ylabel("H_3(x)")
1304 >>> plt.show()
1306 """
1307 if n < 0:
1308 raise ValueError("n must be nonnegative.")
1310 if n == 0:
1311 n1 = n + 1
1312 else:
1313 n1 = n
1314 x, w = roots_hermite(n1)
1315 def wfunc(x):
1316 return exp(-x * x)
1317 if n == 0:
1318 x, w = [], []
1319 hn = 2**n * _gam(n + 1) * sqrt(pi)
1320 kn = 2**n
1321 p = orthopoly1d(x, w, hn, kn, wfunc, (-inf, inf), monic,
1322 lambda x: _ufuncs.eval_hermite(n, x))
1323 return p
1325# Hermite 2 He_n(x)
1328def roots_hermitenorm(n, mu=False):
1329 r"""Gauss-Hermite (statistician's) quadrature.
1331 Compute the sample points and weights for Gauss-Hermite
1332 quadrature. The sample points are the roots of the nth degree
1333 Hermite polynomial, :math:`He_n(x)`. These sample points and
1334 weights correctly integrate polynomials of degree :math:`2n - 1`
1335 or less over the interval :math:`[-\infty, \infty]` with weight
1336 function :math:`w(x) = e^{-x^2/2}`. See 22.2.15 in [AS]_ for more
1337 details.
1339 Parameters
1340 ----------
1341 n : int
1342 quadrature order
1343 mu : bool, optional
1344 If True, return the sum of the weights, optional.
1346 Returns
1347 -------
1348 x : ndarray
1349 Sample points
1350 w : ndarray
1351 Weights
1352 mu : float
1353 Sum of the weights
1355 Notes
1356 -----
1357 For small n up to 150 a modified version of the Golub-Welsch
1358 algorithm is used. Nodes are computed from the eigenvalue
1359 problem and improved by one step of a Newton iteration.
1360 The weights are computed from the well-known analytical formula.
1362 For n larger than 150 an optimal asymptotic algorithm is used
1363 which computes nodes and weights in a numerical stable manner.
1364 The algorithm has linear runtime making computation for very
1365 large n (several thousand or more) feasible.
1367 See Also
1368 --------
1369 scipy.integrate.quadrature
1370 scipy.integrate.fixed_quad
1371 numpy.polynomial.hermite_e.hermegauss
1373 References
1374 ----------
1375 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1376 Handbook of Mathematical Functions with Formulas,
1377 Graphs, and Mathematical Tables. New York: Dover, 1972.
1379 """
1380 m = int(n)
1381 if n < 1 or n != m:
1382 raise ValueError("n must be a positive integer.")
1384 mu0 = np.sqrt(2.0*np.pi)
1385 if n <= 150:
1386 def an_func(k):
1387 return 0.0 * k
1388 def bn_func(k):
1389 return np.sqrt(k)
1390 f = _ufuncs.eval_hermitenorm
1391 def df(n, x):
1392 return n * _ufuncs.eval_hermitenorm(n - 1, x)
1393 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu)
1394 else:
1395 nodes, weights = _roots_hermite_asy(m)
1396 # Transform
1397 nodes *= sqrt(2)
1398 weights *= sqrt(2)
1399 if mu:
1400 return nodes, weights, mu0
1401 else:
1402 return nodes, weights
1405def hermitenorm(n, monic=False):
1406 r"""Normalized (probabilist's) Hermite polynomial.
1408 Defined by
1410 .. math::
1412 He_n(x) = (-1)^ne^{x^2/2}\frac{d^n}{dx^n}e^{-x^2/2};
1414 :math:`He_n` is a polynomial of degree :math:`n`.
1416 Parameters
1417 ----------
1418 n : int
1419 Degree of the polynomial.
1420 monic : bool, optional
1421 If `True`, scale the leading coefficient to be 1. Default is
1422 `False`.
1424 Returns
1425 -------
1426 He : orthopoly1d
1427 Hermite polynomial.
1429 Notes
1430 -----
1432 The polynomials :math:`He_n` are orthogonal over :math:`(-\infty,
1433 \infty)` with weight function :math:`e^{-x^2/2}`.
1435 """
1436 if n < 0:
1437 raise ValueError("n must be nonnegative.")
1439 if n == 0:
1440 n1 = n + 1
1441 else:
1442 n1 = n
1443 x, w = roots_hermitenorm(n1)
1444 def wfunc(x):
1445 return exp(-x * x / 2.0)
1446 if n == 0:
1447 x, w = [], []
1448 hn = sqrt(2 * pi) * _gam(n + 1)
1449 kn = 1.0
1450 p = orthopoly1d(x, w, hn, kn, wfunc=wfunc, limits=(-inf, inf), monic=monic,
1451 eval_func=lambda x: _ufuncs.eval_hermitenorm(n, x))
1452 return p
1454# The remainder of the polynomials can be derived from the ones above.
1456# Ultraspherical (Gegenbauer) C^(alpha)_n(x)
1459def roots_gegenbauer(n, alpha, mu=False):
1460 r"""Gauss-Gegenbauer quadrature.
1462 Compute the sample points and weights for Gauss-Gegenbauer
1463 quadrature. The sample points are the roots of the nth degree
1464 Gegenbauer polynomial, :math:`C^{\alpha}_n(x)`. These sample
1465 points and weights correctly integrate polynomials of degree
1466 :math:`2n - 1` or less over the interval :math:`[-1, 1]` with
1467 weight function :math:`w(x) = (1 - x^2)^{\alpha - 1/2}`. See
1468 22.2.3 in [AS]_ for more details.
1470 Parameters
1471 ----------
1472 n : int
1473 quadrature order
1474 alpha : float
1475 alpha must be > -0.5
1476 mu : bool, optional
1477 If True, return the sum of the weights, optional.
1479 Returns
1480 -------
1481 x : ndarray
1482 Sample points
1483 w : ndarray
1484 Weights
1485 mu : float
1486 Sum of the weights
1488 See Also
1489 --------
1490 scipy.integrate.quadrature
1491 scipy.integrate.fixed_quad
1493 References
1494 ----------
1495 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1496 Handbook of Mathematical Functions with Formulas,
1497 Graphs, and Mathematical Tables. New York: Dover, 1972.
1499 """
1500 m = int(n)
1501 if n < 1 or n != m:
1502 raise ValueError("n must be a positive integer.")
1503 if alpha < -0.5:
1504 raise ValueError("alpha must be greater than -0.5.")
1505 elif alpha == 0.0:
1506 # C(n,0,x) == 0 uniformly, however, as alpha->0, C(n,alpha,x)->T(n,x)
1507 # strictly, we should just error out here, since the roots are not
1508 # really defined, but we used to return something useful, so let's
1509 # keep doing so.
1510 return roots_chebyt(n, mu)
1512 if alpha <= 170:
1513 mu0 = (np.sqrt(np.pi) * _ufuncs.gamma(alpha + 0.5)) \
1514 / _ufuncs.gamma(alpha + 1)
1515 else:
1516 # For large alpha we use a Taylor series expansion around inf,
1517 # expressed as a 6th order polynomial of a^-1 and using Horner's
1518 # method to minimize computation and maximize precision
1519 inv_alpha = 1. / alpha
1520 coeffs = np.array([0.000207186, -0.00152206, -0.000640869,
1521 0.00488281, 0.0078125, -0.125, 1.])
1522 mu0 = coeffs[0]
1523 for term in range(1, len(coeffs)):
1524 mu0 = mu0 * inv_alpha + coeffs[term]
1525 mu0 = mu0 * np.sqrt(np.pi / alpha)
1526 def an_func(k):
1527 return 0.0 * k
1528 def bn_func(k):
1529 return np.sqrt(k * (k + 2 * alpha - 1) / (4 * (k + alpha) * (k + alpha - 1)))
1530 def f(n, x):
1531 return _ufuncs.eval_gegenbauer(n, alpha, x)
1532 def df(n, x):
1533 return (-n * x * _ufuncs.eval_gegenbauer(n, alpha, x) + (n + 2 * alpha - 1) * _ufuncs.eval_gegenbauer(n - 1, alpha, x)) / (1 - x ** 2)
1534 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu)
1537def gegenbauer(n, alpha, monic=False):
1538 r"""Gegenbauer (ultraspherical) polynomial.
1540 Defined to be the solution of
1542 .. math::
1543 (1 - x^2)\frac{d^2}{dx^2}C_n^{(\alpha)}
1544 - (2\alpha + 1)x\frac{d}{dx}C_n^{(\alpha)}
1545 + n(n + 2\alpha)C_n^{(\alpha)} = 0
1547 for :math:`\alpha > -1/2`; :math:`C_n^{(\alpha)}` is a polynomial
1548 of degree :math:`n`.
1550 Parameters
1551 ----------
1552 n : int
1553 Degree of the polynomial.
1554 alpha : float
1555 Parameter, must be greater than -0.5.
1556 monic : bool, optional
1557 If `True`, scale the leading coefficient to be 1. Default is
1558 `False`.
1560 Returns
1561 -------
1562 C : orthopoly1d
1563 Gegenbauer polynomial.
1565 Notes
1566 -----
1567 The polynomials :math:`C_n^{(\alpha)}` are orthogonal over
1568 :math:`[-1,1]` with weight function :math:`(1 - x^2)^{(\alpha -
1569 1/2)}`.
1571 Examples
1572 --------
1573 >>> import numpy as np
1574 >>> from scipy import special
1575 >>> import matplotlib.pyplot as plt
1577 We can initialize a variable ``p`` as a Gegenbauer polynomial using the
1578 `gegenbauer` function and evaluate at a point ``x = 1``.
1580 >>> p = special.gegenbauer(3, 0.5, monic=False)
1581 >>> p
1582 poly1d([ 2.5, 0. , -1.5, 0. ])
1583 >>> p(1)
1584 1.0
1586 To evaluate ``p`` at various points ``x`` in the interval ``(-3, 3)``,
1587 simply pass an array ``x`` to ``p`` as follows:
1589 >>> x = np.linspace(-3, 3, 400)
1590 >>> y = p(x)
1592 We can then visualize ``x, y`` using `matplotlib.pyplot`.
1594 >>> fig, ax = plt.subplots()
1595 >>> ax.plot(x, y)
1596 >>> ax.set_title("Gegenbauer (ultraspherical) polynomial of degree 3")
1597 >>> ax.set_xlabel("x")
1598 >>> ax.set_ylabel("G_3(x)")
1599 >>> plt.show()
1601 """
1602 base = jacobi(n, alpha - 0.5, alpha - 0.5, monic=monic)
1603 if monic:
1604 return base
1605 # Abrahmowitz and Stegan 22.5.20
1606 factor = (_gam(2*alpha + n) * _gam(alpha + 0.5) /
1607 _gam(2*alpha) / _gam(alpha + 0.5 + n))
1608 base._scale(factor)
1609 base.__dict__['_eval_func'] = lambda x: _ufuncs.eval_gegenbauer(float(n),
1610 alpha, x)
1611 return base
1613# Chebyshev of the first kind: T_n(x) =
1614# n! sqrt(pi) / _gam(n+1./2)* P^(-1/2,-1/2)_n(x)
1615# Computed anew.
1618def roots_chebyt(n, mu=False):
1619 r"""Gauss-Chebyshev (first kind) quadrature.
1621 Computes the sample points and weights for Gauss-Chebyshev
1622 quadrature. The sample points are the roots of the nth degree
1623 Chebyshev polynomial of the first kind, :math:`T_n(x)`. These
1624 sample points and weights correctly integrate polynomials of
1625 degree :math:`2n - 1` or less over the interval :math:`[-1, 1]`
1626 with weight function :math:`w(x) = 1/\sqrt{1 - x^2}`. See 22.2.4
1627 in [AS]_ for more details.
1629 Parameters
1630 ----------
1631 n : int
1632 quadrature order
1633 mu : bool, optional
1634 If True, return the sum of the weights, optional.
1636 Returns
1637 -------
1638 x : ndarray
1639 Sample points
1640 w : ndarray
1641 Weights
1642 mu : float
1643 Sum of the weights
1645 See Also
1646 --------
1647 scipy.integrate.quadrature
1648 scipy.integrate.fixed_quad
1649 numpy.polynomial.chebyshev.chebgauss
1651 References
1652 ----------
1653 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1654 Handbook of Mathematical Functions with Formulas,
1655 Graphs, and Mathematical Tables. New York: Dover, 1972.
1657 """
1658 m = int(n)
1659 if n < 1 or n != m:
1660 raise ValueError('n must be a positive integer.')
1661 x = _ufuncs._sinpi(np.arange(-m + 1, m, 2) / (2*m))
1662 w = np.full_like(x, pi/m)
1663 if mu:
1664 return x, w, pi
1665 else:
1666 return x, w
1669def chebyt(n, monic=False):
1670 r"""Chebyshev polynomial of the first kind.
1672 Defined to be the solution of
1674 .. math::
1675 (1 - x^2)\frac{d^2}{dx^2}T_n - x\frac{d}{dx}T_n + n^2T_n = 0;
1677 :math:`T_n` is a polynomial of degree :math:`n`.
1679 Parameters
1680 ----------
1681 n : int
1682 Degree of the polynomial.
1683 monic : bool, optional
1684 If `True`, scale the leading coefficient to be 1. Default is
1685 `False`.
1687 Returns
1688 -------
1689 T : orthopoly1d
1690 Chebyshev polynomial of the first kind.
1692 Notes
1693 -----
1694 The polynomials :math:`T_n` are orthogonal over :math:`[-1, 1]`
1695 with weight function :math:`(1 - x^2)^{-1/2}`.
1697 See Also
1698 --------
1699 chebyu : Chebyshev polynomial of the second kind.
1701 References
1702 ----------
1703 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1704 Handbook of Mathematical Functions with Formulas,
1705 Graphs, and Mathematical Tables. New York: Dover, 1972.
1707 Examples
1708 --------
1709 Chebyshev polynomials of the first kind of order :math:`n` can
1710 be obtained as the determinant of specific :math:`n \times n`
1711 matrices. As an example we can check how the points obtained from
1712 the determinant of the following :math:`3 \times 3` matrix
1713 lay exacty on :math:`T_3`:
1715 >>> import numpy as np
1716 >>> import matplotlib.pyplot as plt
1717 >>> from scipy.linalg import det
1718 >>> from scipy.special import chebyt
1719 >>> x = np.arange(-1.0, 1.0, 0.01)
1720 >>> fig, ax = plt.subplots()
1721 >>> ax.set_ylim(-2.0, 2.0)
1722 >>> ax.set_title(r'Chebyshev polynomial $T_3$')
1723 >>> ax.plot(x, chebyt(3)(x), label=rf'$T_3$')
1724 >>> for p in np.arange(-1.0, 1.0, 0.1):
1725 ... ax.plot(p,
1726 ... det(np.array([[p, 1, 0], [1, 2*p, 1], [0, 1, 2*p]])),
1727 ... 'rx')
1728 >>> plt.legend(loc='best')
1729 >>> plt.show()
1731 They are also related to the Jacobi Polynomials
1732 :math:`P_n^{(-0.5, -0.5)}` through the relation:
1734 .. math::
1735 P_n^{(-0.5, -0.5)}(x) = \frac{1}{4^n} \binom{2n}{n} T_n(x)
1737 Let's verify it for :math:`n = 3`:
1739 >>> from scipy.special import binom
1740 >>> from scipy.special import jacobi
1741 >>> x = np.arange(-1.0, 1.0, 0.01)
1742 >>> np.allclose(jacobi(3, -0.5, -0.5)(x),
1743 ... 1/64 * binom(6, 3) * chebyt(3)(x))
1744 True
1746 We can plot the Chebyshev polynomials :math:`T_n` for some values
1747 of :math:`n`:
1749 >>> x = np.arange(-1.5, 1.5, 0.01)
1750 >>> fig, ax = plt.subplots()
1751 >>> ax.set_ylim(-4.0, 4.0)
1752 >>> ax.set_title(r'Chebyshev polynomials $T_n$')
1753 >>> for n in np.arange(2,5):
1754 ... ax.plot(x, chebyt(n)(x), label=rf'$T_n={n}$')
1755 >>> plt.legend(loc='best')
1756 >>> plt.show()
1758 """
1759 if n < 0:
1760 raise ValueError("n must be nonnegative.")
1762 def wfunc(x):
1763 return 1.0 / sqrt(1 - x * x)
1764 if n == 0:
1765 return orthopoly1d([], [], pi, 1.0, wfunc, (-1, 1), monic,
1766 lambda x: _ufuncs.eval_chebyt(n, x))
1767 n1 = n
1768 x, w, mu = roots_chebyt(n1, mu=True)
1769 hn = pi / 2
1770 kn = 2**(n - 1)
1771 p = orthopoly1d(x, w, hn, kn, wfunc, (-1, 1), monic,
1772 lambda x: _ufuncs.eval_chebyt(n, x))
1773 return p
1775# Chebyshev of the second kind
1776# U_n(x) = (n+1)! sqrt(pi) / (2*_gam(n+3./2)) * P^(1/2,1/2)_n(x)
1779def roots_chebyu(n, mu=False):
1780 r"""Gauss-Chebyshev (second kind) quadrature.
1782 Computes the sample points and weights for Gauss-Chebyshev
1783 quadrature. The sample points are the roots of the nth degree
1784 Chebyshev polynomial of the second kind, :math:`U_n(x)`. These
1785 sample points and weights correctly integrate polynomials of
1786 degree :math:`2n - 1` or less over the interval :math:`[-1, 1]`
1787 with weight function :math:`w(x) = \sqrt{1 - x^2}`. See 22.2.5 in
1788 [AS]_ for details.
1790 Parameters
1791 ----------
1792 n : int
1793 quadrature order
1794 mu : bool, optional
1795 If True, return the sum of the weights, optional.
1797 Returns
1798 -------
1799 x : ndarray
1800 Sample points
1801 w : ndarray
1802 Weights
1803 mu : float
1804 Sum of the weights
1806 See Also
1807 --------
1808 scipy.integrate.quadrature
1809 scipy.integrate.fixed_quad
1811 References
1812 ----------
1813 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1814 Handbook of Mathematical Functions with Formulas,
1815 Graphs, and Mathematical Tables. New York: Dover, 1972.
1817 """
1818 m = int(n)
1819 if n < 1 or n != m:
1820 raise ValueError('n must be a positive integer.')
1821 t = np.arange(m, 0, -1) * pi / (m + 1)
1822 x = np.cos(t)
1823 w = pi * np.sin(t)**2 / (m + 1)
1824 if mu:
1825 return x, w, pi / 2
1826 else:
1827 return x, w
1830def chebyu(n, monic=False):
1831 r"""Chebyshev polynomial of the second kind.
1833 Defined to be the solution of
1835 .. math::
1836 (1 - x^2)\frac{d^2}{dx^2}U_n - 3x\frac{d}{dx}U_n
1837 + n(n + 2)U_n = 0;
1839 :math:`U_n` is a polynomial of degree :math:`n`.
1841 Parameters
1842 ----------
1843 n : int
1844 Degree of the polynomial.
1845 monic : bool, optional
1846 If `True`, scale the leading coefficient to be 1. Default is
1847 `False`.
1849 Returns
1850 -------
1851 U : orthopoly1d
1852 Chebyshev polynomial of the second kind.
1854 Notes
1855 -----
1856 The polynomials :math:`U_n` are orthogonal over :math:`[-1, 1]`
1857 with weight function :math:`(1 - x^2)^{1/2}`.
1859 See Also
1860 --------
1861 chebyt : Chebyshev polynomial of the first kind.
1863 References
1864 ----------
1865 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1866 Handbook of Mathematical Functions with Formulas,
1867 Graphs, and Mathematical Tables. New York: Dover, 1972.
1869 Examples
1870 --------
1871 Chebyshev polynomials of the second kind of order :math:`n` can
1872 be obtained as the determinant of specific :math:`n \times n`
1873 matrices. As an example we can check how the points obtained from
1874 the determinant of the following :math:`3 \times 3` matrix
1875 lay exacty on :math:`U_3`:
1877 >>> import numpy as np
1878 >>> import matplotlib.pyplot as plt
1879 >>> from scipy.linalg import det
1880 >>> from scipy.special import chebyu
1881 >>> x = np.arange(-1.0, 1.0, 0.01)
1882 >>> fig, ax = plt.subplots()
1883 >>> ax.set_ylim(-2.0, 2.0)
1884 >>> ax.set_title(r'Chebyshev polynomial $U_3$')
1885 >>> ax.plot(x, chebyu(3)(x), label=rf'$U_3$')
1886 >>> for p in np.arange(-1.0, 1.0, 0.1):
1887 ... ax.plot(p,
1888 ... det(np.array([[2*p, 1, 0], [1, 2*p, 1], [0, 1, 2*p]])),
1889 ... 'rx')
1890 >>> plt.legend(loc='best')
1891 >>> plt.show()
1893 They satisfy the recurrence relation:
1895 .. math::
1896 U_{2n-1}(x) = 2 T_n(x)U_{n-1}(x)
1898 where the :math:`T_n` are the Chebyshev polynomial of the first kind.
1899 Let's verify it for :math:`n = 2`:
1901 >>> from scipy.special import chebyt
1902 >>> x = np.arange(-1.0, 1.0, 0.01)
1903 >>> np.allclose(chebyu(3)(x), 2 * chebyt(2)(x) * chebyu(1)(x))
1904 True
1906 We can plot the Chebyshev polynomials :math:`U_n` for some values
1907 of :math:`n`:
1909 >>> x = np.arange(-1.0, 1.0, 0.01)
1910 >>> fig, ax = plt.subplots()
1911 >>> ax.set_ylim(-1.5, 1.5)
1912 >>> ax.set_title(r'Chebyshev polynomials $U_n$')
1913 >>> for n in np.arange(1,5):
1914 ... ax.plot(x, chebyu(n)(x), label=rf'$U_n={n}$')
1915 >>> plt.legend(loc='best')
1916 >>> plt.show()
1918 """
1919 base = jacobi(n, 0.5, 0.5, monic=monic)
1920 if monic:
1921 return base
1922 factor = sqrt(pi) / 2.0 * _gam(n + 2) / _gam(n + 1.5)
1923 base._scale(factor)
1924 return base
1926# Chebyshev of the first kind C_n(x)
1929def roots_chebyc(n, mu=False):
1930 r"""Gauss-Chebyshev (first kind) quadrature.
1932 Compute the sample points and weights for Gauss-Chebyshev
1933 quadrature. The sample points are the roots of the nth degree
1934 Chebyshev polynomial of the first kind, :math:`C_n(x)`. These
1935 sample points and weights correctly integrate polynomials of
1936 degree :math:`2n - 1` or less over the interval :math:`[-2, 2]`
1937 with weight function :math:`w(x) = 1 / \sqrt{1 - (x/2)^2}`. See
1938 22.2.6 in [AS]_ for more details.
1940 Parameters
1941 ----------
1942 n : int
1943 quadrature order
1944 mu : bool, optional
1945 If True, return the sum of the weights, optional.
1947 Returns
1948 -------
1949 x : ndarray
1950 Sample points
1951 w : ndarray
1952 Weights
1953 mu : float
1954 Sum of the weights
1956 See Also
1957 --------
1958 scipy.integrate.quadrature
1959 scipy.integrate.fixed_quad
1961 References
1962 ----------
1963 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
1964 Handbook of Mathematical Functions with Formulas,
1965 Graphs, and Mathematical Tables. New York: Dover, 1972.
1967 """
1968 x, w, m = roots_chebyt(n, True)
1969 x *= 2
1970 w *= 2
1971 m *= 2
1972 if mu:
1973 return x, w, m
1974 else:
1975 return x, w
1978def chebyc(n, monic=False):
1979 r"""Chebyshev polynomial of the first kind on :math:`[-2, 2]`.
1981 Defined as :math:`C_n(x) = 2T_n(x/2)`, where :math:`T_n` is the
1982 nth Chebychev polynomial of the first kind.
1984 Parameters
1985 ----------
1986 n : int
1987 Degree of the polynomial.
1988 monic : bool, optional
1989 If `True`, scale the leading coefficient to be 1. Default is
1990 `False`.
1992 Returns
1993 -------
1994 C : orthopoly1d
1995 Chebyshev polynomial of the first kind on :math:`[-2, 2]`.
1997 Notes
1998 -----
1999 The polynomials :math:`C_n(x)` are orthogonal over :math:`[-2, 2]`
2000 with weight function :math:`1/\sqrt{1 - (x/2)^2}`.
2002 See Also
2003 --------
2004 chebyt : Chebyshev polynomial of the first kind.
2006 References
2007 ----------
2008 .. [1] Abramowitz and Stegun, "Handbook of Mathematical Functions"
2009 Section 22. National Bureau of Standards, 1972.
2011 """
2012 if n < 0:
2013 raise ValueError("n must be nonnegative.")
2015 if n == 0:
2016 n1 = n + 1
2017 else:
2018 n1 = n
2019 x, w = roots_chebyc(n1)
2020 if n == 0:
2021 x, w = [], []
2022 hn = 4 * pi * ((n == 0) + 1)
2023 kn = 1.0
2024 p = orthopoly1d(x, w, hn, kn,
2025 wfunc=lambda x: 1.0 / sqrt(1 - x * x / 4.0),
2026 limits=(-2, 2), monic=monic)
2027 if not monic:
2028 p._scale(2.0 / p(2))
2029 p.__dict__['_eval_func'] = lambda x: _ufuncs.eval_chebyc(n, x)
2030 return p
2032# Chebyshev of the second kind S_n(x)
2035def roots_chebys(n, mu=False):
2036 r"""Gauss-Chebyshev (second kind) quadrature.
2038 Compute the sample points and weights for Gauss-Chebyshev
2039 quadrature. The sample points are the roots of the nth degree
2040 Chebyshev polynomial of the second kind, :math:`S_n(x)`. These
2041 sample points and weights correctly integrate polynomials of
2042 degree :math:`2n - 1` or less over the interval :math:`[-2, 2]`
2043 with weight function :math:`w(x) = \sqrt{1 - (x/2)^2}`. See 22.2.7
2044 in [AS]_ for more details.
2046 Parameters
2047 ----------
2048 n : int
2049 quadrature order
2050 mu : bool, optional
2051 If True, return the sum of the weights, optional.
2053 Returns
2054 -------
2055 x : ndarray
2056 Sample points
2057 w : ndarray
2058 Weights
2059 mu : float
2060 Sum of the weights
2062 See Also
2063 --------
2064 scipy.integrate.quadrature
2065 scipy.integrate.fixed_quad
2067 References
2068 ----------
2069 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2070 Handbook of Mathematical Functions with Formulas,
2071 Graphs, and Mathematical Tables. New York: Dover, 1972.
2073 """
2074 x, w, m = roots_chebyu(n, True)
2075 x *= 2
2076 w *= 2
2077 m *= 2
2078 if mu:
2079 return x, w, m
2080 else:
2081 return x, w
2084def chebys(n, monic=False):
2085 r"""Chebyshev polynomial of the second kind on :math:`[-2, 2]`.
2087 Defined as :math:`S_n(x) = U_n(x/2)` where :math:`U_n` is the
2088 nth Chebychev polynomial of the second kind.
2090 Parameters
2091 ----------
2092 n : int
2093 Degree of the polynomial.
2094 monic : bool, optional
2095 If `True`, scale the leading coefficient to be 1. Default is
2096 `False`.
2098 Returns
2099 -------
2100 S : orthopoly1d
2101 Chebyshev polynomial of the second kind on :math:`[-2, 2]`.
2103 Notes
2104 -----
2105 The polynomials :math:`S_n(x)` are orthogonal over :math:`[-2, 2]`
2106 with weight function :math:`\sqrt{1 - (x/2)}^2`.
2108 See Also
2109 --------
2110 chebyu : Chebyshev polynomial of the second kind
2112 References
2113 ----------
2114 .. [1] Abramowitz and Stegun, "Handbook of Mathematical Functions"
2115 Section 22. National Bureau of Standards, 1972.
2117 """
2118 if n < 0:
2119 raise ValueError("n must be nonnegative.")
2121 if n == 0:
2122 n1 = n + 1
2123 else:
2124 n1 = n
2125 x, w = roots_chebys(n1)
2126 if n == 0:
2127 x, w = [], []
2128 hn = pi
2129 kn = 1.0
2130 p = orthopoly1d(x, w, hn, kn,
2131 wfunc=lambda x: sqrt(1 - x * x / 4.0),
2132 limits=(-2, 2), monic=monic)
2133 if not monic:
2134 factor = (n + 1.0) / p(2)
2135 p._scale(factor)
2136 p.__dict__['_eval_func'] = lambda x: _ufuncs.eval_chebys(n, x)
2137 return p
2139# Shifted Chebyshev of the first kind T^*_n(x)
2142def roots_sh_chebyt(n, mu=False):
2143 r"""Gauss-Chebyshev (first kind, shifted) quadrature.
2145 Compute the sample points and weights for Gauss-Chebyshev
2146 quadrature. The sample points are the roots of the nth degree
2147 shifted Chebyshev polynomial of the first kind, :math:`T_n(x)`.
2148 These sample points and weights correctly integrate polynomials of
2149 degree :math:`2n - 1` or less over the interval :math:`[0, 1]`
2150 with weight function :math:`w(x) = 1/\sqrt{x - x^2}`. See 22.2.8
2151 in [AS]_ for more details.
2153 Parameters
2154 ----------
2155 n : int
2156 quadrature order
2157 mu : bool, optional
2158 If True, return the sum of the weights, optional.
2160 Returns
2161 -------
2162 x : ndarray
2163 Sample points
2164 w : ndarray
2165 Weights
2166 mu : float
2167 Sum of the weights
2169 See Also
2170 --------
2171 scipy.integrate.quadrature
2172 scipy.integrate.fixed_quad
2174 References
2175 ----------
2176 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2177 Handbook of Mathematical Functions with Formulas,
2178 Graphs, and Mathematical Tables. New York: Dover, 1972.
2180 """
2181 xw = roots_chebyt(n, mu)
2182 return ((xw[0] + 1) / 2,) + xw[1:]
2185def sh_chebyt(n, monic=False):
2186 r"""Shifted Chebyshev polynomial of the first kind.
2188 Defined as :math:`T^*_n(x) = T_n(2x - 1)` for :math:`T_n` the nth
2189 Chebyshev polynomial of the first kind.
2191 Parameters
2192 ----------
2193 n : int
2194 Degree of the polynomial.
2195 monic : bool, optional
2196 If `True`, scale the leading coefficient to be 1. Default is
2197 `False`.
2199 Returns
2200 -------
2201 T : orthopoly1d
2202 Shifted Chebyshev polynomial of the first kind.
2204 Notes
2205 -----
2206 The polynomials :math:`T^*_n` are orthogonal over :math:`[0, 1]`
2207 with weight function :math:`(x - x^2)^{-1/2}`.
2209 """
2210 base = sh_jacobi(n, 0.0, 0.5, monic=monic)
2211 if monic:
2212 return base
2213 if n > 0:
2214 factor = 4**n / 2.0
2215 else:
2216 factor = 1.0
2217 base._scale(factor)
2218 return base
2221# Shifted Chebyshev of the second kind U^*_n(x)
2222def roots_sh_chebyu(n, mu=False):
2223 r"""Gauss-Chebyshev (second kind, shifted) quadrature.
2225 Computes the sample points and weights for Gauss-Chebyshev
2226 quadrature. The sample points are the roots of the nth degree
2227 shifted Chebyshev polynomial of the second kind, :math:`U_n(x)`.
2228 These sample points and weights correctly integrate polynomials of
2229 degree :math:`2n - 1` or less over the interval :math:`[0, 1]`
2230 with weight function :math:`w(x) = \sqrt{x - x^2}`. See 22.2.9 in
2231 [AS]_ for more details.
2233 Parameters
2234 ----------
2235 n : int
2236 quadrature order
2237 mu : bool, optional
2238 If True, return the sum of the weights, optional.
2240 Returns
2241 -------
2242 x : ndarray
2243 Sample points
2244 w : ndarray
2245 Weights
2246 mu : float
2247 Sum of the weights
2249 See Also
2250 --------
2251 scipy.integrate.quadrature
2252 scipy.integrate.fixed_quad
2254 References
2255 ----------
2256 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2257 Handbook of Mathematical Functions with Formulas,
2258 Graphs, and Mathematical Tables. New York: Dover, 1972.
2260 """
2261 x, w, m = roots_chebyu(n, True)
2262 x = (x + 1) / 2
2263 m_us = _ufuncs.beta(1.5, 1.5)
2264 w *= m_us / m
2265 if mu:
2266 return x, w, m_us
2267 else:
2268 return x, w
2271def sh_chebyu(n, monic=False):
2272 r"""Shifted Chebyshev polynomial of the second kind.
2274 Defined as :math:`U^*_n(x) = U_n(2x - 1)` for :math:`U_n` the nth
2275 Chebyshev polynomial of the second kind.
2277 Parameters
2278 ----------
2279 n : int
2280 Degree of the polynomial.
2281 monic : bool, optional
2282 If `True`, scale the leading coefficient to be 1. Default is
2283 `False`.
2285 Returns
2286 -------
2287 U : orthopoly1d
2288 Shifted Chebyshev polynomial of the second kind.
2290 Notes
2291 -----
2292 The polynomials :math:`U^*_n` are orthogonal over :math:`[0, 1]`
2293 with weight function :math:`(x - x^2)^{1/2}`.
2295 """
2296 base = sh_jacobi(n, 2.0, 1.5, monic=monic)
2297 if monic:
2298 return base
2299 factor = 4**n
2300 base._scale(factor)
2301 return base
2303# Legendre
2306def roots_legendre(n, mu=False):
2307 r"""Gauss-Legendre quadrature.
2309 Compute the sample points and weights for Gauss-Legendre
2310 quadrature [GL]_. The sample points are the roots of the nth degree
2311 Legendre polynomial :math:`P_n(x)`. These sample points and
2312 weights correctly integrate polynomials of degree :math:`2n - 1`
2313 or less over the interval :math:`[-1, 1]` with weight function
2314 :math:`w(x) = 1`. See 2.2.10 in [AS]_ for more details.
2316 Parameters
2317 ----------
2318 n : int
2319 quadrature order
2320 mu : bool, optional
2321 If True, return the sum of the weights, optional.
2323 Returns
2324 -------
2325 x : ndarray
2326 Sample points
2327 w : ndarray
2328 Weights
2329 mu : float
2330 Sum of the weights
2332 See Also
2333 --------
2334 scipy.integrate.quadrature
2335 scipy.integrate.fixed_quad
2336 numpy.polynomial.legendre.leggauss
2338 References
2339 ----------
2340 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2341 Handbook of Mathematical Functions with Formulas,
2342 Graphs, and Mathematical Tables. New York: Dover, 1972.
2343 .. [GL] Gauss-Legendre quadrature, Wikipedia,
2344 https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature
2346 Examples
2347 --------
2348 >>> import numpy as np
2349 >>> from scipy.special import roots_legendre, eval_legendre
2350 >>> roots, weights = roots_legendre(9)
2352 ``roots`` holds the roots, and ``weights`` holds the weights for
2353 Gauss-Legendre quadrature.
2355 >>> roots
2356 array([-0.96816024, -0.83603111, -0.61337143, -0.32425342, 0. ,
2357 0.32425342, 0.61337143, 0.83603111, 0.96816024])
2358 >>> weights
2359 array([0.08127439, 0.18064816, 0.2606107 , 0.31234708, 0.33023936,
2360 0.31234708, 0.2606107 , 0.18064816, 0.08127439])
2362 Verify that we have the roots by evaluating the degree 9 Legendre
2363 polynomial at ``roots``. All the values are approximately zero:
2365 >>> eval_legendre(9, roots)
2366 array([-8.88178420e-16, -2.22044605e-16, 1.11022302e-16, 1.11022302e-16,
2367 0.00000000e+00, -5.55111512e-17, -1.94289029e-16, 1.38777878e-16,
2368 -8.32667268e-17])
2370 Here we'll show how the above values can be used to estimate the
2371 integral from 1 to 2 of f(t) = t + 1/t with Gauss-Legendre
2372 quadrature [GL]_. First define the function and the integration
2373 limits.
2375 >>> def f(t):
2376 ... return t + 1/t
2377 ...
2378 >>> a = 1
2379 >>> b = 2
2381 We'll use ``integral(f(t), t=a, t=b)`` to denote the definite integral
2382 of f from t=a to t=b. The sample points in ``roots`` are from the
2383 interval [-1, 1], so we'll rewrite the integral with the simple change
2384 of variable::
2386 x = 2/(b - a) * t - (a + b)/(b - a)
2388 with inverse::
2390 t = (b - a)/2 * x + (a + 2)/2
2392 Then::
2394 integral(f(t), a, b) =
2395 (b - a)/2 * integral(f((b-a)/2*x + (a+b)/2), x=-1, x=1)
2397 We can approximate the latter integral with the values returned
2398 by `roots_legendre`.
2400 Map the roots computed above from [-1, 1] to [a, b].
2402 >>> t = (b - a)/2 * roots + (a + b)/2
2404 Approximate the integral as the weighted sum of the function values.
2406 >>> (b - a)/2 * f(t).dot(weights)
2407 2.1931471805599276
2409 Compare that to the exact result, which is 3/2 + log(2):
2411 >>> 1.5 + np.log(2)
2412 2.1931471805599454
2414 """
2415 m = int(n)
2416 if n < 1 or n != m:
2417 raise ValueError("n must be a positive integer.")
2419 mu0 = 2.0
2420 def an_func(k):
2421 return 0.0 * k
2422 def bn_func(k):
2423 return k * np.sqrt(1.0 / (4 * k * k - 1))
2424 f = _ufuncs.eval_legendre
2425 def df(n, x):
2426 return (-n * x * _ufuncs.eval_legendre(n, x) + n * _ufuncs.eval_legendre(n - 1, x)) / (1 - x ** 2)
2427 return _gen_roots_and_weights(m, mu0, an_func, bn_func, f, df, True, mu)
2430def legendre(n, monic=False):
2431 r"""Legendre polynomial.
2433 Defined to be the solution of
2435 .. math::
2436 \frac{d}{dx}\left[(1 - x^2)\frac{d}{dx}P_n(x)\right]
2437 + n(n + 1)P_n(x) = 0;
2439 :math:`P_n(x)` is a polynomial of degree :math:`n`.
2441 Parameters
2442 ----------
2443 n : int
2444 Degree of the polynomial.
2445 monic : bool, optional
2446 If `True`, scale the leading coefficient to be 1. Default is
2447 `False`.
2449 Returns
2450 -------
2451 P : orthopoly1d
2452 Legendre polynomial.
2454 Notes
2455 -----
2456 The polynomials :math:`P_n` are orthogonal over :math:`[-1, 1]`
2457 with weight function 1.
2459 Examples
2460 --------
2461 Generate the 3rd-order Legendre polynomial 1/2*(5x^3 + 0x^2 - 3x + 0):
2463 >>> from scipy.special import legendre
2464 >>> legendre(3)
2465 poly1d([ 2.5, 0. , -1.5, 0. ])
2467 """
2468 if n < 0:
2469 raise ValueError("n must be nonnegative.")
2471 if n == 0:
2472 n1 = n + 1
2473 else:
2474 n1 = n
2475 x, w = roots_legendre(n1)
2476 if n == 0:
2477 x, w = [], []
2478 hn = 2.0 / (2 * n + 1)
2479 kn = _gam(2 * n + 1) / _gam(n + 1)**2 / 2.0**n
2480 p = orthopoly1d(x, w, hn, kn, wfunc=lambda x: 1.0, limits=(-1, 1),
2481 monic=monic,
2482 eval_func=lambda x: _ufuncs.eval_legendre(n, x))
2483 return p
2485# Shifted Legendre P^*_n(x)
2488def roots_sh_legendre(n, mu=False):
2489 r"""Gauss-Legendre (shifted) quadrature.
2491 Compute the sample points and weights for Gauss-Legendre
2492 quadrature. The sample points are the roots of the nth degree
2493 shifted Legendre polynomial :math:`P^*_n(x)`. These sample points
2494 and weights correctly integrate polynomials of degree :math:`2n -
2495 1` or less over the interval :math:`[0, 1]` with weight function
2496 :math:`w(x) = 1.0`. See 2.2.11 in [AS]_ for details.
2498 Parameters
2499 ----------
2500 n : int
2501 quadrature order
2502 mu : bool, optional
2503 If True, return the sum of the weights, optional.
2505 Returns
2506 -------
2507 x : ndarray
2508 Sample points
2509 w : ndarray
2510 Weights
2511 mu : float
2512 Sum of the weights
2514 See Also
2515 --------
2516 scipy.integrate.quadrature
2517 scipy.integrate.fixed_quad
2519 References
2520 ----------
2521 .. [AS] Milton Abramowitz and Irene A. Stegun, eds.
2522 Handbook of Mathematical Functions with Formulas,
2523 Graphs, and Mathematical Tables. New York: Dover, 1972.
2525 """
2526 x, w = roots_legendre(n)
2527 x = (x + 1) / 2
2528 w /= 2
2529 if mu:
2530 return x, w, 1.0
2531 else:
2532 return x, w
2535def sh_legendre(n, monic=False):
2536 r"""Shifted Legendre polynomial.
2538 Defined as :math:`P^*_n(x) = P_n(2x - 1)` for :math:`P_n` the nth
2539 Legendre polynomial.
2541 Parameters
2542 ----------
2543 n : int
2544 Degree of the polynomial.
2545 monic : bool, optional
2546 If `True`, scale the leading coefficient to be 1. Default is
2547 `False`.
2549 Returns
2550 -------
2551 P : orthopoly1d
2552 Shifted Legendre polynomial.
2554 Notes
2555 -----
2556 The polynomials :math:`P^*_n` are orthogonal over :math:`[0, 1]`
2557 with weight function 1.
2559 """
2560 if n < 0:
2561 raise ValueError("n must be nonnegative.")
2563 def wfunc(x):
2564 return 0.0 * x + 1.0
2565 if n == 0:
2566 return orthopoly1d([], [], 1.0, 1.0, wfunc, (0, 1), monic,
2567 lambda x: _ufuncs.eval_sh_legendre(n, x))
2568 x, w = roots_sh_legendre(n)
2569 hn = 1.0 / (2 * n + 1.0)
2570 kn = _gam(2 * n + 1) / _gam(n + 1)**2
2571 p = orthopoly1d(x, w, hn, kn, wfunc, limits=(0, 1), monic=monic,
2572 eval_func=lambda x: _ufuncs.eval_sh_legendre(n, x))
2573 return p
2576# Make the old root function names an alias for the new ones
2577_modattrs = globals()
2578for newfun, oldfun in _rootfuns_map.items():
2579 _modattrs[oldfun] = _modattrs[newfun]
2580 __all__.append(oldfun)