# Ticket #9706: orthogonal_polys.py

File orthogonal_polys.py, 28.1 KB (added by , 11 years ago) |
---|

Line | |
---|---|

1 | r""" |

2 | Orthogonal Polynomials |

3 | |

4 | This module wraps some of the orthogonal/special functions in the |

5 | Maxima package "orthopoly". This package was written by Barton |

6 | Willis of the University of Nebraska at Kearney. It is released |

7 | under the terms of the General Public License (GPL). Send |

8 | Maxima-related bug reports and comments on this module to |

9 | willisb@unk.edu. In your report, please include Maxima and specfun |

10 | version information. |

11 | |

12 | |

13 | - The Chebyshev polynomial of the first kind arises as a solution |

14 | to the differential equation |

15 | |

16 | .. math:: |

17 | |

18 | (1-x^2)\,y'' - x\,y' + n^2\,y = 0 |

19 | |

20 | |

21 | and those of the second kind as a solution to |

22 | |

23 | .. math:: |

24 | |

25 | (1-x^2)\,y'' - 3x\,y' + n(n+2)\,y = 0. |

26 | |

27 | |

28 | The Chebyshev polynomials of the first kind are defined by the |

29 | recurrence relation |

30 | |

31 | .. math:: |

32 | |

33 | T_0(x) = 1 \, T_1(x) = x \, T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x). \, |

34 | |

35 | |

36 | The Chebyshev polynomials of the second kind are defined by the |

37 | recurrence relation |

38 | |

39 | .. math:: |

40 | |

41 | U_0(x) = 1 \, U_1(x) = 2x \, U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x). \, |

42 | |

43 | |

44 | |

45 | For integers `m,n`, they satisfy the orthogonality |

46 | relations |

47 | |

48 | .. math:: |

49 | |

50 | \int_{-1}^1 T_n(x)T_m(x)\,\frac{dx}{\sqrt{1-x^2}} =\left\{ \begin{matrix} 0 &: n\ne m~~~~~\\ \pi &: n=m=0\\ \pi/2 &: n=m\ne 0 \end{matrix} \right. |

51 | |

52 | |

53 | and |

54 | |

55 | |

56 | .. math:: |

57 | |

58 | \int_{-1}^1 U_n(x)U_m(x)\sqrt{1-x^2}\,dx =\frac{\pi}{2}\delta_{m,n}. |

59 | |

60 | |

61 | |

62 | They are named after Pafnuty Chebyshev (alternative |

63 | transliterations: Tchebyshef or Tschebyscheff). |

64 | |

65 | - The Hermite polynomials are defined either by |

66 | |

67 | .. math:: |

68 | |

69 | H_n(x)=(-1)^n e^{x^2/2}\frac{d^n}{dx^n}e^{-x^2/2} |

70 | |

71 | |

72 | (the "probabilists' Hermite polynomials"), or by |

73 | |

74 | |

75 | .. math:: |

76 | |

77 | H_n(x)=(-1)^n e^{x^2}\frac{d^n}{dx^n}e^{-x^2} |

78 | |

79 | |

80 | (the "physicists' Hermite polynomials"). Sage (via Maxima) |

81 | implements the latter flavor. These satisfy the orthogonality |

82 | relation |

83 | |

84 | .. math:: |

85 | |

86 | \int_{-\infty}^\infty H_n(x)H_m(x)\,e^{-x^2}\,dx ={n!2^n}{\sqrt{\pi}}\delta_{nm} |

87 | |

88 | |

89 | |

90 | They are named in honor of Charles Hermite. |

91 | |

92 | - Each *Legendre polynomial* `P_n(x)` is an `n`-th degree polynomial. |

93 | It may be expressed using Rodrigues' formula: |

94 | |

95 | .. math:: |

96 | |

97 | P_n(x) = (2^n n!)^{-1} {\frac{d^n}{dx^n} } \left[ (x^2 -1)^n \right]. |

98 | |

99 | These are solutions to Legendre's differential equation: |

100 | |

101 | .. math:: |

102 | |

103 | {\frac{d}{dx}} \left[ (1-x^2) {\frac{d}{dx}} P(x) \right] + n(n+1)P(x) = 0. |

104 | |

105 | and satisfy the orthogonality relation |

106 | |

107 | .. math:: |

108 | |

109 | \int_{-1}^{1} P_m(x) P_n(x)\,dx = {\frac{2}{2n + 1}} \delta_{mn} |

110 | |

111 | The *Legendre function of the second kind* `Q_n(x)` is another |

112 | (linearly independent) solution to the Legendre differential equation. |

113 | It is not an "orthogonal polynomial" however. |

114 | |

115 | The associated Legendre functions of the first kind |

116 | `P_\ell^m(x)` can be given in terms of the "usual" |

117 | Legendre polynomials by |

118 | |

119 | .. math:: |

120 | |

121 | \begin{array}{ll} P_\ell^m(x) &= (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}P_\ell(x) \\ &= \frac{(-1)^m}{2^\ell \ell!} (1-x^2)^{m/2}\frac{d^{\ell+m}}{dx^{\ell+m}}(x^2-1)^\ell. \end{array} |

122 | |

123 | |

124 | Assuming `0 \le m \le \ell`, they satisfy the orthogonality |

125 | relation: |

126 | |

127 | .. math:: |

128 | |

129 | \int_{-1}^{1} P_k ^{(m)} P_\ell ^{(m)} dx = \frac{2 (\ell+m)!}{(2\ell+1)(\ell-m)!}\ \delta _{k,\ell}, |

130 | |

131 | |

132 | where `\delta _{k,\ell}` is the Kronecker delta. |

133 | |

134 | The associated Legendre functions of the second kind |

135 | `Q_\ell^m(x)` can be given in terms of the "usual" |

136 | Legendre polynomials by |

137 | |

138 | |

139 | .. math:: |

140 | |

141 | Q_\ell^m(x) = (-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}Q_\ell(x). |

142 | |

143 | |

144 | |

145 | They are named after Adrien-Marie Legendre. |

146 | |

147 | - Laguerre polynomials may be defined by the Rodrigues formula |

148 | |

149 | .. math:: |

150 | |

151 | L_n(x)=\frac{e^x}{n!}\frac{d^n}{dx^n}\left(e^{-x} x^n\right). |

152 | |

153 | |

154 | They are solutions of Laguerre's equation: |

155 | |

156 | |

157 | .. math:: |

158 | |

159 | x\,y'' + (1 - x)\,y' + n\,y = 0\, |

160 | |

161 | and satisfy the orthogonality relation |

162 | |

163 | |

164 | .. math:: |

165 | |

166 | \int_0^\infty L_m(x) L_n(x) e^{-x}\,dx = \delta_{mn}. |

167 | |

168 | |

169 | |

170 | The generalized Laguerre polynomials may be defined by the |

171 | Rodrigues formula: |

172 | |

173 | |

174 | .. math:: |

175 | |

176 | L_n^{(\alpha)}(x) = {\frac{x^{-\alpha} e^x}{n!}}{\frac{d^n}{dx^n}} \left(e^{-x} x^{n+\alpha}\right) . |

177 | |

178 | |

179 | (These are also sometimes called the associated Laguerre |

180 | polynomials.) The simple Laguerre polynomials are recovered from |

181 | the generalized polynomials by setting `\alpha =0`. |

182 | |

183 | They are named after Edmond Laguerre. |

184 | |

185 | - Jacobi polynomials are a class of orthogonal polynomials. They |

186 | are obtained from hypergeometric series in cases where the series |

187 | is in fact finite: |

188 | |

189 | .. math:: |

190 | |

191 | P_n^{(\alpha,\beta)}(z) =\frac{(\alpha+1)_n}{n!} \,_2F_1\left(-n,1+\alpha+\beta+n;\alpha+1;\frac{1-z}{2}\right) , |

192 | |

193 | |

194 | where `()_n` is Pochhammer's symbol (for the rising |

195 | factorial), (Abramowitz and Stegun p561.) and thus have the |

196 | explicit expression |

197 | |

198 | |

199 | .. math:: |

200 | |

201 | P_n^{(\alpha,\beta)} (z) = \frac{\Gamma (\alpha+n+1)}{n!\Gamma (\alpha+\beta+n+1)} \sum_{m=0}^n {n\choose m} \frac{\Gamma (\alpha + \beta + n + m + 1)}{\Gamma (\alpha + m + 1)} \left(\frac{z-1}{2}\right)^m . |

202 | |

203 | |

204 | |

205 | They are named after Carl Jacobi. |

206 | |

207 | - Ultraspherical or Gegenbauer polynomials are given in terms of |

208 | the Jacobi polynomials `P_n^{(\alpha,\beta)}(x)` with |

209 | `\alpha=\beta=a-1/2` by |

210 | |

211 | |

212 | .. math:: |

213 | |

214 | C_n^{(a)}(x)= \frac{\Gamma(a+1/2)}{\Gamma(2a)}\frac{\Gamma(n+2a)}{\Gamma(n+a+1/2)} P_n^{(a-1/2,a-1/2)}(x). |

215 | |

216 | |

217 | They satisfy the orthogonality relation |

218 | |

219 | .. math:: |

220 | |

221 | \int_{-1}^1(1-x^2)^{a-1/2}C_m^{(a)}(x)C_n^{(a)}(x)\, dx =\delta_{mn}2^{1-2a}\pi \frac{\Gamma(n+2a)}{(n+a)\Gamma^2(a)\Gamma(n+1)} , |

222 | |

223 | |

224 | for `a>-1/2`. They are obtained from hypergeometric series |

225 | in cases where the series is in fact finite: |

226 | |

227 | |

228 | .. math:: |

229 | |

230 | C_n^{(a)}(z) =\frac{(2a)^{\underline{n}}}{n!} \,_2F_1\left(-n,2a+n;a+\frac{1}{2};\frac{1-z}{2}\right) |

231 | |

232 | |

233 | where `\underline{n}` is the falling factorial. (See |

234 | Abramowitz and Stegun p561) |

235 | |

236 | They are named for Leopold Gegenbauer (1849-1903). |

237 | |

238 | |

239 | For completeness, the Pochhammer symbol, introduced by Leo August |

240 | Pochhammer, `(x)_n`, is used in the theory of special |

241 | functions to represent the "rising factorial" or "upper factorial" |

242 | |

243 | .. math:: |

244 | |

245 | (x)_n=x(x+1)(x+2)\cdots(x+n-1)=\frac{(x+n-1)!}{(x-1)!}. |

246 | |

247 | |

248 | On the other hand, the "falling factorial" or "lower factorial" is |

249 | |

250 | .. math:: |

251 | |

252 | x^{\underline{n}}=\frac{x!}{(x-n)!} , |

253 | |

254 | |

255 | in the notation of Ronald L. Graham, Donald E. Knuth and Oren |

256 | Patashnik in their book Concrete Mathematics. |

257 | |

258 | .. note:: |

259 | |

260 | The first call of any of these will usually cost a bit extra |

261 | (it loads "specfun", but I'm not sure if that is the real reason). |

262 | The next call is usually faster but not always. |

263 | |

264 | TODO: Implement associated Legendre polynomials and Zernike |

265 | polynomials. (Neither is in Maxima.) |

266 | http://en.wikipedia.org/wiki/Associated_Legendre_polynomials |

267 | http://en.wikipedia.org/wiki/Zernike_polynomials |

268 | |

269 | REFERENCES: |

270 | |

271 | - Abramowitz and Stegun: Handbook of Mathematical Functions, |

272 | http://www.math.sfu.ca/ cbm/aands/ |

273 | |

274 | - http://en.wikipedia.org/wiki/Chebyshev_polynomials |

275 | |

276 | - http://en.wikipedia.org/wiki/Legendre_polynomials |

277 | |

278 | - http://en.wikipedia.org/wiki/Hermite_polynomials |

279 | |

280 | - http://mathworld.wolfram.com/GegenbauerPolynomial.html |

281 | |

282 | - http://en.wikipedia.org/wiki/Jacobi_polynomials |

283 | |

284 | - http://en.wikipedia.org/wiki/Laguerre_polynomia |

285 | |

286 | - http://en.wikipedia.org/wiki/Associated_Legendre_polynomials |

287 | |

288 | |

289 | AUTHORS: |

290 | |

291 | - David Joyner (2006-06) |

292 | - Stefan Reiterer (2010-) |

293 | """ |

