Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/polynomials.py: 24%

1"""

2This module mainly implements special orthogonal polynomials.

4See also functions.combinatorial.numbers which contains some

5combinatorial polynomials.

7"""

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

21from sympy.polys.orthopolys import (chebyshevt_poly, chebyshevu_poly,

22 gegenbauer_poly, hermite_poly, hermite_prob_poly,

23 jacobi_poly, laguerre_poly, legendre_poly)

25_x = Dummy('x')

28class OrthogonalPolynomial(Function):

29 """Base class for orthogonal polynomials.

30 """

32 @classmethod

33 def _eval_at_order(cls, n, x):

34 if n.is_integer and n >= 0:

35 return cls._ortho_poly(int(n), _x).subs(_x, x)

37 def _eval_conjugate(self):

38 return self.func(self.args[0], self.args[1].conjugate())

40#----------------------------------------------------------------------------

41# Jacobi polynomials

42#

45class jacobi(OrthogonalPolynomial):

46 r"""

47 Jacobi polynomial $P_n^{\left(\alpha, \beta\right)}(x)$.

49 Explanation

50 ===========

52 ``jacobi(n, alpha, beta, x)`` gives the $n$th Jacobi polynomial

53 in $x$, $P_n^{\left(\alpha, \beta\right)}(x)$.

55 The Jacobi polynomials are orthogonal on $[-1, 1]$ with respect

56 to the weight $\left(1-x\right)^\alpha \left(1+x\right)^\beta$.

58 Examples

59 ========

61 >>> from sympy import jacobi, S, conjugate, diff

62 >>> from sympy.abc import a, b, n, x

64 >>> jacobi(0, a, b, x)

65 1

66 >>> jacobi(1, a, b, x)

67 a/2 - b/2 + x*(a/2 + b/2 + 1)

68 >>> jacobi(2, a, b, x)

69 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

71 >>> jacobi(n, a, b, x)

72 jacobi(n, a, b, x)

74 >>> jacobi(n, a, a, x)

75 RisingFactorial(a + 1, n)*gegenbauer(n,

76 a + 1/2, x)/RisingFactorial(2*a + 1, n)

78 >>> jacobi(n, 0, 0, x)

79 legendre(n, x)

81 >>> jacobi(n, S(1)/2, S(1)/2, x)

82 RisingFactorial(3/2, n)*chebyshevu(n, x)/factorial(n + 1)

84 >>> jacobi(n, -S(1)/2, -S(1)/2, x)

85 RisingFactorial(1/2, n)*chebyshevt(n, x)/factorial(n)

87 >>> jacobi(n, a, b, -x)

88 (-1)**n*jacobi(n, b, a, x)

90 >>> jacobi(n, a, b, 0)

91 gamma(a + n + 1)*hyper((-b - n, -n), (a + 1,), -1)/(2**n*factorial(n)*gamma(a + 1))

92 >>> jacobi(n, a, b, 1)

93 RisingFactorial(a + 1, n)/factorial(n)

95 >>> conjugate(jacobi(n, a, b, x))

96 jacobi(n, conjugate(a), conjugate(b), conjugate(x))

98 >>> diff(jacobi(n,a,b,x), x)

99 (a/2 + b/2 + n/2 + 1/2)*jacobi(n - 1, a + 1, b + 1, x)

100

101 See Also

102 ========

103

104 gegenbauer,

105 chebyshevt_root, chebyshevu, chebyshevu_root,

106 legendre, assoc_legendre,

107 hermite, hermite_prob,

108 laguerre, assoc_laguerre,

109 sympy.polys.orthopolys.jacobi_poly,

110 sympy.polys.orthopolys.gegenbauer_poly

111 sympy.polys.orthopolys.chebyshevt_poly

112 sympy.polys.orthopolys.chebyshevu_poly

113 sympy.polys.orthopolys.hermite_poly

114 sympy.polys.orthopolys.legendre_poly

115 sympy.polys.orthopolys.laguerre_poly

116

117 References

118 ==========

119

120 .. [1] https://en.wikipedia.org/wiki/Jacobi_polynomials

121 .. [2] https://mathworld.wolfram.com/JacobiPolynomial.html

122 .. [3] https://functions.wolfram.com/Polynomials/JacobiP/

123

124 """

125

126 @classmethod

127 def eval(cls, n, a, b, x):

128 # Simplify to other polynomials

129 # P^{a, a}_n(x)

130 if a == b:

131 if a == Rational(-1, 2):

132 return RisingFactorial(S.Half, n) / factorial(n) * chebyshevt(n, x)

133 elif a.is_zero:

134 return legendre(n, x)

135 elif a == S.Half:

136 return RisingFactorial(3*S.Half, n) / factorial(n + 1) * chebyshevu(n, x)

137 else:

138 return RisingFactorial(a + 1, n) / RisingFactorial(2*a + 1, n) * gegenbauer(n, a + S.Half, x)

139 elif b == -a:

140 # P^{a, -a}_n(x)

141 return gamma(n + a + 1) / gamma(n + 1) * (1 + x)**(a/2) / (1 - x)**(a/2) * assoc_legendre(n, -a, x)

142

143 if not n.is_Number:

144 # Symbolic result P^{a,b}_n(x)

145 # P^{a,b}_n(-x) ---> (-1)**n * P^{b,a}_n(-x)

146 if x.could_extract_minus_sign():

147 return S.NegativeOne**n * jacobi(n, b, a, -x)

148 # We can evaluate for some special values of x

149 if x.is_zero:

150 return (2**(-n) * gamma(a + n + 1) / (gamma(a + 1) * factorial(n)) *

151 hyper([-b - n, -n], [a + 1], -1))

152 if x == S.One:

153 return RisingFactorial(a + 1, n) / factorial(n)

154 elif x is S.Infinity:

155 if n.is_positive:

156 # Make sure a+b+2*n \notin Z

157 if (a + b + 2*n).is_integer:

158 raise ValueError("Error. a + b + 2*n should not be an integer.")

159 return RisingFactorial(a + b + n + 1, n) * S.Infinity

160 else:

161 # n is a given fixed integer, evaluate into polynomial

162 return jacobi_poly(n, a, b, x)

163

164 def fdiff(self, argindex=4):

165 from sympy.concrete.summations import Sum

166 if argindex == 1:

167 # Diff wrt n

168 raise ArgumentIndexError(self, argindex)

169 elif argindex == 2:

170 # Diff wrt a

171 n, a, b, x = self.args

172 k = Dummy("k")

173 f1 = 1 / (a + b + n + k + 1)

174 f2 = ((a + b + 2*k + 1) * RisingFactorial(b + k + 1, n - k) /

175 ((n - k) * RisingFactorial(a + b + k + 1, n - k)))

176 return Sum(f1 * (jacobi(n, a, b, x) + f2*jacobi(k, a, b, x)), (k, 0, n - 1))

177 elif argindex == 3:

178 # Diff wrt b

179 n, a, b, x = self.args

180 k = Dummy("k")

181 f1 = 1 / (a + b + n + k + 1)

182 f2 = (-1)**(n - k) * ((a + b + 2*k + 1) * RisingFactorial(a + k + 1, n - k) /

183 ((n - k) * RisingFactorial(a + b + k + 1, n - k)))

