# public scipy /scipy

### Subversion checkout URL

You can clone with HTTPS or Subversion.

  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 1 """ 2 A collection of functions to find the weights and abscissas for 3 Gaussian Quadrature. 4 5 These calculations are done by finding the eigenvalues of a 6 tridiagonal matrix whose entries are dependent on the coefficients 7 in the recursion formula for the orthogonal polynomials with the 8 corresponding weighting function over the interval. 9 10 Many recursion relations for orthogonal polynomials are given: 11  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 12 .. math:: 13 14 a1n f_{n+1} (x) = (a2n + a3n x ) f_n (x) - a4n f_{n-1} (x)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 15 16 The recursion relation of interest is 17  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 18 .. math:: 19 20 P_{n+1} (x) = (x - A_n) P_n (x) - B_n P_{n-1} (x)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 21  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 22 where :math:P has a different normalization than :math:f.  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 23 24 The coefficients can be found as: 25  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 26 .. math:: 27 28 A_n = -a2n / a3n 29 \\qquad 30 B_n = ( a4n / a3n \\sqrt{h_n-1 / h_n})^2 31 32 where  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 33  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 34 .. math:: 35 36 h_n = \\int_a^b w(x) f_n(x)^2  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 37 38 assume:  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 39 40 .. math:: 41 42 P_0 (x) = 1 43 \\qquad 44 P_{-1} (x) == 0  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 45  65f545d7 » fperez  2009-01-07 Add bibliographical references contributed by Alan Isaac after discus… 46 For the mathematical background, see [golub.welsch-1969-mathcomp]_ and 47 [abramowitz.stegun-1965]_.  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 48  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 49 Functions::  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 50 51 gen_roots_and_weights -- Generic roots and weights.  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 52 j_roots -- Jacobi  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 53 js_roots -- Shifted Jacobi 54 la_roots -- Generalized Laguerre 55 h_roots -- Hermite 56 he_roots -- Hermite (unit-variance) 57 cg_roots -- Ultraspherical (Gegenbauer) 58 t_roots -- Chebyshev of the first kind 59 u_roots -- Chebyshev of the second kind 60 c_roots -- Chebyshev of the first kind ([-2,2] interval) 61 s_roots -- Chebyshev of the second kind ([-2,2] interval) 62 ts_roots -- Shifted Chebyshev of the first kind. 63 us_roots -- Shifted Chebyshev of the second kind. 64 p_roots -- Legendre 65 ps_roots -- Shifted Legendre 66 l_roots -- Laguerre  65f545d7 » fperez  2009-01-07 Add bibliographical references contributed by Alan Isaac after discus… 67 68 69 .. [golub.welsch-1969-mathcomp] 70 Golub, Gene H, and John H Welsch. 1969. Calculation of Gauss 71 Quadrature Rules. *Mathematics of Computation* 23, 221-230+s1--s10. 72 73 .. [abramowitz.stegun-1965] 74 Abramowitz, Milton, and Irene A Stegun. (1965) *Handbook of 75 Mathematical Functions: with Formulas, Graphs, and Mathematical 76 Tables*. Gaithersburg, MD: National Bureau of Standards. 77 http://www.math.sfu.ca/~cbm/aands/  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 78  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 79 """  cdcb380e » stefanv  2009-01-13 Apply documentation patch. 80 # 81 # Author: Travis Oliphant 2000 82 # Updated Sep. 2003 (fixed bugs --- tested to be accurate)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 83  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 84 # Scipy imports. 85 import numpy as np 86 from numpy import all, any, exp, inf, pi, sqrt  27723486 » rkern  2007-04-15 Remove dependency on scipy.linalg in favor of numpy.dual 87 from numpy.dual import eig  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 88 89 # Local imports.  a6270c97 » Travis Oliphant  2005-10-06 Io and special now build. 90 import _cephes as cephes  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 91 _gam = cephes.gamma 92  f791fa78 » pv  2009-10-26 special: add separate routines to evaluate orthogonal polynomial values 93 __all__ = ['legendre', 'chebyt', 'chebyu', 'chebyc', 'chebys', 94 'jacobi', 'laguerre', 'genlaguerre', 'hermite', 'hermitenorm', 95 'gegenbauer', 'sh_legendre', 'sh_chebyt', 'sh_chebyu', 'sh_jacobi',  b6184152 » Travis Oliphant  2009-12-09 Add _roots names to the orthogonal.py file 96 'p_roots', 'ps_roots', 'j_roots', 'js_roots', 'l_roots', 'la_roots', 97 'he_roots', 'ts_roots', 'us_roots', 's_roots', 't_roots', 'u_roots', 98 'c_roots', 'cg_roots', 'h_roots',  f791fa78 » pv  2009-10-26 special: add separate routines to evaluate orthogonal polynomial values 99 'eval_legendre', 'eval_chebyt', 'eval_chebyu', 'eval_chebyc', 100 'eval_chebys', 'eval_jacobi', 'eval_laguerre', 'eval_genlaguerre', 101 'eval_hermite', 'eval_hermitenorm', 'eval_gegenbauer', 102 'eval_sh_legendre', 'eval_sh_chebyt', 'eval_sh_chebyu', 103 'eval_sh_jacobi', 'poch', 'binom'] 104  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 105 def poch(z, m):  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 106 """Pochhammer symbol (z)_m = (z)(z+1)....(z+m-1) = gamma(z+m)/gamma(z)""" 107 return _gam(z+m) / _gam(z)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 108  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 109 class orthopoly1d(np.poly1d):  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 110 def __init__(self, roots, weights=None, hn=1.0, kn=1.0, wfunc=None, limits=None, monic=0,eval_func=None):  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 111 np.poly1d.__init__(self, roots, r=1)  e9c02633 » cookedm  2006-03-14 replace tab indentations found with tabnanny 112 equiv_weights = [weights[k] / wfunc(roots[k]) for k in range(len(roots))]  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 113 self.__dict__['weights'] = np.array(zip(roots,weights,equiv_weights))  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 114 self.__dict__['weight_func'] = wfunc 115 self.__dict__['limits'] = limits 116 mu = sqrt(hn) 117 if monic:  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 118 evf = eval_func  e9554e8d » pv  2009-12-12 BUG: special: fix issues in orthogonal/orthogonal_eval integration 119 if evf: 120 eval_func = lambda x: evf(x)/kn  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 121 mu = mu / abs(kn) 122 kn = 1.0  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 123 self.__dict__['normcoef'] = mu  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 124 self.__dict__['coeffs'] *= kn  b9c4f2d9 » Travis Oliphant  2002-04-13 Fixed uses of eig in scipy to use the new linalg standard where the e… 125  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 126 # Note: eval_func will be discarded on arithmetic 127 self.__dict__['_eval_func'] = eval_func 128 129 def __call__(self, v): 130 if self._eval_func and (isinstance(v, np.ndarray) or np.isscalar(v)): 131 return self._eval_func(v) 132 else: 133 return np.poly1d.__call__(self, v) 134 135 def _scale(self, p): 136 if p == 1.0: 137 return 138 self.__dict__['coeffs'] *= p 139 evf = self.__dict__['_eval_func']  e9554e8d » pv  2009-12-12 BUG: special: fix issues in orthogonal/orthogonal_eval integration 140 if evf: 141 self.__dict__['_eval_func'] = lambda x: evf(x) * p  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 142 self.__dict__['normcoef'] *= p  bde258f8 » pearu  2003-03-28 Fixed strict scipy-dependence of the special package 143  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 144 def gen_roots_and_weights(n, an_func, sqrt_bn_func, mu):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 145 """[x,w] = gen_roots_and_weights(n,an_func,sqrt_bn_func,mu) 146  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 147 Returns the roots (x) of an nth order orthogonal polynomial,  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 148 and weights (w) to use in appropriate Gaussian quadrature with that 149 orthogonal polynomial. 150 151 The polynomials have the recurrence relation 152 P_n+1(x) = (x - A_n) P_n(x) - B_n P_n-1(x) 153 154 an_func(n) should return A_n 155 sqrt_bn_func(n) should return sqrt(B_n) 156 mu ( = h_0 ) is the integral of the weight over the orthogonal interval 157 """  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 158 nn = np.arange(1.0,n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 159 sqrt_bn = sqrt_bn_func(nn)  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 160 an = an_func(np.concatenate(([0], nn)))  a779ecc4 » Jarrod Millman  2007-10-29 ran reindent.py to clean up whitespace 161 x, v = eig((np.diagflat(an) + 162 np.diagflat(sqrt_bn,1) +  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 163 np.diagflat(sqrt_bn,-1)))  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 164 answer = []  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 165 sortind = x.real.argsort() 166 answer.append(x[sortind]) 167 answer.append((mu*v[0]**2)[sortind])  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 168 return answer  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 169 170 # Jacobi Polynomials 1 P^(alpha,beta)_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 171 def j_roots(n, alpha, beta, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 172 """[x,w] = j_roots(n,alpha,beta)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 173 174 Returns the roots (x) of the nth order Jacobi polynomial, P^(alpha,beta)_n(x) 175 and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting 176 function (1-x)**alpha (1+x)**beta with alpha,beta > -1. 177 """  240b4a01 » jmiller  2005-05-09 Modifications to support scipy.special for numarray. 178 if any(alpha <= -1) or any(beta <= -1):  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 179 raise ValueError("alpha and beta must be greater than -1.")  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 180 if n < 1: 181 raise ValueError("n must be positive.")  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 182  bd311699 » rgommers  2011-01-15 TST: silence floating point test noise in scipy.special. 183 olderr = np.seterr(all='ignore') 184 try: 185 (p,q) = (alpha,beta) 186 # from recurrence relations 187 sbn_J = lambda k: 2.0/(2.0*k+p+q)*sqrt((k+p)*(k+q)/(2*k+q+p+1)) * \ 188 (np.where(k==1,1.0,sqrt(k*(k+p+q)/(2.0*k+p+q-1)))) 189 if any(p == q): # XXX any or all??? 190 an_J = lambda k: 0.0*k 191 else: 192 an_J = lambda k: np.where(k==0,(q-p)/(p+q+2.0), 193 (q*q - p*p)/((2.0*k+p+q)*(2.0*k+p+q+2))) 194 g = cephes.gamma 195 mu0 = 2.0**(p+q+1)*g(p+1)*g(q+1)/(g(p+q+2)) 196 val = gen_roots_and_weights(n,an_J,sbn_J,mu0) 197 finally: 198 np.seterr(**olderr) 199  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 200 if mu: 201 return val + [mu0] 202 else: 203 return val 204  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 205 def jacobi(n, alpha, beta, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 206 """Returns the nth order Jacobi polynomial, P^(alpha,beta)_n(x) 207 orthogonal over [-1,1] with weighting function 208 (1-x)**alpha (1+x)**beta with alpha,beta > -1. 209 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 210 if n < 0: 211 raise ValueError("n must be nonnegative.") 212  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 213 wfunc = lambda x: (1-x)**alpha * (1+x)**beta  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 214 if n==0:  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 215 return orthopoly1d([],[],1.0,1.0,wfunc,(-1,1),monic,  e9554e8d » pv  2009-12-12 BUG: special: fix issues in orthogonal/orthogonal_eval integration 216 eval_func=np.ones_like)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 217 x,w,mu = j_roots(n,alpha,beta,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 218 ab1 = alpha+beta+1.0 219 hn = 2**ab1/(2*n+ab1)*_gam(n+alpha+1) 220 hn *= _gam(n+beta+1.0) / _gam(n+1) / _gam(n+ab1)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 221 kn = _gam(2*n+ab1)/2.0**n / _gam(n+1) / _gam(n+ab1) 222 # here kn = coefficient on x^n term  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 223 p = orthopoly1d(x,w,hn,kn,wfunc,(-1,1),monic, 224 lambda x: eval_jacobi(n,alpha,beta,x))  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 225 return p 226 227 # Jacobi Polynomials shifted G_n(p,q,x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 228 def js_roots(n, p1, q1, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 229 """[x,w] = js_roots(n,p,q)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 230 231 Returns the roots (x) of the nth order shifted Jacobi polynomial, G_n(p,q,x), 232 and weights (w) to use in Gaussian Quadrature over [0,1] with weighting 233 function (1-x)**(p-q) x**(q-1) with p-q > -1 and q > 0. 234 """ 235 # from recurrence relation  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 236 if not ( any((p1 - q1) > -1) and any(q1 > 0) ): 237 raise ValueError("(p - q) > -1 and q > 0 please.")  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 238 if n <= 0:  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 239 raise ValueError("n must be positive.")  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 240  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 241 p,q = p1,q1  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 242  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 243 sbn_Js = lambda k: sqrt(np.where(k==1,q*(p-q+1.0)/(p+2.0), \  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 244 k*(k+q-1.0)*(k+p-1.0)*(k+p-q) \ 245 / ((2.0*k+p-2) * (2.0*k+p))))/(2*k+p-1.0)  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 246 an_Js = lambda k: np.where(k==0,q/(p+1.0),(2.0*k*(k+p)+q*(p-1.0)) / ((2.0*k+p+1.0)*(2*k+p-1.0)))  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 247 248 # could also use definition 249 # Gn(p,q,x) = constant_n * P^(p-q,q-1)_n(2x-1) 250 # so roots of Gn(p,q,x) are (roots of P^(p-q,q-1)_n + 1) / 2.0 251 g = _gam  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 252 # integral of weight over interval 253 mu0 = g(q)*g(p-q+1)/g(p+1) 254 val = gen_roots_and_weights(n,an_Js,sbn_Js,mu0) 255 if mu: 256 return val + [mu0] 257 else: 258 return val  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 259 # What code would look like using jacobi polynomial roots 260 #if mu: 261 # [x,w,mut] = j_roots(n,p-q,q-1,mu=1) 262 # return [(x+1)/2.0,w,mu0] 263 #else: 264 # [x,w] = j_roots(n,p-q,q-1,mu=0)  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 265 # return [(x+1)/2.0,w]  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 266 267 def sh_jacobi(n, p, q, monic=0): 268 """Returns the nth order Jacobi polynomial, G_n(p,q,x) 269 orthogonal over [0,1] with weighting function  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 270 (1-x)**(p-q) (x)**(q-1) with p>q-1 and q > 0.  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 271 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 272 if n < 0: 273 raise ValueError("n must be nonnegative.") 274  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 275 wfunc = lambda x: (1.0-x)**(p-q) * (x)**(q-1.)  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 276 if n==0:  e9554e8d » pv  2009-12-12 BUG: special: fix issues in orthogonal/orthogonal_eval integration 277 return orthopoly1d([],[],1.0,1.0,wfunc,(-1,1),monic, 278 eval_func=np.ones_like)  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 279 n1 = n  470c28c7 » Travis Oliphant  2003-09-24 Fixed some segfaulting and other special function tests. 280 x,w,mu0 = js_roots(n1,p,q,mu=1)  31d6178d » Travis Oliphant  2002-04-04 Fixed spheroidal wave functions. Added lu_solve, lu_factor, cho_solve… 281 hn = _gam(n+1)*_gam(n+q)*_gam(n+p)*_gam(n+p-q+1) 282 hn /= (2*n+p)*(_gam(2*n+p)**2)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 283 # kn = 1.0 in standard form so monic is redundant. Kept for compatibility. 284 kn = 1.0  e9554e8d » pv  2009-12-12 BUG: special: fix issues in orthogonal/orthogonal_eval integration 285 pp = orthopoly1d(x,w,hn,kn,wfunc=wfunc,limits=(0,1),monic=monic, 286 eval_func=lambda x: eval_sh_jacobi(n, p, q, x)) 287 return pp  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 288 289 # Generalized Laguerre L^(alpha)_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 290 def la_roots(n, alpha, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 291 """[x,w] = la_roots(n,alpha)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 292 293 Returns the roots (x) of the nth order generalized (associated) Laguerre 294 polynomial, L^(alpha)_n(x), and weights (w) to use in Gaussian quadrature over 295 [0,inf] with weighting function exp(-x) x**alpha with alpha > -1. 296 """  240b4a01 » jmiller  2005-05-09 Modifications to support scipy.special for numarray. 297 if not all(alpha > -1):  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 298 raise ValueError("alpha > -1")  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 299 if n < 1: 300 raise ValueError("n must be positive.") 301  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 302 (p,q) = (alpha,0.0) 303 sbn_La = lambda k: -sqrt(k*(k + p)) # from recurrence relation  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 304 an_La = lambda k: 2*k + p + 1 305 mu0 = cephes.gamma(alpha+1) # integral of weight over interval  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 306 val = gen_roots_and_weights(n,an_La,sbn_La,mu0) 307 if mu: 308 return val + [mu0] 309 else: 310 return val 311  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 312 def genlaguerre(n, alpha, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 313 """Returns the nth order generalized (associated) Laguerre polynomial, 314 L^(alpha)_n(x), orthogonal over [0,inf) with weighting function 315 exp(-x) x**alpha with alpha > -1 316 """  240b4a01 » jmiller  2005-05-09 Modifications to support scipy.special for numarray. 317 if any(alpha <= -1):  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 318 raise ValueError("alpha must be > -1")  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 319 if n < 0: 320 raise ValueError("n must be nonnegative.") 321  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 322 if n==0: n1 = n+1 323 else: n1 = n  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 324 x,w,mu0 = la_roots(n1,alpha,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 325 wfunc = lambda x: exp(-x) * x**alpha 326 if n==0: x,w = [],[] 327 hn = _gam(n+alpha+1)/_gam(n+1)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 328 kn = (-1)**n / _gam(n+1)  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 329 p = orthopoly1d(x,w,hn,kn,wfunc,(0,inf),monic, 330 lambda x: eval_genlaguerre(n,alpha,x))  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 331 return p 332  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 333 # Laguerre L_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 334 def l_roots(n, mu=0):  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 335 """[x,w] = l_roots(n) 336 337 Returns the roots (x) of the nth order Laguerre polynomial, L_n(x), 338 and weights (w) to use in Gaussian Quadrature over [0,inf] with weighting 339 function exp(-x). 340 """ 341 return la_roots(n,0.0,mu=mu) 342  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 343 def laguerre(n, monic=0):  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 344 """Return the nth order Laguerre polynoimal, L_n(x), orthogonal over 345 [0,inf) with weighting function exp(-x) 346 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 347 if n < 0: 348 raise ValueError("n must be nonnegative.") 349  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 350 if n==0: n1 = n+1 351 else: n1 = n 352 x,w,mu0 = l_roots(n1,mu=1) 353 if n==0: x,w = [],[] 354 hn = 1.0 355 kn = (-1)**n / _gam(n+1)  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 356 p = orthopoly1d(x,w,hn,kn,lambda x: exp(-x),(0,inf),monic, 357 lambda x: eval_laguerre(n,x))  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 358 return p  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 359 360  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 361 # Hermite 1 H_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 362 def h_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 363 """[x,w] = h_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 364 365 Returns the roots (x) of the nth order Hermite polynomial, 366 H_n(x), and weights (w) to use in Gaussian Quadrature over 367 [-inf,inf] with weighting function exp(-x**2). 368 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 369 if n < 1: 370 raise ValueError("n must be positive.") 371  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 372 sbn_H = lambda k: sqrt(k/2) # from recurrence relation  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 373 an_H = lambda k: 0*k 374 mu0 = sqrt(pi) # integral of weight over interval  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 375 val = gen_roots_and_weights(n,an_H,sbn_H,mu0) 376 if mu: 377 return val + [mu0] 378 else: 379 return val 380  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 381 def hermite(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 382 """Return the nth order Hermite polynomial, H_n(x), orthogonal over 383 (-inf,inf) with weighting function exp(-x**2) 384 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 385 if n < 0: 386 raise ValueError("n must be nonnegative.") 387  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 388 if n==0: n1 = n+1 389 else: n1 = n  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 390 x,w,mu0 = h_roots(n1,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 391 wfunc = lambda x: exp(-x*x) 392 if n==0: x,w = [],[] 393 hn = 2**n * _gam(n+1)*sqrt(pi)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 394 kn = 2**n  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 395 p = orthopoly1d(x,w,hn,kn,wfunc,(-inf,inf),monic, 396 lambda x: eval_hermite(n,x))  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 397 return p  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 398  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 399 # Hermite 2 He_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 400 def he_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 401 """[x,w] = he_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 402 403 Returns the roots (x) of the nth order Hermite polynomial, 404 He_n(x), and weights (w) to use in Gaussian Quadrature over 405 [-inf,inf] with weighting function exp(-(x/2)**2). 406 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 407 if n < 1: 408 raise ValueError("n must be positive.") 409  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 410 sbn_He = lambda k: sqrt(k) # from recurrence relation  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 411 an_He = lambda k: 0*k 412 mu0 = sqrt(2*pi) # integral of weight over interval  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 413 val = gen_roots_and_weights(n,an_He,sbn_He,mu0) 414 if mu: 415 return val + [mu0] 416 else: 417 return val 418  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 419 def hermitenorm(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 420 """Return the nth order normalized Hermite polynomial, He_n(x), orthogonal 421 over (-inf,inf) with weighting function exp(-(x/2)**2) 422 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 423 if n < 0: 424 raise ValueError("n must be nonnegative.") 425  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 426 if n==0: n1 = n+1 427 else: n1 = n  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 428 x,w,mu0 = he_roots(n1,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 429 wfunc = lambda x: exp(-x*x/4.0) 430 if n==0: x,w = [],[]  31d6178d » Travis Oliphant  2002-04-04 Fixed spheroidal wave functions. Added lu_solve, lu_factor, cho_solve… 431 hn = sqrt(2*pi)*_gam(n+1)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 432 kn = 1.0  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 433 p = orthopoly1d(x,w,hn,kn,wfunc=wfunc,limits=(-inf,inf),monic=monic,  1ab0f969 » pv  2009-12-12 BUG: hermiternorm had eval_func incorrectly set -- add tests for call… 434 eval_func=lambda x: eval_hermitenorm(n,x))  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 435 return p 436 437 ## The remainder of the polynomials can be derived from the ones above. 438 439 # Ultraspherical (Gegenbauer) C^(alpha)_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 440 def cg_roots(n, alpha, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 441 """[x,w] = cg_roots(n,alpha)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 442 443 Returns the roots (x) of the nth order Ultraspherical (Gegenbauer) 444 polynomial, C^(alpha)_n(x), and weights (w) to use in Gaussian Quadrature 445 over [-1,1] with weighting function (1-x**2)**(alpha-1/2) with alpha>-1/2. 446 """  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 447 return j_roots(n,alpha-0.5,alpha-0.5,mu=mu)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 448  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 449 def gegenbauer(n, alpha, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 450 """Return the nth order Gegenbauer (ultraspherical) polynomial, 451 C^(alpha)_n(x), orthogonal over [-1,1] with weighting function 452 (1-x**2)**(alpha-1/2) with alpha > -1/2 453 """  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 454 base = jacobi(n,alpha-0.5,alpha-0.5,monic=monic) 455 if monic: 456 return base 457 # Abrahmowitz and Stegan 22.5.20 458 factor = _gam(2*alpha+n)*_gam(alpha+0.5) / _gam(2*alpha) / _gam(alpha+0.5+n)  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 459 base._scale(factor) 460 return base  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 461 462 # Chebyshev of the first kind: T_n(x) = n! sqrt(pi) / _gam(n+1./2)* P^(-1/2,-1/2)_n(x) 463 # Computed anew.  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 464 def t_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 465 """[x,w] = t_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 466 467 Returns the roots (x) of the nth order Chebyshev (of the first kind) 468 polynomial, T_n(x), and weights (w) to use in Gaussian Quadrature 469 over [-1,1] with weighting function (1-x**2)**(-1/2). 470 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 471 if n < 1: 472 raise ValueError("n must be positive.") 473  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 474 # from recurrence relation  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 475 sbn_J = lambda k: np.where(k==1,sqrt(2)/2.0,0.5)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 476 an_J = lambda k: 0.0*k 477 g = cephes.gamma 478 mu0 = pi 479 val = gen_roots_and_weights(n,an_J,sbn_J,mu0) 480 if mu: 481 return val + [mu0] 482 else: 483 return val  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 484  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 485 def chebyt(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 486 """Return nth order Chebyshev polynomial of first kind, Tn(x). Orthogonal 487 over [-1,1] with weight function (1-x**2)**(-1/2). 488 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 489 if n < 0: 490 raise ValueError("n must be nonnegative.") 491  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 492 wfunc = lambda x: 1.0/sqrt(1-x*x)  2bbb89ec » rkern  2006-09-24 Use modern numpy idioms. 493 if n==0:  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 494 return orthopoly1d([],[],pi,1.0,wfunc,(-1,1),monic, 495 lambda x: eval_chebyt(n,x))  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 496 n1 = n 497 x,w,mu = t_roots(n1,mu=1) 498 hn = pi/2 499 kn = 2**(n-1)  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 500 p = orthopoly1d(x,w,hn,kn,wfunc,(-1,1),monic, 501 lambda x: eval_chebyt(n,x))  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 502 return p  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 503  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 504 # Chebyshev of the second kind 505 # U_n(x) = (n+1)! sqrt(pi) / (2*_gam(n+3./2)) * P^(1/2,1/2)_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 506 def u_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 507 """[x,w] = u_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 508 509 Returns the roots (x) of the nth order Chebyshev (of the second kind) 510 polynomial, U_n(x), and weights (w) to use in Gaussian Quadrature 511 over [-1,1] with weighting function (1-x**2)**1/2. 512 """  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 513 return j_roots(n,0.5,0.5,mu=mu)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 514  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 515 def chebyu(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 516 """Return nth order Chebyshev polynomial of second kind, Un(x). Orthogonal 517 over [-1,1] with weight function (1-x**2)**(1/2). 518 """  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 519 base = jacobi(n,0.5,0.5,monic=monic) 520 if monic: 521 return base 522 factor = sqrt(pi)/2.0*_gam(n+2) / _gam(n+1.5)  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 523 base._scale(factor) 524 return base  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 525 526 # Chebyshev of the first kind C_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 527 def c_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 528 """[x,w] = c_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 529 530 Returns the roots (x) of the nth order Chebyshev (of the first kind) 531 polynomial, C_n(x), and weights (w) to use in Gaussian Quadrature 532 over [-2,2] with weighting function (1-(x/2)**2)**(-1/2). 533 """ 534 if mu:  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 535 [x,w,mu0] = j_roots(n,-0.5,-0.5,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 536 return [x*2,w,mu0] 537 else:  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 538 [x,w] = j_roots(n,-0.5,-0.5,mu=0)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 539 return [x*2,w] 540  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 541 def chebyc(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 542 """Return nth order Chebyshev polynomial of first kind, Cn(x). Orthogonal 543 over [-2,2] with weight function (1-(x/2)**2)**(-1/2). 544 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 545 if n < 0: 546 raise ValueError("n must be nonnegative.") 547  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 548 if n==0: n1 = n+1 549 else: n1 = n  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 550 x,w,mu0 = c_roots(n1,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 551 if n==0: x,w = [],[]  31d6178d » Travis Oliphant  2002-04-04 Fixed spheroidal wave functions. Added lu_solve, lu_factor, cho_solve… 552 hn = 4*pi * ((n==0)+1)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 553 kn = 1.0 554 p = orthopoly1d(x,w,hn,kn,wfunc=lambda x: 1.0/sqrt(1-x*x/4.0),limits=(-2,2),monic=monic)  31d6178d » Travis Oliphant  2002-04-04 Fixed spheroidal wave functions. Added lu_solve, lu_factor, cho_solve… 555 if not monic:  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 556 p._scale(2.0/p(2))  e9554e8d » pv  2009-12-12 BUG: special: fix issues in orthogonal/orthogonal_eval integration 557 p.__dict__['_eval_func'] = lambda x: eval_chebyc(n,x)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 558 return p 559 560 # Chebyshev of the second kind S_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 561 def s_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 562 """[x,w] = s_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 563 564 Returns the roots (x) of the nth order Chebyshev (of the second kind) 565 polynomial, S_n(x), and weights (w) to use in Gaussian Quadrature 566 over [-2,2] with weighting function (1-(x/2)**2)**1/2. 567 """ 568 if mu:  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 569 [x,w,mu0] = j_roots(n,0.5,0.5,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 570 return [x*2,w,mu0] 571 else:  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 572 [x,w] = j_roots(n,0.5,0.5,mu=0)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 573 return [x*2,w] 574  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 575 def chebys(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 576 """Return nth order Chebyshev polynomial of second kind, Sn(x). Orthogonal 577 over [-2,2] with weight function (1-(x/)**2)**(1/2). 578 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 579 if n < 0: 580 raise ValueError("n must be nonnegative.") 581  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 582 if n==0: n1 = n+1 583 else: n1 = n  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 584 x,w,mu0 = s_roots(n1,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 585 if n==0: x,w = [],[]  31d6178d » Travis Oliphant  2002-04-04 Fixed spheroidal wave functions. Added lu_solve, lu_factor, cho_solve… 586 hn = pi  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 587 kn = 1.0  e9554e8d » pv  2009-12-12 BUG: special: fix issues in orthogonal/orthogonal_eval integration 588 p = orthopoly1d(x,w,hn,kn,wfunc=lambda x: sqrt(1-x*x/4.0),limits=(-2,2),monic=monic)  31d6178d » Travis Oliphant  2002-04-04 Fixed spheroidal wave functions. Added lu_solve, lu_factor, cho_solve… 589 if not monic:  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 590 factor = (n+1.0)/p(2) 591 p._scale(factor)  e9554e8d » pv  2009-12-12 BUG: special: fix issues in orthogonal/orthogonal_eval integration 592 p.__dict__['_eval_func'] = lambda x: eval_chebys(n,x)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 593 return p 594 595 # Shifted Chebyshev of the first kind T^*_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 596 def ts_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 597 """[x,w] = ts_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 598 599 Returns the roots (x) of the nth order shifted Chebyshev (of the first kind) 600 polynomial, T^*_n(x), and weights (w) to use in Gaussian Quadrature 601 over [0,1] with weighting function (x-x**2)**(-1/2). 602 """  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 603 return js_roots(n,0.0,0.5,mu=mu)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 604  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 605 def sh_chebyt(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 606 """Return nth order shifted Chebyshev polynomial of first kind, Tn(x). 607 Orthogonal over [0,1] with weight function (x-x**2)**(-1/2). 608 """  dae01c63 » Travis Oliphant  2003-09-26 Fixed orthogonal polynomials and added tests to verify them against e… 609 base = sh_jacobi(n,0.0,0.5,monic=monic)  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 610 if monic: 611 return base  dae01c63 » Travis Oliphant  2003-09-26 Fixed orthogonal polynomials and added tests to verify them against e… 612 if n > 0: 613 factor = 4**n / 2.0 614 else: 615 factor = 1.0  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 616 base._scale(factor) 617 return base  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 618  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 619 620 # Shifted Chebyshev of the second kind U^*_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 621 def us_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 622 """[x,w] = us_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 623 624 Returns the roots (x) of the nth order shifted Chebyshev (of the second kind) 625 polynomial, U^*_n(x), and weights (w) to use in Gaussian Quadrature 626 over [0,1] with weighting function (x-x**2)**1/2. 627 """  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 628 return js_roots(n,2.0,1.5,mu=mu)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 629  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 630 def sh_chebyu(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 631 """Return nth order shifted Chebyshev polynomial of second kind, Un(x). 632 Orthogonal over [0,1] with weight function (x-x**2)**(1/2). 633 """  dae01c63 » Travis Oliphant  2003-09-26 Fixed orthogonal polynomials and added tests to verify them against e… 634 base = sh_jacobi(n,2.0,1.5,monic=monic) 635 if monic: return base 636 factor = 4**n  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 637 base._scale(factor) 638 return base  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 639  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 640 # Legendre  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 641 def p_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 642 """[x,w] = p_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 643 644 Returns the roots (x) of the nth order Legendre polynomial, P_n(x), 645 and weights (w) to use in Gaussian Quadrature over [-1,1] with weighting 646 function 1. 647 """  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 648 return j_roots(n,0.0,0.0,mu=mu)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 649  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 650 def legendre(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 651 """Returns the nth order Legendre polynomial, P_n(x), orthogonal over 652 [-1,1] with weight function 1. 653 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 654 if n < 0: 655 raise ValueError("n must be nonnegative.") 656  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 657 if n==0: n1 = n+1 658 else: n1 = n  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 659 x,w,mu0 = p_roots(n1,mu=1)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 660 if n==0: x,w = [],[] 661 hn = 2.0/(2*n+1)  74f5ebc0 » Travis Oliphant  2003-09-25 Orthogonal polynomial fixes and nan_to_num fixes for hanging linalg r… 662 kn = _gam(2*n+1)/_gam(n+1)**2 / 2.0**n  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 663 p = orthopoly1d(x,w,hn,kn,wfunc=lambda x: 1.0,limits=(-1,1),monic=monic, 664 eval_func=lambda x: eval_legendre(n,x))  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 665 return p 666 667 # Shifted Legendre P^*_n(x)  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 668 def ps_roots(n, mu=0):  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 669 """[x,w] = ps_roots(n)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 670 671 Returns the roots (x) of the nth order shifted Legendre polynomial, P^*_n(x), 672 and weights (w) to use in Gaussian Quadrature over [0,1] with weighting 673 function 1. 674 """  abf463d3 » eric-jones  2002-04-13 converted upper case function names to lower case names to match the … 675 return js_roots(n,1.0,1.0,mu=mu)  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 676  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 677 def sh_legendre(n, monic=0):  ab1d4868 » Travis Oliphant  2002-04-02 Finished special functions. (Added orthogonal polynomials). 678 """Returns the nth order shifted Legendre polynomial, P^*_n(x), orthogonal 679 over [0,1] with weighting function 1. 680 """  62fa2c26 » warren.weckesser  2011-01-29 ENH: special: Don't use plain assert for argument validation. 681 if n < 0: 682 raise ValueError("n must be nonnegative.") 683  dae01c63 » Travis Oliphant  2003-09-26 Fixed orthogonal polynomials and added tests to verify them against e… 684 wfunc = lambda x: 0.0*x + 1.0  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 685 if n==0: return orthopoly1d([],[],1.0,1.0,wfunc,(0,1),monic, 686 lambda x: eval_sh_legendre(n,x))  dae01c63 » Travis Oliphant  2003-09-26 Fixed orthogonal polynomials and added tests to verify them against e… 687 x,w,mu0 = ps_roots(n,mu=1)  31d6178d » Travis Oliphant  2002-04-04 Fixed spheroidal wave functions. Added lu_solve, lu_factor, cho_solve… 688 hn = 1.0/(2*n+1.0)  dae01c63 » Travis Oliphant  2003-09-26 Fixed orthogonal polynomials and added tests to verify them against e… 689 kn = _gam(2*n+1)/_gam(n+1)**2  cb30b3e6 » pv  2009-10-26 special: Hook up orthopoly evaluation routines to orthopoly1d class 690 p = orthopoly1d(x,w,hn,kn,wfunc,limits=(0,1),monic=monic, 691 eval_func=lambda x: eval_sh_legendre(n,x))  ebea7bf2 » cookedm  2006-03-14 Run reindent.py (script distributed with Python) over the source to r… 692 return p  f791fa78 » pv  2009-10-26 special: add separate routines to evaluate orthogonal polynomial values 693 694 #------------------------------------------------------------------------------ 695 # Vectorized functions for evaluation 696 #------------------------------------------------------------------------------ 697 from orthogonal_eval import \ 698 binom, eval_jacobi, eval_sh_jacobi, eval_gegenbauer, eval_chebyt, \ 699 eval_chebyu, eval_chebys, eval_chebyc, eval_sh_chebyt, eval_sh_chebyu, \ 700 eval_legendre, eval_sh_legendre, eval_genlaguerre, eval_laguerre, \ 701 eval_hermite, eval_hermitenorm