294 | |

295 | #***************************************************************************** |

296 | # Copyright (C) 2006 William Stein <wstein@gmail.com> |

297 | # 2006 David Joyner <wdj@usna.edu> |

298 | # |

299 | # Distributed under the terms of the GNU General Public License (GPL) |

300 | # |

301 | # This code is distributed in the hope that it will be useful, |

302 | # but WITHOUT ANY WARRANTY; without even the implied warranty of |

303 | # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

304 | # General Public License for more details. |

305 | # |

306 | # The full text of the GPL is available at: |

307 | # |

308 | # http://www.gnu.org/licenses/ |

309 | #***************************************************************************** |

310 | |

311 | import copy |

312 | import sage.plot.plot |

313 | import sage.interfaces.all |

314 | from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing |

315 | from sage.rings.rational_field import RationalField |

316 | from sage.rings.real_mpfr import RealField |

317 | from sage.misc.sage_eval import sage_eval |

318 | from sage.rings.all import QQ, ZZ, CDF, RDF |

319 | import sage.rings.commutative_ring as commutative_ring |

320 | import sage.rings.ring as ring |

321 | from sage.calculus.calculus import maxima |

322 | from sage.symbolic.function import BuiltinFunction, GinacFunction |

323 | from sage.symbolic.expression import is_Expression |