184 return Sum(f1 * (jacobi(n, a, b, x) + f2*jacobi(k, a, b, x)), (k, 0, n - 1))

185 elif argindex == 4:

186 # Diff wrt x

187 n, a, b, x = self.args

188 return S.Half * (a + b + n + 1) * jacobi(n - 1, a + 1, b + 1, x)

189 else:

190 raise ArgumentIndexError(self, argindex)

191

192 def _eval_rewrite_as_Sum(self, n, a, b, x, **kwargs):

193 from sympy.concrete.summations import Sum

194 # Make sure n \in N

195 if n.is_negative or n.is_integer is False:

196 raise ValueError("Error: n should be a non-negative integer.")

197 k = Dummy("k")

198 kern = (RisingFactorial(-n, k) * RisingFactorial(a + b + n + 1, k) * RisingFactorial(a + k + 1, n - k) /

199 factorial(k) * ((1 - x)/2)**k)

200 return 1 / factorial(n) * Sum(kern, (k, 0, n))

201

202 def _eval_rewrite_as_polynomial(self, n, a, b, x, **kwargs):

203 # This function is just kept for backwards compatibility

204 # but should not be used

205 return self._eval_rewrite_as_Sum(n, a, b, x, **kwargs)

206

207 def _eval_conjugate(self):

208 n, a, b, x = self.args

209 return self.func(n, a.conjugate(), b.conjugate(), x.conjugate())

210

211

212def jacobi_normalized(n, a, b, x):

213 r"""

214 Jacobi polynomial $P_n^{\left(\alpha, \beta\right)}(x)$.

215

216 Explanation

217 ===========

218

219 ``jacobi_normalized(n, alpha, beta, x)`` gives the $n$th

220 Jacobi polynomial in $x$, $P_n^{\left(\alpha, \beta\right)}(x)$.

221

222 The Jacobi polynomials are orthogonal on $[-1, 1]$ with respect

223 to the weight $\left(1-x\right)^\alpha \left(1+x\right)^\beta$.

224

225 This functions returns the polynomials normilzed:

226

227 .. math::

228

229 \int_{-1}^{1}

230 P_m^{\left(\alpha, \beta\right)}(x)

231 P_n^{\left(\alpha, \beta\right)}(x)

232 (1-x)^{\alpha} (1+x)^{\beta} \mathrm{d}x

233 = \delta_{m,n}

234

235 Examples

236 ========

237

238 >>> from sympy import jacobi_normalized

239 >>> from sympy.abc import n,a,b,x

240

241 >>> jacobi_normalized(n, a, b, x)

242 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)))

243

244 Parameters

245 ==========

246

247 n : integer degree of polynomial

249 a : alpha value

251 b : beta value

253 x : symbol

255 See Also

256 ========

257

258 gegenbauer,

259 chebyshevt_root, chebyshevu, chebyshevu_root,

260 legendre, assoc_legendre,

261 hermite, hermite_prob,

262 laguerre, assoc_laguerre,

263 sympy.polys.orthopolys.jacobi_poly,

264 sympy.polys.orthopolys.gegenbauer_poly

265 sympy.polys.orthopolys.chebyshevt_poly

266 sympy.polys.orthopolys.chebyshevu_poly

267 sympy.polys.orthopolys.hermite_poly

268 sympy.polys.orthopolys.legendre_poly

269 sympy.polys.orthopolys.laguerre_poly

270

271 References

272 ==========

273

274 .. [1] https://en.wikipedia.org/wiki/Jacobi_polynomials

275 .. [2] https://mathworld.wolfram.com/JacobiPolynomial.html

276 .. [3] https://functions.wolfram.com/Polynomials/JacobiP/

277

278 """

279 nfactor = (S(2)**(a + b + 1) * (gamma(n + a + 1) * gamma(n + b + 1))

280 / (2*n + a + b + 1) / (factorial(n) * gamma(n + a + b + 1)))

281

282 return jacobi(n, a, b, x) / sqrt(nfactor)

283

284

285#----------------------------------------------------------------------------

286# Gegenbauer polynomials

287#

288

289

290class gegenbauer(OrthogonalPolynomial):

291 r"""

292 Gegenbauer polynomial $C_n^{\left(\alpha\right)}(x)$.

293

294 Explanation

295 ===========

296

297 ``gegenbauer(n, alpha, x)`` gives the $n$th Gegenbauer polynomial

298 in $x$, $C_n^{\left(\alpha\right)}(x)$.

299

300 The Gegenbauer polynomials are orthogonal on $[-1, 1]$ with

301 respect to the weight $\left(1-x^2\right)^{\alpha-\frac{1}{2}}$.

302

303 Examples

304 ========

305

306 >>> from sympy import gegenbauer, conjugate, diff

307 >>> from sympy.abc import n,a,x

308 >>> gegenbauer(0, a, x)

309 1

310 >>> gegenbauer(1, a, x)

311 2*a*x

312 >>> gegenbauer(2, a, x)

313 -a + x**2*(2*a**2 + 2*a)

314 >>> gegenbauer(3, a, x)

315 x**3*(4*a**3/3 + 4*a**2 + 8*a/3) + x*(-2*a**2 - 2*a)

316

317 >>> gegenbauer(n, a, x)

318 gegenbauer(n, a, x)

319 >>> gegenbauer(n, a, -x)

320 (-1)**n*gegenbauer(n, a, x)

321

322 >>> gegenbauer(n, a, 0)

323 2**n*sqrt(pi)*gamma(a + n/2)/(gamma(a)*gamma(1/2 - n/2)*gamma(n + 1))

324 >>> gegenbauer(n, a, 1)

325 gamma(2*a + n)/(gamma(2*a)*gamma(n + 1))

326

327 >>> conjugate(gegenbauer(n, a, x))

328 gegenbauer(n, conjugate(a), conjugate(x))

329

330 >>> diff(gegenbauer(n, a, x), x)

331 2*a*gegenbauer(n - 1, a + 1, x)

332

333 See Also

334 ========

335

336 jacobi,

337 chebyshevt_root, chebyshevu, chebyshevu_root,

338 legendre, assoc_legendre,

339 hermite, hermite_prob,

340 laguerre, assoc_laguerre,

341 sympy.polys.orthopolys.jacobi_poly

342 sympy.polys.orthopolys.gegenbauer_poly

343 sympy.polys.orthopolys.chebyshevt_poly

344 sympy.polys.orthopolys.chebyshevu_poly

345 sympy.polys.orthopolys.hermite_poly

346 sympy.polys.orthopolys.hermite_prob_poly

347 sympy.polys.orthopolys.legendre_poly

348 sympy.polys.orthopolys.laguerre_poly

349

350 References

351 ==========

352

353 .. [1] https://en.wikipedia.org/wiki/Gegenbauer_polynomials

354 .. [2] https://mathworld.wolfram.com/GegenbauerPolynomial.html

355 .. [3] https://functions.wolfram.com/Polynomials/GegenbauerC3/

356

357 """

358

359 @classmethod

360 def eval(cls, n, a, x):

361 # For negative n the polynomials vanish

362 # See https://functions.wolfram.com/Polynomials/GegenbauerC3/03/01/03/0012/

363 if n.is_negative:

364 return S.Zero

365

366 # Some special values for fixed a

367 if a == S.Half:

368 return legendre(n, x)

