# Ticket #9706: orthogonal_polys.4.py

File orthogonal_polys.4.py, 63.7 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, is_inexact |

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, gamma, factorial, abs, binomial |

330 | from sage.functions.trig import sin |

331 | from sage.functions.log import ln |

332 | import sage.symbolic.expression as expression |

333 | from sage.structure.parent import Parent |

334 | from sage.structure.coerce import parent |

335 | |

336 | _done = False |

337 | def _init(): |

338 | """ |

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

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

341 | file should call this function first. |

342 | |

343 | TEST: |

344 | |

345 | The global starts ``False``:: |

346 | |

347 | sage: sage.functions.orthogonal_polys._done |

348 | False |

349 | |

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

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

352 | |

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

354 | sage: chebyshev_T(2,x) |

355 | 2*x^2 - 1 |

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

357 | False |

358 | sage: legendre_Q(1,0.1) |

359 | -0.989966465226892 |

360 | sage: sage.functions.orthogonal_polys._done |

361 | True |

362 | |

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

364 | the representation of the function is different from |

365 | its actual doctest, where a polynomial indeterminate |

366 | ``x`` is used. |

367 | """ |

368 | global _done |

369 | if _done: |

370 | return |

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

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

373 | # instead of just discarding this info! |

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

375 | _done = True |

376 | |

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

378 | |

379 | class OrthogonalPolynomial(BuiltinFunction): |

380 | """ |

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

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

383 | in maxima maxima_name has to be declared. |

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

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

386 | |

387 | """ |

388 | def __init__(self, name, nargs = 2, latex_name = None, conversions = {}): |

389 | try: |

390 | self._maxima_name = conversions['maxima'] |

391 | except KeyError: |

392 | self._maxima_name = None |

393 | |

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

395 | |

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

397 | """ |

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

399 | *args* in Maxima. |

400 | In fact these are thought to be the old wrappers for the orthogonal |

401 | polynomials. These are used when the other evaluation methods fail, |

402 | or are not fast enough. Experiments showed that for the symbolic |

403 | evaluation for larger n maxima is faster, but for small n simply use |

404 | of the recursion formulas is faster. A little switch does the trick... |

405 | |

406 | EXAMPLES:: |

407 | |

408 | sage: chebyshev_T(3,x) |

409 | 4*x^3 - 3*x |

410 | """ |

411 | return None |

412 | |

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

414 | """ |

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

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

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

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

419 | polynomialsin chebyshev form. |

420 | |

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

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

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

424 | #Wall time: 0.00 s |

425 | #49656746733678312490954442369580252421769338391329426325400124999 |

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

427 | #maxima |

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

429 | #Wall time: 0.05 s |

430 | #49656746733678312490954442369580252421769338391329426325400124999 |

431 | |

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

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

434 | #Wall time: 0.00 s |

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

436 | #maxima |

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

438 | #Wall time: 0.77 s |

439 | """ |

440 | raise NotImplementedError( |

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

442 | |

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

444 | """ |

445 | Evals the polynomial explicitly for special values. |

446 | EXAMPLES: |

447 | |

448 | sage: var('n') |

449 | n |

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

451 | (-1)^n |

452 | """ |

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

454 | |

455 | |

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

457 | """ |

458 | |

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

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

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

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

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

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

465 | |

466 | Update: There are cases where the pure Clanshaw algorithm seems faster |

467 | and sometimes maxima seems faster. e.g. for legendre_P clenshaw is faster |

468 | for legendre_Q maxima is faster. |

469 | |

470 | performance: |

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

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

473 | #Wall time: 0.16 s |

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

475 | |

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

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

478 | #Wall time: 0.01 s |

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

480 | |

481 | #time chebyshev_T(50,x) |

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

483 | #Wall time: 0.20 s |

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

485 | |

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

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

488 | #Wall time: 0.30 s |

489 | |

490 | EXAMPLES:: |

491 | sage: chebyshev_T(5,x) |

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

493 | sage: var('n') |

494 | n |

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

496 | (-1)^n |

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

498 | chebyshev_T(-7, x) |

499 | sage: chebyshev_T(3/2,x) |

500 | chebyshev_T(3/2, x) |

501 | |

502 | """ |

503 | |

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

505 | |

506 | if not is_Expression(args[-1]) and is_inexact(args[-1]): |

507 | try: |

508 | return self._evalf_(*args) |

509 | except AttributeError: |

510 | pass |

511 | except mpmath.NoConvergence: |

512 | print "Warning: mpmath returns NoConvergence!" |

513 | print "Switching to clenshaw_method, but it may not be stable!" |

514 | |

515 | |

516 | #A faster check would be nice... |

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

518 | if not is_Expression(args[-1]): |

519 | try: |

520 | return self._evalf_(*args) |

521 | except AttributeError: |

522 | pass |

523 | else: |

524 | return None |

525 | |

526 | if args[0] < 0: |

527 | return None |

528 | |

529 | |

530 | try: |

531 | return self._eval_special_values_(*args) |

532 | except ValueError: |

533 | pass |

534 | |

535 | |

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

537 | |

538 | try: |

539 | return self._clenshaw_method_(*args) |

540 | except NotImplementedError: |

541 | pass |

542 | |

543 | if self._maxima_name is None: |

544 | return None |

545 | else: |

546 | _init() |

547 | try: |

548 | #s = maxima(self._maxima_init_evaled_(*args)) |

549 | #This above is very inefficient! The older |

550 | #methods were much faster... |

551 | return self._maxima_init_evaled_(*args) |

552 | except TypeError: |

553 | return None |

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

555 | return None |

556 | else: |

557 | return s.sage() |

558 | |

559 | |

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

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

562 | #of ortho polys |

563 | |

564 | |

565 | class Func_chebyshev_T(OrthogonalPolynomial): |

566 | |

567 | """ |

568 | Class for the Chebyshev polynomial of the first kind. |

569 | |

570 | REFERENCE: |

571 | |

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

573 | |

574 | EXAMPLES:: |

575 | |

576 | sage: chebyshev_T(3,x) |

577 | 4*x^3 - 3*x |

578 | sage: var('k') |

579 | k |

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

581 | sage: test |

582 | chebyshev_T(k, x) |

583 | """ |

584 | |

585 | def __init__(self): |

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

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

588 | |

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

590 | """ |

591 | Values known for special values of x. |

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

593 | |

594 | EXAMPLES: |

595 | |

596 | sage: var('n') |

597 | n |

598 | sage: chebyshev_T(n,1) |

599 | 1 |

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

601 | (-1)^n |

602 | """ |

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

604 | return 1 |

605 | |

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

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

608 | |

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

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

611 | |

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

613 | |

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

615 | """ |

616 | Evals chebyshev_T |

617 | numerically with mpmath. |

618 | EXAMPLES:: |

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

620 | 2.261953700000000000000e7 |

621 | """ |

622 | try: |

623 | step_parent = kwds['parent'] |

624 | except KeyError: |

625 | step_parent = parent(args[-1]) |

626 | |

627 | try: |

628 | precision = step_parent.prec() |

629 | except AttributeError: |

630 | precision = mpmath.mp.prec |

631 | |

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

633 | |

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

635 | n = args[0] |

636 | x = args[1] |

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

638 | |

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

640 | """ |

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

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

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

644 | """ |

645 | |

646 | k = args[0] |

647 | x = args[1] |

648 | |

649 | if k == 0: |

650 | return 1 |

651 | elif k == 1: |

652 | return x |

653 | else: |

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

655 | #use these explicit formulas instead! |

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

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

658 | #elif 1 < x: |

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

660 | #else: # x < -1 |

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

662 | |

663 | help1 = 1 |

664 | help2 = x |

665 | if is_Expression(x): |

666 | #raise NotImplementedError |

667 | help3 = 0 |

668 | for j in xrange(0,floor(k/2)+1): |

669 | help3 = help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j-1)/factorial(j)/factorial(k-2*j) |