324 | import sage.functions.special |

325 | from sage.functions.special import MaximaFunction, meval |

326 | import numpy,scipy |

327 | import math |

328 | import sage.libs.mpmath.all as mpmath |

329 | |

330 | _done = False |

331 | def _init(): |

332 | """ |

333 | Internal function which checks if Maxima has loaded the |

334 | "orthopoly" package. All functions using this in this |

335 | file should call this function first. |

336 | |

337 | TEST: |

338 | |

339 | The global starts ``False``:: |

340 | |

341 | sage: sage.functions.orthogonal_polys._done |

342 | False |

343 | |

344 | Then after using one of these functions, it changes:: |

345 | |

346 | sage: from sage.functions.orthogonal_polys import chebyshev_T |

347 | sage: chebyshev_T(2,x) |

348 | 2*(x - 1)^2 + 4*x - 3 |

349 | sage: sage.functions.orthogonal_polys._done |

350 | True |

351 | |

352 | Note that because here we use a Pynac variable ``x``, |

353 | the representation of the function is different from |

354 | its actual doctest, where a polynomial indeterminate |

355 | ``x`` is used. |

356 | """ |

357 | global _done |

358 | if _done: |

359 | return |

360 | maxima.eval('load("orthopoly");') |

361 | # TODO -- make it possible to use the intervals returned |

