# Ticket #9706: orthogonal_polys.8.py

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

Line | |
---|---|

1 | # -*- coding: utf-8 -*- |

2 | r""" |

3 | Orthogonal Polynomials |

4 | |

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

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

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

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

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

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

11 | version information. |

12 | ----------------------------------------------------------------- |

13 | Update (2010): The orthogonal polynomials are now selfstanding classes. |

14 | Altough some of the use the old maxima wrappers, most of the evaluation |

15 | is done with the help of direct methods. Also the orthogonal polys |

16 | are handled as symbolic expressions, if the input is not known |

17 | |

18 | |

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

20 | to the differential equation |

21 | |

22 | .. math:: |

23 | |

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

25 | |

26 | |

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

28 | |

29 | .. math:: |

30 | |

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

32 | |

33 | |

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

35 | recurrence relation |

36 | |

37 | .. math:: |

38 | |

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

40 | |

41 | |

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

43 | recurrence relation |

44 | |

45 | .. math:: |

46 | |

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

48 | |

49 | |

50 | |

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

52 | relations |

53 | |

54 | .. math:: |

55 | |

56 | \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. |

57 | |

58 | |

59 | and |

60 | |

61 | |

62 | .. math:: |

63 | |

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

65 | |

66 | |

67 | |

68 | They are named after Pafnuty Chebyshev (alternative |

69 | transliterations: Tchebyshef or Tschebyscheff). |

70 | |

71 | - The Hermite polynomials are defined either by |

72 | |

73 | .. math:: |

74 | |

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

76 | |

77 | |

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

79 | |

80 | |

81 | .. math:: |

82 | |

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

84 | |

85 | |

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

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

88 | relation |

89 | |

90 | .. math:: |

91 | |

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

93 | |

94 | |

95 | |

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

97 | |

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

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

100 | |

101 | .. math:: |

102 | |

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

104 | |

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

106 | |

107 | .. math:: |

108 | |

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

110 | |

111 | and satisfy the orthogonality relation |

112 | |

113 | .. math:: |

114 | |

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

116 | |

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

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

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

120 | |

121 | The associated Legendre functions of the first kind |

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

123 | Legendre polynomials by |

124 | |

125 | .. math:: |

126 | |

127 | \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} |

128 | |

129 | |

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

131 | relation: |

132 | |

133 | .. math:: |

134 | |

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

136 | |

137 | |

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

139 | |

140 | The associated Legendre functions of the second kind |

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

142 | Legendre polynomials by |

143 | |

144 | |

145 | .. math:: |

146 | |

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

148 | |

149 | |

150 | |

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

152 | |

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

154 | |

155 | .. math:: |

156 | |

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

158 | |

159 | |

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

161 | |

162 | |

163 | .. math:: |

164 | |

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

166 | |

167 | and satisfy the orthogonality relation |

168 | |

169 | |

170 | .. math:: |

171 | |

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

173 | |

174 | |

175 | |

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

177 | Rodrigues formula: |

178 | |

179 | |

180 | .. math:: |

181 | |

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

183 | |

184 | |

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

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

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

188 | |

189 | They are named after Edmond Laguerre. |

190 | |

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

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

193 | is in fact finite: |

194 | |

195 | .. math:: |

196 | |

197 | 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) , |

198 | |

199 | |

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

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

202 | explicit expression |

203 | |

204 | |

205 | .. math:: |

206 | |

207 | 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 . |

208 | |

209 | |

210 | |

211 | They are named after Carl Jacobi. |

212 | |

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

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

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

216 | |

217 | |

218 | .. math:: |

219 | |

220 | 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). |

221 | |

222 | |

223 | They satisfy the orthogonality relation |

224 | |

225 | .. math:: |

226 | |

227 | \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)} , |

228 | |

229 | |

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

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

232 | |

233 | |

234 | .. math:: |

235 | |

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

237 | |

238 | |

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

240 | Abramowitz and Stegun p561) |

241 | |

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

243 | |

244 | |

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

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

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

248 | |

249 | .. math:: |

250 | |

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

252 | |

253 | |

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

255 | |

256 | .. math:: |

257 | |

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

259 | |

260 | |

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

262 | Patashnik in their book Concrete Mathematics. |

263 | |

264 | .. note:: |

265 | |

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

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

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

269 | |

270 | TODO: Implement associated Legendre polynomials and Zernike |

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

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

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

274 | |

275 | REFERENCES: |

276 | |

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

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

279 | |

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

281 | |

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

283 | |

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

285 | |

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

287 | |

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

289 | |

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

291 | |

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

293 | |

294 | |

295 | AUTHORS: |

296 | |

297 | - David Joyner (2006-06) |

298 | - Stefan Reiterer (2010-) |

299 | """ |

300 | |

301 | #***************************************************************************** |

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

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

304 | # 2010 Stefan Reiterer <stefan.reiterer@uni-graz.at> |

305 | # |

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

307 | # |

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

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

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

311 | # General Public License for more details. |

312 | # |

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

314 | # |

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

316 | #***************************************************************************** |

317 | |

318 | import copy |

319 | import sage.plot.plot |

320 | import sage.interfaces.all |

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

322 | from sage.rings.rational_field import RationalField |

323 | from sage.rings.real_mpfr import RealField |

324 | from sage.misc.sage_eval import sage_eval |

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

326 | from sage.symbolic.ring import SR |

327 | import sage.rings.commutative_ring as commutative_ring |

328 | import sage.rings.ring as ring |

329 | from sage.calculus.calculus import maxima |

330 | from sage.symbolic.function import BuiltinFunction, GinacFunction, is_inexact |

331 | from sage.symbolic.expression import is_Expression |

332 | import sage.functions.special |

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

334 | import numpy,scipy |

335 | import math |

336 | #import sage.libs.mpmath.all as mpmath #This line does evil things... |

337 | from sage.functions.other import floor, gamma, factorial, abs, binomial |

338 | from sage.functions.other import sqrt, conjugate |

339 | from sage.functions.trig import sin, cos |

340 | from sage.functions.log import ln |

341 | import sage.symbolic.expression as expression |

342 | from sage.structure.parent import Parent |

343 | from sage.structure.coerce import parent |

344 | |

345 | _done = False |

346 | def _init(): |

347 | """ |

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

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

350 | file should call this function first. |

351 | |

352 | TEST: |

353 | |

354 | The global starts ``False``:: |

355 | |

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

357 | False |

358 | |

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

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

361 | |

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

363 | sage: chebyshev_T(2,x) |

364 | 2*x^2 - 1 |

365 | sage: sage.functions.orthogonal_polys._done |

366 | False |

367 | sage: legendre_Q(1,x) |

368 | 1/2*x*log(-(x + 1)/(x - 1)) - 1 |

369 | sage: sage.functions.orthogonal_polys._done |

370 | True |

371 | |

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

373 | the representation of the function is different from |

374 | its actual doctest, where a polynomial indeterminate |

375 | ``x`` is used. |

376 | """ |

377 | global _done |

378 | if _done: |

379 | return |

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

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

382 | # instead of just discarding this info! |

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

384 | _done = True |

385 | |

386 | #load /home/maldun/sage/sage-4.5.2/devel/sage-ortho/sage/functions/orthogonal_polys.py |

387 | |

388 | class OrthogonalPolynomial(BuiltinFunction): |

389 | """ |

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

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

392 | in maxima maxima_name has to be declared. |

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

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

395 | |

396 | """ |

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

398 | try: |

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

400 | except KeyError: |

401 | self._maxima_name = None |

402 | |

403 | BuiltinFunction.__init__(self, name = name, |

404 | nargs = nargs, latex_name = latex_name, conversions = conversions) |

405 | |

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

407 | """ |

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

409 | *args* in Maxima. |

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

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

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

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

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

415 | |

416 | EXAMPLES:: |

417 | |

418 | sage: chebyshev_T(3,x) |

419 | 4*x^3 - 3*x |

420 | """ |

421 | return None |

422 | |

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

424 | """ |

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

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

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

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

429 | polynomialsin chebyshev form. |

430 | |

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

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

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

434 | #Wall time: 0.00 s |

435 | #49656746733678312490954442369580252421769338391329426325400124999 |

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

437 | #maxima |

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

439 | #Wall time: 0.05 s |

440 | #49656746733678312490954442369580252421769338391329426325400124999 |

441 | |

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

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

444 | #Wall time: 0.00 s |

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

446 | #maxima |

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

448 | #Wall time: 0.77 s |

449 | """ |