369 elif a == S.One:

370 return chebyshevu(n, x)

371 elif a == S.NegativeOne:

372 return S.Zero

373

374 if not n.is_Number:

375 # Handle this before the general sign extraction rule

376 if x == S.NegativeOne:

377 if (re(a) > S.Half) == True:

378 return S.ComplexInfinity

379 else:

380 return (cos(S.Pi*(a+n)) * sec(S.Pi*a) * gamma(2*a+n) /

381 (gamma(2*a) * gamma(n+1)))

382

383 # Symbolic result C^a_n(x)

384 # C^a_n(-x) ---> (-1)**n * C^a_n(x)

385 if x.could_extract_minus_sign():

386 return S.NegativeOne**n * gegenbauer(n, a, -x)

387 # We can evaluate for some special values of x

388 if x.is_zero:

389 return (2**n * sqrt(S.Pi) * gamma(a + S.Half*n) /

390 (gamma((1 - n)/2) * gamma(n + 1) * gamma(a)) )

391 if x == S.One:

392 return gamma(2*a + n) / (gamma(2*a) * gamma(n + 1))

393 elif x is S.Infinity:

394 if n.is_positive:

395 return RisingFactorial(a, n) * S.Infinity

396 else:

397 # n is a given fixed integer, evaluate into polynomial

398 return gegenbauer_poly(n, a, x)

399

400 def fdiff(self, argindex=3):

401 from sympy.concrete.summations import Sum

402 if argindex == 1:

403 # Diff wrt n

404 raise ArgumentIndexError(self, argindex)

405 elif argindex == 2:

406 # Diff wrt a

407 n, a, x = self.args

408 k = Dummy("k")

409 factor1 = 2 * (1 + (-1)**(n - k)) * (k + a) / ((k +

410 n + 2*a) * (n - k))

411 factor2 = 2*(k + 1) / ((k + 2*a) * (2*k + 2*a + 1)) + \

412 2 / (k + n + 2*a)

413 kern = factor1*gegenbauer(k, a, x) + factor2*gegenbauer(n, a, x)

414 return Sum(kern, (k, 0, n - 1))

415 elif argindex == 3:

416 # Diff wrt x

417 n, a, x = self.args

418 return 2*a*gegenbauer(n - 1, a + 1, x)

419 else:

420 raise ArgumentIndexError(self, argindex)

421

422 def _eval_rewrite_as_Sum(self, n, a, x, **kwargs):

423 from sympy.concrete.summations import Sum

424 k = Dummy("k")

425 kern = ((-1)**k * RisingFactorial(a, n - k) * (2*x)**(n - 2*k) /

426 (factorial(k) * factorial(n - 2*k)))

427 return Sum(kern, (k, 0, floor(n/2)))

428

429 def _eval_rewrite_as_polynomial(self, n, a, x, **kwargs):

430 # This function is just kept for backwards compatibility

431 # but should not be used

432 return self._eval_rewrite_as_Sum(n, a, x, **kwargs)

433

434 def _eval_conjugate(self):

435 n, a, x = self.args

436 return self.func(n, a.conjugate(), x.conjugate())

437

438#----------------------------------------------------------------------------

439# Chebyshev polynomials of first and second kind

440#

441

442

443class chebyshevt(OrthogonalPolynomial):

444 r"""

445 Chebyshev polynomial of the first kind, $T_n(x)$.

446

447 Explanation

448 ===========

449

450 ``chebyshevt(n, x)`` gives the $n$th Chebyshev polynomial (of the first

451 kind) in $x$, $T_n(x)$.

452

453 The Chebyshev polynomials of the first kind are orthogonal on

454 $[-1, 1]$ with respect to the weight $\frac{1}{\sqrt{1-x^2}}$.

455

456 Examples

457 ========

458

459 >>> from sympy import chebyshevt, diff

460 >>> from sympy.abc import n,x

461 >>> chebyshevt(0, x)

462 1

463 >>> chebyshevt(1, x)

464 x

465 >>> chebyshevt(2, x)

466 2*x**2 - 1

467

468 >>> chebyshevt(n, x)

469 chebyshevt(n, x)

470 >>> chebyshevt(n, -x)

471 (-1)**n*chebyshevt(n, x)

472 >>> chebyshevt(-n, x)

473 chebyshevt(n, x)

474

475 >>> chebyshevt(n, 0)

476 cos(pi*n/2)

477 >>> chebyshevt(n, -1)

478 (-1)**n

479

480 >>> diff(chebyshevt(n, x), x)

481 n*chebyshevu(n - 1, x)

482

483 See Also

484 ========

485

486 jacobi, gegenbauer,

487 chebyshevt_root, chebyshevu, chebyshevu_root,

488 legendre, assoc_legendre,

489 hermite, hermite_prob,

490 laguerre, assoc_laguerre,

491 sympy.polys.orthopolys.jacobi_poly

492 sympy.polys.orthopolys.gegenbauer_poly

493 sympy.polys.orthopolys.chebyshevt_poly

494 sympy.polys.orthopolys.chebyshevu_poly

495 sympy.polys.orthopolys.hermite_poly

496 sympy.polys.orthopolys.hermite_prob_poly

497 sympy.polys.orthopolys.legendre_poly

498 sympy.polys.orthopolys.laguerre_poly

499

500 References

501 ==========

502

503 .. [1] https://en.wikipedia.org/wiki/Chebyshev_polynomial

504 .. [2] https://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html

505 .. [3] https://mathworld.wolfram.com/ChebyshevPolynomialoftheSecondKind.html

506 .. [4] https://functions.wolfram.com/Polynomials/ChebyshevT/

507 .. [5] https://functions.wolfram.com/Polynomials/ChebyshevU/

508

509 """

510

511 _ortho_poly = staticmethod(chebyshevt_poly)

512

513 @classmethod

514 def eval(cls, n, x):

515 if not n.is_Number:

516 # Symbolic result T_n(x)

517 # T_n(-x) ---> (-1)**n * T_n(x)

518 if x.could_extract_minus_sign():

519 return S.NegativeOne**n * chebyshevt(n, -x)

520 # T_{-n}(x) ---> T_n(x)

521 if n.could_extract_minus_sign():

522 return chebyshevt(-n, x)

523 # We can evaluate for some special values of x

524 if x.is_zero:

525 return cos(S.Half * S.Pi * n)

526 if x == S.One:

527 return S.One

528 elif x is S.Infinity:

529 return S.Infinity

530 else:

531 # n is a given fixed integer, evaluate into polynomial

532 if n.is_negative:

533 # T_{-n}(x) == T_n(x)

534 return cls._eval_at_order(-n, x)

535 else:

536 return cls._eval_at_order(n, x)

537

538 def fdiff(self, argindex=2):

539 if argindex == 1:

540 # Diff wrt n

541 raise ArgumentIndexError(self, argindex)

542 elif argindex == 2:

543 # Diff wrt x

544 n, x = self.args

545 return n * chebyshevu(n - 1, x)

546 else:

547 raise ArgumentIndexError(self, argindex)

548

549 def _eval_rewrite_as_Sum(self, n, x, **kwargs):

550 from sympy.concrete.summations import Sum

551 k = Dummy("k")

552 kern = binomial(n, 2*k) * (x**2 - 1)**k * x**(n - 2*k)