362 | # instead of just discarding this info! |

363 | maxima.eval('orthopoly_returns_intervals:false;') |

364 | _done = True |

365 | |

366 | #load /home/maldun/sage/sage-4.5.1/devel/sage-ortho/sage/functions/orthogonal_polys.py |

367 | |

368 | class OrthogonalPolynomial(BuiltinFunction): |

369 | """ |

370 | Base Class for Orthogonal Polynomials. The evaluation as a polynomial |

371 | is done via maxima due performance reasons. Therefore the internal name |

372 | in maxima maxima_name has to be declared. |

373 | Convention: The first argument is always the order of the polynomial, |

374 | he last one is always the value x where the polynomial is evaluated. |

375 | |

376 | """ |

377 | def __init__(self, name, maxima_name, nargs = 2, latex_name = None, conversions = None): |

378 | self._maxima_name = maxima_name |

379 | BuiltinFunction.__init__(self, name = name, nargs = nargs, latex_name = latex_name, conversions = conversions) |

380 | |

381 | def _maxima_init_evaled_(self, *args): |

382 | """ |

383 | Returns a string which represents this function evaluated at |

384 | *args* in Maxima. |

385 | Remark: This function is stolen from the class MaximaFunction, |

386 | from sage.functions.special.py |

387 | |

388 | EXAMPLES:: |

389 | |

390 | sage: chebyshev_T(3,x) |

391 | sage: 4*(x - 1)^3 + 12*(x - 1)^2 + 9*x - 8 |

392 | """ |

393 | args_maxima = [] |

394 | for a in args: |

395 | if isinstance(a, str): |

396 | args_maxima.append(a) |

397 | elif hasattr(a, '_maxima_init_'): |

398 | args_maxima.append(a._maxima_init_()) |

399 | else: |

400 | args_maxima.append(str(a)) |

401 | return "%s(%s)"%(self._maxima_name, ', '.join(args_maxima)) |

402 | |

403 | def _clenshaw_method_(self,*args): |

404 | """ |

405 | The Clenshaw method uses the three term recursion of the polynomial, |

406 | or explicit formulas instead of maxima to evaluate the polynomial |

407 | efficiently, if the x argument is not a symbolic expression. |

408 | The name comes from the Clenshaw algorithm for fast evaluation of |

409 | polynomialsin chebyshev form. |

410 | |

411 | comparison: Maxima and Clenshaw algorithm for non-symbolic evaluation: |

412 | sage: time chebyshev_T(50,10) #clenshaw |

413 | CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s |

414 | Wall time: 0.00 s |

415 | 49656746733678312490954442369580252421769338391329426325400124999 |

416 | sage: time sage.functions.orthogonal_polys.chebyshev_T(50,10) |

417 | #maxima |

418 | CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s |

419 | Wall time: 0.05 s |

420 | 49656746733678312490954442369580252421769338391329426325400124999 |

421 | |

422 | sage: time chebyshev_T(500,10); #clenshaw |

423 | CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s |

424 | Wall time: 0.00 s |

425 | sage: time sage.functions.orthogonal_polys.chebyshev_T(500,10); |

426 | #maxima |

427 | CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s |

428 | Wall time: 0.77 s |

429 | """ |

430 | raise NotImplementedError( |

431 | "No recursive calculation of values implemented (yet)!") |

432 | |

433 | def _eval_special_values_(self,*args): |

434 | """ |

435 | Evals the polynomial explicitly for special values. |

436 | EXAMPLES: |

437 | sage: chebyshev_T(n,-1) |

438 | (-1)^n |

439 | """ |

440 | raise ValueError("No special values known!") |

441 | |

442 | def _eval_(self, *args): |

443 | """ |

444 | |

445 | The symbolic evaluation is done with maxima, because the evaluation of |

446 | the Polynomial representation seems to be quite clever. |

447 | For the fast numerical evaluation an other method should be used... |

448 | Therefore I suggest Clenshaw's algorithm, which uses the rekursion! |

449 | |

450 | EXAMPLES:: |

451 | |

452 | sage: time chebyshev_T(5,x) |

453 | CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s |

454 | Wall time: 0.16 s |

455 | 16*(x - 1)^5 + 80*(x - 1)^4 + 140*(x - 1)^3 + 100*(x - 1)^2 + 25*x - 24 |

456 | |

457 | time chebyshev_T(50,x) |

458 | CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s |

459 | Wall time: 0.20 s |

460 | 562949953421312*(x - 1)^50 + 28147497671065600*(x - 1)^49 +.... |

461 | |

462 | time chebyshev_T(100,x); |

463 | CPU times: user 0.12 s, sys: 0.00 s, total: 0.12 s |

464 | Wall time: 0.30 s |

465 | |

466 | The function also checks for special values, and if |

467 | the order is an integer and in range! |

468 | EXAMPLES:: |

469 | sage: var('n') |

470 | n |

471 | sage: chebyshev_T(n,-1) |

472 | (-1)^n |

473 | sage: chebyshev_T(-7,x) |

474 | ... |

475 | TypeError: Order is not valid! |

476 | sage: chebyshev_T(7.5,x) |

477 | TypeError: Order is not an Integer! |

478 | |

479 | """ |