450 | raise NotImplementedError( |

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

452 | |

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

454 | """ |

455 | Evals the polynomial explicitly for special values. |

456 | EXAMPLES: |

457 | |

458 | sage: var('n') |

459 | n |

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

461 | (-1)^n |

462 | """ |

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

464 | |

465 | |

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

467 | """ |

468 | |

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

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

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

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

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

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

475 | |

476 | performance: |

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

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

479 | #Wall time: 0.16 s |

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

481 | |

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

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

484 | #Wall time: 0.01 s |

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

486 | |

487 | #time chebyshev_T(50,x) |

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

489 | #Wall time: 0.04 s |

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

491 | |

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

493 | #CPU times: user 0.08 s, sys: 0.00 s, total: 0.08 s |

494 | #Wall time: 0.08 s |

495 | |

496 | EXAMPLES:: |

497 | sage: chebyshev_T(5,x) |

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

499 | sage: var('n') |

500 | n |

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

502 | (-1)^n |

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

504 | chebyshev_T(-7, x) |

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

506 | chebyshev_T(3/2, x) |

507 | |

508 | """ |

509 | |

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

511 | |

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

513 | try: |

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

515 | return self._evalf_(*args) |

516 | except AttributeError: |

517 | pass |

518 | except mpmath.NoConvergence: |

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

520 | print "Switching to clenshaw_method, but it \ |

521 | may not be stable!" |

522 | except ValueError: |

523 | pass |

524 | |

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

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

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

528 | try: |

529 | return self._evalf_(*args) |

530 | except AttributeError: |

531 | pass |

532 | else: |

533 | return None |

534 | |

535 | if args[0] < 0: |

536 | return None |

537 | |

538 | |

539 | try: |

540 | return self._eval_special_values_(*args) |

541 | except ValueError: |

542 | pass |

543 | |

544 | |

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

546 | |

547 | try: |

548 | return self._clenshaw_method_(*args) |

549 | except NotImplementedError: |

550 | pass |

551 | |

552 | if self._maxima_name is None: |

553 | return None |

554 | else: |

555 | _init() |

556 | try: |

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

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

559 | #methods were much faster... |

560 | return self._maxima_init_evaled_(*args) |

561 | except TypeError: |

562 | return None |

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

564 | return None |

565 | else: |

566 | return s.sage() |

567 | |

568 | #def _eval_numpy_(self, *args): |

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

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

571 | #of ortho polys |

572 | #Now this only evaluates the array pointwise, and only the first one... |

573 | #if isinstance(arg[0], numpy.ndarray) |

574 | |

575 | |

576 | |

577 | |

578 | class Func_chebyshev_T(OrthogonalPolynomial): |

579 | |

580 | """ |

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

582 | |

583 | REFERENCE: |

584 | |

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

586 | |

587 | EXAMPLES:: |

588 | |

589 | sage: chebyshev_T(3,x) |

590 | 4*x^3 - 3*x |

591 | sage: var('k') |

592 | k |

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

594 | sage: test |

595 | chebyshev_T(k, x) |

596 | """ |

597 | |

598 | def __init__(self): |

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

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

601 | |

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

603 | """ |

604 | Values known for special values of x. |

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

606 | |

607 | EXAMPLES: |

608 | |

609 | sage: var('n') |

610 | n |

611 | sage: chebyshev_T(n,1) |

612 | 1 |

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

614 | (-1)^n |

615 | """ |

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

617 | return 1 |

618 | |

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

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

621 | |

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

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

624 | |

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

626 | |

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

628 | """ |

629 | Evals chebyshev_T |

630 | numerically with mpmath. |

631 | EXAMPLES:: |

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

633 | 2.261953700000000000000e7 |

634 | """ |

635 | try: |

636 | step_parent = kwds['parent'] |

637 | except KeyError: |

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

639 | |

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

641 | |

642 | try: |

643 | precision = step_parent.prec() |

644 | except AttributeError: |

645 | precision = mpmath.mp.prec |

646 | |

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

648 | |

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

650 | n = args[0] |

651 | x = args[1] |

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

653 | |

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

655 | """ |

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

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

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

659 | """ |

660 | |

661 | k = args[0] |

662 | x = args[1] |

663 | |

664 | if k == 0: |

665 | return 1 |

666 | elif k == 1: |

667 | return x |

668 | else: |

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

670 | #use these explicit formulas instead! |

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

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

673 | #elif 1 < x: |

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

675 | #else: # x < -1 |

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

677 | |

678 | help1 = 1 |

679 | help2 = x |

680 | if is_Expression(x): |

681 | #raise NotImplementedError |

682 | help3 = 0 |

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

684 | help3 = \ |

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

686 | help3 = help3*k/2 |

687 | else: |

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

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

690 | help1 = help2 |

691 | help2 = help3 |

692 | |

693 | return help3 |

694 | |

695 | def _eval_numpy_(self, *args): |

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

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

698 | #of ortho polys |

699 | #Now this only evaluates the array pointwise, and only the first one... |

700 | #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability! |

701 | """ |

702 | EXAMPLES:: |

703 | sage: import numpy |

704 | sage: z = numpy.array([1,2]) |

705 | sage: z2 = numpy.array([[1,2],[1,2]]) |

706 | sage: z3 = numpy.array([1,2,3.]) |

707 | sage: chebyshev_T(1,z) |

708 | array([1, 2]) |

709 | sage: chebyshev_T(1,z2) |

710 | array([[1, 2], |

711 | [1, 2]]) |

712 | sage: chebyshev_T(1,z3) |

713 | array([ 1., 2., 3.]) |

714 | |

715 | """ |

716 | if isinstance(args[0],numpy.ndarray): |

717 | raise NotImplementedError("Support for numpy array in \ |

718 | first argument(s) is not supported yet!") |

719 | |

720 | result = numpy.zeros(args[-1].shape).tolist() |

721 | if isinstance(args[-1][0],int): |

722 | for k in xrange(len(args[-1])): |

723 | result[k] = chebyshev_T(args[0],ZZ(args[-1][k])) |

724 | |

725 | if isinstance(args[-1][0],float): |

726 | for k in xrange(len(args[-1])): |

727 | result[k] = chebyshev_T(args[0],RR(args[-1][k])) |

728 | |

729 | if isinstance(args[-1][0],numpy.ndarray): |

730 | for k in xrange(len(args[-1])): |

731 | result[k] = chebyshev_T(args[0],args[-1][k]) |

732 | |

733 | return numpy.array(result) |

734 | |

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

736 | """ |

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

738 | of the second kind chebyshev_U |

739 | EXAMPLES:: |

740 | sage: var('k') |

741 | k |

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

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

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

745 | 12*x^2 - 3 |

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

747 | Traceback (most recent call last): |

748 | ... |

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

750 | |

751 | """ |

752 | diff_param = kwds['diff_param'] |

753 | if diff_param == 0: |

754 | raise NotImplementedError( |

755 | "Derivative w.r.t. to the index is not supported, yet, \ |

756 | and perhaps never will be...") |

757 | else: |

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

759 | |

760 | |

761 | chebyshev_T = Func_chebyshev_T() |

762 | |

763 | class Func_chebyshev_U(OrthogonalPolynomial): |

764 | |

765 | """ |

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

767 | |

768 | REFERENCE: |

769 | |

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

771 | |

772 | EXAMPLES:: |

773 | |

774 | sage: chebyshev_U(3,x) |

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

776 | """ |

777 | |

778 | def __init__(self): |

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

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

781 | |

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

783 | """ |

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

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

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

787 | """ |

788 | k = args[0] |

789 | x = args[1] |

790 | |

791 | if k == 0: |

792 | return 1 |

793 | elif k == 1: |

794 | return 2*x |

795 | else: |

796 | help1 = 1 |

797 | help2 = 2*x |

798 | if is_Expression(x): |

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

800 | help3 = 0 |

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

802 | help3 = \ |

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

804 | |

805 | else: |

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

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

808 | help1 = help2 |

809 | help2 = help3 |

810 | |

811 | return help3 |

812 | |

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

814 | """ |

815 | Uses |

816 | EXAMPLES:: |

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

818 | sage: chebyshev_T(2,x) |

819 | 2*x^2 - 1 |

820 | """ |

821 | n = args[0] |

822 | x = args[1] |

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

824 | |

825 | |

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

827 | """ |

828 | Evals chebyshev_U |

829 | numerically with mpmath. |

830 | EXAMPLES:: |

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

832 | 4.661117900000000000000e7 |

833 | """ |

834 | try: |

835 | step_parent = kwds['parent'] |

836 | except KeyError: |

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

838 | |

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

840 | |

841 | try: |

842 | precision = step_parent.prec() |

843 | except AttributeError: |

844 | precision = mpmath.mp.prec |

845 | |

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

847 | |

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

849 | """ |

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

851 | EXAMPLES: |

852 | |

853 | sage: var('n') |

854 | n |

855 | sage: chebyshev_U(n,1) |

856 | n + 1 |

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

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

859 | """ |

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

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

862 | |

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

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

865 | |

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

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

868 | |

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

870 | |

871 | def _eval_numpy_(self, *args): |

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

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

874 | #of ortho polys |

875 | #Now this only evaluates the array pointwise, and only the first one... |

876 | #if isinstance(arg[0], numpy.ndarray). |

877 | #This is a hack to provide compability! |

878 | """ |

879 | EXAMPLES:: |

880 | sage: import numpy |

881 | sage: z = numpy.array([1,2]) |

882 | sage: z2 = numpy.array([[1,2],[1,2]]) |

883 | sage: z3 = numpy.array([1,2,3.]) |

884 | sage: chebyshev_U(1,z) |

885 | array([2, 4]) |

886 | sage: chebyshev_U(1,z2) |

887 | array([[2, 4], |

888 | [2, 4]]) |

889 | sage: chebyshev_U(1,z3) |

890 | array([ 2., 4., 6.]) |

891 | |

892 | """ |

893 | if isinstance(args[0],numpy.ndarray): |

894 | raise NotImplementedError( |

895 | "Support for numpy array in first argument(s) is not supported yet!") |

896 | |

897 | result = numpy.zeros(args[-1].shape).tolist() |

898 | if isinstance(args[-1][0],int): |

899 | for k in xrange(len(args[-1])): |

900 | result[k] = chebyshev_U(args[0],ZZ(args[-1][k])) |

901 | |

902 | if isinstance(args[-1][0],float): |

903 | for k in xrange(len(args[-1])): |

904 | result[k] = chebyshev_U(args[0],RR(args[-1][k])) |

905 | |

906 | if isinstance(args[-1][0],numpy.ndarray): |

907 | for k in xrange(len(args[-1])): |

908 | result[k] = chebyshev_U(args[0],args[-1][k]) |

909 | |

910 | return numpy.array(result) |

911 | |

912 | |

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

914 | """ |

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

916 | Polynomials of the first and second kind |

917 | |

918 | EXAMPLES:: |

919 | sage: var('k') |

920 | k |

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

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

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

924 | 24*x^2 - 4 |

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

926 | Traceback (most recent call last): |

927 | ... |

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

929 | |

930 | """ |

931 | diff_param = kwds['diff_param'] |

932 | if diff_param == 0: |

933 | raise NotImplementedError( |

934 | "Derivative w.r.t. to the index is not supported, \ |

935 | yet, and perhaps never will be...") |

936 | else: |

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

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

939 | |

940 | chebyshev_U = Func_chebyshev_U() |

941 | |

942 | class Func_gegenbauer(OrthogonalPolynomial): |

943 | """ |

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

945 | |

946 | REFERENCE: |

947 | |

948 | - AS 22.5.27 |

949 | |

950 | EXAMPLES:: |

951 | |

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

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

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

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

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

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

958 | 2*x |

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

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

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

962 | """ |

963 | def __init__(self): |

964 | OrthogonalPolynomial.__init__(self,"gegenbauer",nargs = 3, |

965 | conversions =dict(maxima='ultraspherical',mathematica='GegenbauerC')) |

966 | |

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

968 | """ |

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

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

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

972 | """ |

973 | k = args[0] |

974 | x = args[-1] |

975 | alpha = args[1] |

976 | |

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

978 | alpha_zero = False |

979 | else: |

980 | alpha_zero = True |

981 | |

982 | if k == 0: |

983 | return 1 |

984 | elif k == 1: |

985 | if not alpha_zero: |

986 | return 2*x*alpha |

987 | else: |

988 | return 2*x # It seems that maxima evals this the wrong way! |

989 | #(see A.S. 22.4 (p.777)) |

990 | else: |

991 | help1 = 1 |

992 | if alpha_zero: |

993 | help2 = 2*x |

994 | gamma_alpha = 1 |

995 | else: |

996 | help2 = 2*x*alpha |

997 | gamma_alpha = gamma(alpha) |

998 | |

999 | help3 = 0 |

1000 | if is_Expression(x): |

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

1002 | help3 = \ |

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

1004 | |

1005 | else: |

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

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

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

1009 | help1 = help2 |

1010 | help2 = help3 |

1011 | |

1012 | return help3 |

1013 | |

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

1015 | """ |

1016 | Uses |

1017 | EXAMPLES:: |

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

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

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

1021 | """ |

1022 | n = args[0] |

1023 | a = args[1] |

1024 | x = args[2] |

1025 | return sage_eval(maxima.eval('ultraspherical(%s,%s,x)'%(ZZ(n),a)),\ |

1026 | locals={'x':x}) |

1027 | |

1028 | |

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

1030 | """ |

1031 | Evals gegenbauer |

1032 | numerically with mpmath. |

1033 | EXAMPLES:: |

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

1035 | 5.25360702000000e8 |

1036 | """ |

1037 | if args[1] == 0: |

1038 | raise ValueError("mpmath don't handle alpha = 0 correctly!") |

1039 | |

1040 | try: |

1041 | step_parent = kwds['parent'] |

1042 | except KeyError: |

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

1044 | |

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

1046 | |

1047 | try: |

1048 | precision = step_parent.prec() |

1049 | except AttributeError: |

1050 | precision = mpmath.mp.prec |

1051 | |

1052 | return mpmath.call( |

1053 | mpmath.gegenbauer,args[0],args[1],args[-1],prec = precision) |

1054 | |

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

1056 | """ |

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

1058 | EXAMPLES: |

1059 | |

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

1061 | (n, a) |

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

1063 | legendre_P(n, x) |

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

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

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

1067 | chebyshev_U(n, x) |

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

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

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

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

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

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

1074 | """ |

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

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

1077 | |

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

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

1080 | |

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

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

1083 | |

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

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

1086 | |

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

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

1089 | |

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

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

1092 | gamma(args[1]+args[0]/2)/gamma(args[1])/gamma(args[0]/2) |

1093 | |

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

1095 | |

1096 | def _eval_numpy_(self, *args): |

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

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

1099 | #of ortho polys |

1100 | #Now this only evaluates the array pointwise, and only the first one... |

1101 | #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability! |

1102 | """ |

1103 | EXAMPLES:: |

1104 | sage: import numpy |

1105 | sage: z = numpy.array([1,2]) |

1106 | sage: z2 = numpy.array([[1,2],[1,2]]) |

1107 | sage: z3 = numpy.array([1,2,3.]) |

1108 | sage: gegenbauer(2,0,z) |

1109 | array([1, 7]) |

1110 | sage: gegenbauer(2,1,z) |

1111 | array([ 3, 15]) |

1112 | sage: gegenbauer(2,1,z2) |

1113 | array([[ 3, 15], |

1114 | [ 3, 15]]) |

1115 | sage: gegenbauer(2,1,z3) |

1116 | array([ 3., 15., 35.]) |

1117 | """ |

1118 | |

1119 | if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray): |

1120 | raise NotImplementedError( |

1121 | "Support for numpy array in first argument(s) is not supported yet!") |

1122 | |

1123 | result = numpy.zeros(args[-1].shape).tolist() |

1124 | if isinstance(args[-1][0],int): |

1125 | for k in xrange(len(args[-1])): |

1126 | result[k] = gegenbauer(args[0],args[1],ZZ(args[-1][k])) |

1127 | |

1128 | if isinstance(args[-1][0],float): |

1129 | for k in xrange(len(args[-1])): |

1130 | result[k] = gegenbauer(args[0],args[1],RR(args[-1][k])) |

1131 | |

1132 | if isinstance(args[-1][0],numpy.ndarray): |

1133 | for k in xrange(len(args[-1])): |

1134 | result[k] = gegenbauer(args[0],args[1],args[-1][k]) |

1135 | |

1136 | return numpy.array(result) |

1137 | |

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

1139 | """ |

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

1141 | Polynomials of the first and second kind |

1142 | |

1143 | EXAMPLES:: |

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

1145 | (k, a) |

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

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

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

1149 | 960*x^3 - 240*x |

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

1151 | Traceback (most recent call last): |

1152 | ... |

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

1154 | |

1155 | """ |

1156 | diff_param = kwds['diff_param'] |

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

1158 | raise NotImplementedError( |

1159 | "Derivative w.r.t. to the index is not supported, yet,\ |

1160 | and perhaps never will be...") |

1161 | else: |

1162 | return (-args[0]*args[-1]*gegenbauer(args[0],args[1],args[2])+\ |

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

1164 | |

1165 | gegenbauer = Func_gegenbauer() |

1166 | ultraspherical = Func_gegenbauer() |

1167 | |

1168 | class Func_gen_laguerre(OrthogonalPolynomial): |

1169 | """ |

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

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

1172 | |

1173 | REFERENCE: |

1174 | |

1175 | - table on page 789 in AS. |

1176 | |

1177 | EXAMPLES:: |

1178 | |

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

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

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

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

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

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

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

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

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

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

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

1190 | """ |

1191 | def __init__(self): |

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

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

1194 | |

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

1196 | """ |

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

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

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

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

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

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

1203 | """ |

1204 | k = args[0] |

1205 | a = args[1] |

1206 | x = args[-1] |

1207 | |

1208 | if k == 0: |

1209 | return 1 |

1210 | elif k == 1: |

1211 | return 1-x + a |

1212 | else: |

1213 | help1 = 1 |

1214 | help2 = 1-x + a |

1215 | if is_Expression(x): |

1216 | help3 = 0 |

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

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

1219 | else: |

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

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

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

1223 | help1 = help2 |

1224 | help2 = help3 |

1225 | |

1226 | return help3 |

1227 | |

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

1229 | n = args[0] |

1230 | a = args[1] |

1231 | x = args[-1] |

1232 | |

1233 | from sage.functions.all import sqrt |

1234 | _init() |

1235 | return sage_eval(maxima.eval('gen_laguerre(%s,%s,x)'%(ZZ(n),a)),\ |

1236 | locals={'x':x}) |

1237 | |

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

1239 | """ |

1240 | Evals laguerre polynomial |

1241 | numerically with mpmath. |

1242 | EXAMPLES:: |

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

1244 | 2.66666666666667 |

1245 | """ |

1246 | try: |

1247 | step_parent = kwds['parent'] |

1248 | except KeyError: |

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

1250 | |

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

1252 | |

1253 | try: |

1254 | precision = step_parent.prec() |

1255 | except AttributeError: |

1256 | precision = mpmath.mp.prec |

1257 | |

1258 | return mpmath.call( |

1259 | mpmath.laguerre,args[0],args[1],args[-1],prec = precision) |

1260 | |

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

1262 | """ |

1263 | Special values known. |

1264 | EXAMPLES: |

1265 | |

1266 | sage: var('n') |

1267 | n |

1268 | sage: gen_laguerre(n,1,0) |

1269 | binomial(n + 1, n) |

1270 | """ |

1271 | |

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

1273 | try: |

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

1275 | except TypeError: |

1276 | pass |

1277 | |

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

1279 | |

1280 | def _eval_numpy_(self, *args): |

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

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

1283 | #of ortho polys |

1284 | #Now this only evaluates the array pointwise, and only the first one... |

1285 | #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability! |

1286 | """ |

1287 | EXAMPLES:: |

1288 | sage: import numpy |

1289 | sage: z = numpy.array([1,2]) |

1290 | sage: z2 = numpy.array([[1,2],[1,2]]) |

1291 | sage: z3 = numpy.array([1,2,3.]) |

1292 | sage: gen_laguerre(1,1,z) |

1293 | array([1, 0]) |

1294 | sage: gen_laguerre(1,1,z2) |

1295 | array([[1, 0], |

1296 | [1, 0]]) |

1297 | sage: gen_laguerre(1,1,z3) |

1298 | array([ 1., 0., -1.]) |

1299 | """ |

1300 | |

1301 | if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray): |

1302 | raise NotImplementedError( |

1303 | "Support for numpy array in first argument(s) is not supported yet!") |

1304 | |

1305 | result = numpy.zeros(args[-1].shape).tolist() |

1306 | if isinstance(args[-1][0],int): |

1307 | for k in xrange(len(args[-1])): |

1308 | result[k] = gen_laguerre(args[0],args[1],ZZ(args[-1][k])) |

1309 | |

1310 | if isinstance(args[-1][0],float): |

1311 | for k in xrange(len(args[-1])): |

1312 | result[k] = gen_laguerre(args[0],args[1],RR(args[-1][k])) |

1313 | |

1314 | if isinstance(args[-1][0],numpy.ndarray): |

1315 | for k in xrange(len(args[-1])): |

1316 | result[k] = gen_laguerre(args[0],args[1],args[-1][k]) |

1317 | |

1318 | return numpy.array(result) |

1319 | |

1320 | |

1321 | |

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

1323 | """return the derivative of gen_laguerre in |

1324 | form of the gen. Laguerre Polynomial. |

1325 | EXAMPLES:: |

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

1327 | sage: derivative(gen_laguerre(3,1,x),x) |

1328 | -1/2*x^2 + 4*x - 6 |

1329 | sage: derivative(gen_laguerre(n,1,x),x) |

1330 | -((n + 1)*gen_laguerre(n - 1, 1, x) - n*gen_laguerre(n, 1, x))/x |

1331 | """ |

1332 | diff_param = kwds['diff_param'] |

1333 | if diff_param == 0: |

1334 | raise NotImplementedError( |

1335 | "Derivative w.r.t. to the index is not supported, \ |

1336 | yet, and perhaps never will be...") |

1337 | else: |

1338 | return (args[0]*gen_laguerre(args[0],args[1],args[-1])-(args[0]+args[1])*\ |

1339 | gen_laguerre(args[0]-1,args[1],args[-1]))/args[-1] |

1340 | |

1341 | gen_laguerre = Func_gen_laguerre() |

1342 | |

1343 | class Func_hermite(OrthogonalPolynomial): |

1344 | """ |

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

1346 | |

1347 | REFERENCE: |

1348 | |

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

1350 | |

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

1352 | # <- Here a bug of the base class BuiltinFunction occours! |

1353 | |

1354 | EXAMPLES:: |

1355 | |

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

1357 | sage: hermite(2,x) |

1358 | 4*x^2 - 2 |

1359 | sage: hermite(3,x) |

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

1361 | sage: hermite(3,2) |

1362 | 40 |

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

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

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

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

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

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

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

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

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

1372 | """ |

1373 | |

1374 | def __init__(self): |

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

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

1377 | |

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

1379 | """ |

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

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

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

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

1384 | n = 25 |

1385 | """ |

1386 | k = args[0] |

1387 | x = args[1] |

1388 | |

1389 | if k == 0: |

1390 | return 1 |

1391 | elif k == 1: |

1392 | return 2*x |

1393 | else: |

1394 | help1 = 1 |

1395 | help2 = 2*x |

1396 | if is_Expression(x): |

1397 | help3 = 0 |

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

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

1400 | factorial(k-2*j) |

1401 | help3 = help3*factorial(k) |

1402 | else: |

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

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

1405 | help1 = help2 |

1406 | help2 = help3 |

1407 | |

1408 | return help3 |

1409 | |

1410 | |

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

1412 | """ |

1413 | Evals hermite |

1414 | numerically with mpmath. |

1415 | EXAMPLES:: |

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

1417 | 40.0000000000000000000 |

1418 | """ |

1419 | try: |

1420 | step_parent = kwds['parent'] |

1421 | except KeyError: |

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

1423 | |

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

1425 | |

1426 | try: |

1427 | precision = step_parent.prec() |

1428 | except AttributeError: |

1429 | precision = mpmath.mp.prec |

1430 | |

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

1432 | |

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

1434 | """ |

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

1436 | EXAMPLES: |

1437 | |

1438 | sage: var('n') |

1439 | n |

1440 | sage: hermite(n,0) |

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

1442 | """ |

1443 | |

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

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

1446 | factorial(args[0])/gamma(args[0]/2+1) |

1447 | |

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

1449 | |

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

1451 | """ |

1452 | Old maxima method. |

1453 | """ |

1454 | n = args[0] |

1455 | x = args[1] |

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

1457 | |

1458 | def _eval_numpy_(self, *args): |

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

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

1461 | #of ortho polys |

1462 | #Now this only evaluates the array pointwise, and only the first one... |

1463 | #if isinstance(arg[0], numpy.ndarray). |

1464 | #This is a hack to provide compability! |

1465 | """ |

1466 | EXAMPLES:: |

1467 | sage: import numpy |

1468 | sage: z = numpy.array([1,2]) |

1469 | sage: z2 = numpy.array([[1,2],[1,2]]) |

1470 | sage: z3 = numpy.array([1,2,3.]) |

1471 | sage: hermite(1,z) |

1472 | array([2, 4]) |

1473 | sage: hermite(1,z2) |

1474 | array([[2, 4], |

1475 | [2, 4]]) |

1476 | sage: hermite(1,z3) |

1477 | array([ 2., 4., 6.]) |

1478 | |

1479 | """ |

1480 | if isinstance(args[0],numpy.ndarray): |

1481 | raise NotImplementedError( |

1482 | "Support for numpy array in first argument(s) is not supported yet!") |

1483 | |

1484 | result = numpy.zeros(args[-1].shape).tolist() |

1485 | if isinstance(args[-1][0],int): |

1486 | for k in xrange(len(args[-1])): |

1487 | result[k] = hermite(args[0],ZZ(args[-1][k])) |

1488 | |

1489 | if isinstance(args[-1][0],float): |

1490 | for k in xrange(len(args[-1])): |

1491 | result[k] = hermite(args[0],RR(args[-1][k])) |

1492 | |

1493 | if isinstance(args[-1][0],numpy.ndarray): |

1494 | for k in xrange(len(args[-1])): |

1495 | result[k] = hermite(args[0],args[-1][k]) |

1496 | |

1497 | return numpy.array(result) |

1498 | |

1499 | |

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

1501 | """ |

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

1503 | Polynomials of the first and second kind |

1504 | |

1505 | EXAMPLES:: |

1506 | sage: var('k') |

1507 | k |

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

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

1510 | |

1511 | """ |

1512 | diff_param = kwds['diff_param'] |

1513 | if diff_param == 0: |

1514 | raise NotImplementedError( |

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

1516 | else: |

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

1518 | |

1519 | hermite = Func_hermite() |

1520 | |

1521 | |

1522 | class Func_jacobi_P(OrthogonalPolynomial): |

1523 | r""" |

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

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

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

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

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

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

1530 | |

1531 | REFERENCE: |

1532 | |

1533 | - table on page 789 in AS. |

1534 | |

1535 | EXAMPLES:: |

1536 | |

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

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

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

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

1541 | 5.009999999999998 |

1542 | """ |

1543 | |

1544 | def __init__(self): |

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

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

1547 | |

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

1549 | """ |

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

1551 | or sum formula) |

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

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

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

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

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

1557 | """ |

1558 | k = args[0] |

1559 | x = args[-1] |

1560 | alpha = args[1] |

1561 | beta = args[2] |

1562 | |

1563 | if k == 0: |

1564 | return 1 |

1565 | elif k == 1: |

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

1567 | else: |

1568 | |

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

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

1571 | #optimal for use. |

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

1573 | help2 = 0 |

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

1575 | help2 = help2 + binomial(k,j)*gamma(alpha+beta+k+j+1)/\ |

1576 | gamma(alpha+j+1)*((x-1)/2)**j |

1577 | return help1*help2 |

1578 | else: |

1579 | help1 = 1 |

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

1581 | |

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

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

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

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

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

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

1588 | help3 = help3/a1 |

1589 | help1 = help2 |

1590 | help2 = help3 |

1591 | |

1592 | return help3 |

1593 | |

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

1595 | """ |

1596 | Old maxima method. |

1597 | """ |

1598 | n = args[0] |

1599 | a = args[1] |

1600 | b = args[2] |

1601 | x = args[3] |

1602 | return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)),\ |

1603 | locals={'x':x}) |

1604 | |

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

1606 | """ |

1607 | Evals jacobi_P |

1608 | numerically with mpmath. |

1609 | EXAMPLES:: |

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

1611 | 1.322776620000000000000e8 |

1612 | """ |

1613 | try: |

1614 | step_parent = kwds['parent'] |

1615 | except KeyError: |

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

1617 | |

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

1619 | |

1620 | try: |

1621 | precision = step_parent.prec() |

1622 | except AttributeError: |

1623 | precision = mpmath.mp.prec |

1624 | |

1625 | return mpmath.call(mpmath.jacobi,args[0],args[1],args[2],args[-1],\ |

1626 | prec = precision) |

1627 | |

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

1629 | """ |

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

1631 | EXAMPLES: |

1632 | |

1633 | sage: var('n') |

1634 | n |

1635 | sage: legendre_P(n,1) |

1636 | 1 |

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

1638 | (-1)^n |

1639 | """ |

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

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

1642 | |

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

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

1645 | |

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

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

1648 | |

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

1650 | try: |

1651 | return binomial(2*args[0],args[0])*\ |

1652 | chebyshev_T(args[0],args[-1])/4**args[0] |

1653 | except TypeError: |

1654 | pass |

1655 | |

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

1657 | |

1658 | def _eval_numpy_(self, *args): |

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

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

1661 | #of ortho polys |

1662 | #Now this only evaluates the array pointwise, and only the first one... |

1663 | #if isinstance(arg[0], numpy.ndarray). |

1664 | #This is a hack to provide compability! |

1665 | """ |

1666 | EXAMPLES:: |

1667 | sage: import numpy |

1668 | sage: z = numpy.array([1,2]) |

1669 | sage: z2 = numpy.array([[1,2],[1,2]]) |

1670 | sage: z3 = numpy.array([1,2,3.]) |

1671 | sage: gen_laguerre(1,1,z) |

1672 | array([1, 0]) |

1673 | sage: gen_laguerre(1,1,z2) |

1674 | array([[1, 0], |

1675 | [1, 0]]) |

1676 | sage: gen_laguerre(1,1,z3) |

1677 | array([ 1., 0., -1.]) |

1678 | """ |

1679 | |

1680 | if isinstance(args[0],numpy.ndarray) or \ |

1681 | isinstance(args[1],numpy.ndarray) or \ |

1682 | isinstance(args[2],numpy.ndarray): |

1683 | raise NotImplementedError( |

1684 | "Support for numpy array in first argument(s) is not supported yet!") |

1685 | |

1686 | result = numpy.zeros(args[-1].shape).tolist() |

1687 | if isinstance(args[-1][0],int): |

1688 | for k in xrange(len(args[-1])): |

1689 | result[k] = jacobi_P(args[0],args[1],args[2],ZZ(args[-1][k])) |

1690 | |

1691 | if isinstance(args[-1][0],float): |

1692 | for k in xrange(len(args[-1])): |

1693 | result[k] = jacobi_P(args[0],args[1],args[2],RR(args[-1][k])) |

1694 | |

1695 | if isinstance(args[-1][0],numpy.ndarray): |

1696 | for k in xrange(len(args[-1])): |

1697 | result[k] = jacobi_P(args[0],args[1],args[2],args[-1][k]) |

1698 | |

1699 | return numpy.array(result) |

1700 | |

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

1702 | """ |

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

1704 | |

1705 | EXAMPLES:: |

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

1707 | (k, a, b) |

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

1709 | -(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)) |

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

1711 | 14*x - 7/2 |

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

1713 | Traceback (most recent call last): |

1714 | ... |

1715 | |

1716 | |

1717 | """ |

1718 | diff_param = kwds['diff_param'] |

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

1720 | raise NotImplementedError( |

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

1722 | else: |

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

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

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

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

1727 | |

1728 | |

1729 | jacobi_P = Func_jacobi_P() |

1730 | |

1731 | class Func_laguerre(OrthogonalPolynomial): |

1732 | """ |

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

1734 | |

1735 | REFERENCE: |

1736 | |

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

1738 | |

1739 | EXAMPLES:: |

1740 | |

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

1742 | sage: laguerre(2,x) |

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

1744 | sage: laguerre(3,x) |

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

1746 | sage: laguerre(2,2) |

1747 | -1 |

1748 | """ |

1749 | def __init__(self): |

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

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

1752 | |

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

1754 | """ |

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

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

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

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

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

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

1761 | """ |

1762 | k = args[0] |

1763 | x = args[-1] |

1764 | |

1765 | if k == 0: |

1766 | return 1 |

1767 | elif k == 1: |

1768 | return 1-x |

1769 | else: |

1770 | help1 = 1 |

1771 | help2 = 1-x |

1772 | if is_Expression(x): |

1773 | help3 = 0 |

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

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

1776 | else: |

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

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

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

1780 | help1 = help2 |

1781 | help2 = help3 |

1782 | |

1783 | return help3 |

1784 | |

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

1786 | n = args[0] |

1787 | x = args[1] |

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

1789 | |

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

1791 | """ |

1792 | Evals laguerre polynomial |

1793 | numerically with mpmath. |

1794 | EXAMPLES:: |

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

1796 | 2.66666666666667 |

1797 | """ |

1798 | try: |

1799 | step_parent = kwds['parent'] |

1800 | except KeyError: |

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

1802 | |

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

1804 | |

1805 | try: |

1806 | precision = step_parent.prec() |

1807 | except AttributeError: |

1808 | precision = mpmath.mp.prec |

1809 | |

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

1811 | |

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

1813 | """ |

1814 | Special values known. |

1815 | EXAMPLES: |

1816 | |

1817 | sage: var('n') |

1818 | n |

1819 | sage: laguerre(n,0) |

1820 | 1 |

1821 | """ |

1822 | |

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

1824 | try: |

1825 | return 1 |

1826 | except TypeError: |

1827 | pass |

1828 | |

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

1830 | |

1831 | |

1832 | def _eval_numpy_(self, *args): |

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

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

1835 | #of ortho polys |

1836 | #Now this only evaluates the array pointwise, and only the first one... |

1837 | #if isinstance(arg[0], numpy.ndarray). |

1838 | #This is a hack to provide compability! |

1839 | """ |

1840 | EXAMPLES:: |

1841 | sage: import numpy |

1842 | sage: z = numpy.array([1,2]) |

1843 | sage: z2 = numpy.array([[1,2],[1,2]]) |

1844 | sage: z3 = numpy.array([1,2,3.]) |

1845 | sage: laguerre(1,z) |

1846 | array([ 0, -1]) |

1847 | sage: laguerre(1,z2) |

1848 | array([[ 0, -1], |

1849 | [ 0, -1]]) |

1850 | sage: laguerre(1,z3) |

1851 | array([ 0., -1., -2.]) |

1852 | |

1853 | """ |

1854 | if isinstance(args[0],numpy.ndarray): |

1855 | raise NotImplementedError( |

1856 | "Support for numpy array in first argument(s) is not supported yet!") |

1857 | |

1858 | result = numpy.zeros(args[-1].shape).tolist() |

1859 | if isinstance(args[-1][0],int): |

1860 | for k in xrange(len(args[-1])): |

1861 | result[k] = laguerre(args[0],ZZ(args[-1][k])) |

1862 | |

1863 | if isinstance(args[-1][0],float): |

1864 | for k in xrange(len(args[-1])): |

1865 | result[k] = laguerre(args[0],RR(args[-1][k])) |

1866 | |

1867 | if isinstance(args[-1][0],numpy.ndarray): |

1868 | for k in xrange(len(args[-1])): |

1869 | result[k] = laguerre(args[0],args[-1][k]) |

1870 | |

1871 | return numpy.array(result) |

1872 | |

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

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

1875 | form of the Laguerre Polynomial. |

1876 | EXAMPLES:: |

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

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

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

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

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

1882 | """ |

1883 | diff_param = kwds['diff_param'] |

1884 | if diff_param == 0: |

1885 | raise NotImplementedError( |

1886 | "Derivative w.r.t. to the index is not supported, \ |

1887 | yet, and perhaps never will be...") |

1888 | else: |

1889 | return (args[0]*laguerre(args[0],args[-1])-args[0]*\ |

1890 | laguerre(args[0]-1,args[-1]))/args[1] |

1891 | |

1892 | laguerre = Func_laguerre() |

1893 | |

1894 | class Func_legendre_P(OrthogonalPolynomial): |

1895 | """ |

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

1897 | |

1898 | REFERENCE: |

1899 | |

1900 | - AS 22.5.35 page 779. |

1901 | |

1902 | EXAMPLES:: |

1903 | |

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

1905 | sage: legendre_P(2,t) |

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

1907 | sage: legendre_P(3, 1.1) |

1908 | 1.67750000000000 |

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

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

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

1912 | Traceback (most recent call last): |

1913 | ... |

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

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

1916 | 8 |

1917 | """ |

1918 | def __init__(self): |

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

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

1921 | |

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

1923 | """ |

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

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

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

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

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

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

1930 | """ |

1931 | k = args[0] |

1932 | x = args[-1] |

1933 | |

1934 | if k == 0: |

1935 | return 1 |

1936 | elif k == 1: |

1937 | return x |

1938 | else: |

1939 | help1 = 1 |

1940 | help2 = x |

1941 | if is_Expression(x): |

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

1943 | help1 = ZZ(2**k) #Workarround because of segmentation fault... |

1944 | help3 = 0 |

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

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

1947 | binomial(2*(k-j),k) |

1948 | |

1949 | help3 = help3/help1 |

1950 | else: |

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

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

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

1954 | help1 = help2 |

1955 | help2 = help3 |

1956 | |

1957 | return help3 |

1958 | |

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

1960 | n = args[0] |

1961 | x = args[1] |

1962 | return sage_eval(maxima.eval('legendre_p(%s,x)'%ZZ(n)),\ |

1963 | locals={'x':x}) |

1964 | |

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

1966 | """ |

1967 | Evals legendre_P |

1968 | numerically with mpmath. |

1969 | EXAMPLES:: |

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

1971 | 8.097453000000000000000e6 |

1972 | """ |

1973 | try: |

1974 | step_parent = kwds['parent'] |

1975 | except KeyError: |

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

1977 | |

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

1979 | |

1980 | try: |

1981 | precision = step_parent.prec() |

1982 | except AttributeError: |

1983 | precision = mpmath.mp.prec |

1984 | |

1985 | return mpmath.call( |

1986 | mpmath.legendre,args[0],args[-1],prec = precision) |

1987 | |

1988 | |

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

1990 | """ |

1991 | Special values known. |

1992 | EXAMPLES: |

1993 | |

1994 | sage: var('n') |

1995 | n |

1996 | sage: legendre_P(n,1) |

1997 | 1 |

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

1999 | (-1)^n |

2000 | """ |

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

2002 | return 1 |

2003 | |

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

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

2006 | |

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

2008 | try: |

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

2010 | except TypeError: |

2011 | pass |

2012 | |

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

2014 | |

2015 | def _eval_numpy_(self, *args): |

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

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

2018 | #of ortho polys |

2019 | #Now this only evaluates the array pointwise, and only the first one... |

2020 | #if isinstance(arg[0], numpy.ndarray). |

2021 | #This is a hack to provide compability! |

2022 | """ |

2023 | EXAMPLES:: |

2024 | sage: import numpy |

2025 | sage: z = numpy.array([1,2]) |

2026 | sage: z2 = numpy.array([[1,2],[1,2]]) |

2027 | sage: z3 = numpy.array([1,2,3.]) |

2028 | sage: legendre_P(1,z) |

2029 | array([1, 2]) |

2030 | sage: legendre_P(1,z2) |

2031 | array([[1, 2], |

2032 | [1, 2]]) |

2033 | sage: legendre_P(1,z3) |

2034 | array([ 1., 2., 3.]) |

2035 | |

2036 | """ |

2037 | if isinstance(args[0],numpy.ndarray): |

2038 | raise NotImplementedError( |

2039 | "Support for numpy array in first argument(s) is not supported yet!") |

2040 | |

2041 | result = numpy.zeros(args[-1].shape).tolist() |

2042 | if isinstance(args[-1][0],int): |

2043 | for k in xrange(len(args[-1])): |

2044 | result[k] = legendre_P(args[0],ZZ(args[-1][k])) |

2045 | |

2046 | if isinstance(args[-1][0],float): |

2047 | for k in xrange(len(args[-1])): |

2048 | result[k] = legendre_P(args[0],RR(args[-1][k])) |

2049 | |

2050 | if isinstance(args[-1][0],numpy.ndarray): |

2051 | for k in xrange(len(args[-1])): |

2052 | result[k] = legendre_P(args[0],args[-1][k]) |

2053 | |

2054 | return numpy.array(result) |

2055 | |

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

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

2058 | form of the Legendre Polynomial. |

2059 | EXAMPLES:: |

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

2061 | sage: derivative(legendre_P(n,x),x) |

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

2063 | sage: derivative(legendre_P(3,x),x) |

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

2065 | """ |

2066 | diff_param = kwds['diff_param'] |

2067 | if diff_param == 0: |

2068 | raise NotImplementedError( |

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

2070 | else: |

2071 | return (args[0]*legendre_P(args[0]-1,args[-1])-args[0]*args[-1]*\ |

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

2073 | |

2074 | |

2075 | legendre_P = Func_legendre_P() |

2076 | |

2077 | class Func_legendre_Q(OrthogonalPolynomial): |

2078 | """ |

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

2080 | |

2081 | REFERENCE: |

2082 | |

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

2084 | |

2085 | EXAMPLES:: |

2086 | |

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

2088 | sage: legendre_Q(2, t) |

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

2090 | sage: legendre_Q(3, 0.5) |

2091 | -0.198654771479482 |

2092 | sage: legendre_Q(4, 2) |

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

2094 | sage: legendre_Q(4, 2.0) |

2095 | 0.00116107583162041 + 86.9828465962674*I |

2096 | |

2097 | """ |

2098 | def __init__(self): |

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

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

2101 | |

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

2103 | """ |

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

2105 | """ |

2106 | n = args[0] |

2107 | x = args[1] |

2108 | return sage_eval(maxima.eval('legendre_q(%s,x)'%ZZ(n)),\ |

2109 | locals={'x':x}) |

2110 | |

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

2112 | """ |

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

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

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

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

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

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

2119 | """ |

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

2121 | |

2122 | k = args[0] |

2123 | x = args[-1] |

2124 | |

2125 | if k == 0: |

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

2127 | elif k == 1: |

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

2129 | else: |

2130 | if is_Expression(x): |

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

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

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

2134 | else: |

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

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

2137 | |

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

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

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

2141 | help1 = help2 |

2142 | help2 = help3 |

2143 | |

2144 | return help3 |

2145 | |

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

2147 | """ |

2148 | Special values known. |

2149 | EXAMPLES: |

2150 | |

2151 | sage: var('n') |

2152 | n |

2153 | sage: legendre_Q(n,0) |

2154 | -1/2*sqrt(pi)*gamma(1/2*n + 1/2)*sin(1/2*pi*n)/gamma(1/2*n + 1) |

2155 | """ |

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

2157 | return NaN |

2158 | |

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

2160 | return NaN |

2161 | |

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

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

2164 | try: |

2165 | return -(sqrt(SR.pi()))/2*sin(SR.pi()/2*args[0])*\ |

2166 | gamma((args[0]+1)/2)/gamma(args[0]/2 + 1) |

2167 | except TypeError: |

2168 | pass |

2169 | else: |

2170 | return -(sqrt(math.pi))/2*sin(math.pi/2*args[0])*\ |

2171 | gamma((args[0]+1)/2)/gamma(args[0]/2. + 1) |

2172 | |

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

2174 | |

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

2176 | """ |

2177 | Evals legendre_Q |

2178 | numerically with mpmath. |

2179 | EXAMPLES:: |

2180 | sage: legendre_Q(10,3).n(75) |

2181 | 2.079454941572578263731e-9 + 1.271944942879431601408e7*I |

2182 | """ |

2183 | |

2184 | try: |

2185 | step_parent = kwds['parent'] |

2186 | except KeyError: |

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

2188 | |

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

2190 | |

2191 | try: |

2192 | precision = step_parent.prec() |

2193 | except AttributeError: |

2194 | precision = mpmath.mp.prec |

2195 | |

2196 | return conjugate( |

2197 | mpmath.call(mpmath.legenq,args[0],0,args[-1],prec = precision)) |

2198 | #it seems that mpmath uses here a different branch of the logarithm |

2199 | |

2200 | def _eval_numpy_(self, *args): |

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

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

2203 | #of ortho polys |

2204 | #Now this only evaluates the array pointwise, and only the first one... |

2205 | #if isinstance(arg[0], numpy.ndarray). This is a hack to provide compability! |

2206 | """ |

2207 | EXAMPLES:: |

2208 | sage: import numpy |

2209 | sage: z = numpy.array([1,2]) |

2210 | sage: z2 = numpy.array([[1,2],[1,2]]) |

2211 | sage: z3 = numpy.array([1,2,3.]) |

2212 | sage: legendre_Q(1,z/5.) |

2213 | array([-0.95945349, -0.83054043]) |

2214 | sage: legendre_Q(1,z2/5.) |

2215 | array([[-0.95945349, -0.83054043], |

2216 | [-0.95945349, -0.83054043]]) |

2217 | sage: legendre_Q(1,z3/5.) |

2218 | array([-0.95945349, -0.83054043, -0.58411169]) |

2219 | |

2220 | """ |

2221 | if isinstance(args[0],numpy.ndarray): |

2222 | raise NotImplementedError( |

2223 | "Support for numpy array in first argument(s) is not supported yet!") |

2224 | |

2225 | result = numpy.zeros(args[-1].shape).tolist() |

2226 | if isinstance(args[-1][0],int): |

2227 | for k in xrange(len(args[-1])): |

2228 | result[k] = legendre_Q(args[0],ZZ(args[-1][k])) |

2229 | |

2230 | if isinstance(args[-1][0],float): |

2231 | for k in xrange(len(args[-1])): |

2232 | result[k] = legendre_Q(args[0],RR(args[-1][k])) |

2233 | |

2234 | if isinstance(args[-1][0],numpy.ndarray): |

2235 | for k in xrange(len(args[-1])): |

2236 | result[k] = legendre_Q(args[0],args[-1][k]) |

2237 | |

2238 | return numpy.array(result) |

2239 | |

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

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

2242 | form of the Legendre Function. |

2243 | EXAMPLES:: |

2244 | n = var('n') |

2245 | derivative(legendre_Q(n,x),x) |

2246 | (n*x*legendre_Q(n, x) - n*legendre_Q(n - 1, x))/(x^2 - 1) |

2247 | sage: derivative(legendre_Q(3,x),x) |

2248 | 5/4*(x - 1)*(1/(x - 1) - (x + 1)/(x - 1)^2)*x^3/(x + 1) + 15/4*x^2*log(-(x + 1)/(x - 1)) - 3/4*(x - 1)*(1/(x - 1) - (x + 1)/(x - 1)^2)*x/(x + 1) - 5*x - 3/4*log(-(x + 1)/(x - 1)) |

2249 | """ |

2250 | diff_param = kwds['diff_param'] |

2251 | if diff_param == 0: |

2252 | raise NotImplementedError( |

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

2254 | else: |

2255 | return (args[0]*args[-1]*legendre_Q(args[0],args[-1])-args[0]*\ |

2256 | legendre_Q(args[0]-1,args[-1]))/(args[-1]**2-1) |

2257 | |

2258 | legendre_Q = Func_legendre_Q() |

2259 | |

2260 | |

2261 | class Func_gen_legendre_P(OrthogonalPolynomial): |

2262 | |

2263 | def __init__(self): |

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

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

2266 | |

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

2268 | """ |

2269 | Evals gen_legendre_P |

2270 | numerically with mpmath. |

2271 | EXAMPLES:: |

2272 | sage: gen_legendre_P(10,2,3).n(75) |

2273 | -7.194963600000000000000e8 |

2274 | """ |

2275 | |

2276 | try: |

2277 | step_parent = kwds['parent'] |

2278 | except KeyError: |

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

2280 | |

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

2282 | |

2283 | try: |

2284 | precision = step_parent.prec() |

2285 | except AttributeError: |

2286 | precision = mpmath.mp.prec |

2287 | |

2288 | return mpmath.call( |

2289 | mpmath.legenp,args[0],args[1],args[-1],prec = precision) |

2290 | |

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

2292 | """ |

2293 | Special values known. |

2294 | EXAMPLES: |

2295 | |

2296 | sage: n, m = var('n m') |

2297 | sage: gen_legendre_P(n,m,0) |

2298 | 2^m*gamma(1/2*m + 1/2*n + 1/2)*cos(1/2*(m + n)*pi)/(sqrt(pi)*gamma(-1/2*m + 1/2*n + 1)) |

2299 | """ |

2300 | |

2301 | if args[1] == 0: |

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

2303 | |

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

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

2306 | try: |

2307 | return cos(SR.pi()/2*(args[0]+args[1]))/(sqrt(SR.pi()))*\ |

2308 | gamma((args[0]+args[1]+1)/2)/\ |

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

2310 | except TypeError: |

2311 | pass |

2312 | else: |

2313 | return cos(math.pi/2*(args[0]+args[1]))/(sqrt(math.pi))*\ |

2314 | gamma((args[0]+args[1]+1)/2)/\ |

2315 | gamma((args[0]-args[1])/2. + 1)*2**args[1] |

2316 | |

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

2318 | |

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

2320 | n = args[0] |

2321 | m = args[1] |

2322 | x = args[2] |

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

2324 | return None |

2325 | |

2326 | from sage.functions.all import sqrt |

2327 | |

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

2329 | return sage_eval(maxima.eval('assoc_legendre_p(%s,%s,x)'\ |

2330 | %(ZZ(n),ZZ(m))), locals={'x':x}) |

2331 | else: |

2332 | return sqrt(1-x**2)*(((n-m+1)*x*gen_legendre_P(n,m-1,x)-\ |

2333 | (n+m-1)*gen_legendre_P(n-1,m-1,x))/(1-x**2)) |

2334 | |

2335 | |

2336 | def _eval_numpy_(self, *args): |

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

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

2339 | #of ortho polys |

2340 | #Now this only evaluates the array pointwise, and only the first one... |

2341 | #if isinstance(arg[0], numpy.ndarray). |

2342 | #This is a hack to provide compability! |

2343 | """ |

2344 | EXAMPLES:: |

2345 | sage: import numpy |

2346 | sage: z = numpy.array([1,2]) |

2347 | sage: z2 = numpy.array([[1,2],[1,2]]) |

2348 | sage: z3 = numpy.array([1,2,3.]) |

2349 | sage: gen_legendre_P(1,1,z/5.) |

2350 | array([-0.9797959 , -0.91651514]) |

2351 | sage: gen_legendre_P(1,1,z2/5.) |

2352 | array([[-0.9797959 , -0.91651514], |

2353 | [-0.9797959 , -0.91651514]]) |

2354 | sage: gen_legendre_P(1,1,z3/5.) |

2355 | array([-0.9797959 , -0.91651514, -0.8 ]) |

2356 | """ |

2357 | |

2358 | if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray): |

2359 | raise NotImplementedError( |

2360 | "Support for numpy array in first argument(s) is not supported yet!") |

2361 | |

2362 | result = numpy.zeros(args[-1].shape).tolist() |

2363 | if isinstance(args[-1][0],int): |

2364 | for k in xrange(len(args[-1])): |

2365 | result[k] = gen_legendre_P(args[0],args[1],ZZ(args[-1][k])) |

2366 | |

2367 | if isinstance(args[-1][0],float): |

2368 | for k in xrange(len(args[-1])): |

2369 | result[k] = gen_legendre_P(args[0],args[1],RR(args[-1][k])) |

2370 | |

2371 | if isinstance(args[-1][0],numpy.ndarray): |

2372 | for k in xrange(len(args[-1])): |

2373 | result[k] = gen_legendre_P(args[0],args[1],args[-1][k]) |

2374 | |

2375 | return numpy.array(result) |

2376 | |

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

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

2379 | form of the Legendre Function. |

2380 | EXAMPLES:: |

2381 | sage: n,m = var('n m') |

2382 | sage: derivative(gen_legendre_P(n,m,x),x) |

2383 | (n*x*gen_legendre_P(n, m, x) - (m + n)*gen_legendre_P(n - 1, m, x))/(x^2 - 1) |

2384 | sage: derivative(gen_legendre_P(3,1,x),x) |

2385 | -15*sqrt(-x^2 + 1)*x + 3/2*(5*(x - 1)^2 + 10*x - 6)*x/sqrt(-x^2 + 1) |

2386 | """ |

2387 | diff_param = kwds['diff_param'] |

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

2389 | raise NotImplementedError( |

2390 | "Derivative w.r.t. to the index is not supported,\ |

2391 | yet, and perhaps never will be...") |

2392 | else: |

2393 | return (args[0]*args[-1]*gen_legendre_P(args[0],args[1],args[-1])\ |

2394 | -(args[0]+args[1])*\ |

2395 | gen_legendre_P(args[0]-1,args[1],args[-1]))/(args[-1]**2-1) |

2396 | |

2397 | gen_legendre_P = Func_gen_legendre_P() |

2398 | |

2399 | class Func_gen_legendre_Q(OrthogonalPolynomial): |

2400 | |

2401 | def __init__(self): |

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

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

2404 | |

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

2406 | """ |

2407 | Evals gen_legendre_Q |

2408 | numerically with mpmath. |

2409 | EXAMPLES:: |

2410 | sage: gen_legendre_Q(10,2,3).n(75) |

2411 | -2.773909528741569374688e-7 - 1.130182239430298584113e9*I |

2412 | """ |

2413 | |

2414 | try: |

2415 | step_parent = kwds['parent'] |

2416 | except KeyError: |

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

2418 | |

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

2420 | |

2421 | try: |

2422 | precision = step_parent.prec() |

2423 | except AttributeError: |

2424 | precision = mpmath.mp.prec |

2425 | |

2426 | return mpmath.call( |

2427 | mpmath.legenq,args[0],args[1],args[-1],prec = precision) |

2428 | |

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

2430 | """ |

2431 | Special values known. |

2432 | EXAMPLES: |

2433 | |

2434 | sage: n, m = var('n m') |

2435 | sage: gen_legendre_Q(n,m,0) |

2436 | -sqrt(pi)*2^(m - 1)*gamma(1/2*m + 1/2*n + 1/2)*sin(1/2*(m + n)*pi)/gamma(-1/2*m + 1/2*n + 1) |

2437 | """ |

2438 | |

2439 | if args[1] == 0: |

2440 | return legendre_Q(args[0],args[-1]) |

2441 | |

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

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

2444 | try: |

2445 | return -(sqrt(SR.pi()))*sin(SR.pi()/2*(args[0]+args[1]))\ |

2446 | *gamma((args[0]+args[1]+1)/2)/\ |

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

2448 | except TypeError: |

2449 | pass |

2450 | else: |

2451 | return -(sqrt(math.pi))/2*sin(math.pi/2*(args[0]+args[1]))\ |

2452 | *gamma((args[0]+args[1]+1)/2)/\ |

2453 | gamma((args[0]-args[1])/2. + 1)*2**args[1] |

2454 | |

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

2456 | |

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

2458 | n = args[0] |

2459 | m = args[1] |

2460 | x = args[2] |

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

2462 | return None |

2463 | |

2464 | from sage.functions.all import sqrt |

2465 | |

2466 | if m <= n: |

2467 | return sage_eval(maxima.eval('assoc_legendre_q(%s,%s,x)'\ |

2468 | %(ZZ(n),ZZ(m))), locals={'x':x}) |

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

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

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

2472 | else: |

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

2474 | if m == n + 1: |

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

2476 | else: |

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

2478 | else: |

2479 | return ((n-m+1)*x*gen_legendre_Q(n,m-1,x)-\ |

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

2481 | |

2482 | |

2483 | def _eval_numpy_(self, *args): |

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

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

2486 | #of ortho polys |

2487 | #Now this only evaluates the array pointwise, and only the first one... |

2488 | #if isinstance(arg[0], numpy.ndarray). |

2489 | #This is a hack to provide compability! |

2490 | """ |

2491 | EXAMPLES:: |

2492 | sage: import numpy |

2493 | sage: z = numpy.array([1,2]) |

2494 | sage: z2 = numpy.array([[1,2],[1,2]]) |

2495 | sage: z3 = numpy.array([1,2,3.]) |

2496 | sage: gen_legendre_Q(1,1,z/5.) |

2497 | array([-0.40276067, -0.82471644]) |

2498 | sage: gen_legendre_Q(1,1,z2/5.) |

2499 | array([[-0.40276067, -0.82471644], |

2500 | [-0.40276067, -0.82471644]]) |

2501 | sage: gen_legendre_Q(1,1,z3/5.) |

2502 | array([-0.40276067, -0.82471644, -1.30451774]) |

2503 | """ |

2504 | |

2505 | if isinstance(args[0],numpy.ndarray) or isinstance(args[1],numpy.ndarray): |

2506 | raise NotImplementedError( |

2507 | "Support for numpy array in first argument(s) is not supported yet!") |

2508 | |

2509 | result = numpy.zeros(args[-1].shape).tolist() |

2510 | if isinstance(args[-1][0],int): |

2511 | for k in xrange(len(args[-1])): |

2512 | result[k] = gen_legendre_Q(args[0],args[1],ZZ(args[-1][k])) |

2513 | |

2514 | if isinstance(args[-1][0],float): |

2515 | for k in xrange(len(args[-1])): |

2516 | result[k] = gen_legendre_Q(args[0],args[1],RR(args[-1][k])) |

2517 | |

2518 | if isinstance(args[-1][0],numpy.ndarray): |

2519 | for k in xrange(len(args[-1])): |

2520 | result[k] = gen_legendre_Q(args[0],args[1],args[-1][k]) |

2521 | |

2522 | return numpy.array(result) |

2523 | |

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

2525 | """return the derivative of gen_legendre_Q in |

2526 | form of the Legendre Function. |

2527 | EXAMPLES:: |

2528 | sage: n,m = var('n m') |

2529 | sage: derivative(gen_legendre_Q(n,m,x),x) |

2530 | (n*x*gen_legendre_Q(n, m, x) - (m + n)*gen_legendre_Q(n - 1, m, x))/(x^2 - 1) |

2531 | sage: derivative(gen_legendre_Q(0,1,x),x) |

2532 | -x/(-x^2 + 1)^(3/2) |

2533 | """ |

2534 | diff_param = kwds['diff_param'] |

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

2536 | raise NotImplementedError("Derivative w.r.t. to the index \ |

2537 | is not supported, yet, and perhaps never will be...") |

2538 | else: |

2539 | return (args[0]*args[-1]*gen_legendre_Q(args[0],args[1],args[-1])\ |

2540 | -(args[0]+args[1])*\ |

2541 | gen_legendre_Q(args[0]-1,args[1],args[-1]))/(args[-1]**2-1) |

2542 | |

2543 | |

2544 | gen_legendre_Q = Func_gen_legendre_Q() |

2545 | |

2546 | |

2547 | |

2548 | |

2549 | |

2550 | |

2551 | |

2552 | |

2553 | |

2554 | |

2555 | |

2556 |