553 return Sum(kern, (k, 0, floor(n/2)))

554

555 def _eval_rewrite_as_polynomial(self, n, x, **kwargs):

556 # This function is just kept for backwards compatibility

557 # but should not be used

558 return self._eval_rewrite_as_Sum(n, x, **kwargs)

559

560

561class chebyshevu(OrthogonalPolynomial):

562 r"""

563 Chebyshev polynomial of the second kind, $U_n(x)$.

564

565 Explanation

566 ===========

567

568 ``chebyshevu(n, x)`` gives the $n$th Chebyshev polynomial of the second

569 kind in x, $U_n(x)$.

570

571 The Chebyshev polynomials of the second kind are orthogonal on

572 $[-1, 1]$ with respect to the weight $\sqrt{1-x^2}$.

573

574 Examples

575 ========

576

577 >>> from sympy import chebyshevu, diff

578 >>> from sympy.abc import n,x

579 >>> chebyshevu(0, x)

580 1

581 >>> chebyshevu(1, x)

582 2*x

583 >>> chebyshevu(2, x)

584 4*x**2 - 1

585

586 >>> chebyshevu(n, x)

587 chebyshevu(n, x)

588 >>> chebyshevu(n, -x)

589 (-1)**n*chebyshevu(n, x)

590 >>> chebyshevu(-n, x)

591 -chebyshevu(n - 2, x)

592

593 >>> chebyshevu(n, 0)

594 cos(pi*n/2)

595 >>> chebyshevu(n, 1)

596 n + 1

597

598 >>> diff(chebyshevu(n, x), x)

599 (-x*chebyshevu(n, x) + (n + 1)*chebyshevt(n + 1, x))/(x**2 - 1)

600

601 See Also

602 ========

603

604 jacobi, gegenbauer,

605 chebyshevt, chebyshevt_root, chebyshevu_root,

606 legendre, assoc_legendre,

607 hermite, hermite_prob,

608 laguerre, assoc_laguerre,

609 sympy.polys.orthopolys.jacobi_poly

610 sympy.polys.orthopolys.gegenbauer_poly

611 sympy.polys.orthopolys.chebyshevt_poly

612 sympy.polys.orthopolys.chebyshevu_poly

613 sympy.polys.orthopolys.hermite_poly

614 sympy.polys.orthopolys.hermite_prob_poly

615 sympy.polys.orthopolys.legendre_poly

616 sympy.polys.orthopolys.laguerre_poly

617

618 References

619 ==========

620

621 .. [1] https://en.wikipedia.org/wiki/Chebyshev_polynomial

622 .. [2] https://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html

623 .. [3] https://mathworld.wolfram.com/ChebyshevPolynomialoftheSecondKind.html

624 .. [4] https://functions.wolfram.com/Polynomials/ChebyshevT/

625 .. [5] https://functions.wolfram.com/Polynomials/ChebyshevU/

626

627 """

628

629 _ortho_poly = staticmethod(chebyshevu_poly)

630

631 @classmethod

632 def eval(cls, n, x):

633 if not n.is_Number:

634 # Symbolic result U_n(x)

635 # U_n(-x) ---> (-1)**n * U_n(x)

636 if x.could_extract_minus_sign():

637 return S.NegativeOne**n * chebyshevu(n, -x)

638 # U_{-n}(x) ---> -U_{n-2}(x)

639 if n.could_extract_minus_sign():

640 if n == S.NegativeOne:

641 # n can not be -1 here

642 return S.Zero

643 elif not (-n - 2).could_extract_minus_sign():

644 return -chebyshevu(-n - 2, x)

645 # We can evaluate for some special values of x

646 if x.is_zero:

647 return cos(S.Half * S.Pi * n)

648 if x == S.One:

649 return S.One + n

650 elif x is S.Infinity:

651 return S.Infinity

652 else:

653 # n is a given fixed integer, evaluate into polynomial

654 if n.is_negative:

655 # U_{-n}(x) ---> -U_{n-2}(x)

656 if n == S.NegativeOne:

657 return S.Zero

658 else:

659 return -cls._eval_at_order(-n - 2, x)

660 else:

661 return cls._eval_at_order(n, x)

662

663 def fdiff(self, argindex=2):

664 if argindex == 1:

665 # Diff wrt n

666 raise ArgumentIndexError(self, argindex)

667 elif argindex == 2:

668 # Diff wrt x

669 n, x = self.args

670 return ((n + 1) * chebyshevt(n + 1, x) - x * chebyshevu(n, x)) / (x**2 - 1)

671 else:

672 raise ArgumentIndexError(self, argindex)

673

674 def _eval_rewrite_as_Sum(self, n, x, **kwargs):

675 from sympy.concrete.summations import Sum

676 k = Dummy("k")

677 kern = S.NegativeOne**k * factorial(

678 n - k) * (2*x)**(n - 2*k) / (factorial(k) * factorial(n - 2*k))

679 return Sum(kern, (k, 0, floor(n/2)))

680

681 def _eval_rewrite_as_polynomial(self, n, x, **kwargs):

682 # This function is just kept for backwards compatibility

683 # but should not be used

684 return self._eval_rewrite_as_Sum(n, x, **kwargs)

685

686

687class chebyshevt_root(Function):

688 r"""

689 ``chebyshev_root(n, k)`` returns the $k$th root (indexed from zero) of

690 the $n$th Chebyshev polynomial of the first kind; that is, if

691 $0 \le k < n$, ``chebyshevt(n, chebyshevt_root(n, k)) == 0``.

692

693 Examples

694 ========

695

696 >>> from sympy import chebyshevt, chebyshevt_root

697 >>> chebyshevt_root(3, 2)

698 -sqrt(3)/2

699 >>> chebyshevt(3, chebyshevt_root(3, 2))

700 0

701

702 See Also

703 ========

704

705 jacobi, gegenbauer,

706 chebyshevt, chebyshevu, chebyshevu_root,

707 legendre, assoc_legendre,

708 hermite, hermite_prob,

709 laguerre, assoc_laguerre,

710 sympy.polys.orthopolys.jacobi_poly

711 sympy.polys.orthopolys.gegenbauer_poly

712 sympy.polys.orthopolys.chebyshevt_poly

713 sympy.polys.orthopolys.chebyshevu_poly

714 sympy.polys.orthopolys.hermite_poly

715 sympy.polys.orthopolys.hermite_prob_poly

716 sympy.polys.orthopolys.legendre_poly

717 sympy.polys.orthopolys.laguerre_poly

718 """

719

720 @classmethod

721 def eval(cls, n, k):

722 if not ((0 <= k) and (k < n)):

723 raise ValueError("must have 0 <= k < n, "

724 "got k = %s and n = %s" % (k, n))

725 return cos(S.Pi*(2*k + 1)/(2*n))

726

727

728class chebyshevu_root(Function):