480 | if not is_Expression(args[0]): |

481 | if args[0] != floor(args[0]): |

482 | raise TypeError("Order is not an Integer!") |

483 | |

484 | if args[0] < -1: |

485 | raise TypeError("Order is not valid!") |

486 | |

487 | try: |

488 | return self._eval_special_values_(*args) |

489 | except ValueError: |

490 | pass |

491 | |

492 | if not (is_Expression(args[0]) or is_Expression(args[-1])): |

493 | try: |

494 | return self._clenshaw_method_(*args) |

495 | except NotImplementedError: |

496 | pass |

497 | |

498 | _init() |

499 | try: |

500 | s = maxima(self._maxima_init_evaled_(*args)) |

501 | except TypeError: |

502 | return None |

503 | if self._maxima_name in repr(s): |

504 | return None |

505 | else: |

506 | return s.sage() |

507 | |

508 | |

509 | #TODO: Numerical Evals with help of the new scipy version!!!! |

510 | #Reason scipy 0.8 supports stable and fast numerical evaluation |

511 | #of ortho polys |

512 | |

513 | |

514 | class Func_chebyshev_T(OrthogonalPolynomial): |

515 | |

516 | """ |

517 | Class for the Chebyshev polynomial of the first kind for integers |

518 | `n>-1`. |

519 | |

520 | REFERENCE: |

521 | |

522 | - AS 22.5.31 page 778 and AS 6.1.22 page 256. |

523 | |

524 | EXAMPLES:: |

525 | |

526 | sage: chebyshev_T(3,x) |

527 | sage: 4*(x - 1)^3 + 12*(x - 1)^2 + 9*x - 8 |

528 | sage: var('k') |

529 | k |

530 | sage: test = chebyshev_T(k,x) |

531 | sage: test |

532 | chebyshev_T(k, x) |

533 | """ |

534 | |

535 | def __init__(self): |

536 | OrthogonalPolynomial.__init__(self,"chebyshev_T","chebyshev_t",nargs = 2, |

537 | conversions =dict(maxima='chebyshev_t',mathematica='ChebyshevT')) |

538 | |

539 | def _eval_special_values_(self,*args): |

540 | """ |

541 | Values known at the boundary. |

542 | EXAMPLES: |

543 | sage: chebyshev_T(n,1) |

544 | 1 |

545 | sage: chebyshev_T(n,-1) |

546 | (-1)^n |

547 | """ |

548 | if args[-1] == 1: |

549 | return 1 |

550 | |

551 | if args[-1] == -1: |

552 | return (-1)**args[0] |

553 | |

554 | raise ValueError("Value not found!") |

555 | |

556 | def _clenshaw_method_(self,*args): |

557 | """ |

558 | Clenshaw method for chebyshev_T (means use recursions in this case) |

559 | This is much faster for numerical evaluation than maxima! |

560 | """ |

561 | |

562 | k = args[0] |

563 | x = args[1] |

564 | |

565 | if k == -1: |

566 | return 0 |

567 | elif k == 0: |

568 | return 1 |

569 | elif k == 1: |

570 | return x |

571 | else: |

572 | #TODO: When evaluation of Symbolic Expressions works better |

573 | #use these explicit formulas instead! |

574 | #if -1 <= x <= 1: |

575 | # return cos(k*acos(x)) |

576 | #elif 1 < x: |

577 | # return cosh(k*acosh(x)) |

578 | #else: # x < -1 |

579 | # return (-1)**(k%2)*cosh(k*acosh(-x)) |

580 | |

581 | help1 = 1 |

582 | help2 = x |

583 | for j in xrange(0,k-1): |

584 | help3 = 2*x*help2 - help1 |

585 | help1 = help2 |

586 | help2 = help3 |

587 | |

588 | return help3 |

589 | |

590 | def _evalf_(self, *args,**kwds): |

591 | """ |

592 | Evals chebyshev_T numerically with mpmath. |

593 | EXAMPLES:: |

594 | chebyshev_T(10,3).n(75) |

595 | 2.261953700000000000000e7 |

596 | """ |

597 | parent = kwds['parent'] |

598 | return mpmath.call(mpmath.chebyt,args[0],args[-1],prec = parent.prec()) |

599 | |

600 | def _derivative_(self, *args, **kwds): |

601 | """ |

602 | Returns the derivative of chebyshev_T in form of the chebyshev Polynomial |

603 | of the second kind chebyshev_U |

604 | EXAMPLES:: |

605 | sage: derivative(chebyshev_T(k,x),x) |

606 | k*chebyshev_U(k - 1, x) |

607 | sage: derivative(chebyshev_T(3,x),x) |

608 | 12*(x - 1)^2 + 24*x - 15 |

609 | sage: derivative(chebyshev_T(k,x),k) |

610 | ... |

611 | NotImplementedError: Derivative w.r.t. to the index is not |

612 | supported, yet, and perhaps never will be... |

613 | |

614 | """ |

