# Ticket #9706: orthogonal_polys.2.py

File orthogonal_polys.2.py, 33.2 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 | from sage.functions.other import floor |

330 | import sage.symbolic.expression as expression |

331 | |

332 | _done = False |

333 | def _init(): |

334 | """ |

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

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

337 | file should call this function first. |

338 | |

339 | TEST: |

340 | |

341 | The global starts ``False``:: |

342 | |

343 | sage: sage.functions.orthogonal_polys._done |

344 | False |

345 | |

346 | Then after using one of these functions, it changes:: (The value is now |

347 | False for chebyshev_T because chebyshev_T uses clenshaw method instead...) |

348 | |

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

350 | sage: chebyshev_T(2,x) |

351 | 2*x^2 - 1 |

352 | sage: sage.functions.orthogonal_polys._done |

353 | False |

354 | sage: gegenbauer(1,2,3) |

355 | 12 |

356 | sage: sage.functions.orthogonal_polys._done |

357 | True |

358 | |

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

360 | the representation of the function is different from |

361 | its actual doctest, where a polynomial indeterminate |

362 | ``x`` is used. |

363 | """ |

364 | global _done |

365 | if _done: |

366 | return |

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

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

369 | # instead of just discarding this info! |

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

371 | _done = True |

372 | |

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

374 | |

375 | class OrthogonalPolynomial(BuiltinFunction): |

376 | """ |

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

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

379 | in maxima maxima_name has to be declared. |

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

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

382 | |

383 | """ |

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

385 | self._maxima_name = maxima_name |

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

387 | |

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

389 | """ |

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

391 | *args* in Maxima. |

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

393 | from sage.functions.special.py |

394 | Comment: The usefulness is in question.... |

395 | |

396 | EXAMPLES:: |

397 | |

398 | sage: chebyshev_T(3,x) |

399 | 4*x^3 - 3*x |

400 | """ |

401 | args_maxima = [] |

402 | for a in args: |

403 | if isinstance(a, str): |

404 | args_maxima.append(a) |

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

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

407 | else: |

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

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

410 | |

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

412 | """ |

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

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

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

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

417 | polynomialsin chebyshev form. |

418 | |

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

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

421 | #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s |

422 | #Wall time: 0.00 s |

423 | #49656746733678312490954442369580252421769338391329426325400124999 |

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

425 | #maxima |

426 | #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s |

427 | #Wall time: 0.05 s |

428 | #49656746733678312490954442369580252421769338391329426325400124999 |

429 | |

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

431 | #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s |

432 | #Wall time: 0.00 s |

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

434 | #maxima |

435 | #CPU times: user 0.11 s, sys: 0.00 s, total: 0.11 s |

436 | #Wall time: 0.77 s |

437 | """ |