729 r"""

730 ``chebyshevu_root(n, k)`` returns the $k$th root (indexed from zero) of the

731 $n$th Chebyshev polynomial of the second kind; that is, if $0 \le k < n$,

732 ``chebyshevu(n, chebyshevu_root(n, k)) == 0``.

733

734 Examples

735 ========

736

737 >>> from sympy import chebyshevu, chebyshevu_root

738 >>> chebyshevu_root(3, 2)

739 -sqrt(2)/2

740 >>> chebyshevu(3, chebyshevu_root(3, 2))

741 0

742

743 See Also

744 ========

745

746 chebyshevt, chebyshevt_root, chebyshevu,

747 legendre, assoc_legendre,

748 hermite, hermite_prob,

749 laguerre, assoc_laguerre,

750 sympy.polys.orthopolys.jacobi_poly

751 sympy.polys.orthopolys.gegenbauer_poly

752 sympy.polys.orthopolys.chebyshevt_poly

753 sympy.polys.orthopolys.chebyshevu_poly

754 sympy.polys.orthopolys.hermite_poly

755 sympy.polys.orthopolys.hermite_prob_poly

756 sympy.polys.orthopolys.legendre_poly

757 sympy.polys.orthopolys.laguerre_poly

758 """

759

760

761 @classmethod

762 def eval(cls, n, k):

763 if not ((0 <= k) and (k < n)):

764 raise ValueError("must have 0 <= k < n, "

765 "got k = %s and n = %s" % (k, n))

766 return cos(S.Pi*(k + 1)/(n + 1))

767

768#----------------------------------------------------------------------------

769# Legendre polynomials and Associated Legendre polynomials

770#

771

772

773class legendre(OrthogonalPolynomial):

774 r"""

775 ``legendre(n, x)`` gives the $n$th Legendre polynomial of $x$, $P_n(x)$

776

777 Explanation

778 ===========

779

780 The Legendre polynomials are orthogonal on $[-1, 1]$ with respect to

781 the constant weight 1. They satisfy $P_n(1) = 1$ for all $n$; further,

782 $P_n$ is odd for odd $n$ and even for even $n$.

783

784 Examples

785 ========

786

787 >>> from sympy import legendre, diff

788 >>> from sympy.abc import x, n

789 >>> legendre(0, x)

790 1

791 >>> legendre(1, x)

792 x

793 >>> legendre(2, x)

794 3*x**2/2 - 1/2

795 >>> legendre(n, x)

796 legendre(n, x)

797 >>> diff(legendre(n,x), x)

798 n*(x*legendre(n, x) - legendre(n - 1, x))/(x**2 - 1)

799

800 See Also

801 ========

802

803 jacobi, gegenbauer,

804 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root,

805 assoc_legendre,

806 hermite, hermite_prob,

807 laguerre, assoc_laguerre,

808 sympy.polys.orthopolys.jacobi_poly

809 sympy.polys.orthopolys.gegenbauer_poly

810 sympy.polys.orthopolys.chebyshevt_poly

811 sympy.polys.orthopolys.chebyshevu_poly

812 sympy.polys.orthopolys.hermite_poly

813 sympy.polys.orthopolys.hermite_prob_poly

814 sympy.polys.orthopolys.legendre_poly

815 sympy.polys.orthopolys.laguerre_poly

816

817 References

818 ==========

819

820 .. [1] https://en.wikipedia.org/wiki/Legendre_polynomial

821 .. [2] https://mathworld.wolfram.com/LegendrePolynomial.html

822 .. [3] https://functions.wolfram.com/Polynomials/LegendreP/

823 .. [4] https://functions.wolfram.com/Polynomials/LegendreP2/

824

825 """

826

827 _ortho_poly = staticmethod(legendre_poly)

828

829 @classmethod

830 def eval(cls, n, x):

831 if not n.is_Number:

832 # Symbolic result L_n(x)

833 # L_n(-x) ---> (-1)**n * L_n(x)

834 if x.could_extract_minus_sign():

835 return S.NegativeOne**n * legendre(n, -x)

836 # L_{-n}(x) ---> L_{n-1}(x)

837 if n.could_extract_minus_sign() and not(-n - 1).could_extract_minus_sign():

838 return legendre(-n - S.One, x)

839 # We can evaluate for some special values of x

840 if x.is_zero:

841 return sqrt(S.Pi)/(gamma(S.Half - n/2)*gamma(S.One + n/2))

842 elif x == S.One:

843 return S.One

844 elif x is S.Infinity:

845 return S.Infinity

846 else:

847 # n is a given fixed integer, evaluate into polynomial;

848 # L_{-n}(x) ---> L_{n-1}(x)

849 if n.is_negative:

850 n = -n - S.One

851 return cls._eval_at_order(n, x)

852

853 def fdiff(self, argindex=2):

854 if argindex == 1:

855 # Diff wrt n

856 raise ArgumentIndexError(self, argindex)

857 elif argindex == 2:

858 # Diff wrt x

859 # Find better formula, this is unsuitable for x = +/-1

860 # https://www.autodiff.org/ad16/Oral/Buecker_Legendre.pdf says

861 # at x = 1:

862 # n*(n + 1)/2 , m = 0

863 # oo , m = 1

864 # -(n-1)*n*(n+1)*(n+2)/4 , m = 2

865 # 0 , m = 3, 4, ..., n

866 #

867 # at x = -1

868 # (-1)**(n+1)*n*(n + 1)/2 , m = 0

869 # (-1)**n*oo , m = 1

870 # (-1)**n*(n-1)*n*(n+1)*(n+2)/4 , m = 2

871 # 0 , m = 3, 4, ..., n

872 n, x = self.args

873 return n/(x**2 - 1)*(x*legendre(n, x) - legendre(n - 1, x))

874 else:

875 raise ArgumentIndexError(self, argindex)

876

877 def _eval_rewrite_as_Sum(self, n, x, **kwargs):

878 from sympy.concrete.summations import Sum

879 k = Dummy("k")

880 kern = S.NegativeOne**k*binomial(n, k)**2*((1 + x)/2)**(n - k)*((1 - x)/2)**k

881 return Sum(kern, (k, 0, n))

882

883 def _eval_rewrite_as_polynomial(self, n, x, **kwargs):

884 # This function is just kept for backwards compatibility

885 # but should not be used

886 return self._eval_rewrite_as_Sum(n, x, **kwargs)

887

888

889class assoc_legendre(Function):

890 r"""

891 ``assoc_legendre(n, m, x)`` gives $P_n^m(x)$, where $n$ and $m$ are

892 the degree and order or an expression which is related to the nth

893 order Legendre polynomial, $P_n(x)$ in the following manner:

894

895 .. math::

896 P_n^m(x) = (-1)^m (1 - x^2)^{\frac{m}{2}}

897 \frac{\mathrm{d}^m P_n(x)}{\mathrm{d} x^m}

898

899 Explanation

900 ===========

901

902 Associated Legendre polynomials are orthogonal on $[-1, 1]$ with:

903

904 - weight $= 1$ for the same $m$ and different $n$.

905 - weight $= \frac{1}{1-x^2}$ for the same $n$ and different $m$.

906

907 Examples

908 ========

909

910 >>> from sympy import assoc_legendre

911 >>> from sympy.abc import x, m, n

912 >>> assoc_legendre(0,0, x)

913 1

914 >>> assoc_legendre(1,0, x)

915 x

916 >>> assoc_legendre(1,1, x)

917 -sqrt(1 - x**2)

918 >>> assoc_legendre(n,m,x)

919 assoc_legendre(n, m, x)

920

921 See Also

922 ========

923

924 jacobi, gegenbauer,

925 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root,

926 legendre,

927 hermite, hermite_prob,

928 laguerre, assoc_laguerre,

929 sympy.polys.orthopolys.jacobi_poly

930 sympy.polys.orthopolys.gegenbauer_poly

931 sympy.polys.orthopolys.chebyshevt_poly

932 sympy.polys.orthopolys.chebyshevu_poly

933 sympy.polys.orthopolys.hermite_poly

934 sympy.polys.orthopolys.hermite_prob_poly

935 sympy.polys.orthopolys.legendre_poly

936 sympy.polys.orthopolys.laguerre_poly

937

938 References

939 ==========

940

941 .. [1] https://en.wikipedia.org/wiki/Associated_Legendre_polynomials

942 .. [2] https://mathworld.wolfram.com/LegendrePolynomial.html

943 .. [3] https://functions.wolfram.com/Polynomials/LegendreP/

944 .. [4] https://functions.wolfram.com/Polynomials/LegendreP2/

945

946 """