615 | diff_param = kwds['diff_param'] |

616 | if diff_param == 0: #args[0] is args[diff_param]: |

617 | raise NotImplementedError( |

618 | "Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...") |

619 | else: |

620 | return args[0]*chebyshev_U(args[0]-1,args[1]) |

621 | |

622 | chebyshev_T = Func_chebyshev_T() |

623 | |

624 | class Func_chebyshev_U(OrthogonalPolynomial): |

625 | |

626 | """ |

627 | Class for the Chebyshev polynomial of the second kind for integers `n>=-1`. |

628 | |

629 | REFERENCE: |

630 | |

631 | - AS, 22.8.3 page 783 and AS 6.1.22 page 256. |

632 | |

633 | EXAMPLES:: |

634 | |

635 | sage: chebyshev_U(3,x) |

636 | 8*(x - 1)^3 + 24*(x - 1)^2 + 20*x - 16 |

637 | """ |

638 | |

639 | def __init__(self): |

640 | OrthogonalPolynomial.__init__(self,"chebyshev_U","chebyshev_u", |

641 | nargs = 2,conversions =dict(maxima='chebyshev_u',mathematica='ChebyshevU')) |

642 | |

643 | def _clenshaw_method_(self,*args): |

644 | """ |

645 | Clenshaw method for chebyshev_T (means use the recursion...) |

646 | This is much faster for numerical evaluation than maxima! |

647 | """ |

648 | |

649 | k = args[0] |

650 | x = args[1] |

651 | |

652 | if k == -1: |

653 | return 0 |

654 | elif k == 0: |

655 | return 1 |

656 | elif k == 1: |

657 | return 2*x |

658 | else: |

659 | help1 = 1 |

660 | help2 = 2*x |

661 | for j in xrange(0,k-1): |

662 | help3 = 2*x*help2 - help1 |

663 | help1 = help2 |

664 | help2 = help3 |

665 | |

666 | return help3 |

667 | |

668 | def _eval_special_values_(self,*args): |

669 | """ |

670 | Values known at the boundary. |

671 | EXAMPLES: |

672 | sage: chebyshev_U(n,1) |

673 | n + 1 |

674 | sage: chebyshev_U(n,-1) |

675 | (n + 1)*(-1)^n |

676 | """ |

677 | if args[-1] == 1: |

678 | return (args[0]+1) |

679 | |

680 | if args[-1] == -1: |

681 | return (-1)**args[0]*(args[0]+1) |

682 | |

683 | raise ValueError("Value not found") |

684 | |

685 | def _evalf_(self, *args,**kwds): |

686 | """ |

687 | Evals chebyshev_U numerically with mpmath. |

688 | EXAMPLES:: |

689 | chebyshev_U(10,3).n(75) |

690 | 4.661117900000000000000e7 |

691 | """ |

692 | parent = kwds['parent'] |

693 | return mpmath.call(mpmath.chebyu,args[0],args[-1],prec = parent.prec()) |

694 | |

695 | def _derivative_(self, *args, **kwds): |

696 | """ |

697 | Returns the derivative of chebyshev_U in form of the chebyshev |

698 | Polynomials of the first and second kind |

699 | EXAMPLES:: |

700 | sage: derivative(chebyshev_U(k,x),x) |

701 | ((k + 1)*chebyshev_T(k + 1, x) - x*chebyshev_U(k, x))/(x^2 - 1) |

702 | sage: derivative(chebyshev_U(3,x),x) |

703 | 24*(x - 1)^2 + 48*x - 28 |

704 | sage: derivative(chebyshev_U(k,x),k) |

705 | ... |

706 | NotImplementedError: Derivative w.r.t. to the index is not |

707 | supported, yet, and perhaps never will be... |

708 | |

709 | """ |

710 | diff_param = kwds['diff_param'] |

711 | if diff_param == 0: #args[0] is args[diff_param]: |

712 | raise NotImplementedError( |

713 | "Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...") |

714 | else: |

715 | return ((args[0]+1)*chebyshev_T(args[0]+1,args[1])-args[1]* |

716 | chebyshev_U(args[0],args[1]))/(args[1]**2-1) |

717 | |

718 | chebyshev_U = Func_chebyshev_U() |

719 | |

720 | def gen_laguerre(n,a,x): |

721 | """ |

722 | Returns the generalized Laguerre polynomial for integers `n > -1`. |

723 | Typically, a = 1/2 or a = -1/2. |

724 | |

725 | REFERENCE: |

726 | |

727 | - table on page 789 in AS. |

728 | |

729 | EXAMPLES:: |

730 | |

731 | sage: x = PolynomialRing(QQ, 'x').gen() |

732 | sage: gen_laguerre(2,1,x) |

733 | 1/2*x^2 - 3*x + 3 |

734 | sage: gen_laguerre(2,1/2,x) |

735 | 1/2*x^2 - 5/2*x + 15/8 |

736 | sage: gen_laguerre(2,-1/2,x) |

737 | 1/2*x^2 - 3/2*x + 3/8 |

738 | sage: gen_laguerre(2,0,x) |

739 | 1/2*x^2 - 2*x + 1 |

740 | sage: gen_laguerre(3,0,x) |

741 | -1/6*x^3 + 3/2*x^2 - 3*x + 1 |

742 | """ |