670 | help3 = help3*k/2 |

671 | else: |

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

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

674 | help1 = help2 |

675 | help2 = help3 |

676 | |

677 | return help3 |

678 | |

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

680 | """ |

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

682 | of the second kind chebyshev_U |

683 | EXAMPLES:: |

684 | sage: var('k') |

685 | k |

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

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

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

689 | 12*x^2 - 3 |

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

691 | Traceback (most recent call last): |

692 | ... |

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

694 | |

695 | """ |

696 | diff_param = kwds['diff_param'] |

697 | if diff_param == 0: |

698 | raise NotImplementedError( |

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

700 | else: |

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

702 | |

703 | chebyshev_T = Func_chebyshev_T() |

704 | |

705 | class Func_chebyshev_U(OrthogonalPolynomial): |

706 | |

707 | """ |

708 | Class for the Chebyshev polynomial of the second kind. |

709 | |

710 | REFERENCE: |

711 | |

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

713 | |

714 | EXAMPLES:: |

715 | |

716 | sage: chebyshev_U(3,x) |

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

718 | """ |

719 | |

720 | def __init__(self): |

721 | OrthogonalPolynomial.__init__(self,"chebyshev_U", |

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

723 | |

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

725 | """ |

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

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

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

729 | """ |

730 | k = args[0] |

731 | x = args[1] |

732 | |

733 | if k == 0: |

734 | return 1 |

735 | elif k == 1: |

736 | return 2*x |

737 | else: |

738 | help1 = 1 |

739 | help2 = 2*x |

740 | if is_Expression(x): |

741 | #raise NotImplementedError("Maxima is faster here!") |

742 | help3 = 0 |

743 | for j in xrange(0,floor(k/2)+1): |

744 | help3 = help3 + (-1)**j*(2*x)**(k-2*j)*factorial(k-j)/factorial(j)/factorial(k-2*j) |

745 | #help3 = help3*k/2 |

746 | else: |

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

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

749 | help1 = help2 |

750 | help2 = help3 |

751 | |

752 | return help3 |

753 | |

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

755 | """ |

756 | Uses |

757 | EXAMPLES:: |

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

759 | sage: chebyshev_T(2,x) |

760 | 2*x^2 - 1 |

761 | """ |

762 | n = args[0] |

763 | x = args[1] |

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

765 | |

766 | |

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

768 | """ |

769 | Evals chebyshev_U |

770 | numerically with mpmath. |

771 | EXAMPLES:: |

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

773 | 4.661117900000000000000e7 |

774 | """ |

775 | try: |

776 | step_parent = kwds['parent'] |

777 | except KeyError: |

778 | step_parent = parent(args[-1]) |

779 | |

780 | try: |

781 | precision = step_parent.prec() |

782 | except AttributeError: |

783 | precision = mpmath.mp.prec |

784 | |

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

786 | |

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

788 | """ |

789 | Special values known. A.S. 22.4 (p.777). |

790 | EXAMPLES: |

791 | |

792 | sage: var('n') |

793 | n |

794 | sage: chebyshev_U(n,1) |

795 | n + 1 |

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

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

798 | """ |

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

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

801 | |

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

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

804 | |

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

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

807 | |

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

809 | |

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

811 | """ |

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

813 | Polynomials of the first and second kind |

814 | |

815 | EXAMPLES:: |

816 | sage: var('k') |

817 | k |

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

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

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

821 | 24*x^2 - 4 |

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

823 | Traceback (most recent call last): |

824 | ... |

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

826 | |

827 | """ |

828 | diff_param = kwds['diff_param'] |

829 | if diff_param == 0: |

830 | raise NotImplementedError( |

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

832 | else: |

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

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

835 | |

836 | chebyshev_U = Func_chebyshev_U() |

837 | |

838 | class Func_gegenbauer(OrthogonalPolynomial): |

839 | """ |

840 | Returns the ultraspherical (or Gegenbauer) polynomial. |

841 | |

842 | REFERENCE: |

843 | |

844 | - AS 22.5.27 |

845 | |

846 | EXAMPLES:: |

847 | |

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

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

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

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

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

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

854 | 2*x |

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

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

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

858 | """ |

859 | def __init__(self): |

860 | OrthogonalPolynomial.__init__(self,"gegenbauer",nargs = 3,conversions =dict(maxima='ultraspherical',mathematica='GegenbauerC')) |

861 | |

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

863 | """ |

864 | Clenshaw method for gegenbauer poly (means use the recursion...) |

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

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

867 | """ |

868 | k = args[0] |

869 | x = args[-1] |

870 | alpha = args[1] |

871 | |

872 | if is_Expression(alpha) or abs(alpha) > numpy.finfo(float).eps: |

873 | alpha_zero = False |

874 | else: |

875 | alpha_zero = True |

876 | |

877 | if k == 0: |

878 | return 1 |

879 | elif k == 1: |

880 | if not alpha_zero: |

881 | return 2*x*alpha |

882 | else: |

883 | return 2*x # It seems that maxima evals this the wrong way! (see A.S. 22.4 (p.777)) |

884 | else: |

885 | help1 = 1 |

886 | if alpha_zero: |

887 | help2 = 2*x |

888 | else: |

889 | help2 = 2*x*alpha |

890 | |

891 | help3 = 0 |

892 | if is_Expression(x): |

893 | for j in xrange(0,floor(k/2)+1): |

894 | help3 = help3 + (-1)**j*gamma(alpha+k-j)/factorial(j)/factorial(k-2*j)*(2*x)**(k-2*j) |

895 | else: |

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

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

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

899 | help1 = help2 |

900 | help2 = help3 |

901 | |

902 | return help3 |

903 | |

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

905 | """ |

906 | Uses |

907 | EXAMPLES:: |

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

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

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

911 | """ |

912 | n = args[0] |

913 | a = args[1] |

914 | x = args[2] |

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

916 | |

917 | |

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

919 | """ |

920 | Evals gegenbauer |

921 | numerically with mpmath. |

922 | EXAMPLES:: |

923 | sage: gegenbauer(10,2,3.).n(54) |

924 | 5.25360702000000e8 |

925 | """ |

926 | try: |

927 | step_parent = kwds['parent'] |

928 | except KeyError: |

929 | step_parent = parent(args[-1]) |

930 | |

931 | try: |

932 | precision = step_parent.prec() |

933 | except AttributeError: |

934 | precision = mpmath.mp.prec |

935 | |

936 | return mpmath.call(mpmath.gegenbauer,args[0],args[1],args[-1],prec = precision) |

937 | |

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

939 | """ |

940 | Special values known. A.S. 22.4 (p.777) |

941 | EXAMPLES: |

942 | |

943 | sage: var('n a') |

944 | (n, a) |

945 | sage: gegenbauer(n,1/2,x) |

946 | legendre_P(n, x) |

947 | sage: gegenbauer(n,0,x) |

948 | 1/2*n*chebyshev_T(n, x) |

949 | sage: gegenbauer(n,1,x) |

950 | chebyshev_U(n, x) |

951 | sage: gegenbauer(n,a,1) |

952 | binomial(2*a + n - 1, n) |

953 | sage: gegenbauer(n,a,-1) |

954 | (-1)^n*binomial(2*a + n - 1, n) |

955 | sage: gegenbauer(n,a,0) |

956 | 1/2*((-1)^n + 1)*(-1)^(1/2*n)*gamma(a + 1/2*n)/(gamma(a)*gamma(1/2*n)) |

957 | """ |

958 | if args[1] == 0 and args[0] != 0: |

959 | return args[0]*chebyshev_T(args[0],args[-1])/2 |

960 | |

961 | if args[1] == 0.5: |

962 | return legendre_P(args[0],args[-1]) |

963 | |

964 | if args[1] == 1: |

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

966 | |

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

968 | return binomial(args[0] + 2*args[1] - 1,args[0]) |

969 | |

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

971 | return (-1)**args[0]*binomial(args[0] + 2*args[1] - 1,args[0]) |

972 | |

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

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

975 | |

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

977 | |

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

979 | """ |

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

981 | Polynomials of the first and second kind |

982 | |

983 | EXAMPLES:: |

984 | sage: var('k a') |

985 | (k, a) |

986 | sage: derivative(gegenbauer(k,a,x),x) |

987 | (k*x*gegenbauer(k, a, x) - (2*a + k - 1)*gegenbauer(k - 1, a, x))/(x^2 - 1) |

988 | sage: derivative(gegenbauer(4,3,x),x) |

989 | 1920*x^3 - 480*x |

990 | sage: derivative(gegenbauer(k,a,x),a) |

991 | Traceback (most recent call last): |

992 | ... |

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

994 | |

995 | """ |

996 | diff_param = kwds['diff_param'] |

997 | if diff_param in [0,1]: |

998 | raise NotImplementedError( |

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

1000 | else: |

1001 | return (-args[0]*args[-1]*gegenbauer(args[0],args[1],args[2])+ (args[0] + 2*args[1]-1)*gegenbauer(args[0]-1,args[1],args[2]))/(1-args[-1]**2) |

1002 | |

1003 | gegenbauer = Func_gegenbauer() |

1004 | |

1005 | class Func_gen_laguerre(OrthogonalPolynomial): |

1006 | """ |

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

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

1009 | |

1010 | REFERENCE: |

1011 | |

1012 | - table on page 789 in AS. |

1013 | |

1014 | EXAMPLES:: |

1015 | |

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

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

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

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

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

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

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

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

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

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

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

1027 | """ |

1028 | def __init__(self): |

1029 | OrthogonalPolynomial.__init__(self,"gen_laguerre",nargs = 3, |

1030 | conversions =dict(maxima='gen_laguerre',mathematica='LaguerreL')) |

1031 | |

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

1033 | """ |

1034 | Clenshaw method for gen_laguerre polynomial (means use the recursion...) |

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

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

1037 | Warning: The clanshaw method for the laguerre Polynomials |

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

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

1040 | """ |

1041 | k = args[0] |

1042 | a = args[1] |

1043 | x = args[-1] |

1044 | |

1045 | if k == 0: |

1046 | return 1 |

1047 | elif k == 1: |

1048 | return 1-x + a |

1049 | else: |

1050 | help1 = 1 |

1051 | help2 = 1-x + a |

1052 | if is_Expression(x): |

1053 | help3 = 0 |

1054 | for j in xrange(0,k+1): |

1055 | help3 = help3 + (-1)**j*binomial(k+a,k-j)/factorial(j)*x**j |

1056 | else: |

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

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

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

1060 | help1 = help2 |

1061 | help2 = help3 |

1062 | |

1063 | return help3 |

1064 | |

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

1066 | n = args[0] |

1067 | a = args[1] |

1068 | x = args[-1] |

1069 | |

1070 | from sage.functions.all import sqrt |

1071 | _init() |

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

1073 | |

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

1075 | """ |

1076 | Evals laguerre polynomial |

1077 | numerically with mpmath. |

1078 | EXAMPLES:: |

1079 | sage: laguerre(3,5.).n(53) |

1080 | 2.66666666666667 |

1081 | """ |

1082 | try: |

1083 | step_parent = kwds['parent'] |

1084 | except KeyError: |

1085 | step_parent = parent(args[-1]) |

1086 | |

1087 | try: |

1088 | precision = step_parent.prec() |

1089 | except AttributeError: |

1090 | precision = mpmath.mp.prec |

1091 | |

1092 | return mpmath.call(mpmath.laguerre,args[0],args[1],args[-1],prec = precision) |

1093 | |

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

1095 | """ |

1096 | Special values known. |

1097 | EXAMPLES: |

1098 | |

1099 | sage: var('n') |

1100 | n |

1101 | sage: laguerre(n,0) |

1102 | 1 |

1103 | """ |

1104 | |

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

1106 | try: |

1107 | return binomial(args[0]+args[1],args[0]) |

1108 | except TypeError: |

1109 | pass |

1110 | |

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

1112 | |

1113 | |

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

1115 | """return the derivative of laguerre in |

1116 | form of the Laguerre Polynomial. |

1117 | EXAMPLES:: |

1118 | sage: n = var('n') |

1119 | sage: derivative(laguerre(3,x),x) |

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

1121 | sage: derivative(laguerre(n,x),x) |

1122 | -(n*laguerre(n - 1, x) - n*laguerre(n, x))/x |

1123 | """ |

1124 | diff_param = kwds['diff_param'] |

1125 | if diff_param == 0: |

1126 | raise NotImplementedError( |

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

1128 | else: |

1129 | return (args[0]*laguerre(args[0],args[1],args[-1])-(args[0]+args[1])*laguerre(args[0]-1,args[1],args[-1]))/args[-1] |

1130 | |

1131 | gen_laguerre = Func_gen_laguerre() |

1132 | |

1133 | class Func_hermite(OrthogonalPolynomial): |

1134 | """ |

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

1136 | |

1137 | REFERENCE: |

1138 | |

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

1140 | |

1141 | #sage: S.<y> = PolynomialRing(RR)<- Here a bug of the base class BuiltinFunction occours! |

1142 | |

1143 | EXAMPLES:: |

1144 | |

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

1146 | sage: hermite(2,x) |

1147 | 4*x^2 - 2 |

1148 | sage: hermite(3,x) |

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

1150 | sage: hermite(3,2) |

1151 | 40 |

1152 | sage: y = var('y') |

1153 | sage: hermite(3,y).polynomial(PolynomialRing(RR, 'y')) |

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

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

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

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

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

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

1160 | 64*w^3 - 24*w |

1161 | """ |

1162 | |

1163 | def __init__(self): |

1164 | OrthogonalPolynomial.__init__(self,"hermite",nargs = 2, |

1165 | conversions =dict(maxima='hermite',mathematica='HermiteH')) |

1166 | |

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

1168 | """ |

1169 | Clenshaw method for hermite polynomial (means use the recursion...) |

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

1171 | For the symbolic evaluation, maxima seems to be quite fast. |

1172 | The break even point between the recursion and Maxima is about |

1173 | n = 25 |

1174 | """ |

1175 | k = args[0] |

1176 | x = args[1] |

1177 | |

1178 | if k == 0: |

1179 | return 1 |

1180 | elif k == 1: |

1181 | return 2*x |

1182 | else: |

1183 | help1 = 1 |

1184 | help2 = 2*x |

1185 | if is_Expression(x): |

1186 | help3 = 0 |

1187 | for j in xrange(0,floor(k/2)+1): |

1188 | help3 = help3 + (-1)**j*(2*x)**(k-2*j)/factorial(j)/factorial(k-2*j) |

1189 | help3 = help3*factorial(k) |

1190 | else: |

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

1192 | help3 = 2*x*help2 - 2*j*help1 |

1193 | help1 = help2 |

1194 | help2 = help3 |

1195 | |

1196 | return help3 |

1197 | |

1198 | |

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

1200 | """ |

1201 | Evals hermite |

1202 | numerically with mpmath. |

1203 | EXAMPLES:: |

1204 | sage: hermite(3,2.).n(74) |

1205 | 40.0000000000000000000 |

1206 | """ |

1207 | try: |

1208 | step_parent = kwds['parent'] |

1209 | except KeyError: |

1210 | step_parent = parent(args[-1]) |

1211 | |

1212 | try: |

1213 | precision = step_parent.prec() |

1214 | except AttributeError: |

1215 | precision = mpmath.mp.prec |

1216 | |

1217 | return mpmath.call(mpmath.hermite,args[0],args[-1],prec = precision) |

1218 | |

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

1220 | """ |

1221 | Special values known. A.S. 22.4 (p.777) |

1222 | EXAMPLES: |

1223 | |

1224 | sage: var('n') |

1225 | n |

1226 | sage: hermite(n,0) |

1227 | ((-1)^n + 1)*(-1)^(1/2*n)*factorial(n)/gamma(1/2*n + 1) |

1228 | """ |

1229 | |

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

1231 | return (1+(-1)**args[0])*(-1)**(args[0]/2)*factorial(args[0])/gamma(args[0]/2+1) |

1232 | |

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

1234 | |

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

1236 | """ |

1237 | Old maxima method. |

1238 | """ |

1239 | n = args[0] |

1240 | x = args[1] |

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

1242 | |

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

1244 | """ |

1245 | Returns the derivative of the hermite polynomial in form of the chebyshev |

1246 | Polynomials of the first and second kind |

1247 | |

1248 | EXAMPLES:: |

1249 | sage: var('k') |

1250 | k |

1251 | sage: derivative(hermite(k,x),x) |

1252 | 2*k*hermite(k - 1, x) |

1253 | |

1254 | """ |

1255 | diff_param = kwds['diff_param'] |

1256 | if diff_param == 0: |

1257 | raise NotImplementedError( |

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

1259 | else: |

1260 | return 2*args[0]*hermite(args[0]-1,args[1]) |

1261 | |

1262 | hermite = Func_hermite() |

1263 | |

1264 | |

1265 | class Func_jacobi_P(OrthogonalPolynomial): |

1266 | r""" |

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

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

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

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

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

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

1273 | |

1274 | REFERENCE: |

1275 | |

1276 | - table on page 789 in AS. |

1277 | |

1278 | EXAMPLES:: |

1279 | |

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

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

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

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

1284 | 5.009999999999998 |

1285 | """ |

1286 | |

1287 | def __init__(self): |

1288 | OrthogonalPolynomial.__init__(self,"jacobi_P",nargs = 4, |

1289 | conversions =dict(maxima='jacobi_p',mathematica='JacobiP')) |

1290 | |

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

1292 | """ |

1293 | Clenshaw method for jacobi_P (means use the recursion, |

1294 | or sum formula) |

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

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

1297 | Warning: The clanshaw method for the Jacobi Polynomials |

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

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

1300 | """ |

1301 | k = args[0] |

1302 | x = args[-1] |

1303 | alpha = args[1] |

1304 | beta = args[2] |

1305 | |

1306 | if k == 0: |

1307 | return 1 |

1308 | elif k == 1: |

1309 | return (alpha-beta + (alpha+beta+2)*x)/2 |

1310 | else: |

1311 | |

1312 | if is_Expression(x) or is_Expression(alpha) or is_Expression(beta): |

1313 | #Here we use the sum formula of jacobi_P it seems this is rather |

1314 | #optimal for use. |

1315 | help1 = gamma(alpha+k+1)/factorial(k)/gamma(alpha+beta+k+1) |

1316 | help2 = 0 |

1317 | for j in xrange(0,k+1): |

1318 | help2 = help2 + binomial(k,j)*gamma(alpha+beta+k+j+1)/gamma(alpha+j+1)*((x-1)/2)**j |

1319 | return help1*help2 |

1320 | else: |

1321 | help1 = 1 |

1322 | help2 = (alpha-beta + (alpha+beta+2)*x)/2 |

1323 | |

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

1325 | a1 = 2*(j+1)*(j+alpha+beta+1)*(2*j+alpha+beta) |

1326 | a2 = (2*j+alpha+beta+1)*(alpha**2-beta**2) |

1327 | a3 = (2*j+alpha+beta)*(2*j+alpha+beta+1)*(2*j+alpha+beta+2) |

1328 | a4 = 2*(j+alpha)*(j+beta)*(2*j+alpha+beta+2) |

1329 | help3 = (a2+a3*x)*help2 - a4*help1 |

1330 | help3 = help3/a1 |

1331 | help1 = help2 |

1332 | help2 = help3 |

1333 | |

1334 | return help3 |

1335 | |

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

1337 | """ |

1338 | Old maxima method. |

1339 | """ |

1340 | n = args[0] |

1341 | a = args[1] |

1342 | b = args[2] |

1343 | x = args[3] |

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

1345 | |

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

1347 | """ |

1348 | Evals jacobi_P |

1349 | numerically with mpmath. |

1350 | EXAMPLES:: |

1351 | sage: jacobi_P(10,2,3,3).n(75) |

1352 | 1.322776620000000000000e8 |

1353 | """ |

1354 | try: |

1355 | step_parent = kwds['parent'] |

1356 | except KeyError: |

1357 | step_parent = parent(args[-1]) |

1358 | |

1359 | try: |

1360 | precision = step_parent.prec() |

1361 | except AttributeError: |

1362 | precision = mpmath.mp.prec |

1363 | |

1364 | return mpmath.call(mpmath.jacobi,args[0],args[1],args[2],args[-1],prec = precision) |

1365 | |

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

1367 | """ |

1368 | Special values known. A.S. 22.4 (p.777) |

1369 | EXAMPLES: |

1370 | |

1371 | sage: var('n') |

1372 | n |

1373 | sage: legendre_P(n,1) |

1374 | 1 |

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

1376 | (-1)^n |

1377 | """ |

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

1379 | return binomial(args[0]+args[1],args[0]) |

1380 | |

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

1382 | return (-1)**args[0]*binomial(args[0]+args[1],args[0]) |

1383 | |

1384 | if args[1] == 0 and args[2] == 0: |

1385 | return legendre_P(args[0],args[-1]) |

1386 | |

1387 | if args[1] == -0.5 and args[2] == -0.5: |

1388 | try: |

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

1390 | except TypeError: |

1391 | pass |

1392 | |

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

1394 | |

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

1396 | """ |

1397 | Returns the derivative of jacobi_P in form of jacobi_polynomials |

1398 | |

1399 | EXAMPLES:: |

1400 | sage: var('k a b') |

1401 | (k, a, b) |

1402 | sage: derivative(jacobi_P(k,a,b,x),x) |

1403 | -(2*(b + k)*(a + k)*jacobi_P(k - 1, a, b, x) - ((a + b + 2*k)*x - a + b)*k*jacobi_P(k, a, b, x))/((x^2 - 1)*(a + b + 2*k)) |

1404 | sage: derivative(jacobi_P(2,1,3,x),x) |

1405 | 14*x - 7/2 |

1406 | sage: derivative(jacobi_P(k,a,b,x),a) |

1407 | Traceback (most recent call last): |

1408 | ... |

1409 | |

1410 | |

1411 | """ |

1412 | diff_param = kwds['diff_param'] |

1413 | if diff_param in [0,1,2]: |

1414 | raise NotImplementedError( |

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

1416 | else: |

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

1418 | jacobi_P(args[0],args[1],args[2],args[-1])+ 2*(args[0]+args[1])* |

1419 | (args[0]+args[2])*jacobi_P(args[0]-1,args[1],args[2],args[-1]))/(2*args[0]+args[1]+args[2])/(1-args[-1]**2) |

1420 | |

1421 | |

1422 | jacobi_P = Func_jacobi_P() |

1423 | |

1424 | class Func_laguerre(OrthogonalPolynomial): |

1425 | """ |

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

1427 | |

1428 | REFERENCE: |

1429 | |

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

1431 | |

1432 | EXAMPLES:: |

1433 | |

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

1435 | sage: laguerre(2,x) |

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

1437 | sage: laguerre(3,x) |

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

1439 | sage: laguerre(2,2) |

1440 | -1 |

1441 | """ |

1442 | def __init__(self): |

1443 | OrthogonalPolynomial.__init__(self,"laguerre",nargs = 2, |

1444 | conversions =dict(maxima='laguerre',mathematica='LaguerreL')) |

1445 | |

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

1447 | """ |

1448 | Clenshaw method for laguerre polynomial (means use the recursion...) |

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

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

1451 | Warning: The clanshaw method for the laguerre Polynomials |

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

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

1454 | """ |

1455 | k = args[0] |

1456 | x = args[-1] |

1457 | |

1458 | if k == 0: |

1459 | return 1 |

1460 | elif k == 1: |

1461 | return 1-x |

1462 | else: |

1463 | help1 = 1 |

1464 | help2 = 1-x |

1465 | if is_Expression(x): |

1466 | help3 = 0 |

1467 | for j in xrange(0,k+1): |

1468 | help3 = help3 + (-1)**j*binomial(k,k-j)/factorial(j)*x**j |

1469 | else: |

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

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

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

1473 | help1 = help2 |

1474 | help2 = help3 |

1475 | |

1476 | return help3 |

1477 | |

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

1479 | n = args[0] |

1480 | x = args[1] |

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

1482 | |

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

1484 | """ |

1485 | Evals laguerre polynomial |

1486 | numerically with mpmath. |

1487 | EXAMPLES:: |

1488 | sage: laguerre(3,5.).n(53) |

1489 | 2.66666666666667 |

1490 | """ |

1491 | try: |

1492 | step_parent = kwds['parent'] |

1493 | except KeyError: |

1494 | step_parent = parent(args[-1]) |

1495 | |

1496 | try: |

1497 | precision = step_parent.prec() |

1498 | except AttributeError: |

1499 | precision = mpmath.mp.prec |

1500 | |

1501 | return mpmath.call(mpmath.laguerre,args[0],0,args[-1],prec = precision) |

1502 | |

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

1504 | """ |

1505 | Special values known. |

1506 | EXAMPLES: |

1507 | |

1508 | sage: var('n') |

1509 | n |

1510 | sage: laguerre(n,0) |

1511 | 1 |

1512 | """ |

1513 | |

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

1515 | try: |

1516 | return 1 |

1517 | except TypeError: |

1518 | pass |

1519 | |

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

1521 | |

1522 | |

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

1524 | """return the derivative of laguerre in |

1525 | form of the Laguerre Polynomial. |

1526 | EXAMPLES:: |

1527 | sage: n = var('n') |

1528 | sage: derivative(laguerre(3,x),x) |

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

1530 | sage: derivative(laguerre(n,x),x) |

1531 | -(n*laguerre(n - 1, x) - n*laguerre(n, x))/x |

1532 | """ |

1533 | diff_param = kwds['diff_param'] |

1534 | if diff_param == 0: |

1535 | raise NotImplementedError( |

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

1537 | else: |

1538 | return (args[0]*laguerre(args[0],args[-1])-args[0]*laguerre(args[0]-1,args[-1]))/args[1] |

1539 | |

1540 | laguerre = Func_laguerre() |

1541 | |

1542 | class Func_legendre_P(OrthogonalPolynomial): |

1543 | """ |

1544 | Returns the Legendre polynomial of the first kind.. |

1545 | |

1546 | REFERENCE: |

1547 | |

1548 | - AS 22.5.35 page 779. |

1549 | |

1550 | EXAMPLES:: |

1551 | |

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

1553 | sage: legendre_P(2,t) |

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

1555 | sage: legendre_P(3, 1.1) |

1556 | 1.67750000000000 |

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

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

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

1560 | Traceback (most recent call last): |

1561 | ... |

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

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

1564 | 8 |

1565 | """ |

1566 | def __init__(self): |

1567 | OrthogonalPolynomial.__init__(self,"legendre_P",nargs = 2, |

1568 | conversions =dict(maxima='legendre_p',mathematica='LegendreP')) |

1569 | |

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

1571 | """ |

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

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

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

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

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

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

1578 | """ |

1579 | k = args[0] |

1580 | x = args[-1] |

1581 | |

1582 | if k == 0: |

1583 | return 1 |

1584 | elif k == 1: |

1585 | return x |

1586 | else: |

1587 | help1 = 1 |

1588 | help2 = x |

1589 | if is_Expression(x): |

1590 | #raise NotImplementedError("Maxima is faster here...") |

1591 | help3 = 0 |

1592 | for j in xrange(0,floor(k/2)+1): |

1593 | help3 = help3 + (-1)**j*x**(k-2*j)*binomial(k,j)*binomial(2*k-2*j,k) |

1594 | help3 = help3/2**k |

1595 | |

1596 | else: |

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

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

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

1600 | help1 = help2 |

1601 | help2 = help3 |

1602 | |

1603 | return help3 |

1604 | |

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

1606 | n = args[0] |

1607 | x = args[1] |

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

1609 | |

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

1611 | """ |

1612 | Evals legendre_P |

1613 | numerically with mpmath. |

1614 | EXAMPLES:: |

1615 | sage: legendre_P(10,3).n(75) |

1616 | 8.097453000000000000000e6 |

1617 | """ |

1618 | try: |

1619 | step_parent = kwds['parent'] |

1620 | except KeyError: |

1621 | step_parent = parent(args[-1]) |

1622 | |

1623 | try: |

1624 | precision = step_parent.prec() |

1625 | except AttributeError: |

1626 | precision = mpmath.mp.prec |

1627 | |

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

1629 | |

1630 | |

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

1632 | """ |

1633 | Special values known. |

1634 | EXAMPLES: |

1635 | |

1636 | sage: var('n') |

1637 | n |

1638 | sage: legendre_P(n,1) |

1639 | 1 |

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

1641 | (-1)^n |

1642 | """ |

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

1644 | return 1 |

1645 | |

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

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

1648 | |

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

1650 | try: |

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

1652 | except TypeError: |

1653 | pass |

1654 | |

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

1656 | |

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

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

1659 | form of the Legendre Polynomial. |

1660 | EXAMPLES:: |

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

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

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

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

1665 | """ |

1666 | diff_param = kwds['diff_param'] |

1667 | if diff_param == 0: |

1668 | raise NotImplementedError( |

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

1670 | else: |

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

1672 | |

1673 | |

1674 | legendre_P = Func_legendre_P() |

1675 | |

1676 | class Func_legendre_Q(OrthogonalPolynomial): |

1677 | """ |

1678 | Return the chebyshev function of the second kind. |

1679 | |

1680 | REFERENCE: |

1681 | |

1682 | - A.S. 8 (p. 332) |

1683 | |

1684 | EXAMPLES:: |

1685 | |

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

1687 | sage: legendre_Q(2, t) |

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

1689 | sage: legendre_Q(3, 0.5) |

1690 | -0.198654771479482 |

1691 | sage: legendre_Q(4, 2) |

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

1693 | sage: legendre_Q(4, 2.0) |

1694 | 0.00116107583162324 + 86.9828465962674*I |

1695 | |

1696 | """ |

1697 | def __init__(self): |

1698 | OrthogonalPolynomial.__init__(self,"legendre_Q",nargs = 2, |

1699 | conversions =dict(maxima='legendre_q',mathematica='LegendreQ')) |

1700 | |

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

1702 | """ |

1703 | Maxima seems just fine for legendre Q. So we use it here! |

1704 | """ |

1705 | n = args[0] |

1706 | x = args[1] |

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

1708 | |

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

1710 | """ |

1711 | Clenshaw method for legendre_q (means use the recursion...) |

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

1713 | See A.S. 8.5.3 (p. 334) for details for the recurions. |

1714 | Warning: The clanshaw method for the Legendre fUNCTIONS |

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

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

1717 | """ |

1718 | raise NotImplementedError("Function not ready yet...") |

1719 | |

1720 | k = args[0] |

1721 | x = args[-1] |

1722 | |

1723 | if k == 0: |

1724 | return ln((1+x)/(1-x))/2 |

1725 | elif k == 1: |

1726 | return x/2*ln((1+x)/(1-x))-1 |

1727 | else: |

1728 | if is_Expression(x): |

1729 | raise NotImplementedError("Maxima works fine here!") |

1730 | #it seems that the old method just works fine here... |

1731 | #raise NotImplementedError("clenshaw does not work well...") |

1732 | else: |

1733 | help1 = ln((1+x)/(1-x))/2 |

1734 | help2 = x/2*ln((1+x)/(1-x))-1 |

1735 | |

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

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

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

1739 | help1 = help2 |

1740 | help2 = help3 |

1741 | |

1742 | return help3 |

1743 | |

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

1745 | """ |

1746 | Special values known. |

1747 | EXAMPLES: |

1748 | |

1749 | sage: var('n') |

1750 | n |

1751 | sage: legendre_P(n,1) |

1752 | 1 |

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

1754 | (-1)^n |

1755 | """ |

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

1757 | return NaN |

1758 | |

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

1760 | return NaN |

1761 | |

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

1763 | if is_Expression(args[0]): |

1764 | try: |

1765 | return -(SR.pi()**(QQ(1)/QQ(2)))/2*sin(SR.pi()/2*args[0])*gamma((args[0]+1)/2)*gamma(args[0]/2 + 1) |

1766 | except TypeError: |

1767 | pass |

1768 | else: |

1769 | return -(math.pi**(QQ(1)/QQ(2)))/2*sin(math.pi/2*args[0])*gamma((args[0]+1)/2)*gamma(args[0]/2. + 1) |

1770 | |

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

1772 | |

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

1774 | """ |

1775 | Evals legendre_Q |

1776 | numerically with mpmath. |

1777 | Function does not work yet... |

1778 | EXAMPLES:: |

1779 | sage: legendre_P(10,3).n(75) |

1780 | 8.097453000000000000000e6 |

1781 | """ |

1782 | raise AttributeError("mpmath function does not work correctly!") |

1783 | |

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

1785 | """return the derivative of legendre_Q in |

1786 | form of the Legendre Function. |

1787 | EXAMPLES:: |

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

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

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

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

1792 | """ |

1793 | diff_param = kwds['diff_param'] |

1794 | if diff_param == 0: |

1795 | raise NotImplementedError( |

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

1797 | else: |

1798 | return (gen_legendre_Q(args[0],1,args[-1]))/sqrt(1-args[-1]**2) |

1799 | |

1800 | legendre_Q = Func_legendre_Q() |

1801 | |

1802 | def legendre_Q_old(n,x): |

1803 | """ |

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

1805 | `n>-1`. (old version) |

1806 | |

1807 | Computed using Maxima. |

1808 | |

1809 | EXAMPLES:: |

1810 | |

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

1812 | sage: legendre_Q(2, t) |

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

1814 | sage: legendre_Q(3, 0.5) |

1815 | -0.198654771479482 |

1816 | sage: legendre_Q(4, 2) |

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

1818 | sage: legendre_Q(4, 2.0) |

1819 | 0.00116107583162324 + 86.9828465962674*I |

1820 | """ |

1821 | _init() |

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

1823 | |

1824 | |

1825 | |

1826 | def gen_laguerre_old(n,a,x): |

1827 | """ |

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

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

1830 | |

1831 | REFERENCE: |

1832 | |

1833 | - table on page 789 in AS. |

1834 | |

1835 | EXAMPLES:: |

1836 | |

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

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

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

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

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

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

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

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

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

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

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

1848 | """ |

1849 | from sage.functions.all import sqrt |

1850 | _init() |

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

1852 | |

1853 | class Func_gen_legendre_P(OrthogonalPolynomial): |

1854 | |

1855 | def __init__(self): |

1856 | OrthogonalPolynomial.__init__(self,"gen_legendre_P",nargs = 3, |

1857 | conversions =dict(maxima='assoc_legendre_p',mathematica='LegendreP')) |

1858 | |

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

1860 | n = args[0] |

1861 | m = args[1] |

1862 | x = args[2] |

1863 | if is_Expression(n) or is_Expression(m): |

1864 | return None |

1865 | |

1866 | from sage.functions.all import sqrt |

1867 | |

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

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

1870 | else: |

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

1872 | |

1873 | |

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

1875 | """return the derivative of gen_legendre_P in |

1876 | form of the Legendre Function. |

1877 | EXAMPLES:: |

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

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

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

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

1882 | """ |

1883 | diff_param = kwds['diff_param'] |

1884 | if diff_param == 0: |

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

1886 | else: |

1887 | return ((-args[1]*args[-1]*gen_legendre_P(args[0],args[1],args[-1]))+(args*[0]+args[1])*(args[0]-args[1]-1)*sqrt(args[0]^2-1)*gen_legendre_P(args[0],args[1]-1,args[-1]))/(args[0]**2-1) |

1888 | |

1889 | gen_legendre_P = Func_gen_legendre_P() |

1890 | |

1891 | class Func_gen_legendre_Q(OrthogonalPolynomial): |

1892 | |

1893 | def __init__(self): |

1894 | OrthogonalPolynomial.__init__(self,"gen_legendre_Q",nargs = 3, |

1895 | conversions =dict(maxima='assoc_legendre_q',mathematica='LegendreQ')) |

1896 | |

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

1898 | n = args[0] |

1899 | m = args[1] |

1900 | x = args[2] |

1901 | if is_Expression(n) or is_Expression(m): |

1902 | return None |

1903 | |

1904 | from sage.functions.all import sqrt |

1905 | |

1906 | if m <= n: |

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

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

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

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

1911 | else: |

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

1913 | if m == n + 1: |

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

1915 | else: |

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

1917 | else: |

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

1919 | |

1920 | |

1921 | #def _derivative_(self,*args,**kwds): |

1922 | # """return the derivative of gen_legendre_P in |

1923 | # form of the Legendre Function. |

1924 | # EXAMPLES:: |

1925 | # derivative(legendre_P(n,x),x) |

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

1927 | # derivative(legendre_P(3,x),x) |

1928 | # 15/2*x^2 - 3/2 |

1929 | # """ |

1930 | # diff_param = kwds['diff_param'] |

1931 | # if diff_param == 0: |

1932 | # raise NotImplementedError("Derivative w.r.t. to the index is not supported, yet, and perhaps never will be...") |

1933 | # else: |

1934 | # return ((-args[1]*args[-1]*gen_legendre_P(args[0],args[1],args[-1]))+(args*[0]+args[1])*(args[0]-args[1]-1)*sqrt(args[0]^2-1)*gen_legendre_P(args[0],args[1]-1,args[-1]))/(args[0]**2-1) |

1935 | |

1936 | gen_legendre_Q = Func_gen_legendre_Q() |

1937 | |

1938 | def gen_legendre_P_old(n,m,x): |

1939 | r""" |

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

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

1942 | |

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

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

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

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

1947 | functions. |

1948 | |

1949 | REFERENCE: |

1950 | |

1951 | - Gradshteyn and Ryzhik 8.706 page 1000. |

1952 | |

1953 | EXAMPLES:: |

1954 | |

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

1956 | sage: expand(gen_legendre_P(2, 0, t)) |

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

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

1959 | 0 |

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

1961 | -3/2*sqrt(-t^2 + 1)*(5*(t - 1)^2 + 10*t - 6) |

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

1963 | 15*sqrt(-t^2 + 1)*((t^2 - 1)*(7*(t - 1)^2 + 14*t - 8)*t - 6*(t^2 - 1)*t)/(t^2 - 1) |

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

1965 | -16695*sqrt(2) |

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

1967 | -583.562373654533*I |

1968 | """ |

1969 | from sage.functions.all import sqrt |

1970 | _init() |

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

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

1973 | else: |

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

1975 | |

1976 | def gen_legendre_Q_old(n,m,x): |

1977 | """ |

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

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

1980 | |

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

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

1983 | 1. |

1984 | |

1985 | EXAMPLES:: |

1986 | |

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

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

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

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

1991 | 0 |

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

1993 | 2.49185259170895 |

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

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

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

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

1998 | """ |

1999 | from sage.functions.all import sqrt |

2000 | if m <= n: |

2001 | _init() |

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

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

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

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

2006 | else: |

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

2008 | if m == n + 1: |

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

2010 | else: |

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

2012 | else: |

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

2014 | |

2015 | |

2016 | def jacobi_P_old(n,a,b,x): |

2017 | r""" |

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

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

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

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

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

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

2024 | |

2025 | REFERENCE: |

2026 | |

2027 | - table on page 789 in AS. |

2028 | |

2029 | EXAMPLES:: |

2030 | |

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

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

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

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

2035 | 5.009999999999998 |

2036 | """ |

2037 | _init() |

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

2039 | |

2040 | def laguerre_old(n,x): |

2041 | """ |

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

2043 | |

2044 | REFERENCE: |

2045 | |

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

2047 | |

2048 | EXAMPLES:: |

2049 | |

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

2051 | sage: laguerre(2,x) |

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

2053 | sage: laguerre(3,x) |

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

2055 | sage: laguerre(2,2) |

2056 | -1 |

2057 | """ |

2058 | _init() |

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

2060 | |

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

2062 | """ |

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

2064 | `n > -1`. |

2065 | |

2066 | Computed using Maxima. |

2067 | |

2068 | REFERENCE: |

2069 | |

2070 | - AS 22.5.27 |

2071 | |

2072 | EXAMPLES:: |

2073 | |

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

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

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

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

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

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

2080 | 2*x |

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

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

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

2084 | """ |

2085 | _init() |

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

2087 | |

2088 |