947

948 @classmethod

949 def _eval_at_order(cls, n, m):

950 P = legendre_poly(n, _x, polys=True).diff((_x, m))

951 return S.NegativeOne**m * (1 - _x**2)**Rational(m, 2) * P.as_expr()

952

953 @classmethod

954 def eval(cls, n, m, x):

955 if m.could_extract_minus_sign():

956 # P^{-m}_n ---> F * P^m_n

957 return S.NegativeOne**(-m) * (factorial(m + n)/factorial(n - m)) * assoc_legendre(n, -m, x)

958 if m == 0:

959 # P^0_n ---> L_n

960 return legendre(n, x)

961 if x == 0:

962 return 2**m*sqrt(S.Pi) / (gamma((1 - m - n)/2)*gamma(1 - (m - n)/2))

963 if n.is_Number and m.is_Number and n.is_integer and m.is_integer:

964 if n.is_negative:

965 raise ValueError("%s : 1st index must be nonnegative integer (got %r)" % (cls, n))

966 if abs(m) > n:

967 raise ValueError("%s : abs('2nd index') must be <= '1st index' (got %r, %r)" % (cls, n, m))

968 return cls._eval_at_order(int(n), abs(int(m))).subs(_x, x)

969

970 def fdiff(self, argindex=3):

971 if argindex == 1:

972 # Diff wrt n

973 raise ArgumentIndexError(self, argindex)

974 elif argindex == 2:

975 # Diff wrt m

976 raise ArgumentIndexError(self, argindex)

977 elif argindex == 3:

978 # Diff wrt x

979 # Find better formula, this is unsuitable for x = 1

980 n, m, x = self.args

981 return 1/(x**2 - 1)*(x*n*assoc_legendre(n, m, x) - (m + n)*assoc_legendre(n - 1, m, x))

982 else:

983 raise ArgumentIndexError(self, argindex)

984

985 def _eval_rewrite_as_Sum(self, n, m, x, **kwargs):

986 from sympy.concrete.summations import Sum

987 k = Dummy("k")

988 kern = factorial(2*n - 2*k)/(2**n*factorial(n - k)*factorial(

989 k)*factorial(n - 2*k - m))*S.NegativeOne**k*x**(n - m - 2*k)

990 return (1 - x**2)**(m/2) * Sum(kern, (k, 0, floor((n - m)*S.Half)))

991

992 def _eval_rewrite_as_polynomial(self, n, m, x, **kwargs):

993 # This function is just kept for backwards compatibility

994 # but should not be used

995 return self._eval_rewrite_as_Sum(n, m, x, **kwargs)

996

997 def _eval_conjugate(self):

998 n, m, x = self.args

999 return self.func(n, m.conjugate(), x.conjugate())

1000

1001#----------------------------------------------------------------------------

1002# Hermite polynomials

1003#

1004

1005

1006class hermite(OrthogonalPolynomial):

1007 r"""

1008 ``hermite(n, x)`` gives the $n$th Hermite polynomial in $x$, $H_n(x)$.

1009

1010 Explanation

1011 ===========

1012

1013 The Hermite polynomials are orthogonal on $(-\infty, \infty)$

1014 with respect to the weight $\exp\left(-x^2\right)$.

1015

1016 Examples

1017 ========

1018

1019 >>> from sympy import hermite, diff

1020 >>> from sympy.abc import x, n

1021 >>> hermite(0, x)

1022 1

1023 >>> hermite(1, x)

1024 2*x

1025 >>> hermite(2, x)

1026 4*x**2 - 2

1027 >>> hermite(n, x)

1028 hermite(n, x)

1029 >>> diff(hermite(n,x), x)

1030 2*n*hermite(n - 1, x)

1031 >>> hermite(n, -x)

1032 (-1)**n*hermite(n, x)

1033

1034 See Also

1035 ========

1036

1037 jacobi, gegenbauer,

1038 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root,

1039 legendre, assoc_legendre,

1040 hermite_prob,

1041 laguerre, assoc_laguerre,

1042 sympy.polys.orthopolys.jacobi_poly

1043 sympy.polys.orthopolys.gegenbauer_poly

1044 sympy.polys.orthopolys.chebyshevt_poly

1045 sympy.polys.orthopolys.chebyshevu_poly

1046 sympy.polys.orthopolys.hermite_poly

1047 sympy.polys.orthopolys.hermite_prob_poly

1048 sympy.polys.orthopolys.legendre_poly

1049 sympy.polys.orthopolys.laguerre_poly

1050

1051 References

1052 ==========

1053

1054 .. [1] https://en.wikipedia.org/wiki/Hermite_polynomial

1055 .. [2] https://mathworld.wolfram.com/HermitePolynomial.html

1056 .. [3] https://functions.wolfram.com/Polynomials/HermiteH/

1057

1058 """

1059

1060 _ortho_poly = staticmethod(hermite_poly)

1061

1062 @classmethod

1063 def eval(cls, n, x):

1064 if not n.is_Number:

1065 # Symbolic result H_n(x)

1066 # H_n(-x) ---> (-1)**n * H_n(x)

1067 if x.could_extract_minus_sign():

1068 return S.NegativeOne**n * hermite(n, -x)

1069 # We can evaluate for some special values of x

1070 if x.is_zero:

1071 return 2**n * sqrt(S.Pi) / gamma((S.One - n)/2)

1072 elif x is S.Infinity:

1073 return S.Infinity

1074 else:

1075 # n is a given fixed integer, evaluate into polynomial

1076 if n.is_negative:

1077 raise ValueError(

1078 "The index n must be nonnegative integer (got %r)" % n)

1079 else:

1080 return cls._eval_at_order(n, x)

1081

1082 def fdiff(self, argindex=2):

1083 if argindex == 1:

1084 # Diff wrt n

1085 raise ArgumentIndexError(self, argindex)

1086 elif argindex == 2:

1087 # Diff wrt x

1088 n, x = self.args

1089 return 2*n*hermite(n - 1, x)

1090 else:

1091 raise ArgumentIndexError(self, argindex)

1092

1093 def _eval_rewrite_as_Sum(self, n, x, **kwargs):

1094 from sympy.concrete.summations import Sum

1095 k = Dummy("k")

1096 kern = S.NegativeOne**k / (factorial(k)*factorial(n - 2*k)) * (2*x)**(n - 2*k)