438 | raise NotImplementedError( |

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

440 | |

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

442 | """ |

443 | Evals the polynomial explicitly for special values. |

444 | EXAMPLES: |

445 | |

446 | sage: var('n') |

447 | n |

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

449 | (-1)^n |

450 | """ |

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

452 | |

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

454 | """ |

455 | |

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

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

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

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

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

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

462 | |

463 | Update: I take back my argument about maxima... |

464 | it seems that just adding the expand command to |

465 | the for loop in the Clanshaw algorithm speeds up |

466 | the symbolic evaluation enormously. |

467 | Now the maxima option stays in the class, because |

468 | it could be good for more complicated polynomials. |

469 | But perhaps I will throw it out later... |

470 | |

471 | performance: |

472 | #sage: time chebyshev_T(5,x) #maxima |

473 | #CPU times: user 0.00 s, sys: 0.00 s, total: 0.00 s |

474 | #Wall time: 0.16 s |

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

476 | |

477 | #sage: time chebyshev_T(5,x) #clenshaw |

478 | #CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s |

479 | #Wall time: 0.01 s |

480 | #16*x^5 - 20*x^3 + 5*x |

481 | |

482 | #time chebyshev_T(50,x) |

483 | #CPU times: user 0.03 s, sys: 0.00 s, total: 0.03 s |

484 | #Wall time: 0.20 s |

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

486 | |

487 | #time chebyshev_T(100,x); |

488 | #CPU times: user 0.12 s, sys: 0.00 s, total: 0.12 s |

489 | #Wall time: 0.30 s |

490 | |

491 | EXAMPLES:: |

492 | sage: chebyshev_T(5,x) |

493 | 16*x^5 - 20*x^3 + 5*x |

494 | sage: var('n') |

495 | n |

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

497 | (-1)^n |

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

499 | 0 |

500 | sage: chebyshev_T(7.5,x) |

501 | Traceback (most recent call last): |

502 | ... |

503 | TypeError: Order is not an Integer! |

504 | |

505 | |

506 | |

507 | |

508 | """ |

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

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

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

512 | |

513 | if args[0] < 0: |

514 | return 0 #raise TypeError("Order is not valid!") |

515 | |

516 | try: |

517 | return self._eval_special_values_(*args) |

518 | except ValueError: |

519 | pass |

520 | |

521 | |

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

523 | try: |

524 | return self._clenshaw_method_(*args) |

525 | except NotImplementedError: |

526 | pass |

527 | |

528 | _init() |

529 | try: |

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

531 | except TypeError: |

532 | return None |

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

534 | return None |

535 | else: |

536 | return s.sage() |

537 | |

538 | |

539 | #TODO: numpy_eval with help of the new scipy version!!!! |

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

541 | #of ortho polys |

542 | |

543 | |

544 | class Func_chebyshev_T(OrthogonalPolynomial): |

545 | |

546 | """ |

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

548 | `n>-1`. |

549 | |

550 | REFERENCE: |

551 | |

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

553 | |

554 | EXAMPLES:: |

555 | |

556 | sage: chebyshev_T(3,x) |

557 | 4*x^3 - 3*x |

558 | sage: var('k') |

559 | k |

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

561 | sage: test |

562 | chebyshev_T(k, x) |

563 | """ |

564 | |

565 | def __init__(self): |

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

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

568 | |

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

570 | """ |

571 | Values known for special values of x. |

572 | For details see A.S. 22.4 (p. 777) |

573 | |

574 | EXAMPLES: |

575 | |

576 | sage: var('n') |

577 | n |

578 | sage: chebyshev_T(n,1) |

579 | 1 |

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

581 | (-1)^n |

582 | """ |

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

584 | return 1 |

585 | |

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

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

588 | |

589 | if (args[-1] == 0): |

590 | return (1+(-1)**args[0])*(-1)**(args[0]/2)/2 |

591 | |

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

593 | |

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

595 | """ |

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

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

598 | See A.S. 227 (p. 782) for details for the recurions |

599 | """ |

600 | |

601 | k = args[0] |

602 | x = args[1] |

603 | |

604 | if k == 0: |

605 | return 1 |

606 | elif k == 1: |

607 | return x |

608 | else: |

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

610 | #use these explicit formulas instead! |

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

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

613 | #elif 1 < x: |

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

615 | #else: # x < -1 |

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

617 | |

618 | help1 = 1 |

619 | help2 = x |

620 | if is_Expression(x): |

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

622 | help3 = expression.Expression.expand(2*x*help2 - help1) |

623 | help1 = help2 |

624 | help2 = help3 |

625 | else: |

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

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

628 | help1 = help2 |

629 | help2 = help3 |

630 | |

631 | return help3 |

632 | |

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

634 | """ |

635 | Evals chebyshev_T numerically with mpmath. |

636 | EXAMPLES:: |

637 | sage: chebyshev_T(10,3).n(75) |

638 | 2.261953700000000000000e7 |

639 | """ |

640 | parent = kwds['parent'] |

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

642 | |

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

644 | """ |

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

646 | of the second kind chebyshev_U |

647 | EXAMPLES:: |

648 | sage: var('k') |

649 | k |

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

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

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

653 | 12*x^2 - 3 |

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

655 | Traceback (most recent call last): |

656 | ... |

657 | NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be... |

658 | |

659 | """ |

660 | diff_param = kwds['diff_param'] |

661 | if diff_param == 0: |

662 | raise NotImplementedError( |

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

664 | else: |

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

666 | |

667 | chebyshev_T = Func_chebyshev_T() |

668 | |

669 | class Func_chebyshev_U(OrthogonalPolynomial): |

670 | |

671 | """ |

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

673 | |

674 | REFERENCE: |

675 | |

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

677 | |

678 | EXAMPLES:: |

679 | |

680 | sage: chebyshev_U(3,x) |

681 | 8*x^3 - 4*x |

682 | """ |

683 | |

684 | def __init__(self): |

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

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

687 | |

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

689 | """ |

690 | Clenshaw method for chebyshev_U (means use the recursion...) |

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

692 | See A.S. 227 (p. 782) for details for the recurions |

693 | """ |

694 | |

695 | k = args[0] |

696 | x = args[1] |

697 | |

698 | if k == 0: |

699 | return 1 |

700 | elif k == 1: |

701 | return 2*x |

702 | else: |

703 | help1 = 1 |

704 | help2 = 2*x |

705 | if is_Expression(x): |

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

707 | help3 = expression.Expression.expand(2*x*help2 - help1) |

708 | help1 = help2 |

709 | help2 = help3 |

710 | else: |

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

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

713 | help1 = help2 |

714 | help2 = help3 |

715 | |

716 | return help3 |

717 | |

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

719 | """ |

720 | Values known at the boundary. |

721 | EXAMPLES: |

722 | |

723 | sage: var('n') |

724 | n |

725 | sage: chebyshev_U(n,1) |

726 | n + 1 |

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

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

729 | """ |

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

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

732 | |

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

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

735 | |

736 | if (args[-1] == 0): |

737 | return (1+(-1)**args[0])*(-1)**(args[0]/2)/2 |

738 | |

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

740 | |

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

742 | """ |

743 | Evals chebyshev_U numerically with mpmath. |

744 | EXAMPLES:: |

745 | sage: chebyshev_U(10,3).n(75) |

746 | 4.661117900000000000000e7 |

747 | """ |

748 | parent = kwds['parent'] |

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

750 | |

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

752 | """ |

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

754 | Polynomials of the first and second kind |

755 | |

756 | EXAMPLES:: |

757 | sage: var('k') |

758 | k |

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

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

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

762 | 24*x^2 - 4 |

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

764 | Traceback (most recent call last): |

765 | ... |

766 | NotImplementedError: Derivative w.r.t. to the index is not supported, yet, and perhaps never will be... |

767 | |

768 | """ |

769 | diff_param = kwds['diff_param'] |

770 | if diff_param == 0: |

771 | raise NotImplementedError( |

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

773 | else: |

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

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

776 | |

777 | chebyshev_U = Func_chebyshev_U() |

778 | |

779 | class Func_legendre_P(OrthogonalPolynomial): |

780 | """ |

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

782 | `n >= -1`. |

783 | |

784 | REFERENCE: |

785 | |

786 | - AS 22.5.35 page 779. |

787 | |

788 | EXAMPLES:: |

789 | |

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

791 | sage: legendre_P(2,t) |

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

793 | sage: legendre_P(3, 1.1) |

794 | 1.67750000000000 |

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

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

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

798 | Traceback (most recent call last): |

799 | ... |

800 | TypeError: cannot coerce arguments: no canonical coercion from Full MatrixSpace of 2 by 2 dense matrices over Integer Ring to Symbolic Ring |

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

802 | 8 |

803 | """ |

804 | def __init__(self): |

805 | OrthogonalPolynomial.__init__(self,"legendre_P","legendre_p",nargs = 2,conversions =dict(maxima='legendre_p',mathematica='LegendreP')) |

806 | |

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

808 | """ |

809 | Clenshaw method for legendre_P (means use the recursion...) |

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

811 | See A.S. 227 (p. 782) for details for the recurions. |

812 | Warning: The clanshaw method for the Legendre Polynomials |

813 | should only used for exact data types, when high orders are |

814 | used, due to weak instabilities of the recursion! |

815 | """ |

816 | k = args[0] |

817 | x = args[-1] |

818 | |

819 | if k == 0: |

820 | return 1 |

821 | elif k == 1: |

822 | return x |

823 | else: |

824 | help1 = 1 |

825 | help2 = x |

826 | if is_Expression(x): |

827 | for j in xrange(1,k): |

828 | help3 = (2*j+1)*x*help2 - j*help1 |

829 | help3 = expression.Expression.expand(help3/(j+1)) |

830 | help1 = help2 |

831 | help2 = help3 |

832 | else: |

833 | for j in xrange(1,k): |

834 | help3 = (2*j+1)*x*help2 - j*help1 |

835 | help3 = help3/(j+1) |

836 | help1 = help2 |

837 | help2 = help3 |

838 | |

839 | return help3 |

840 | |

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

842 | """ |

843 | Special values known. |

844 | EXAMPLES: |

845 | |

846 | sage: var('n') |

847 | n |

848 | sage: legendre_P(n,1) |

849 | 1 |

850 | sage: legendre_P(n,-1) |

851 | (-1)^n |

852 | """ |

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

854 | return 1 |

855 | |

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

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

858 | |

859 | if (args[-1] == 0): |

860 | try: |

861 | return (1+(-1)**args[0])/2*binomial(args[0],args[0]/2)/4**(args[0]/2) |

862 | except TypeError: |

863 | pass |

864 | |

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

866 | |

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

868 | """ |

869 | Evals Legendre_P numerically with mpmath. |

870 | EXAMPLES:: |

871 | sage: legendre_P(10,5).n(53) |

872 | 1.60047267700000e9 |

873 | """ |

874 | parent = kwds['parent'] |

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

876 | |

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

878 | """return the derivative of legendre_P in |

879 | form of the Legendre Polynomial. |

880 | EXAMPLES:: |

881 | derivative(legendre_P(n,x),x) |

882 | -(n*x*legendre_P(n, x) - n*legendre_P(n - 1, x))/(x - 1)^2 |

883 | derivative(legendre_P(3,x),x) |

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

885 | """ |

886 | diff_param = kwds['diff_param'] |

887 | if diff_param == 0: |

888 | raise NotImplementedError( |

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

890 | else: |

891 | return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*legendre_P(args[0],args[-1]))/(1-args[-1])**2 |

892 | |

893 | |

894 | legendre_P = Func_legendre_P() |

895 | |

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

897 | """ |

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

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

900 | |

901 | REFERENCE: |

902 | |

903 | - table on page 789 in AS. |

904 | |

905 | EXAMPLES:: |

906 | |

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

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

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

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

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

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

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

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

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

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

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

918 | """ |

919 | from sage.functions.all import sqrt |

920 | _init() |

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

922 | |

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

924 | r""" |

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

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

927 | |

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

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

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

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

932 | functions. |

933 | |

934 | REFERENCE: |

935 | |

936 | - Gradshteyn and Ryzhik 8.706 page 1000. |

937 | |

938 | EXAMPLES:: |

939 | |

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

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

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

943 | sage: gen_legendre_P(2, 0, t) - legendre_P(2, t) |

944 | 0 |

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

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

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

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

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

950 | -16695*sqrt(2) |

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

952 | -583.562373654533*I |

953 | """ |

954 | from sage.functions.all import sqrt |

955 | _init() |

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

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

958 | else: |

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

960 | |

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

962 | """ |

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

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

965 | |

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

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

968 | 1. |

969 | |

970 | EXAMPLES:: |

971 | |

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

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

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

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

976 | 0 |

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

978 | 2.49185259170895 |

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

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

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

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

983 | """ |

984 | from sage.functions.all import sqrt |

985 | if m <= n: |

986 | _init() |

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

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

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

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

991 | else: |

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

993 | if m == n + 1: |

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

995 | else: |

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

997 | else: |

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

999 | |

1000 | def hermite(n,x): |

1001 | """ |

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

1003 | |

1004 | REFERENCE: |

1005 | |

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

1007 | |

1008 | EXAMPLES:: |

1009 | |

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

1011 | sage: hermite(2,x) |

1012 | 4*x^2 - 2 |

1013 | sage: hermite(3,x) |

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

1015 | sage: hermite(3,2) |

1016 | 40 |

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

1018 | sage: hermite(3,y) |

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

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

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

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

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

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

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

1026 | """ |

1027 | _init() |

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

1029 | |

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

1031 | r""" |

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

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

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

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

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

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

1038 | |

1039 | REFERENCE: |

1040 | |

1041 | - table on page 789 in AS. |

1042 | |

1043 | EXAMPLES:: |

1044 | |

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

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

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

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

1049 | 5.009999999999998 |

1050 | """ |

1051 | _init() |

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

1053 | |

1054 | def laguerre(n,x): |

1055 | """ |

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

1057 | |

1058 | REFERENCE: |

1059 | |

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

1061 | |

1062 | EXAMPLES:: |

1063 | |

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

1065 | sage: laguerre(2,x) |

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

1067 | sage: laguerre(3,x) |

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

1069 | sage: laguerre(2,2) |

1070 | -1 |

1071 | """ |

1072 | _init() |

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

1074 | |

1075 | |

1076 | def legendre_Q(n,x): |

1077 | """ |

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

1079 | `n>-1`. |

1080 | |

1081 | Computed using Maxima. |

1082 | |

1083 | EXAMPLES:: |

1084 | |

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

1086 | sage: legendre_Q(2, t) |

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

1088 | sage: legendre_Q(3, 0.5) |

1089 | -0.198654771479482 |

1090 | sage: legendre_Q(4, 2) |

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

1092 | sage: legendre_Q(4, 2.0) |

1093 | 0.00116107583162324 + 86.9828465962674*I |

1094 | """ |

1095 | _init() |

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

1097 | |

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

1099 | """ |

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

1101 | `n > -1`. |

1102 | |

1103 | Computed using Maxima. |

1104 | |

1105 | REFERENCE: |

1106 | |

1107 | - AS 22.5.27 |

1108 | |

1109 | EXAMPLES:: |

1110 | |

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

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

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

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

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

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

1117 | 2*x |

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

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

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

1121 | """ |

1122 | _init() |

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

1124 | |

1125 | gegenbauer = ultraspherical |