1""" 2This module mainly implements special orthogonal polynomials. 3 4See also functions.combinatorial.numbers which contains some 5combinatorial polynomials. 6 7""" 8 9from sympy.core import Rational 10from sympy.core.function import Function, ArgumentIndexError 11from sympy.core.singleton import S 12from sympy.core.symbol import Dummy 13from sympy.functions.combinatorial.factorials import binomial, factorial, RisingFactorial 14from sympy.functions.elementary.complexes import re 15from sympy.functions.elementary.exponential import exp 16from sympy.functions.elementary.integers import floor 17from sympy.functions.elementary.miscellaneous import sqrt 18from sympy.functions.elementary.trigonometric import cos, sec 19from sympy.functions.special.gamma_functions import gamma 20from sympy.functions.special.hyper import hyper 21 22from sympy.polys.orthopolys import ( 23 jacobi_poly, 24 gegenbauer_poly, 25 chebyshevt_poly, 26 chebyshevu_poly, 27 laguerre_poly, 28 hermite_poly, 29 legendre_poly 30) 31 32_x = Dummy('x') 33 34 35class OrthogonalPolynomial(Function): 36 """Base class for orthogonal polynomials. 37 """ 38 39 @classmethod 40 def _eval_at_order(cls, n, x): 41 if n.is_integer and n >= 0: 42 return cls._ortho_poly(int(n), _x).subs(_x, x) 43 44 def _eval_conjugate(self): 45 return self.func(self.args[0], self.args[1].conjugate()) 46 47#---------------------------------------------------------------------------- 48# Jacobi polynomials 49# 50 51 52class jacobi(OrthogonalPolynomial): 53 r""" 54 Jacobi polynomial $P_n^{\left(\alpha, \beta\right)}(x)$. 55 56 Explanation 57 =========== 58 59 ``jacobi(n, alpha, beta, x)`` gives the nth Jacobi polynomial 60 in x, $P_n^{\left(\alpha, \beta\right)}(x)$. 61 62 The Jacobi polynomials are orthogonal on $[-1, 1]$ with respect 63 to the weight $\left(1-x\right)^\alpha \left(1+x\right)^\beta$. 64 65 Examples 66 ======== 67 68 >>> from sympy import jacobi, S, conjugate, diff 69 >>> from sympy.abc import a, b, n, x 70 71 >>> jacobi(0, a, b, x) 72 1 73 >>> jacobi(1, a, b, x) 74 a/2 - b/2 + x*(a/2 + b/2 + 1) 75 >>> jacobi(2, a, b, x) 76 a**2/8 - a*b/4 - a/8 + b**2/8 - b/8 + x**2*(a**2/8 + a*b/4 + 7*a/8 + b**2/8 + 7*b/8 + 3/2) + x*(a**2/4 + 3*a/4 - b**2/4 - 3*b/4) - 1/2 77 78 >>> jacobi(n, a, b, x) 79 jacobi(n, a, b, x) 80 81 >>> jacobi(n, a, a, x) 82 RisingFactorial(a + 1, n)*gegenbauer(n, 83 a + 1/2, x)/RisingFactorial(2*a + 1, n) 84 85 >>> jacobi(n, 0, 0, x) 86 legendre(n, x) 87 88 >>> jacobi(n, S(1)/2, S(1)/2, x) 89 RisingFactorial(3/2, n)*chebyshevu(n, x)/factorial(n + 1) 90 91 >>> jacobi(n, -S(1)/2, -S(1)/2, x) 92 RisingFactorial(1/2, n)*chebyshevt(n, x)/factorial(n) 93 94 >>> jacobi(n, a, b, -x) 95 (-1)**n*jacobi(n, b, a, x) 96 97 >>> jacobi(n, a, b, 0) 98 gamma(a + n + 1)*hyper((-b - n, -n), (a + 1,), -1)/(2**n*factorial(n)*gamma(a + 1)) 99 >>> jacobi(n, a, b, 1) 100 RisingFactorial(a + 1, n)/factorial(n) 101 102 >>> conjugate(jacobi(n, a, b, x)) 103 jacobi(n, conjugate(a), conjugate(b), conjugate(x)) 104 105 >>> diff(jacobi(n,a,b,x), x) 106 (a/2 + b/2 + n/2 + 1/2)*jacobi(n - 1, a + 1, b + 1, x) 107 108 See Also 109 ======== 110 111 gegenbauer, 112 chebyshevt_root, chebyshevu, chebyshevu_root, 113 legendre, assoc_legendre, 114 hermite, 115 laguerre, assoc_laguerre, 116 sympy.polys.orthopolys.jacobi_poly, 117 sympy.polys.orthopolys.gegenbauer_poly 118 sympy.polys.orthopolys.chebyshevt_poly 119 sympy.polys.orthopolys.chebyshevu_poly 120 sympy.polys.orthopolys.hermite_poly 121 sympy.polys.orthopolys.legendre_poly 122 sympy.polys.orthopolys.laguerre_poly 123 124 References 125 ========== 126 127 .. [1] https://en.wikipedia.org/wiki/Jacobi_polynomials 128 .. [2] http://mathworld.wolfram.com/JacobiPolynomial.html 129 .. [3] http://functions.wolfram.com/Polynomials/JacobiP/ 130 131 """ 132 133 @classmethod 134 def eval(cls, n, a, b, x): 135 # Simplify to other polynomials 136 # P^{a, a}_n(x) 137 if a == b: 138 if a == Rational(-1, 2): 139 return RisingFactorial(S.Half, n) / factorial(n) * chebyshevt(n, x) 140 elif a.is_zero: 141 return legendre(n, x) 142 elif a == S.Half: 143 return RisingFactorial(3*S.Half, n) / factorial(n + 1) * chebyshevu(n, x) 144 else: 145 return RisingFactorial(a + 1, n) / RisingFactorial(2*a + 1, n) * gegenbauer(n, a + S.Half, x) 146 elif b == -a: 147 # P^{a, -a}_n(x) 148 return gamma(n + a + 1) / gamma(n + 1) * (1 + x)**(a/2) / (1 - x)**(a/2) * assoc_legendre(n, -a, x) 149 150 if not n.is_Number: 151 # Symbolic result P^{a,b}_n(x) 152 # P^{a,b}_n(-x) ---> (-1)**n * P^{b,a}_n(-x) 153 if x.could_extract_minus_sign(): 154 return S.NegativeOne**n * jacobi(n, b, a, -x) 155 # We can evaluate for some special values of x 156 if x.is_zero: 157 return (2**(-n) * gamma(a + n + 1) / (gamma(a + 1) * factorial(n)) * 158 hyper([-b - n, -n], [a + 1], -1)) 159 if x == S.One: 160 return RisingFactorial(a + 1, n) / factorial(n) 161 elif x is S.Infinity: 162 if n.is_positive: 163 # Make sure a+b+2*n \notin Z 164 if (a + b + 2*n).is_integer: 165 raise ValueError("Error. a + b + 2*n should not be an integer.") 166 return RisingFactorial(a + b + n + 1, n) * S.Infinity 167 else: 168 # n is a given fixed integer, evaluate into polynomial 169 return jacobi_poly(n, a, b, x) 170 171 def fdiff(self, argindex=4): 172 from sympy import Sum 173 if argindex == 1: 174 # Diff wrt n 175 raise ArgumentIndexError(self, argindex) 176 elif argindex == 2: 177 # Diff wrt a 178 n, a, b, x = self.args 179 k = Dummy("k") 180 f1 = 1 / (a + b + n + k + 1) 181 f2 = ((a + b + 2*k + 1) * RisingFactorial(b + k + 1, n - k) / 182 ((n - k) * RisingFactorial(a + b + k + 1, n - k))) 183 return Sum(f1 * (jacobi(n, a, b, x) + f2*jacobi(k, a, b, x)), (k, 0, n - 1)) 184 elif argindex == 3: 185 # Diff wrt b 186 n, a, b, x = self.args 187 k = Dummy("k") 188 f1 = 1 / (a + b + n + k + 1) 189 f2 = (-1)**(n - k) * ((a + b + 2*k + 1) * RisingFactorial(a + k + 1, n - k) / 190 ((n - k) * RisingFactorial(a + b + k + 1, n - k))) 191 return Sum(f1 * (jacobi(n, a, b, x) + f2*jacobi(k, a, b, x)), (k, 0, n - 1)) 192 elif argindex == 4: 193 # Diff wrt x 194 n, a, b, x = self.args 195 return S.Half * (a + b + n + 1) * jacobi(n - 1, a + 1, b + 1, x) 196 else: 197 raise ArgumentIndexError(self, argindex) 198 199 def _eval_rewrite_as_polynomial(self, n, a, b, x, **kwargs): 200 from sympy import Sum 201 # Make sure n \in N 202 if n.is_negative or n.is_integer is False: 203 raise ValueError("Error: n should be a non-negative integer.") 204 k = Dummy("k") 205 kern = (RisingFactorial(-n, k) * RisingFactorial(a + b + n + 1, k) * RisingFactorial(a + k + 1, n - k) / 206 factorial(k) * ((1 - x)/2)**k) 207 return 1 / factorial(n) * Sum(kern, (k, 0, n)) 208 209 def _eval_conjugate(self): 210 n, a, b, x = self.args 211 return self.func(n, a.conjugate(), b.conjugate(), x.conjugate()) 212 213 214def jacobi_normalized(n, a, b, x): 215 r""" 216 Jacobi polynomial $P_n^{\left(\alpha, \beta\right)}(x)$. 217 218 Explanation 219 =========== 220 221 ``jacobi_normalized(n, alpha, beta, x)`` gives the nth 222 Jacobi polynomial in *x*, $P_n^{\left(\alpha, \beta\right)}(x)$. 223 224 The Jacobi polynomials are orthogonal on $[-1, 1]$ with respect 225 to the weight $\left(1-x\right)^\alpha \left(1+x\right)^\beta$. 226 227 This functions returns the polynomials normilzed: 228 229 .. math:: 230 231 \int_{-1}^{1} 232 P_m^{\left(\alpha, \beta\right)}(x) 233 P_n^{\left(\alpha, \beta\right)}(x) 234 (1-x)^{\alpha} (1+x)^{\beta} \mathrm{d}x 235 = \delta_{m,n} 236 237 Examples 238 ======== 239 240 >>> from sympy import jacobi_normalized 241 >>> from sympy.abc import n,a,b,x 242 243 >>> jacobi_normalized(n, a, b, x) 244 jacobi(n, a, b, x)/sqrt(2**(a + b + 1)*gamma(a + n + 1)*gamma(b + n + 1)/((a + b + 2*n + 1)*factorial(n)*gamma(a + b + n + 1))) 245 246 Parameters 247 ========== 248 249 n : integer degree of polynomial 250 251 a : alpha value 252 253 b : beta value 254 255 x : symbol 256 257 See Also 258 ======== 259 260 gegenbauer, 261 chebyshevt_root, chebyshevu, chebyshevu_root, 262 legendre, assoc_legendre, 263 hermite, 264 laguerre, assoc_laguerre, 265 sympy.polys.orthopolys.jacobi_poly, 266 sympy.polys.orthopolys.gegenbauer_poly 267 sympy.polys.orthopolys.chebyshevt_poly 268 sympy.polys.orthopolys.chebyshevu_poly 269 sympy.polys.orthopolys.hermite_poly 270 sympy.polys.orthopolys.legendre_poly 271 sympy.polys.orthopolys.laguerre_poly 272 273 References 274 ========== 275 276 .. [1] https://en.wikipedia.org/wiki/Jacobi_polynomials 277 .. [2] http://mathworld.wolfram.com/JacobiPolynomial.html 278 .. [3] http://functions.wolfram.com/Polynomials/JacobiP/ 279 280 """ 281 nfactor = (S(2)**(a + b + 1) * (gamma(n + a + 1) * gamma(n + b + 1)) 282 / (2*n + a + b + 1) / (factorial(n) * gamma(n + a + b + 1))) 283 284 return jacobi(n, a, b, x) / sqrt(nfactor) 285 286 287#---------------------------------------------------------------------------- 288# Gegenbauer polynomials 289# 290 291 292class gegenbauer(OrthogonalPolynomial): 293 r""" 294 Gegenbauer polynomial $C_n^{\left(\alpha\right)}(x)$. 295 296 Explanation 297 =========== 298 299 ``gegenbauer(n, alpha, x)`` gives the nth Gegenbauer polynomial 300 in x, $C_n^{\left(\alpha\right)}(x)$. 301 302 The Gegenbauer polynomials are orthogonal on $[-1, 1]$ with 303 respect to the weight $\left(1-x^2\right)^{\alpha-\frac{1}{2}}$. 304 305 Examples 306 ======== 307 308 >>> from sympy import gegenbauer, conjugate, diff 309 >>> from sympy.abc import n,a,x 310 >>> gegenbauer(0, a, x) 311 1 312 >>> gegenbauer(1, a, x) 313 2*a*x 314 >>> gegenbauer(2, a, x) 315 -a + x**2*(2*a**2 + 2*a) 316 >>> gegenbauer(3, a, x) 317 x**3*(4*a**3/3 + 4*a**2 + 8*a/3) + x*(-2*a**2 - 2*a) 318 319 >>> gegenbauer(n, a, x) 320 gegenbauer(n, a, x) 321 >>> gegenbauer(n, a, -x) 322 (-1)**n*gegenbauer(n, a, x) 323 324 >>> gegenbauer(n, a, 0) 325 2**n*sqrt(pi)*gamma(a + n/2)/(gamma(a)*gamma(1/2 - n/2)*gamma(n + 1)) 326 >>> gegenbauer(n, a, 1) 327 gamma(2*a + n)/(gamma(2*a)*gamma(n + 1)) 328 329 >>> conjugate(gegenbauer(n, a, x)) 330 gegenbauer(n, conjugate(a), conjugate(x)) 331 332 >>> diff(gegenbauer(n, a, x), x) 333 2*a*gegenbauer(n - 1, a + 1, x) 334 335 See Also 336 ======== 337 338 jacobi, 339 chebyshevt_root, chebyshevu, chebyshevu_root, 340 legendre, assoc_legendre, 341 hermite, 342 laguerre, assoc_laguerre, 343 sympy.polys.orthopolys.jacobi_poly 344 sympy.polys.orthopolys.gegenbauer_poly 345 sympy.polys.orthopolys.chebyshevt_poly 346 sympy.polys.orthopolys.chebyshevu_poly 347 sympy.polys.orthopolys.hermite_poly 348 sympy.polys.orthopolys.legendre_poly 349 sympy.polys.orthopolys.laguerre_poly 350 351 References 352 ========== 353 354 .. [1] https://en.wikipedia.org/wiki/Gegenbauer_polynomials 355 .. [2] http://mathworld.wolfram.com/GegenbauerPolynomial.html 356 .. [3] http://functions.wolfram.com/Polynomials/GegenbauerC3/ 357 358 """ 359 360 @classmethod 361 def eval(cls, n, a, x): 362 # For negative n the polynomials vanish 363 # See http://functions.wolfram.com/Polynomials/GegenbauerC3/03/01/03/0012/ 364 if n.is_negative: 365 return S.Zero 366 367 # Some special values for fixed a 368 if a == S.Half: 369 return legendre(n, x) 370 elif a == S.One: 371 return chebyshevu(n, x) 372 elif a == S.NegativeOne: 373 return S.Zero 374 375 if not n.is_Number: 376 # Handle this before the general sign extraction rule 377 if x == S.NegativeOne: 378 if (re(a) > S.Half) == True: 379 return S.ComplexInfinity 380 else: 381 return (cos(S.Pi*(a+n)) * sec(S.Pi*a) * gamma(2*a+n) / 382 (gamma(2*a) * gamma(n+1))) 383 384 # Symbolic result C^a_n(x) 385 # C^a_n(-x) ---> (-1)**n * C^a_n(x) 386 if x.could_extract_minus_sign(): 387 return S.NegativeOne**n * gegenbauer(n, a, -x) 388 # We can evaluate for some special values of x 389 if x.is_zero: 390 return (2**n * sqrt(S.Pi) * gamma(a + S.Half*n) / 391 (gamma((1 - n)/2) * gamma(n + 1) * gamma(a)) ) 392 if x == S.One: 393 return gamma(2*a + n) / (gamma(2*a) * gamma(n + 1)) 394 elif x is S.Infinity: 395 if n.is_positive: 396 return RisingFactorial(a, n) * S.Infinity 397 else: 398 # n is a given fixed integer, evaluate into polynomial 399 return gegenbauer_poly(n, a, x) 400 401 def fdiff(self, argindex=3): 402 from sympy import Sum 403 if argindex == 1: 404 # Diff wrt n 405 raise ArgumentIndexError(self, argindex) 406 elif argindex == 2: 407 # Diff wrt a 408 n, a, x = self.args 409 k = Dummy("k") 410 factor1 = 2 * (1 + (-1)**(n - k)) * (k + a) / ((k + 411 n + 2*a) * (n - k)) 412 factor2 = 2*(k + 1) / ((k + 2*a) * (2*k + 2*a + 1)) + \ 413 2 / (k + n + 2*a) 414 kern = factor1*gegenbauer(k, a, x) + factor2*gegenbauer(n, a, x) 415 return Sum(kern, (k, 0, n - 1)) 416 elif argindex == 3: 417 # Diff wrt x 418 n, a, x = self.args 419 return 2*a*gegenbauer(n - 1, a + 1, x) 420 else: 421 raise ArgumentIndexError(self, argindex) 422 423 def _eval_rewrite_as_polynomial(self, n, a, x, **kwargs): 424 from sympy import Sum 425 k = Dummy("k") 426 kern = ((-1)**k * RisingFactorial(a, n - k) * (2*x)**(n - 2*k) / 427 (factorial(k) * factorial(n - 2*k))) 428 return Sum(kern, (k, 0, floor(n/2))) 429 430 def _eval_conjugate(self): 431 n, a, x = self.args 432 return self.func(n, a.conjugate(), x.conjugate()) 433 434#---------------------------------------------------------------------------- 435# Chebyshev polynomials of first and second kind 436# 437 438 439class chebyshevt(OrthogonalPolynomial): 440 r""" 441 Chebyshev polynomial of the first kind, $T_n(x)$. 442 443 Explanation 444 =========== 445 446 ``chebyshevt(n, x)`` gives the nth Chebyshev polynomial (of the first 447 kind) in x, $T_n(x)$. 448 449 The Chebyshev polynomials of the first kind are orthogonal on 450 $[-1, 1]$ with respect to the weight $\frac{1}{\sqrt{1-x^2}}$. 451 452 Examples 453 ======== 454 455 >>> from sympy import chebyshevt, diff 456 >>> from sympy.abc import n,x 457 >>> chebyshevt(0, x) 458 1 459 >>> chebyshevt(1, x) 460 x 461 >>> chebyshevt(2, x) 462 2*x**2 - 1 463 464 >>> chebyshevt(n, x) 465 chebyshevt(n, x) 466 >>> chebyshevt(n, -x) 467 (-1)**n*chebyshevt(n, x) 468 >>> chebyshevt(-n, x) 469 chebyshevt(n, x) 470 471 >>> chebyshevt(n, 0) 472 cos(pi*n/2) 473 >>> chebyshevt(n, -1) 474 (-1)**n 475 476 >>> diff(chebyshevt(n, x), x) 477 n*chebyshevu(n - 1, x) 478 479 See Also 480 ======== 481 482 jacobi, gegenbauer, 483 chebyshevt_root, chebyshevu, chebyshevu_root, 484 legendre, assoc_legendre, 485 hermite, 486 laguerre, assoc_laguerre, 487 sympy.polys.orthopolys.jacobi_poly 488 sympy.polys.orthopolys.gegenbauer_poly 489 sympy.polys.orthopolys.chebyshevt_poly 490 sympy.polys.orthopolys.chebyshevu_poly 491 sympy.polys.orthopolys.hermite_poly 492 sympy.polys.orthopolys.legendre_poly 493 sympy.polys.orthopolys.laguerre_poly 494 495 References 496 ========== 497 498 .. [1] https://en.wikipedia.org/wiki/Chebyshev_polynomial 499 .. [2] http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html 500 .. [3] http://mathworld.wolfram.com/ChebyshevPolynomialoftheSecondKind.html 501 .. [4] http://functions.wolfram.com/Polynomials/ChebyshevT/ 502 .. [5] http://functions.wolfram.com/Polynomials/ChebyshevU/ 503 504 """ 505 506 _ortho_poly = staticmethod(chebyshevt_poly) 507 508 @classmethod 509 def eval(cls, n, x): 510 if not n.is_Number: 511 # Symbolic result T_n(x) 512 # T_n(-x) ---> (-1)**n * T_n(x) 513 if x.could_extract_minus_sign(): 514 return S.NegativeOne**n * chebyshevt(n, -x) 515 # T_{-n}(x) ---> T_n(x) 516 if n.could_extract_minus_sign(): 517 return chebyshevt(-n, x) 518 # We can evaluate for some special values of x 519 if x.is_zero: 520 return cos(S.Half * S.Pi * n) 521 if x == S.One: 522 return S.One 523 elif x is S.Infinity: 524 return S.Infinity 525 else: 526 # n is a given fixed integer, evaluate into polynomial 527 if n.is_negative: 528 # T_{-n}(x) == T_n(x) 529 return cls._eval_at_order(-n, x) 530 else: 531 return cls._eval_at_order(n, x) 532 533 def fdiff(self, argindex=2): 534 if argindex == 1: 535 # Diff wrt n 536 raise ArgumentIndexError(self, argindex) 537 elif argindex == 2: 538 # Diff wrt x 539 n, x = self.args 540 return n * chebyshevu(n - 1, x) 541 else: 542 raise ArgumentIndexError(self, argindex) 543 544 def _eval_rewrite_as_polynomial(self, n, x, **kwargs): 545 from sympy import Sum 546 k = Dummy("k") 547 kern = binomial(n, 2*k) * (x**2 - 1)**k * x**(n - 2*k) 548 return Sum(kern, (k, 0, floor(n/2))) 549 550 551class chebyshevu(OrthogonalPolynomial): 552 r""" 553 Chebyshev polynomial of the second kind, $U_n(x)$. 554 555 Explanation 556 =========== 557 558 ``chebyshevu(n, x)`` gives the nth Chebyshev polynomial of the second 559 kind in x, $U_n(x)$. 560 561 The Chebyshev polynomials of the second kind are orthogonal on 562 $[-1, 1]$ with respect to the weight $\sqrt{1-x^2}$. 563 564 Examples 565 ======== 566 567 >>> from sympy import chebyshevu, diff 568 >>> from sympy.abc import n,x 569 >>> chebyshevu(0, x) 570 1 571 >>> chebyshevu(1, x) 572 2*x 573 >>> chebyshevu(2, x) 574 4*x**2 - 1 575 576 >>> chebyshevu(n, x) 577 chebyshevu(n, x) 578 >>> chebyshevu(n, -x) 579 (-1)**n*chebyshevu(n, x) 580 >>> chebyshevu(-n, x) 581 -chebyshevu(n - 2, x) 582 583 >>> chebyshevu(n, 0) 584 cos(pi*n/2) 585 >>> chebyshevu(n, 1) 586 n + 1 587 588 >>> diff(chebyshevu(n, x), x) 589 (-x*chebyshevu(n, x) + (n + 1)*chebyshevt(n + 1, x))/(x**2 - 1) 590 591 See Also 592 ======== 593 594 jacobi, gegenbauer, 595 chebyshevt, chebyshevt_root, chebyshevu_root, 596 legendre, assoc_legendre, 597 hermite, 598 laguerre, assoc_laguerre, 599 sympy.polys.orthopolys.jacobi_poly 600 sympy.polys.orthopolys.gegenbauer_poly 601 sympy.polys.orthopolys.chebyshevt_poly 602 sympy.polys.orthopolys.chebyshevu_poly 603 sympy.polys.orthopolys.hermite_poly 604 sympy.polys.orthopolys.legendre_poly 605 sympy.polys.orthopolys.laguerre_poly 606 607 References 608 ========== 609 610 .. [1] https://en.wikipedia.org/wiki/Chebyshev_polynomial 611 .. [2] http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html 612 .. [3] http://mathworld.wolfram.com/ChebyshevPolynomialoftheSecondKind.html 613 .. [4] http://functions.wolfram.com/Polynomials/ChebyshevT/ 614 .. [5] http://functions.wolfram.com/Polynomials/ChebyshevU/ 615 616 """ 617 618 _ortho_poly = staticmethod(chebyshevu_poly) 619 620 @classmethod 621 def eval(cls, n, x): 622 if not n.is_Number: 623 # Symbolic result U_n(x) 624 # U_n(-x) ---> (-1)**n * U_n(x) 625 if x.could_extract_minus_sign(): 626 return S.NegativeOne**n * chebyshevu(n, -x) 627 # U_{-n}(x) ---> -U_{n-2}(x) 628 if n.could_extract_minus_sign(): 629 if n == S.NegativeOne: 630 # n can not be -1 here 631 return S.Zero 632 elif not (-n - 2).could_extract_minus_sign(): 633 return -chebyshevu(-n - 2, x) 634 # We can evaluate for some special values of x 635 if x.is_zero: 636 return cos(S.Half * S.Pi * n) 637 if x == S.One: 638 return S.One + n 639 elif x is S.Infinity: 640 return S.Infinity 641 else: 642 # n is a given fixed integer, evaluate into polynomial 643 if n.is_negative: 644 # U_{-n}(x) ---> -U_{n-2}(x) 645 if n == S.NegativeOne: 646 return S.Zero 647 else: 648 return -cls._eval_at_order(-n - 2, x) 649 else: 650 return cls._eval_at_order(n, x) 651 652 def fdiff(self, argindex=2): 653 if argindex == 1: 654 # Diff wrt n 655 raise ArgumentIndexError(self, argindex) 656 elif argindex == 2: 657 # Diff wrt x 658 n, x = self.args 659 return ((n + 1) * chebyshevt(n + 1, x) - x * chebyshevu(n, x)) / (x**2 - 1) 660 else: 661 raise ArgumentIndexError(self, argindex) 662 663 def _eval_rewrite_as_polynomial(self, n, x, **kwargs): 664 from sympy import Sum 665 k = Dummy("k") 666 kern = S.NegativeOne**k * factorial( 667 n - k) * (2*x)**(n - 2*k) / (factorial(k) * factorial(n - 2*k)) 668 return Sum(kern, (k, 0, floor(n/2))) 669 670 671class chebyshevt_root(Function): 672 r""" 673 ``chebyshev_root(n, k)`` returns the kth root (indexed from zero) of 674 the nth Chebyshev polynomial of the first kind; that is, if 675 0 <= k < n, ``chebyshevt(n, chebyshevt_root(n, k)) == 0``. 676 677 Examples 678 ======== 679 680 >>> from sympy import chebyshevt, chebyshevt_root 681 >>> chebyshevt_root(3, 2) 682 -sqrt(3)/2 683 >>> chebyshevt(3, chebyshevt_root(3, 2)) 684 0 685 686 See Also 687 ======== 688 689 jacobi, gegenbauer, 690 chebyshevt, chebyshevu, chebyshevu_root, 691 legendre, assoc_legendre, 692 hermite, 693 laguerre, assoc_laguerre, 694 sympy.polys.orthopolys.jacobi_poly 695 sympy.polys.orthopolys.gegenbauer_poly 696 sympy.polys.orthopolys.chebyshevt_poly 697 sympy.polys.orthopolys.chebyshevu_poly 698 sympy.polys.orthopolys.hermite_poly 699 sympy.polys.orthopolys.legendre_poly 700 sympy.polys.orthopolys.laguerre_poly 701 """ 702 703 @classmethod 704 def eval(cls, n, k): 705 if not ((0 <= k) and (k < n)): 706 raise ValueError("must have 0 <= k < n, " 707 "got k = %s and n = %s" % (k, n)) 708 return cos(S.Pi*(2*k + 1)/(2*n)) 709 710 711class chebyshevu_root(Function): 712 r""" 713 ``chebyshevu_root(n, k)`` returns the kth root (indexed from zero) of the 714 nth Chebyshev polynomial of the second kind; that is, if 0 <= k < n, 715 ``chebyshevu(n, chebyshevu_root(n, k)) == 0``. 716 717 Examples 718 ======== 719 720 >>> from sympy import chebyshevu, chebyshevu_root 721 >>> chebyshevu_root(3, 2) 722 -sqrt(2)/2 723 >>> chebyshevu(3, chebyshevu_root(3, 2)) 724 0 725 726 See Also 727 ======== 728 729 chebyshevt, chebyshevt_root, chebyshevu, 730 legendre, assoc_legendre, 731 hermite, 732 laguerre, assoc_laguerre, 733 sympy.polys.orthopolys.jacobi_poly 734 sympy.polys.orthopolys.gegenbauer_poly 735 sympy.polys.orthopolys.chebyshevt_poly 736 sympy.polys.orthopolys.chebyshevu_poly 737 sympy.polys.orthopolys.hermite_poly 738 sympy.polys.orthopolys.legendre_poly 739 sympy.polys.orthopolys.laguerre_poly 740 """ 741 742 743 @classmethod 744 def eval(cls, n, k): 745 if not ((0 <= k) and (k < n)): 746 raise ValueError("must have 0 <= k < n, " 747 "got k = %s and n = %s" % (k, n)) 748 return cos(S.Pi*(k + 1)/(n + 1)) 749 750#---------------------------------------------------------------------------- 751# Legendre polynomials and Associated Legendre polynomials 752# 753 754 755class legendre(OrthogonalPolynomial): 756 r""" 757 ``legendre(n, x)`` gives the nth Legendre polynomial of x, $P_n(x)$ 758 759 Explanation 760 =========== 761 762 The Legendre polynomials are orthogonal on [-1, 1] with respect to 763 the constant weight 1. They satisfy $P_n(1) = 1$ for all n; further, 764 $P_n$ is odd for odd n and even for even n. 765 766 Examples 767 ======== 768 769 >>> from sympy import legendre, diff 770 >>> from sympy.abc import x, n 771 >>> legendre(0, x) 772 1 773 >>> legendre(1, x) 774 x 775 >>> legendre(2, x) 776 3*x**2/2 - 1/2 777 >>> legendre(n, x) 778 legendre(n, x) 779 >>> diff(legendre(n,x), x) 780 n*(x*legendre(n, x) - legendre(n - 1, x))/(x**2 - 1) 781 782 See Also 783 ======== 784 785 jacobi, gegenbauer, 786 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, 787 assoc_legendre, 788 hermite, 789 laguerre, assoc_laguerre, 790 sympy.polys.orthopolys.jacobi_poly 791 sympy.polys.orthopolys.gegenbauer_poly 792 sympy.polys.orthopolys.chebyshevt_poly 793 sympy.polys.orthopolys.chebyshevu_poly 794 sympy.polys.orthopolys.hermite_poly 795 sympy.polys.orthopolys.legendre_poly 796 sympy.polys.orthopolys.laguerre_poly 797 798 References 799 ========== 800 801 .. [1] https://en.wikipedia.org/wiki/Legendre_polynomial 802 .. [2] http://mathworld.wolfram.com/LegendrePolynomial.html 803 .. [3] http://functions.wolfram.com/Polynomials/LegendreP/ 804 .. [4] http://functions.wolfram.com/Polynomials/LegendreP2/ 805 806 """ 807 808 _ortho_poly = staticmethod(legendre_poly) 809 810 @classmethod 811 def eval(cls, n, x): 812 if not n.is_Number: 813 # Symbolic result L_n(x) 814 # L_n(-x) ---> (-1)**n * L_n(x) 815 if x.could_extract_minus_sign(): 816 return S.NegativeOne**n * legendre(n, -x) 817 # L_{-n}(x) ---> L_{n-1}(x) 818 if n.could_extract_minus_sign() and not(-n - 1).could_extract_minus_sign(): 819 return legendre(-n - S.One, x) 820 # We can evaluate for some special values of x 821 if x.is_zero: 822 return sqrt(S.Pi)/(gamma(S.Half - n/2)*gamma(S.One + n/2)) 823 elif x == S.One: 824 return S.One 825 elif x is S.Infinity: 826 return S.Infinity 827 else: 828 # n is a given fixed integer, evaluate into polynomial; 829 # L_{-n}(x) ---> L_{n-1}(x) 830 if n.is_negative: 831 n = -n - S.One 832 return cls._eval_at_order(n, x) 833 834 def fdiff(self, argindex=2): 835 if argindex == 1: 836 # Diff wrt n 837 raise ArgumentIndexError(self, argindex) 838 elif argindex == 2: 839 # Diff wrt x 840 # Find better formula, this is unsuitable for x = +/-1 841 # http://www.autodiff.org/ad16/Oral/Buecker_Legendre.pdf says 842 # at x = 1: 843 # n*(n + 1)/2 , m = 0 844 # oo , m = 1 845 # -(n-1)*n*(n+1)*(n+2)/4 , m = 2 846 # 0 , m = 3, 4, ..., n 847 # 848 # at x = -1 849 # (-1)**(n+1)*n*(n + 1)/2 , m = 0 850 # (-1)**n*oo , m = 1 851 # (-1)**n*(n-1)*n*(n+1)*(n+2)/4 , m = 2 852 # 0 , m = 3, 4, ..., n 853 n, x = self.args 854 return n/(x**2 - 1)*(x*legendre(n, x) - legendre(n - 1, x)) 855 else: 856 raise ArgumentIndexError(self, argindex) 857 858 def _eval_rewrite_as_polynomial(self, n, x, **kwargs): 859 from sympy import Sum 860 k = Dummy("k") 861 kern = (-1)**k*binomial(n, k)**2*((1 + x)/2)**(n - k)*((1 - x)/2)**k 862 return Sum(kern, (k, 0, n)) 863 864 865class assoc_legendre(Function): 866 r""" 867 ``assoc_legendre(n, m, x)`` gives $P_n^m(x)$, where n and m are 868 the degree and order or an expression which is related to the nth 869 order Legendre polynomial, $P_n(x)$ in the following manner: 870 871 .. math:: 872 P_n^m(x) = (-1)^m (1 - x^2)^{\frac{m}{2}} 873 \frac{\mathrm{d}^m P_n(x)}{\mathrm{d} x^m} 874 875 Explanation 876 =========== 877 878 Associated Legendre polynomials are orthogonal on [-1, 1] with: 879 880 - weight = 1 for the same m, and different n. 881 - weight = 1/(1-x**2) for the same n, and different m. 882 883 Examples 884 ======== 885 886 >>> from sympy import assoc_legendre 887 >>> from sympy.abc import x, m, n 888 >>> assoc_legendre(0,0, x) 889 1 890 >>> assoc_legendre(1,0, x) 891 x 892 >>> assoc_legendre(1,1, x) 893 -sqrt(1 - x**2) 894 >>> assoc_legendre(n,m,x) 895 assoc_legendre(n, m, x) 896 897 See Also 898 ======== 899 900 jacobi, gegenbauer, 901 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, 902 legendre, 903 hermite, 904 laguerre, assoc_laguerre, 905 sympy.polys.orthopolys.jacobi_poly 906 sympy.polys.orthopolys.gegenbauer_poly 907 sympy.polys.orthopolys.chebyshevt_poly 908 sympy.polys.orthopolys.chebyshevu_poly 909 sympy.polys.orthopolys.hermite_poly 910 sympy.polys.orthopolys.legendre_poly 911 sympy.polys.orthopolys.laguerre_poly 912 913 References 914 ========== 915 916 .. [1] https://en.wikipedia.org/wiki/Associated_Legendre_polynomials 917 .. [2] http://mathworld.wolfram.com/LegendrePolynomial.html 918 .. [3] http://functions.wolfram.com/Polynomials/LegendreP/ 919 .. [4] http://functions.wolfram.com/Polynomials/LegendreP2/ 920 921 """ 922 923 @classmethod 924 def _eval_at_order(cls, n, m): 925 P = legendre_poly(n, _x, polys=True).diff((_x, m)) 926 return (-1)**m * (1 - _x**2)**Rational(m, 2) * P.as_expr() 927 928 @classmethod 929 def eval(cls, n, m, x): 930 if m.could_extract_minus_sign(): 931 # P^{-m}_n ---> F * P^m_n 932 return S.NegativeOne**(-m) * (factorial(m + n)/factorial(n - m)) * assoc_legendre(n, -m, x) 933 if m == 0: 934 # P^0_n ---> L_n 935 return legendre(n, x) 936 if x == 0: 937 return 2**m*sqrt(S.Pi) / (gamma((1 - m - n)/2)*gamma(1 - (m - n)/2)) 938 if n.is_Number and m.is_Number and n.is_integer and m.is_integer: 939 if n.is_negative: 940 raise ValueError("%s : 1st index must be nonnegative integer (got %r)" % (cls, n)) 941 if abs(m) > n: 942 raise ValueError("%s : abs('2nd index') must be <= '1st index' (got %r, %r)" % (cls, n, m)) 943 return cls._eval_at_order(int(n), abs(int(m))).subs(_x, x) 944 945 def fdiff(self, argindex=3): 946 if argindex == 1: 947 # Diff wrt n 948 raise ArgumentIndexError(self, argindex) 949 elif argindex == 2: 950 # Diff wrt m 951 raise ArgumentIndexError(self, argindex) 952 elif argindex == 3: 953 # Diff wrt x 954 # Find better formula, this is unsuitable for x = 1 955 n, m, x = self.args 956 return 1/(x**2 - 1)*(x*n*assoc_legendre(n, m, x) - (m + n)*assoc_legendre(n - 1, m, x)) 957 else: 958 raise ArgumentIndexError(self, argindex) 959 960 def _eval_rewrite_as_polynomial(self, n, m, x, **kwargs): 961 from sympy import Sum 962 k = Dummy("k") 963 kern = factorial(2*n - 2*k)/(2**n*factorial(n - k)*factorial( 964 k)*factorial(n - 2*k - m))*(-1)**k*x**(n - m - 2*k) 965 return (1 - x**2)**(m/2) * Sum(kern, (k, 0, floor((n - m)*S.Half))) 966 967 def _eval_conjugate(self): 968 n, m, x = self.args 969 return self.func(n, m.conjugate(), x.conjugate()) 970 971#---------------------------------------------------------------------------- 972# Hermite polynomials 973# 974 975 976class hermite(OrthogonalPolynomial): 977 r""" 978 ``hermite(n, x)`` gives the nth Hermite polynomial in x, $H_n(x)$ 979 980 Explanation 981 =========== 982 983 The Hermite polynomials are orthogonal on $(-\infty, \infty)$ 984 with respect to the weight $\exp\left(-x^2\right)$. 985 986 Examples 987 ======== 988 989 >>> from sympy import hermite, diff 990 >>> from sympy.abc import x, n 991 >>> hermite(0, x) 992 1 993 >>> hermite(1, x) 994 2*x 995 >>> hermite(2, x) 996 4*x**2 - 2 997 >>> hermite(n, x) 998 hermite(n, x) 999 >>> diff(hermite(n,x), x) 1000 2*n*hermite(n - 1, x) 1001 >>> hermite(n, -x) 1002 (-1)**n*hermite(n, x) 1003 1004 See Also 1005 ======== 1006 1007 jacobi, gegenbauer, 1008 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, 1009 legendre, assoc_legendre, 1010 laguerre, assoc_laguerre, 1011 sympy.polys.orthopolys.jacobi_poly 1012 sympy.polys.orthopolys.gegenbauer_poly 1013 sympy.polys.orthopolys.chebyshevt_poly 1014 sympy.polys.orthopolys.chebyshevu_poly 1015 sympy.polys.orthopolys.hermite_poly 1016 sympy.polys.orthopolys.legendre_poly 1017 sympy.polys.orthopolys.laguerre_poly 1018 1019 References 1020 ========== 1021 1022 .. [1] https://en.wikipedia.org/wiki/Hermite_polynomial 1023 .. [2] http://mathworld.wolfram.com/HermitePolynomial.html 1024 .. [3] http://functions.wolfram.com/Polynomials/HermiteH/ 1025 1026 """ 1027 1028 _ortho_poly = staticmethod(hermite_poly) 1029 1030 @classmethod 1031 def eval(cls, n, x): 1032 if not n.is_Number: 1033 # Symbolic result H_n(x) 1034 # H_n(-x) ---> (-1)**n * H_n(x) 1035 if x.could_extract_minus_sign(): 1036 return S.NegativeOne**n * hermite(n, -x) 1037 # We can evaluate for some special values of x 1038 if x.is_zero: 1039 return 2**n * sqrt(S.Pi) / gamma((S.One - n)/2) 1040 elif x is S.Infinity: 1041 return S.Infinity 1042 else: 1043 # n is a given fixed integer, evaluate into polynomial 1044 if n.is_negative: 1045 raise ValueError( 1046 "The index n must be nonnegative integer (got %r)" % n) 1047 else: 1048 return cls._eval_at_order(n, x) 1049 1050 def fdiff(self, argindex=2): 1051 if argindex == 1: 1052 # Diff wrt n 1053 raise ArgumentIndexError(self, argindex) 1054 elif argindex == 2: 1055 # Diff wrt x 1056 n, x = self.args 1057 return 2*n*hermite(n - 1, x) 1058 else: 1059 raise ArgumentIndexError(self, argindex) 1060 1061 def _eval_rewrite_as_polynomial(self, n, x, **kwargs): 1062 from sympy import Sum 1063 k = Dummy("k") 1064 kern = (-1)**k / (factorial(k)*factorial(n - 2*k)) * (2*x)**(n - 2*k) 1065 return factorial(n)*Sum(kern, (k, 0, floor(n/2))) 1066 1067#---------------------------------------------------------------------------- 1068# Laguerre polynomials 1069# 1070 1071 1072class laguerre(OrthogonalPolynomial): 1073 r""" 1074 Returns the nth Laguerre polynomial in x, $L_n(x)$. 1075 1076 Examples 1077 ======== 1078 1079 >>> from sympy import laguerre, diff 1080 >>> from sympy.abc import x, n 1081 >>> laguerre(0, x) 1082 1 1083 >>> laguerre(1, x) 1084 1 - x 1085 >>> laguerre(2, x) 1086 x**2/2 - 2*x + 1 1087 >>> laguerre(3, x) 1088 -x**3/6 + 3*x**2/2 - 3*x + 1 1089 1090 >>> laguerre(n, x) 1091 laguerre(n, x) 1092 1093 >>> diff(laguerre(n, x), x) 1094 -assoc_laguerre(n - 1, 1, x) 1095 1096 Parameters 1097 ========== 1098 1099 n : int 1100 Degree of Laguerre polynomial. Must be ``n >= 0``. 1101 1102 See Also 1103 ======== 1104 1105 jacobi, gegenbauer, 1106 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, 1107 legendre, assoc_legendre, 1108 hermite, 1109 assoc_laguerre, 1110 sympy.polys.orthopolys.jacobi_poly 1111 sympy.polys.orthopolys.gegenbauer_poly 1112 sympy.polys.orthopolys.chebyshevt_poly 1113 sympy.polys.orthopolys.chebyshevu_poly 1114 sympy.polys.orthopolys.hermite_poly 1115 sympy.polys.orthopolys.legendre_poly 1116 sympy.polys.orthopolys.laguerre_poly 1117 1118 References 1119 ========== 1120 1121 .. [1] https://en.wikipedia.org/wiki/Laguerre_polynomial 1122 .. [2] http://mathworld.wolfram.com/LaguerrePolynomial.html 1123 .. [3] http://functions.wolfram.com/Polynomials/LaguerreL/ 1124 .. [4] http://functions.wolfram.com/Polynomials/LaguerreL3/ 1125 1126 """ 1127 1128 _ortho_poly = staticmethod(laguerre_poly) 1129 1130 @classmethod 1131 def eval(cls, n, x): 1132 if n.is_integer is False: 1133 raise ValueError("Error: n should be an integer.") 1134 if not n.is_Number: 1135 # Symbolic result L_n(x) 1136 # L_{n}(-x) ---> exp(-x) * L_{-n-1}(x) 1137 # L_{-n}(x) ---> exp(x) * L_{n-1}(-x) 1138 if n.could_extract_minus_sign() and not(-n - 1).could_extract_minus_sign(): 1139 return exp(x)*laguerre(-n - 1, -x) 1140 # We can evaluate for some special values of x 1141 if x.is_zero: 1142 return S.One 1143 elif x is S.NegativeInfinity: 1144 return S.Infinity 1145 elif x is S.Infinity: 1146 return S.NegativeOne**n * S.Infinity 1147 else: 1148 if n.is_negative: 1149 return exp(x)*laguerre(-n - 1, -x) 1150 else: 1151 return cls._eval_at_order(n, x) 1152 1153 def fdiff(self, argindex=2): 1154 if argindex == 1: 1155 # Diff wrt n 1156 raise ArgumentIndexError(self, argindex) 1157 elif argindex == 2: 1158 # Diff wrt x 1159 n, x = self.args 1160 return -assoc_laguerre(n - 1, 1, x) 1161 else: 1162 raise ArgumentIndexError(self, argindex) 1163 1164 def _eval_rewrite_as_polynomial(self, n, x, **kwargs): 1165 from sympy import Sum 1166 # Make sure n \in N_0 1167 if n.is_negative: 1168 return exp(x) * self._eval_rewrite_as_polynomial(-n - 1, -x, **kwargs) 1169 if n.is_integer is False: 1170 raise ValueError("Error: n should be an integer.") 1171 k = Dummy("k") 1172 kern = RisingFactorial(-n, k) / factorial(k)**2 * x**k 1173 return Sum(kern, (k, 0, n)) 1174 1175 1176class assoc_laguerre(OrthogonalPolynomial): 1177 r""" 1178 Returns the nth generalized Laguerre polynomial in x, $L_n(x)$. 1179 1180 Examples 1181 ======== 1182 1183 >>> from sympy import assoc_laguerre, diff 1184 >>> from sympy.abc import x, n, a 1185 >>> assoc_laguerre(0, a, x) 1186 1 1187 >>> assoc_laguerre(1, a, x) 1188 a - x + 1 1189 >>> assoc_laguerre(2, a, x) 1190 a**2/2 + 3*a/2 + x**2/2 + x*(-a - 2) + 1 1191 >>> assoc_laguerre(3, a, x) 1192 a**3/6 + a**2 + 11*a/6 - x**3/6 + x**2*(a/2 + 3/2) + 1193 x*(-a**2/2 - 5*a/2 - 3) + 1 1194 1195 >>> assoc_laguerre(n, a, 0) 1196 binomial(a + n, a) 1197 1198 >>> assoc_laguerre(n, a, x) 1199 assoc_laguerre(n, a, x) 1200 1201 >>> assoc_laguerre(n, 0, x) 1202 laguerre(n, x) 1203 1204 >>> diff(assoc_laguerre(n, a, x), x) 1205 -assoc_laguerre(n - 1, a + 1, x) 1206 1207 >>> diff(assoc_laguerre(n, a, x), a) 1208 Sum(assoc_laguerre(_k, a, x)/(-a + n), (_k, 0, n - 1)) 1209 1210 Parameters 1211 ========== 1212 1213 n : int 1214 Degree of Laguerre polynomial. Must be ``n >= 0``. 1215 1216 alpha : Expr 1217 Arbitrary expression. For ``alpha=0`` regular Laguerre 1218 polynomials will be generated. 1219 1220 See Also 1221 ======== 1222 1223 jacobi, gegenbauer, 1224 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root, 1225 legendre, assoc_legendre, 1226 hermite, 1227 laguerre, 1228 sympy.polys.orthopolys.jacobi_poly 1229 sympy.polys.orthopolys.gegenbauer_poly 1230 sympy.polys.orthopolys.chebyshevt_poly 1231 sympy.polys.orthopolys.chebyshevu_poly 1232 sympy.polys.orthopolys.hermite_poly 1233 sympy.polys.orthopolys.legendre_poly 1234 sympy.polys.orthopolys.laguerre_poly 1235 1236 References 1237 ========== 1238 1239 .. [1] https://en.wikipedia.org/wiki/Laguerre_polynomial#Generalized_Laguerre_polynomials 1240 .. [2] http://mathworld.wolfram.com/AssociatedLaguerrePolynomial.html 1241 .. [3] http://functions.wolfram.com/Polynomials/LaguerreL/ 1242 .. [4] http://functions.wolfram.com/Polynomials/LaguerreL3/ 1243 1244 """ 1245 1246 @classmethod 1247 def eval(cls, n, alpha, x): 1248 # L_{n}^{0}(x) ---> L_{n}(x) 1249 if alpha.is_zero: 1250 return laguerre(n, x) 1251 1252 if not n.is_Number: 1253 # We can evaluate for some special values of x 1254 if x.is_zero: 1255 return binomial(n + alpha, alpha) 1256 elif x is S.Infinity and n > 0: 1257 return S.NegativeOne**n * S.Infinity 1258 elif x is S.NegativeInfinity and n > 0: 1259 return S.Infinity 1260 else: 1261 # n is a given fixed integer, evaluate into polynomial 1262 if n.is_negative: 1263 raise ValueError( 1264 "The index n must be nonnegative integer (got %r)" % n) 1265 else: 1266 return laguerre_poly(n, x, alpha) 1267 1268 def fdiff(self, argindex=3): 1269 from sympy import Sum 1270 if argindex == 1: 1271 # Diff wrt n 1272 raise ArgumentIndexError(self, argindex) 1273 elif argindex == 2: 1274 # Diff wrt alpha 1275 n, alpha, x = self.args 1276 k = Dummy("k") 1277 return Sum(assoc_laguerre(k, alpha, x) / (n - alpha), (k, 0, n - 1)) 1278 elif argindex == 3: 1279 # Diff wrt x 1280 n, alpha, x = self.args 1281 return -assoc_laguerre(n - 1, alpha + 1, x) 1282 else: 1283 raise ArgumentIndexError(self, argindex) 1284 1285 def _eval_rewrite_as_polynomial(self, n, alpha, x, **kwargs): 1286 from sympy import Sum 1287 # Make sure n \in N_0 1288 if n.is_negative or n.is_integer is False: 1289 raise ValueError("Error: n should be a non-negative integer.") 1290 k = Dummy("k") 1291 kern = RisingFactorial( 1292 -n, k) / (gamma(k + alpha + 1) * factorial(k)) * x**k 1293 return gamma(n + alpha + 1) / factorial(n) * Sum(kern, (k, 0, n)) 1294 1295 def _eval_conjugate(self): 1296 n, alpha, x = self.args 1297 return self.func(n, alpha.conjugate(), x.conjugate()) 1298