1097 return factorial(n)*Sum(kern, (k, 0, floor(n/2)))

1098

1099 def _eval_rewrite_as_polynomial(self, n, x, **kwargs):

1100 # This function is just kept for backwards compatibility

1101 # but should not be used

1102 return self._eval_rewrite_as_Sum(n, x, **kwargs)

1103

1104 def _eval_rewrite_as_hermite_prob(self, n, x, **kwargs):

1105 return sqrt(2)**n * hermite_prob(n, x*sqrt(2))

1106

1107

1108class hermite_prob(OrthogonalPolynomial):

1109 r"""

1110 ``hermite_prob(n, x)`` gives the $n$th probabilist's Hermite polynomial

1111 in $x$, $He_n(x)$.

1112

1113 Explanation

1114 ===========

1115

1116 The probabilist's Hermite polynomials are orthogonal on $(-\infty, \infty)$

1117 with respect to the weight $\exp\left(-\frac{x^2}{2}\right)$. They are monic

1118 polynomials, related to the plain Hermite polynomials (:py:class:`~.hermite`) by

1119

1120 .. math :: He_n(x) = 2^{-n/2} H_n(x/\sqrt{2})

1121

1122 Examples

1123 ========

1124

1125 >>> from sympy import hermite_prob, diff, I

1126 >>> from sympy.abc import x, n

1127 >>> hermite_prob(1, x)

1128 x

1129 >>> hermite_prob(5, x)

1130 x**5 - 10*x**3 + 15*x

1131 >>> diff(hermite_prob(n,x), x)

1132 n*hermite_prob(n - 1, x)

1133 >>> hermite_prob(n, -x)

1134 (-1)**n*hermite_prob(n, x)

1135

1136 The sum of absolute values of coefficients of $He_n(x)$ is the number of

1137 matchings in the complete graph $K_n$ or telephone number, A000085 in the OEIS:

1138

1139 >>> [hermite_prob(n,I) / I**n for n in range(11)]

1140 [1, 1, 2, 4, 10, 26, 76, 232, 764, 2620, 9496]

1141

1142 See Also

1143 ========

1144

1145 jacobi, gegenbauer,

1146 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root,

1147 legendre, assoc_legendre,

1148 hermite,

1149 laguerre, assoc_laguerre,

1150 sympy.polys.orthopolys.jacobi_poly

1151 sympy.polys.orthopolys.gegenbauer_poly

1152 sympy.polys.orthopolys.chebyshevt_poly

1153 sympy.polys.orthopolys.chebyshevu_poly

1154 sympy.polys.orthopolys.hermite_poly

1155 sympy.polys.orthopolys.hermite_prob_poly

1156 sympy.polys.orthopolys.legendre_poly

1157 sympy.polys.orthopolys.laguerre_poly

1158

1159 References

1160 ==========

1161

1162 .. [1] https://en.wikipedia.org/wiki/Hermite_polynomial

1163 .. [2] https://mathworld.wolfram.com/HermitePolynomial.html

1164 """

1165

1166 _ortho_poly = staticmethod(hermite_prob_poly)

1167

1168 @classmethod

1169 def eval(cls, n, x):

1170 if not n.is_Number:

1171 if x.could_extract_minus_sign():

1172 return S.NegativeOne**n * hermite_prob(n, -x)

1173 if x.is_zero:

1174 return sqrt(S.Pi) / gamma((S.One-n) / 2)

1175 elif x is S.Infinity:

1176 return S.Infinity

1177 else:

1178 if n.is_negative:

1179 ValueError("n must be a nonnegative integer, not %r" % n)

1180 else:

1181 return cls._eval_at_order(n, x)

1182

1183 def fdiff(self, argindex=2):

1184 if argindex == 2:

1185 n, x = self.args

1186 return n*hermite_prob(n-1, x)

1187 else:

1188 raise ArgumentIndexError(self, argindex)

1189

1190 def _eval_rewrite_as_Sum(self, n, x, **kwargs):

1191 from sympy.concrete.summations import Sum

1192 k = Dummy("k")

1193 kern = (-S.Half)**k * x**(n-2*k) / (factorial(k) * factorial(n-2*k))

1194 return factorial(n)*Sum(kern, (k, 0, floor(n/2)))

1195

1196 def _eval_rewrite_as_polynomial(self, n, x, **kwargs):

1197 # This function is just kept for backwards compatibility

1198 # but should not be used

1199 return self._eval_rewrite_as_Sum(n, x, **kwargs)

1200

1201 def _eval_rewrite_as_hermite(self, n, x, **kwargs):

1202 return sqrt(2)**(-n) * hermite(n, x/sqrt(2))

1203

1204

1205#----------------------------------------------------------------------------

1206# Laguerre polynomials

1207#

1208

1209

1210class laguerre(OrthogonalPolynomial):

1211 r"""

1212 Returns the $n$th Laguerre polynomial in $x$, $L_n(x)$.

1213

1214 Examples

1215 ========

1216

1217 >>> from sympy import laguerre, diff

1218 >>> from sympy.abc import x, n

1219 >>> laguerre(0, x)

1220 1

1221 >>> laguerre(1, x)

1222 1 - x

1223 >>> laguerre(2, x)

1224 x**2/2 - 2*x + 1

1225 >>> laguerre(3, x)

1226 -x**3/6 + 3*x**2/2 - 3*x + 1

1227

1228 >>> laguerre(n, x)

1229 laguerre(n, x)

1230

1231 >>> diff(laguerre(n, x), x)

1232 -assoc_laguerre(n - 1, 1, x)

1233

1234 Parameters

1235 ==========

1236

1237 n : int

1238 Degree of Laguerre polynomial. Must be `n \ge 0`.

1239

1240 See Also

1241 ========

1242

1243 jacobi, gegenbauer,

1244 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root,

1245 legendre, assoc_legendre,

1246 hermite, hermite_prob,

1247 assoc_laguerre,

1248 sympy.polys.orthopolys.jacobi_poly

1249 sympy.polys.orthopolys.gegenbauer_poly

1250 sympy.polys.orthopolys.chebyshevt_poly

1251 sympy.polys.orthopolys.chebyshevu_poly

1252 sympy.polys.orthopolys.hermite_poly

1253 sympy.polys.orthopolys.hermite_prob_poly

1254 sympy.polys.orthopolys.legendre_poly

1255 sympy.polys.orthopolys.laguerre_poly

1256

1257 References

1258 ==========

1259

1260 .. [1] https://en.wikipedia.org/wiki/Laguerre_polynomial

1261 .. [2] https://mathworld.wolfram.com/LaguerrePolynomial.html

1262 .. [3] https://functions.wolfram.com/Polynomials/LaguerreL/

1263 .. [4] https://functions.wolfram.com/Polynomials/LaguerreL3/

1264

1265 """

1266

1267 _ortho_poly = staticmethod(laguerre_poly)

1268

1269 @classmethod

1270 def eval(cls, n, x):

1271 if n.is_integer is False:

1272 raise ValueError("Error: n should be an integer.")

1273 if not n.is_Number:

1274 # Symbolic result L_n(x)

1275 # L_{n}(-x) ---> exp(-x) * L_{-n-1}(x)

1276 # L_{-n}(x) ---> exp(x) * L_{n-1}(-x)