743 | from sage.functions.all import sqrt |

744 | _init() |

745 | return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)), locals={'x':x}) |

746 | |

747 | def gen_legendre_P(n,m,x): |

748 | r""" |

749 | Returns the generalized (or associated) Legendre function of the |

750 | first kind for integers `n > -1, m > -1`. |

751 | |

752 | The awkward code for when m is odd and 1 results from the fact that |

753 | Maxima is happy with, for example, `(1 - t^2)^3/2`, but |

754 | Sage is not. For these cases the function is computed from the |

755 | (m-1)-case using one of the recursions satisfied by the Legendre |

756 | functions. |

757 | |

758 | REFERENCE: |

759 | |

760 | - Gradshteyn and Ryzhik 8.706 page 1000. |

761 | |

762 | EXAMPLES:: |

763 | |

764 | sage: P.<t> = QQ[] |

765 | sage: gen_legendre_P(2, 0, t) |

766 | 3/2*t^2 - 1/2 |

767 | sage: gen_legendre_P(2, 0, t) == legendre_P(2, t) |

768 | True |

769 | sage: gen_legendre_P(3, 1, t) |

770 | -3/2*sqrt(-t^2 + 1)*(5*t^2 - 1) |

771 | sage: gen_legendre_P(4, 3, t) |

772 | 105*sqrt(-t^2 + 1)*(t^2 - 1)*t |

773 | sage: gen_legendre_P(7, 3, I).expand() |

774 | -16695*sqrt(2) |

775 | sage: gen_legendre_P(4, 1, 2.5) |

776 | -583.562373654533*I |

777 | """ |

778 | from sage.functions.all import sqrt |

779 | _init() |

780 | if m.mod(2).is_zero() or m.is_one(): |

781 | return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x}) |

782 | else: |

783 | return sqrt(1-x**2)*(((n-m+1)*x*gen_legendre_P(n,m-1,x)-(n+m-1)*gen_legendre_P(n-1,m-1,x))/(1-x**2)) |

784 | |

785 | def gen_legendre_Q(n,m,x): |

786 | """ |

787 | Returns the generalized (or associated) Legendre function of the |

788 | second kind for integers `n>-1`, `m>-1`. |

789 | |

790 | Maxima restricts m = n. Hence the cases m n are computed using the |

791 | same recursion used for gen_legendre_P(n,m,x) when m is odd and |

792 | 1. |

793 | |

794 | EXAMPLES:: |

795 | |

796 | sage: P.<t> = QQ[] |

797 | sage: gen_legendre_Q(2,0,t) |

798 | 3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1)) |

799 | sage: gen_legendre_Q(2,0,t) - legendre_Q(2, t) |

800 | 0 |

801 | sage: gen_legendre_Q(3,1,0.5) |

802 | 2.49185259170895 |

803 | sage: gen_legendre_Q(0, 1, x) |

804 | -1/sqrt(-x^2 + 1) |

805 | sage: gen_legendre_Q(2, 4, x).factor() |

806 | 48*x/((x - 1)^2*(x + 1)^2) |

807 | """ |

808 | from sage.functions.all import sqrt |

809 | if m <= n: |

810 | _init() |

811 | return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'%(ZZ(n),ZZ(m))), locals={'x':x}) |

812 | if m == n + 1 or n == 0: |

813 | if m.mod(2).is_zero(): |

814 | denom = (1 - x**2)**(m/2) |

815 | else: |

816 | denom = sqrt(1 - x**2)*(1 - x**2)**((m-1)/2) |

817 | if m == n + 1: |

818 | return (-1)**m*(m-1).factorial()*2**n/denom |

819 | else: |

820 | return (-1)**m*(m-1).factorial()*((x+1)**m - (x-1)**m)/(2*denom) |

821 | else: |

822 | return ((n-m+1)*x*gen_legendre_Q(n,m-1,x)-(n+m-1)*gen_legendre_Q(n-1,m-1,x))/sqrt(1-x**2) |

823 | |

824 | def hermite(n,x): |

825 | """ |

826 | Returns the Hermite polynomial for integers `n > -1`. |

827 | |

828 | REFERENCE: |

829 | |

830 | - AS 22.5.40 and 22.5.41, page 779. |

831 | |

832 | EXAMPLES:: |

833 | |

834 | sage: x = PolynomialRing(QQ, 'x').gen() |

835 | sage: hermite(2,x) |

836 | 4*x^2 - 2 |

837 | sage: hermite(3,x) |

838 | 8*x^3 - 12*x |

839 | sage: hermite(3,2) |

840 | 40 |

841 | sage: S.<y> = PolynomialRing(RR) |

842 | sage: hermite(3,y) |

843 | 8.00000000000000*y^3 - 12.0000000000000*y |

844 | sage: R.<x,y> = QQ[] |

845 | sage: hermite(3,y^2) |

846 | 8*y^6 - 12*y^2 |

847 | sage: w = var('w') |

848 | sage: hermite(3,2*w) |

849 | 8*(8*w^2 - 3)*w |

850 | """ |