1277 if n.could_extract_minus_sign() and not(-n - 1).could_extract_minus_sign():

1278 return exp(x)*laguerre(-n - 1, -x)

1279 # We can evaluate for some special values of x

1280 if x.is_zero:

1281 return S.One

1282 elif x is S.NegativeInfinity:

1283 return S.Infinity

1284 elif x is S.Infinity:

1285 return S.NegativeOne**n * S.Infinity

1286 else:

1287 if n.is_negative:

1288 return exp(x)*laguerre(-n - 1, -x)

1289 else:

1290 return cls._eval_at_order(n, x)

1291

1292 def fdiff(self, argindex=2):

1293 if argindex == 1:

1294 # Diff wrt n

1295 raise ArgumentIndexError(self, argindex)

1296 elif argindex == 2:

1297 # Diff wrt x

1298 n, x = self.args

1299 return -assoc_laguerre(n - 1, 1, x)

1300 else:

1301 raise ArgumentIndexError(self, argindex)

1302

1303 def _eval_rewrite_as_Sum(self, n, x, **kwargs):

1304 from sympy.concrete.summations import Sum

1305 # Make sure n \in N_0

1306 if n.is_negative:

1307 return exp(x) * self._eval_rewrite_as_Sum(-n - 1, -x, **kwargs)

1308 if n.is_integer is False:

1309 raise ValueError("Error: n should be an integer.")

1310 k = Dummy("k")

1311 kern = RisingFactorial(-n, k) / factorial(k)**2 * x**k

1312 return Sum(kern, (k, 0, n))

1313

1314 def _eval_rewrite_as_polynomial(self, n, x, **kwargs):

1315 # This function is just kept for backwards compatibility

1316 # but should not be used

1317 return self._eval_rewrite_as_Sum(n, x, **kwargs)

1318

1319

1320class assoc_laguerre(OrthogonalPolynomial):

1321 r"""

1322 Returns the $n$th generalized Laguerre polynomial in $x$, $L_n(x)$.

1323

1324 Examples

1325 ========

1326

1327 >>> from sympy import assoc_laguerre, diff

1328 >>> from sympy.abc import x, n, a

1329 >>> assoc_laguerre(0, a, x)

1330 1

1331 >>> assoc_laguerre(1, a, x)

1332 a - x + 1

1333 >>> assoc_laguerre(2, a, x)

1334 a**2/2 + 3*a/2 + x**2/2 + x*(-a - 2) + 1

1335 >>> assoc_laguerre(3, a, x)

1336 a**3/6 + a**2 + 11*a/6 - x**3/6 + x**2*(a/2 + 3/2) +

1337 x*(-a**2/2 - 5*a/2 - 3) + 1

1338

1339 >>> assoc_laguerre(n, a, 0)

1340 binomial(a + n, a)

1341

1342 >>> assoc_laguerre(n, a, x)

1343 assoc_laguerre(n, a, x)

1344

1345 >>> assoc_laguerre(n, 0, x)

1346 laguerre(n, x)

1347

1348 >>> diff(assoc_laguerre(n, a, x), x)

1349 -assoc_laguerre(n - 1, a + 1, x)

1350

1351 >>> diff(assoc_laguerre(n, a, x), a)

1352 Sum(assoc_laguerre(_k, a, x)/(-a + n), (_k, 0, n - 1))

1353

1354 Parameters

1355 ==========

1356

1357 n : int

1358 Degree of Laguerre polynomial. Must be `n \ge 0`.

1359

1360 alpha : Expr

1361 Arbitrary expression. For ``alpha=0`` regular Laguerre

1362 polynomials will be generated.

1363

1364 See Also

1365 ========

1366

1367 jacobi, gegenbauer,

1368 chebyshevt, chebyshevt_root, chebyshevu, chebyshevu_root,

1369 legendre, assoc_legendre,

1370 hermite, hermite_prob,

1371 laguerre,

1372 sympy.polys.orthopolys.jacobi_poly

1373 sympy.polys.orthopolys.gegenbauer_poly

1374 sympy.polys.orthopolys.chebyshevt_poly

1375 sympy.polys.orthopolys.chebyshevu_poly

1376 sympy.polys.orthopolys.hermite_poly

1377 sympy.polys.orthopolys.hermite_prob_poly

1378 sympy.polys.orthopolys.legendre_poly

1379 sympy.polys.orthopolys.laguerre_poly

1380

1381 References

1382 ==========

1383

1384 .. [1] https://en.wikipedia.org/wiki/Laguerre_polynomial#Generalized_Laguerre_polynomials

1385 .. [2] https://mathworld.wolfram.com/AssociatedLaguerrePolynomial.html

1386 .. [3] https://functions.wolfram.com/Polynomials/LaguerreL/

1387 .. [4] https://functions.wolfram.com/Polynomials/LaguerreL3/

1388

1389 """

1390

1391 @classmethod

1392 def eval(cls, n, alpha, x):

1393 # L_{n}^{0}(x) ---> L_{n}(x)

1394 if alpha.is_zero:

1395 return laguerre(n, x)

1396

1397 if not n.is_Number:

1398 # We can evaluate for some special values of x

1399 if x.is_zero:

1400 return binomial(n + alpha, alpha)

1401 elif x is S.Infinity and n > 0:

1402 return S.NegativeOne**n * S.Infinity

1403 elif x is S.NegativeInfinity and n > 0:

1404 return S.Infinity

1405 else:

1406 # n is a given fixed integer, evaluate into polynomial

1407 if n.is_negative:

1408 raise ValueError(

1409 "The index n must be nonnegative integer (got %r)" % n)

1410 else:

1411 return laguerre_poly(n, x, alpha)

1412

1413 def fdiff(self, argindex=3):

1414 from sympy.concrete.summations import Sum

1415 if argindex == 1:

1416 # Diff wrt n

1417 raise ArgumentIndexError(self, argindex)

1418 elif argindex == 2:

1419 # Diff wrt alpha

1420 n, alpha, x = self.args

1421 k = Dummy("k")

1422 return Sum(assoc_laguerre(k, alpha, x) / (n - alpha), (k, 0, n - 1))

1423 elif argindex == 3:

1424 # Diff wrt x

1425 n, alpha, x = self.args

1426 return -assoc_laguerre(n - 1, alpha + 1, x)

1427 else:

1428 raise ArgumentIndexError(self, argindex)

1429

1430 def _eval_rewrite_as_Sum(self, n, alpha, x, **kwargs):

1431 from sympy.concrete.summations import Sum

1432 # Make sure n \in N_0

1433 if n.is_negative or n.is_integer is False:

1434 raise ValueError("Error: n should be a non-negative integer.")

1435 k = Dummy("k")

1436 kern = RisingFactorial(

1437 -n, k) / (gamma(k + alpha + 1) * factorial(k)) * x**k

1438 return gamma(n + alpha + 1) / factorial(n) * Sum(kern, (k, 0, n))

1439

1440 def _eval_rewrite_as_polynomial(self, n, alpha, x, **kwargs):

1441 # This function is just kept for backwards compatibility

1442 # but should not be used

1443 return self._eval_rewrite_as_Sum(n, alpha, x, **kwargs)

1444

1445 def _eval_conjugate(self):

1446 n, alpha, x = self.args

1447 return self.func(n, alpha.conjugate(), x.conjugate())