851 | _init() |

852 | return sage_eval(maxima.eval('hermite(%s,x)'%ZZ(n)), locals={'x':x}) |

853 | |

854 | def jacobi_P(n,a,b,x): |

855 | r""" |

856 | Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for |

857 | integers `n > -1` and a and b symbolic or `a > -1` |

858 | and `b > -1`. The Jacobi polynomials are actually defined |

859 | for all a and b. However, the Jacobi polynomial weight |

860 | `(1-x)^a(1+x)^b` isn't integrable for `a \leq -1` |

861 | or `b \leq -1`. |

862 | |

863 | REFERENCE: |

864 | |

865 | - table on page 789 in AS. |

866 | |

867 | EXAMPLES:: |

868 | |

869 | sage: x = PolynomialRing(QQ, 'x').gen() |

870 | sage: jacobi_P(2,0,0,x) |

871 | 3/2*x^2 - 1/2 |

872 | sage: jacobi_P(2,1,2,1.2) # random output of low order bits |

873 | 5.009999999999998 |

874 | """ |

875 | _init() |

876 | return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x}) |

877 | |

878 | def laguerre(n,x): |

879 | """ |

880 | Returns the Laguerre polynomial for integers `n > -1`. |

881 | |

882 | REFERENCE: |

883 | |

884 | - AS 22.5.16, page 778 and AS page 789. |

885 | |

886 | EXAMPLES:: |

887 | |

888 | sage: x = PolynomialRing(QQ, 'x').gen() |

889 | sage: laguerre(2,x) |

890 | 1/2*x^2 - 2*x + 1 |

891 | sage: laguerre(3,x) |

892 | -1/6*x^3 + 3/2*x^2 - 3*x + 1 |

893 | sage: laguerre(2,2) |

894 | -1 |

895 | """ |

896 | _init() |

897 | return sage_eval(maxima.eval('laguerre(%s,x)'%ZZ(n)), locals={'x':x}) |

898 | |

899 | def legendre_P(n,x): |

900 | """ |

901 | Returns the Legendre polynomial of the first kind for integers |

902 | `n > -1`. |

903 | |

904 | REFERENCE: |

905 | |

906 | - AS 22.5.35 page 779. |

907 | |

908 | EXAMPLES:: |

909 | |

910 | sage: P.<t> = QQ[] |

911 | sage: legendre_P(2,t) |

912 | 3/2*t^2 - 1/2 |

913 | sage: legendre_P(3, 1.1) |

914 | 1.67750000000000 |

915 | sage: legendre_P(3, 1 + I) |

916 | 7/2*I - 13/2 |

917 | sage: legendre_P(3, MatrixSpace(ZZ, 2)([1, 2, -4, 7])) |

918 | [-179 242] |

919 | [-484 547] |

920 | sage: legendre_P(3, GF(11)(5)) |

921 | 8 |

922 | """ |

923 | _init() |

924 | return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)), locals={'x':x}) |

925 | |

926 | def legendre_Q(n,x): |

927 | """ |

928 | Returns the Legendre function of the second kind for integers |

929 | `n>-1`. |

930 | |

931 | Computed using Maxima. |

932 | |

933 | EXAMPLES:: |

934 | |

935 | sage: P.<t> = QQ[] |

936 | sage: legendre_Q(2, t) |

937 | 3/4*t^2*log(-(t + 1)/(t - 1)) - 3/2*t - 1/4*log(-(t + 1)/(t - 1)) |

938 | sage: legendre_Q(3, 0.5) |

939 | -0.198654771479482 |

940 | sage: legendre_Q(4, 2) |

941 | 443/16*I*pi + 443/16*log(3) - 365/12 |

942 | sage: legendre_Q(4, 2.0) |

943 | 0.00116107583162324 + 86.9828465962674*I |

944 | """ |

945 | _init() |

946 | return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)), locals={'x':x}) |

947 | |

948 | def ultraspherical(n,a,x): |

949 | """ |

950 | Returns the ultraspherical (or Gegenbauer) polynomial for integers |

951 | `n > -1`. |

952 | |

953 | Computed using Maxima. |

954 | |

955 | REFERENCE: |

956 | |

957 | - AS 22.5.27 |

958 | |

959 | EXAMPLES:: |

960 | |

961 | sage: x = PolynomialRing(QQ, 'x').gen() |

962 | sage: ultraspherical(2,3/2,x) |

963 | 15/2*x^2 - 3/2 |

964 | sage: ultraspherical(2,1/2,x) |

965 | 3/2*x^2 - 1/2 |

966 | sage: ultraspherical(1,1,x) |

967 | 2*x |

968 | sage: t = PolynomialRing(RationalField(),"t").gen() |

969 | sage: gegenbauer(3,2,t) |

970 | 32*t^3 - 12*t |

971 | """ |

972 | _init() |

973 | return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)), locals={'x':x}) |

974 | |

975 | gegenbauer = ultraspherical |