Ticket #4102: trac_symbolic_bessel_v3.patch
File trac_symbolic_bessel_v3.patch, 56.5 KB (added by , 9 years ago) |
---|
-
doc/en/reference/functions.rst
# HG changeset patch # User Benjamin Jones <benjaminfjones@gmail.com> # Date 1356654134 28800 # Node ID eb9efa4564e2b4f095c2abc0796e81fed6ec4510 # Parent 3df9be5aadd2e4159e69bdf911d38784cb6fa130 Trac 4102: Implement symbolic Bessel functions diff --git a/doc/en/reference/functions.rst b/doc/en/reference/functions.rst
a b 11 11 sage/functions/orthogonal_polys 12 12 sage/functions/other 13 13 sage/functions/special 14 sage/functions/bessel 14 15 sage/functions/exp_integral 15 16 sage/functions/wigner 16 17 sage/functions/generalized -
sage/functions/all.py
diff --git a/sage/functions/all.py b/sage/functions/all.py
a b 26 26 27 27 from transcendental import (zeta, zetaderiv, zeta_symmetric, dickman_rho) 28 28 29 from special import (bessel_I, bessel_J, bessel_K, bessel_Y, 30 hypergeometric_U, Bessel, 29 from bessel import (bessel_I, bessel_J, bessel_K, bessel_Y, Bessel) 30 31 from special import (hypergeometric_U, 31 32 spherical_bessel_J, spherical_bessel_Y, 32 33 spherical_hankel1, spherical_hankel2, 33 34 spherical_harmonic, jacobi, … … 36 37 elliptic_f, elliptic_ec, elliptic_eu, 37 38 elliptic_kc, elliptic_pi, elliptic_j, 38 39 airy_ai, airy_bi) 39 40 40 41 from orthogonal_polys import (chebyshev_T, 41 42 chebyshev_U, 42 43 gen_laguerre, -
new file sage/functions/bessel.py
diff --git a/sage/functions/bessel.py b/sage/functions/bessel.py new file mode 100644
- + 1 r""" 2 Bessel Functions 3 4 This module provides symbolic Bessel Functions. These functions use the 5 `mpmath Library`_ for numerical evaluation. 6 7 - Bessel functions, first defined by the Swiss mathematician 8 Daniel Bernoulli and named after Friedrich Bessel, are canonical 9 solutions y(x) of Bessel's differential equation: 10 11 .. math:: 12 13 x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \alpha^2)y = 0, 14 15 for an arbitrary real number `\alpha` (the order). 16 17 - In this module, `J_\alpha` denotes the unqiue solution of Bessel's equation 18 which is non-singular at `x = 0`. This function is known as the Bessel 19 Function of the First Kind. This function also arrises as a special case 20 of the hypergeometric function `{}_0F_1`: 21 22 .. math:: 23 24 J_\alpha(x) = \frac{x^n}{2^\alpha \Gamma(\alpha + 1)} {}_0F_1(\alpha + 1, -\frac{x^2}{4}). 25 26 - The second linearly independent solution to Bessel's equation (which is 27 singular at `x=0`) is denoted by `Y_\alpha` and is called the Bessel Function 28 of the Second Kind: 29 30 .. math:: 31 32 Y_\alpha(x) = \frac{ J_\alpha(x) \cos(\pi \alpha) - J_{-\alpha}(x)}{\sin(\pi \alpha)}. 33 34 - There are also two commonly used combinations of the Bessel J and Y Functions. 35 The Bessel I Function, or the Modified Bessel Function of the First Kind, is 36 defined by: 37 38 .. math:: 39 40 I_\alpha(x) = i^{-\alpha} J_\alpha(ix). 41 42 The Bessel K Function, or the Modified Bessel Function of the Second Kind, 43 is deinfed by: 44 45 .. math:: 46 47 K_\alpha(x) = \frac{\pi}{2} \cdot \frac{I_{-\alpha}(x) - I_n(x)}{\sin(\pi \alpha)}. 48 49 We should note here that the above formulas for Bessel Y and K functions 50 should be understood as limits when `\alpha` is an integer. 51 52 - It follows from Bessel's differential equation that the derivative of 53 `J_n(x)` with respect to `x` is: 54 55 ..math:: 56 57 \frac{d}{dx} J_n(x) = \frac{1}{x^n} \left\(x^n J_{n-1}(x) - n x^{n-1} J_n(z) \right\) 58 59 - Another important formulation of the two linearly independent 60 solutions to Bessel's equation are the Hankel functions 61 `H_\alpha^{(1)}(x)` and `H_\alpha^{(2)}(x)`, 62 defined by: 63 64 .. math:: 65 66 H_\alpha^{(1)}(x) = J_\alpha(x) + i Y_\alpha(x) 67 68 .. math:: 69 70 H_\alpha^{(2)}(x) = J_\alpha(x) - i Y_\alpha(x) 71 72 where `i` is the imaginary unit (and `J_*` and 73 `Y_*` are the usual J- and Y-Bessel functions). These 74 linear combinations are also known as Bessel functions of the third 75 kind; they are two linearly independent solutions of Bessel's 76 differential equation. They are named for Hermann Hankel. 77 78 EXAMPLES: 79 80 Evaluate the Bessel J function symbolically and numerically:: 81 82 sage: bessel_J(0, x) 83 bessel_J(0, x) 84 sage: bessel_J(0, 0) 85 bessel_J(0, 0) 86 sage: bessel_J(0, x).diff(x) 87 1/2*bessel_J(-1, x) - 1/2*bessel_J(1, x) 88 89 sage: N(bessel_J(0, 0), digits = 20) 90 1.0000000000000000000 91 sage: find_root(bessel_J(0,x), 0, 5) 92 2.404825557695773 93 94 Plot the Bessel J function:: 95 96 sage: f(x) = Bessel(0)(x); f 97 x |--> bessel_J(0, x) 98 sage: plot(f, (x, 1, 10)) 99 100 Visualize the Bessel Y function on the complex plane:: 101 102 sage: complex_plot(bessel_Y(0, x), (-5, 5), (-5, 5)) 103 104 Evaluate a combination of Bessel functions:: 105 106 sage: f(x) = bessel_J(1, x) - bessel_Y(0, x) 107 sage: f(pi) 108 bessel_J(1, pi) - bessel_Y(0, pi) 109 sage: f(pi).n() 110 -0.0437509653365599 111 sage: f(pi).n(digits=50) 112 -0.043750965336559909054985168023342675387737118378169 113 114 Symbolically solve a second order differential equation with initial 115 conditions `y(1) = a` and `y'(1) = b` in terms of Bessel functions:: 116 117 sage: y = function('y', x) 118 sage: a, b = var('a, b') 119 sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0 120 sage: f = desolve(diffeq, y, [1, a, b]); f 121 (a*bessel_Y(1, 1) + b*bessel_Y(0, 1))*bessel_J(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) - (a*bessel_J(1, 1) + b*bessel_J(0, 1))*bessel_Y(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) 122 123 For more examples, see the docstring for :meth:`Bessel`. 124 125 AUTHORS: 126 127 - Benjamin Jones (2012-12-27): initial version. 128 129 - Some of the documentation here has been adapted from David Joyner's 130 original documentation of Sage's special functions module (2006). 131 132 REFERENCES: 133 134 - Abramowitz and Stegun: Handbook of Mathematical Functions, 135 http://www.math.sfu.ca/~cbm/aands/ 136 137 - http://en.wikipedia.org/wiki/Bessel_function 138 139 - mpmath Library `Bessel Functions`_ 140 141 .. _`mpmath Library`: http://code.google.com/p/mpmath/ 142 .. _`Bessel Functions`: http://mpmath.googlecode.com/svn/trunk/doc/build/functions/bessel.html 143 144 """ 145 146 #***************************************************************************** 147 # Copyright (C) 2012 Benjamin Jones <benjaminfjones@gmail.com> 148 # 149 # Distributed under the terms of the GNU General Public License (GPL) 150 # 151 # This code is distributed in the hope that it will be useful, 152 # but WITHOUT ANY WARRANTY; without even the implied warranty of 153 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 154 # General Public License for more details. 155 # 156 # The full text of the GPL is available at: 157 # 158 # http://www.gnu.org/licenses/ 159 #***************************************************************************** 160 161 from sage.rings.all import RR, Integer 162 from sage.misc.latex import latex 163 from sage.symbolic.function import BuiltinFunction, is_inexact 164 from sage.symbolic.expression import Expression 165 from sage.structure.coerce import parent 166 import sage.structure.element 167 from sage.calculus.calculus import maxima 168 from sage.libs.mpmath import utils as mpmath_utils 169 from sage.functions.other import sqrt 170 from sage.functions.log import exp 171 from sage.functions.hyperbolic import sinh, cosh 172 from sage.symbolic.constants import pi 173 174 class Function_Bessel_J(BuiltinFunction): 175 r""" 176 The Bessel J Function, denoted by bessel_J(`\alpha`, x) or `J_\alpha(x)`. 177 As a Taylor series about `x=0` it is equal to: 178 179 .. math:: 180 181 J_\alpha(x) = \sum_{k=0}^\infty \frac{(-1)^k}{k! \Gamma(k+\alpha+1)} (\frac{x}{2})^{2k+\alpha} 182 183 For integer orders `\alpha = n` there is an integral representation: 184 185 .. math:: 186 187 J_n(x) = \frac{1}{\pi} \int_0^\pi \cos(n t - x \sin(t)) \; dt 188 189 This function also arrises as a special case of the hypergeometric 190 function `{}_0F_1`: 191 192 .. math:: 193 194 J_\alpha(x) = \frac{x^n}{2^\alpha \Gamma(\alpha + 1)} {}_0F_1(\alpha + 1, -\frac{x^2}{4}). 195 196 EXAMPLES:: 197 198 sage: bessel_J(1, x) 199 bessel_J(1, x) 200 sage: bessel_J(1.0, 1.0) 201 0.440050585744933 202 sage: n = var('n') 203 sage: bessel_J(n, x) 204 bessel_J(n, x) 205 sage: bessel_J(2, I).n(digits=30) 206 -0.135747669767038281182852569995 207 208 Examples of symbolic manipulation:: 209 210 sage: a = bessel_J(pi, bessel_J(1, I)) 211 sage: N(a, digits=20) 212 0.00059023706363796717363 - 0.0026098820470081958110*I 213 214 sage: f = bessel_J(2, x) 215 sage: f.diff(x) 216 1/2*bessel_J(1, x) - 1/2*bessel_J(3, x) 217 218 Currently, integration is not supported (directly) since we cannot yet convert 219 hypergeometric functions to and from Maxima:: 220 221 sage: f = bessel_J(2, x) 222 sage: f.integrate(x) 223 Traceback (most recent call last): 224 ... 225 TypeError: cannot coerce arguments: no canonical coercion from <type 'list'> to Symbolic Ring 226 227 sage: m = maxima(bessel_J(2, x)) 228 sage: m.integrate(x) 229 hypergeometric([3/2],[5/2,3],-x^2/4)*x^3/24 230 231 Visualization:: 232 233 sage: plot(bessel_J(1,x), (x,0,5), color='blue') 234 sage: complex_plot(bessel_J(1, x), (-5, 5), (-5, 5)) # long time 235 236 ALGORITHM: 237 238 Numerical evaluation is handled by the mpmath library. Symbolics are 239 handled by a combination of Maxima and Sage (Ginac/Pynac). 240 """ 241 def __init__(self): 242 """ 243 See the docstring for :meth:`Function_Bessel_J`. 244 245 EXAMPLES:: 246 247 sage: sage.functions.bessel.Function_Bessel_J() 248 bessel_J 249 """ 250 BuiltinFunction.__init__(self, "bessel_J", nargs=2, 251 conversions=dict(maxima='bessel_j', mathematica='BesselJ')) 252 253 def _eval_(self, n, x): 254 """ 255 EXAMPLES:: 256 257 sage: a, b = var('a, b') 258 sage: bessel_J(a, b) 259 bessel_J(a, b) 260 sage: bessel_J(1.0, 1.0) 261 0.440050585744933 262 """ 263 if not isinstance(n, Expression) and not isinstance(x, Expression) and \ 264 (is_inexact(n) or is_inexact(x)): 265 coercion_model = sage.structure.element.get_coercion_model() 266 n, x = coercion_model.canonical_coercion(n, x) 267 return self._evalf_(n, x, parent(n)) 268 269 return None # leaves the expression unevaluated 270 271 def _evalf_(self, n, x, parent=None): 272 """ 273 EXAMPLES:: 274 275 sage: bessel_J(0.0, 1.0) 276 0.765197686557966 277 sage: bessel_J(0, 1).n(digits=20) 278 0.76519768655796655145 279 """ 280 import mpmath 281 return mpmath_utils.call(mpmath.besselj, n, x, parent=parent) 282 283 def _derivative_(self, n, x, diff_param=None): 284 """ 285 Return the derivative of the Bessel J function. 286 287 EXAMPLES:: 288 289 sage: f(z) = bessel_J(10, z) 290 sage: derivative(f, z) 291 z |--> 1/2*bessel_J(9, z) - 1/2*bessel_J(11, z) 292 """ 293 return (bessel_J(n-1, x) - bessel_J(n+1, x))/Integer(2) 294 295 def _print_latex_(self, n, z): 296 """ 297 Custom _print_latex_ method. 298 299 EXAMPLES:: 300 301 sage: latex(bessel_J(1, x)) 302 \operatorname{J_{1}}(x) 303 """ 304 return r"\operatorname{J_{%s}}(%s)" % (latex(n), latex(z)) 305 306 bessel_J = Function_Bessel_J() 307 308 class Function_Bessel_Y(BuiltinFunction): 309 r""" 310 The Bessel Y function. 311 312 DEFINITION: 313 314 .. math: 315 316 bessel\_Y(n, z) = Y_n(z) = ... 317 318 The derivative with respect to `z` is: 319 320 ..math:: 321 322 \frac{d}{dz} Y_n(z) = \frac{1}{z^n} \left\(z^n Y_{n-1}(z) - n z^{n-1} Y_n(z) \right\) 323 324 EXAMPLES:: 325 326 sage: bessel_Y(1, x) 327 bessel_Y(1, x) 328 sage: bessel_Y(1.0, 1.0) 329 -0.781212821300289 330 sage: n = var('n') 331 sage: bessel_Y(n, x) 332 bessel_Y(n, x) 333 sage: bessel_Y(2, I).n() 334 1.03440456978312 - 0.135747669767038*I 335 336 Examples of symbolic manipulation:: 337 338 sage: a = bessel_Y(pi, bessel_Y(1, I)); a 339 bessel_Y(pi, bessel_Y(1, I)) 340 sage: N(a, digits=20) 341 4.2059146571791095708 + 21.307914215321993526*I 342 343 sage: f = bessel_Y(2, x) 344 sage: f.diff(x) 345 1/2*bessel_Y(1, x) - 1/2*bessel_Y(3, x) 346 347 High preceision and complex valued inputs (fixes Trac issue #4230):: 348 349 sage: bessel_Y(0, 1).n(128) 350 0.088256964215676957982926766023515162828 351 sage: bessel_Y(0, RealField(200)(1)) 352 0.088256964215676957982926766023515162827817523090675546711044 353 sage: bessel_Y(0, ComplexField(200)(0.5+I)) 354 0.077763160184438051408593468823822434235010300228009867784073 + 1.0142336049916069152644677682828326441579314239591288411739*I 355 356 Visualization:: 357 358 sage: plot(bessel_Y(1,x), (x,0,5), color='blue') 359 sage: complex_plot(bessel_Y(1, x), (-5, 5), (-5, 5)) # long time 360 361 ALGORITHM: 362 363 Numerical evaluation is handled by the mpmath library. Symbolics are 364 handled by a combination of Maxima and Sage (Ginac/Pynac). 365 """ 366 def __init__(self): 367 """ 368 See the docstring for :meth:`Function_Bessel_Y`. 369 370 EXAMPLES:: 371 372 sage: sage.functions.bessel.Function_Bessel_Y()(0, x) 373 bessel_Y(0, x) 374 """ 375 BuiltinFunction.__init__(self, "bessel_Y", nargs=2, 376 conversions=dict(maxima='bessel_y', mathematica='BesselY')) 377 378 def _eval_(self, n, x): 379 """ 380 EXAMPLES:: 381 382 sage: a,b = var('a, b') 383 sage: bessel_Y(a, b) 384 bessel_Y(a, b) 385 sage: bessel_Y(0, 1).n(128) 386 0.088256964215676957982926766023515162828 387 """ 388 if not isinstance(n, Expression) and not isinstance(x, Expression) and \ 389 (is_inexact(n) or is_inexact(x)): 390 coercion_model = sage.structure.element.get_coercion_model() 391 n, x = coercion_model.canonical_coercion(n, x) 392 return self._evalf_(n, x, parent(n)) 393 394 return None # leaves the expression unevaluated 395 396 def _evalf_(self, n, x, parent=None): 397 """ 398 EXAMPLES:: 399 400 sage: bessel_Y(1.0+2*I, 3.0+4*I) 401 0.699410324467538 + 0.228917940896421*I 402 sage: bessel_Y(0, 1).n(256) 403 0.08825696421567695798292676602351516282781752309067554671104384761199978932351 404 """ 405 import mpmath 406 return mpmath_utils.call(mpmath.bessely, n, x, parent=parent) 407 408 def _derivative_(self, n, x, diff_param=None): 409 """ 410 Return the derivative of the Bessel Y function. 411 412 EXAMPLES:: 413 414 sage: f(x) = bessel_Y(10, x) 415 sage: derivative(f, x) 416 x |--> 1/2*bessel_Y(9, x) - 1/2*bessel_Y(11, x) 417 """ 418 return (bessel_Y(n-1, x) - bessel_Y(n+1, x))/Integer(2) 419 420 def _print_latex_(self, n, z): 421 """ 422 Custom _print_latex_ method. 423 424 EXAMPLES:: 425 426 sage: latex(bessel_Y(1, x)) 427 \operatorname{Y_{1}}(x) 428 """ 429 return r"\operatorname{Y_{%s}}(%s)" % (latex(n), latex(z)) 430 431 bessel_Y = Function_Bessel_Y() 432 433 class Function_Bessel_I(BuiltinFunction): 434 r""" 435 The Bessel I function, or the Modified Bessel Function of the First Kind. 436 437 DEFINITION: 438 439 .. math:: 440 441 \operatorname{bessel\_I}(\alpha, x) = I_\alpha(x) = i^{-\alpha} J_\alpha(ix) 442 443 EXAMPLES:: 444 445 sage: bessel_I(1, x) 446 bessel_I(1, x) 447 sage: bessel_I(1.0, 1.0) 448 0.565159103992485 449 sage: n = var('n') 450 sage: bessel_I(n, x) 451 bessel_I(n, x) 452 sage: bessel_I(2, I).n() 453 -0.114903484931900 454 455 Examples of symbolic manipulation:: 456 457 sage: a = bessel_I(pi, bessel_I(1, I)) 458 sage: N(a, digits=20) 459 0.00026073272117205890528 - 0.0011528954889080572266*I 460 461 sage: f = bessel_I(2, x) 462 sage: f.diff(x) 463 1/2*bessel_I(1, x) + 1/2*bessel_I(3, x) 464 465 Special identities that bessel_I satisfies:: 466 467 sage: bessel_I(1/2, x) 468 sqrt(1/(pi*x))*sqrt(2)*sinh(x) 469 sage: eq = bessel_I(1/2, x) == bessel_I(0.5, x) 470 sage: eq.test_relation() 471 True 472 sage: bessel_I(-1/2, x) 473 sqrt(1/(pi*x))*sqrt(2)*cosh(x) 474 sage: eq = bessel_I(-1/2, x) == bessel_I(-0.5, x) 475 sage: eq.test_relation() 476 True 477 478 High preceision and complex valued inputs:: 479 480 sage: bessel_I(0, 1).n(128) 481 1.2660658777520083355982446252147175376 482 sage: bessel_I(0, RealField(200)(1)) 483 1.2660658777520083355982446252147175376076703113549622068081 484 sage: bessel_I(0, ComplexField(200)(0.5+I)) 485 0.80644357583493619472428518415019222845373366024179916785502 + 0.22686958987911161141397453401487525043310874687430711021434*I 486 487 Visualization:: 488 489 sage: plot(bessel_I(1,x), (x,0,5), color='blue') 490 sage: complex_plot(bessel_I(1, x), (-5, 5), (-5, 5)) # long time 491 492 ALGORITHM: 493 494 Numerical evaluation is handled by the mpmath library. Symbolics are 495 handled by a combination of Maxima and Sage (Ginac/Pynac). 496 """ 497 def __init__(self): 498 """ 499 See the docstring for :meth:`Function_Bessel_I`. 500 501 EXAMPLES:: 502 503 sage: bessel_I(1,x) 504 bessel_I(1, x) 505 """ 506 BuiltinFunction.__init__(self, "bessel_I", nargs=2, 507 conversions=dict(maxima='bessel_i', mathematica='BesselI')) 508 509 def _eval_(self, n, x): 510 """ 511 EXAMPLES:: 512 513 sage: y=var('y') 514 sage: bessel_I(y,x) 515 bessel_I(y, x) 516 sage: bessel_I(0.0, 1.0) 517 1.26606587775201 518 sage: bessel_I(1/2, 1) 519 sqrt(2)*sinh(1)/sqrt(pi) 520 sage: bessel_I(-1/2, pi) 521 sqrt(2)*cosh(pi)/pi 522 """ 523 if not isinstance(n, Expression) and not isinstance(x, Expression) and \ 524 (is_inexact(n) or is_inexact(x)): 525 coercion_model = sage.structure.element.get_coercion_model() 526 n, x = coercion_model.canonical_coercion(n, x) 527 return self._evalf_(n, x, parent(n)) 528 529 # special identities 530 if n == Integer(1)/Integer(2): 531 return sqrt(2/(pi*x))*sinh(x) 532 elif n == -Integer(1)/Integer(2): 533 return sqrt(2/(pi*x))*cosh(x) 534 535 return None # leaves the expression unevaluated 536 537 def _evalf_(self, n, x, parent=None): 538 """ 539 EXAMPLES:: 540 541 sage: bessel_I(1,3).n(digits=20) 542 3.9533702174026093965 543 """ 544 import mpmath 545 return mpmath_utils.call(mpmath.besseli, n, x, parent=parent) 546 547 def _derivative_(self, n, x, diff_param=None): 548 """ 549 Return the derivative of the Bessel I function `I_n(x)` with repect 550 to `x`. 551 552 EXAMPLES:: 553 554 sage: f(z) = bessel_I(10, x) 555 sage: derivative(f, x) 556 z |--> 1/2*bessel_I(9, x) + 1/2*bessel_I(11, x) 557 """ 558 return (bessel_I(n-1, x) + bessel_I(n+1, x))/Integer(2) 559 560 def _print_latex_(self, n, z): 561 """ 562 Custom _print_latex_ method. 563 564 EXAMPLES:: 565 566 sage: latex(bessel_I(1, x)) 567 \operatorname{I_{1}}(x) 568 """ 569 return r"\operatorname{I_{%s}}(%s)" % (latex(n), latex(z)) 570 571 bessel_I = Function_Bessel_I() 572 573 class Function_Bessel_K(BuiltinFunction): 574 r""" 575 The Bessel K function. 576 577 DEFINITION: 578 579 .. math: 580 581 \operatorname{bessel\_K}(\alpha, x) = K_\alpha(x) = ... 582 583 EXAMPLES:: 584 585 sage: bessel_K(1, x) 586 bessel_K(1, x) 587 sage: bessel_K(1.0, 1.0) 588 0.601907230197235 589 sage: n = var('n') 590 sage: bessel_K(n, x) 591 bessel_K(n, x) 592 sage: bessel_K(2, I).n() 593 -2.59288617549120 + 0.180489972066962*I 594 595 Examples of symbolic manipulation:: 596 597 sage: a = bessel_K(pi, bessel_K(1, I)); a 598 bessel_K(pi, bessel_K(1, I)) 599 sage: N(a, digits=20) 600 3.8507583115005220157 + 0.068528298579883425792*I 601 602 sage: f = bessel_K(2, x) 603 sage: f.diff(x) 604 1/2*bessel_K(1, x) + 1/2*bessel_K(3, x) 605 606 sage: bessel_K(1/2, x) 607 bessel_K(1/2, x) 608 sage: bessel_K(1/2, -1) 609 bessel_K(1/2, -1) 610 sage: bessel_K(1/2, 1) 611 sqrt(pi)*sqrt(1/2)*e^(-1) 612 613 High preceision and complex valued inputs:: 614 615 sage: bessel_K(0, 1).n(128) 616 0.42102443824070833333562737921260903614 617 sage: bessel_K(0, RealField(200)(1)) 618 0.42102443824070833333562737921260903613621974822666047229897 619 sage: bessel_K(0, ComplexField(200)(0.5+I)) 620 0.058365979093103864080375311643360048144715516692187818271179 - 0.67645499731334483535184142196073004335768129348518210260256*I 621 622 Visualization:: 623 624 sage: plot(bessel_K(1,x), (x,0,5), color='blue') 625 sage: complex_plot(bessel_K(1, x), (-5, 5), (-5, 5)) # long time 626 627 ALGORITHM: 628 629 Numerical evaluation is handled by the mpmath library. Symbolics are 630 handled by a combination of Maxima and Sage (Ginac/Pynac). 631 632 TESTS: 633 634 Verfify that Trac #3426 is fixed: 635 636 The Bessel K function can be evaluated numerically at complex orders:: 637 638 sage: bessel_K(10 * I, 10).n() 639 9.82415743819925e-8 640 641 For a fixed imaginary order and increasin, real, second component the 642 value of Bessel K is exponentially decaying:: 643 644 sage: for x in [10, 20, 50, 100, 200]: print bessel_K(5*I, x).n() 645 5.27812176514912e-6 646 3.11005908421801e-10 647 2.66182488515423e-23 - 8.59622057747552e-58*I 648 4.11189776828337e-45 - 1.01494840019482e-80*I 649 1.15159692553603e-88 - 6.75787862113718e-125*I 650 """ 651 def __init__(self): 652 """ 653 See the docstring for :meth:`Function_Bessel_K`. 654 655 EXAMPLES:: 656 657 sage: sage.functions.bessel.Function_Bessel_K() 658 bessel_K 659 """ 660 BuiltinFunction.__init__(self, "bessel_K", nargs=2, 661 conversions=dict(maxima='bessel_k', mathematica='BesselK')) 662 663 def _eval_(self, n, x): 664 """ 665 EXAMPLES:: 666 667 sage: bessel_K(1,0) 668 bessel_K(1, 0) 669 sage: bessel_K(1.0, 0.0) 670 +infinity 671 sage: bessel_K(-1, 1).n(128) 672 0.60190723019723457473754000153561733926 673 """ 674 if not isinstance(n, Expression) and not isinstance(x, Expression) and \ 675 (is_inexact(n) or is_inexact(x)): 676 coercion_model = sage.structure.element.get_coercion_model() 677 n, x = coercion_model.canonical_coercion(n, x) 678 return self._evalf_(n, x, parent(n)) 679 680 # special identity 681 if n == Integer(1)/Integer(2) and x > 0: 682 return sqrt(pi/2)*exp(-x)*x**(Integer(1)/Integer(2)) 683 684 return None # leaves the expression unevaluated 685 686 def _evalf_(self, n, x, parent=None): 687 """ 688 EXAMPLES:: 689 690 sage: bessel_K(0.0, 1.0) 691 0.421024438240708 692 sage: bessel_K(0, RealField(128)(1)) 693 0.42102443824070833333562737921260903614 694 """ 695 import mpmath 696 return mpmath_utils.call(mpmath.besselk, n, x, parent=parent) 697 698 def _derivative_(self, n, x, diff_param=None): 699 """ 700 Return the derivative of the Bessel K function. 701 702 EXAMPLES:: 703 704 sage: f(x) = bessel_K(10, x) 705 sage: derivative(f, x) 706 x |--> 1/2*bessel_K(9, x) + 1/2*bessel_K(11, x) 707 """ 708 return (bessel_K(n-1, x) + bessel_K(n+1, x))/Integer(2) 709 710 def _print_latex_(self, n, z): 711 """ 712 Custom _print_latex_ method. 713 714 EXAMPLES:: 715 716 sage: latex(bessel_K(1, x)) 717 \operatorname{K_{1}}(x) 718 """ 719 return r"\operatorname{K_{%s}}(%s)" % (latex(n), latex(z)) 720 721 bessel_K = Function_Bessel_K() 722 723 724 # dictionary used in Bessel 725 bessel_type_dict = {'I': bessel_I, 'J': bessel_J, 'K': bessel_K, 'Y': bessel_Y} 726 727 def Bessel(*args, **kwds): 728 """ 729 A factory function that produces symbolic I, J, K, and Y Bessel functions. There 730 are several ways to call this function: 731 732 - ``Bessel(order, type)`` 733 - ``Bessel(order)`` -- type defaults to 'J' 734 - ``Bessel(order, typ=T)`` 735 - ``Bessel(typ=T)`` -- order is unspecified, this is a 2-parameter function 736 - ``Bessel()`` -- order is unspecified, type is 'J' 737 738 where `order` can be any integer and T must be one of the strings 'I', 'J', 739 'K', or 'Y'. 740 741 See the EXAMPLES below. 742 743 EXAMPLES: 744 745 Construction of Bessel functions with various orders and types:: 746 747 sage: Bessel() 748 bessel_J 749 sage: Bessel(1)(x) 750 bessel_J(1, x) 751 sage: Bessel(1, 'Y')(x) 752 bessel_Y(1, x) 753 sage: Bessel(-2, 'Y')(x) 754 bessel_Y(-2, x) 755 sage: Bessel(typ='K') 756 bessel_K 757 sage: Bessel(0, typ='I')(x) 758 bessel_I(0, x) 759 760 Evaluation:: 761 762 sage: f = Bessel(1) 763 sage: f(3.0) 764 0.339058958525936 765 sage: f(3) 766 bessel_J(1, 3) 767 sage: f(3).n(digits=50) 768 0.33905895852593645892551459720647889697308041819801 769 770 sage: g = Bessel(typ='J') 771 sage: g(1,3) 772 bessel_J(1, 3) 773 sage: g(2, 3+I).n() 774 0.634160370148554 + 0.0253384000032695*I 775 sage: abs(numerical_integral(1/pi*cos(3*sin(x)), 0.0, pi)[0] - Bessel(0, 'J')(3.0)) < 1e-15 776 True 777 778 Symbolic calculus:: 779 780 sage: f(x) = Bessel(0, 'J')(x) 781 sage: derivative(f, x) 782 x |--> 1/2*bessel_J(-1, x) - 1/2*bessel_J(1, x) 783 sage: derivative(f, x, x) 784 x |--> 1/4*bessel_J(-2, x) - 1/2*bessel_J(0, x) + 1/4*bessel_J(2, x) 785 786 Verify that `J_0` satisfies Bessel's differential equation numerically 787 using the `test_relation()` method:: 788 789 sage: y = bessel_J(0, x) 790 sage: diffeq = x^2*derivative(y,x,x) + x*derivative(y,x) + x^2*y == 0 791 sage: diffeq.test_relation(proof=False) 792 True 793 794 Conversion to other systems:: 795 796 sage: x,y = var('x,y') 797 sage: f = maxima(Bessel(typ='K')(x,y)) 798 sage: f.derivative('x') 799 %pi*csc(%pi*x)*('diff(bessel_i(-x,y),x,1)-'diff(bessel_i(x,y),x,1))/2-%pi*bessel_k(x,y)*cot(%pi*x) 800 sage: f.derivative('y') 801 -(bessel_k(x+1,y)+bessel_k(x-1,y))/2 802 803 Compute the particular solution to Bessel's Differential Equation that 804 satisfies `y(1) = 1` and `y'(1) = 1`, then verify the initial conditions 805 and plot it:: 806 807 sage: y = function('y', x) 808 sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0 809 sage: f = desolve(diffeq, y, [1, 1, 1]); f 810 (bessel_J(0, 1) + bessel_J(1, 1))*bessel_Y(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1)) - (bessel_Y(0, 1) + bessel_Y(1, 1))*bessel_J(0, x)/(bessel_Y(0, 1)*bessel_J(1, 1) - bessel_Y(1, 1)*bessel_J(0, 1)) 811 812 sage: f.subs(x=1).n() # numerical verification 813 1.00000000000000 814 sage: fp = f.diff(x) 815 sage: fp.subs(x=1).n() 816 1.00000000000000 817 818 sage: f.subs(x=1).simplify_full() # symbolic verification 819 1 820 sage: fp = f.diff(x) 821 sage: fp.subs(x=1).simplify_full() 822 1 823 824 sage: plot(f, (x,0,5)) 825 826 Plotting:: 827 828 sage: f(x) = Bessel(0)(x); f 829 x |--> bessel_J(0, x) 830 sage: plot(f, (x, 1, 10)) 831 832 sage: plot([ Bessel(i, 'J') for i in range(5) ], 2, 10) 833 834 sage: G = Graphics() 835 sage: G += sum([ plot(Bessel(i), 0, 4*pi, rgbcolor=hue(sin(pi*i/10))) for i in range(5) ]) 836 sage: show(G) 837 838 A recreation of Abramowitz and Stegun Figure 9.1:: 839 840 sage: G = plot(Bessel(0, 'J'), 0, 15, color='black') 841 sage: G += plot(Bessel(0, 'Y'), 0, 15, color='black') 842 sage: G += plot(Bessel(1, 'J'), 0, 15, color='black', linestyle='dotted') 843 sage: G += plot(Bessel(1, 'Y'), 0, 15, color='black', linestyle='dotted') 844 sage: show(G, ymin=-1, ymax=1) 845 846 """ 847 # Determine the order and type of function from the arguments and keywords. 848 # These are recored in local variables: _type, _order, _system, _nargs. 849 _type = None 850 if len(args) == 0: # no order specified 851 _order = None 852 _nargs = 2 853 elif len(args) == 1: # order is specified 854 _order = args[0] 855 _nargs = 1 856 elif len(args) == 2: # both order and type are positional arguments 857 _order = args[0] 858 _type = args[1] 859 _nargs = 1 860 else: 861 raise TypeError("at most two positional arguments may be specified, " 862 +"see the docstring for Bessel") 863 # check for type inconsistency 864 if _type is not None and 'typ' in kwds and _type != kwds['typ']: 865 raise ValueError("inconsistent types given") 866 # record the function type 867 if _type is None: 868 if 'typ' in kwds: 869 _type = kwds['typ'] 870 else: 871 _type = 'J' 872 if not (_type in ['I', 'J', 'K', 'Y']): 873 raise ValueError("type must be one of I, J, K, Y") 874 # record the numerical evaluation system 875 if 'algorithm' in kwds: 876 _system = kwds['algorithm'] 877 else: 878 _system = 'mpmath' 879 880 # return the function 881 # TODO remove assertions 882 assert _type in ['I', 'J', 'K', 'Y'] 883 assert _order is None or _order in RR 884 assert _nargs == 1 or _nargs == 2 885 assert _system == 'mpmath' 886 887 # TODO what to do with _system? 888 _f = bessel_type_dict[_type] 889 if _nargs == 1: 890 return lambda x: _f(_order, x) 891 else: 892 return _f -
sage/functions/special.py
diff --git a/sage/functions/special.py b/sage/functions/special.py
a b 26 26 Toy. It is placed under the terms of the General Public License 27 27 (GPL) that governs the distribution of Maxima. 28 28 29 The (usual) Bessel functions and Airy functions are part of the30 standard Maxima package. Some Bessel functions also are implemented31 in PARI. (Caution: The PARI versions are sometimes different than32 the Maxima version.) For example, the J-Bessel function33 `J_\nu (z)` can be computed using either Maxima or PARI,34 depending on an optional variable you pass to bessel_J.35 36 29 Next, we summarize some of the properties of the functions 37 30 implemented here. 38 31 39 40 - Bessel functions, first defined by the Swiss mathematician41 Daniel Bernoulli and named after Friedrich Bessel, are canonical42 solutions y(x) of Bessel's differential equation:43 44 45 .. math::46 47 x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \alpha^2)y = 0,48 49 50 for an arbitrary real number `\alpha` (the order).51 52 - Another important formulation of the two linearly independent53 solutions to Bessel's equation are the Hankel functions54 `H_\alpha^{(1)}(x)` and `H_\alpha^{(2)}(x)`,55 defined by:56 57 58 .. math::59 60 H_\alpha^{(1)}(x) = J_\alpha(x) + i Y_\alpha(x)61 62 63 64 .. math::65 66 H_\alpha^{(2)}(x) = J_\alpha(x) - i Y_\alpha(x)67 68 69 where `i` is the imaginary unit (and `J_*` and70 `Y_*` are the usual J- and Y-Bessel functions). These71 linear combinations are also known as Bessel functions of the third72 kind; they are two linearly independent solutions of Bessel's73 differential equation. They are named for Hermann Hankel.74 75 32 - Airy function The function `Ai(x)` and the related 76 33 function `Bi(x)`, which is also called an Airy function, 77 34 are solutions to the differential equation 78 35 79 36 80 37 .. math:: 81 38 82 y'' - xy = 0, 39 y'' - xy = 0, 83 40 84 41 known as the Airy equation. They belong to the class of 'Bessel functions of 85 42 fractional order'. The initial conditions … … 93 50 94 51 - Spherical harmonics: Laplace's equation in spherical coordinates 95 52 is: 96 53 97 54 .. math:: 98 55 99 {\frac{1}{r^2}}{\frac{\partial}{\partial r}} \left(r^2 {\frac{\partial f}{\partial r}}\right) + {\frac{1}{r^2}\sin\theta}{\frac{\partial}{\partial \theta}} \left(\sin\theta {\frac{\partial f}{\partial \theta}}\right) + {\frac{1}{r^2\sin^2\theta}}{\frac{\partial^2 f}{\partial \varphi^2}} = 0. 56 {\frac{1}{r^2}}{\frac{\partial}{\partial r}} \left(r^2 {\frac{\partial f}{\partial r}}\right) + {\frac{1}{r^2}\sin\theta}{\frac{\partial}{\partial \theta}} \left(\sin\theta {\frac{\partial f}{\partial \theta}}\right) + {\frac{1}{r^2\sin^2\theta}}{\frac{\partial^2 f}{\partial \varphi^2}} = 0. 100 57 101 58 102 59 Note that the spherical coordinates `\theta` and … … 107 64 108 65 The general solution which remains finite towards infinity is a 109 66 linear combination of functions of the form 110 67 111 68 .. math:: 112 69 113 r^{-1-\ell} \cos (m \varphi) P_\ell^m (\cos{\theta} ) 70 r^{-1-\ell} \cos (m \varphi) P_\ell^m (\cos{\theta} ) 114 71 115 72 116 73 and 117 74 118 75 .. math:: 119 76 120 r^{-1-\ell} \sin (m \varphi) P_\ell^m (\cos{\theta} ) 77 r^{-1-\ell} \sin (m \varphi) P_\ell^m (\cos{\theta} ) 121 78 122 79 123 80 where `P_\ell^m` are the associated Legendre polynomials, … … 126 83 solutions with integer parameters `\ell \ge 0` and 127 84 `- \ell\leq m\leq \ell`, can be written as linear 128 85 combinations of: 129 86 130 87 .. math:: 131 88 132 U_{\ell,m}(r,\theta , \varphi ) = r^{-1-\ell} Y_\ell^m( \theta , \varphi ) 89 U_{\ell,m}(r,\theta , \varphi ) = r^{-1-\ell} Y_\ell^m( \theta , \varphi ) 133 90 134 91 135 92 where the functions `Y` are the spherical harmonic 136 93 functions with parameters `\ell`, `m`, which can be 137 94 written as: 138 95 139 96 .. math:: 140 97 141 Y_\ell^m( \theta , \varphi ) = \sqrt{{\frac{(2\ell+1)}{4\pi}}{\frac{(\ell-m)!}{(\ell+m)!}}} \cdot e^{i m \varphi } \cdot P_\ell^m ( \cos{\theta} ) . 98 Y_\ell^m( \theta , \varphi ) = \sqrt{{\frac{(2\ell+1)}{4\pi}}{\frac{(\ell-m)!}{(\ell+m)!}}} \cdot e^{i m \varphi } \cdot P_\ell^m ( \cos{\theta} ) . 142 99 143 100 144 101 145 102 The spherical harmonics obey the normalisation condition 146 103 147 104 148 105 .. math:: 149 106 150 \int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi} Y_\ell^mY_{\ell'}^{m'*}\,d\Omega =\delta_{\ell\ell'}\delta_{mm'}\quad\quad d\Omega =\sin\theta\,d\varphi\,d\theta . 107 \int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi} Y_\ell^mY_{\ell'}^{m'*}\,d\Omega =\delta_{\ell\ell'}\delta_{mm'}\quad\quad d\Omega =\sin\theta\,d\varphi\,d\theta . 151 108 152 109 153 110 154 111 - When solving for separable solutions of Laplace's equation in 155 112 spherical coordinates, the radial equation has the form: 156 113 157 114 .. math:: 158 115 159 x^2 \frac{d^2 y}{dx^2} + 2x \frac{dy}{dx} + [x^2 - n(n+1)]y = 0. 116 x^2 \frac{d^2 y}{dx^2} + 2x \frac{dy}{dx} + [x^2 - n(n+1)]y = 0. 160 117 161 118 162 119 The spherical Bessel functions `j_n` and `y_n`, 163 120 are two linearly independent solutions to this equation. They are 164 121 related to the ordinary Bessel functions `J_n` and 165 122 `Y_n` by: 166 123 167 124 .. math:: 168 125 169 j_n(x) = \sqrt{\frac{\pi}{2x}} J_{n+1/2}(x), 126 j_n(x) = \sqrt{\frac{\pi}{2x}} J_{n+1/2}(x), 170 127 171 128 172 129 173 130 .. math:: 174 131 175 y_n(x) = \sqrt{\frac{\pi}{2x}} Y_{n+1/2}(x) = (-1)^{n+1} \sqrt{\frac{\pi}{2x}} J_{-n-1/2}(x). 132 y_n(x) = \sqrt{\frac{\pi}{2x}} Y_{n+1/2}(x) = (-1)^{n+1} \sqrt{\frac{\pi}{2x}} J_{-n-1/2}(x). 176 133 177 134 178 135 … … 180 137 `y = U(a,b,x)` is defined to be the solution to Kummer's 181 138 differential equation 182 139 183 140 184 141 .. math:: 185 142 186 xy'' + (b-x)y' - ay = 0, 143 xy'' + (b-x)y' - ay = 0, 187 144 188 145 which satisfies `U(a,b,x) \sim x^{-a}`, as 189 146 `x\rightarrow \infty`. (There is a linearly independent … … 213 170 doubly-periodic, meromorphic functions satisfying the following 214 171 three properties: 215 172 216 173 217 174 #. There is a simple zero at the corner p, and a simple pole at the 218 175 corner q. 219 176 … … 232 189 is 1. 233 190 234 191 235 We can write 236 192 We can write 193 237 194 .. math:: 238 195 239 pq(x)=\frac{pr(x)}{qr(x)} 196 pq(x)=\frac{pr(x)}{qr(x)} 240 197 241 198 242 199 where `p`, `q`, and `r` are any of the … … 244 201 the understanding that `ss=cc=dd=nn=1`. 245 202 246 203 Let 247 204 248 205 .. math:: 249 206 250 u=\int_0^\phi \frac{d\theta} {\sqrt {1-m \sin^2 \theta}} 207 u=\int_0^\phi \frac{d\theta} {\sqrt {1-m \sin^2 \theta}} 251 208 252 209 253 210 Then the *Jacobi elliptic function* `sn(u)` is given by 254 211 255 212 .. math:: 256 213 257 {sn}\; u = \sin \phi 214 {sn}\; u = \sin \phi 258 215 259 216 and `cn(u)` is given by 260 217 261 218 .. math:: 262 219 263 {cn}\; u = \cos \phi 220 {cn}\; u = \cos \phi 264 221 265 222 and 266 223 267 224 .. math:: 268 225 269 {dn}\; u = \sqrt {1-m\sin^2 \phi}. 226 {dn}\; u = \sqrt {1-m\sin^2 \phi}. 270 227 271 228 To emphasize the dependence on `m`, one can write 272 229 `sn(u,m)` for example (and similarly for `cn` and … … 276 233 solutions to the following nonlinear ordinary differential 277 234 equations: 278 235 279 236 280 237 - `\mathrm{sn}\,(x;k)` solves the differential equations 281 238 282 239 .. math:: 283 240 284 241 \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1+k^2) y - 2 k^2 y^3 = 0, 285 242 286 243 and 287 244 288 245 .. math:: … … 291 248 292 249 - `\mathrm{cn}\,(x;k)` solves the differential equations 293 250 294 251 295 252 .. math:: 296 297 \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1-2k^2) y + 2 k^2 y^3 = 0, 253 254 \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} + (1-2k^2) y + 2 k^2 y^3 = 0, 298 255 299 256 and `\left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2 = (1-y^2) (1-k^2 + k^2 y^2)`. 300 257 301 258 - `\mathrm{dn}\,(x;k)` solves the differential equations 302 259 303 260 .. math:: 304 261 305 \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} - (2 - k^2) y + 2 y^3 = 0, 262 \frac{\mathrm{d}^2 y}{\mathrm{d}x^2} - (2 - k^2) y + 2 y^3 = 0, 306 263 307 264 and `\left(\frac{\mathrm{d} y}{\mathrm{d}x}\right)^2= y^2 (1 - k^2 - y^2)`. 308 265 … … 320 277 `dn(x, 1) = \mbox{\rm sech} x`. 321 278 322 279 - The incomplete elliptic integrals (of the first kind, etc.) are: 323 280 324 281 .. math:: 325 326 \begin{array}{c} \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2}}\, dx,\\ \displaystyle\int_0^\phi \sqrt{1 - m\sin(x)^2}\, dx,\\ \displaystyle\int_0^\phi \frac{\sqrt{1-mt^2}}{\sqrt(1 - t^2)}\, dx,\\ \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2\sqrt{1 - n\sin(x)^2}}}\, dx, \end{array} 282 283 \begin{array}{c} \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2}}\, dx,\\ \displaystyle\int_0^\phi \sqrt{1 - m\sin(x)^2}\, dx,\\ \displaystyle\int_0^\phi \frac{\sqrt{1-mt^2}}{\sqrt(1 - t^2)}\, dx,\\ \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2\sqrt{1 - n\sin(x)^2}}}\, dx, \end{array} 327 284 328 285 and the complete ones are obtained by taking `\phi =\pi/2`. 329 286 330 331 287 REFERENCES: 332 288 333 289 - Abramowitz and Stegun: Handbook of Mathematical Functions, 334 290 http://www.math.sfu.ca/~cbm/aands/ 335 291 336 - http://en.wikipedia.org/wiki/Bessel_function337 338 292 - http://en.wikipedia.org/wiki/Airy_function 339 293 340 294 - http://en.wikipedia.org/wiki/Spherical_harmonics … … 349 303 - Online Encyclopedia of Special Function 350 304 http://algo.inria.fr/esf/index.html 351 305 352 TODO: Resolve weird bug in commented out code in hypergeometric_U353 below.354 355 306 AUTHORS: 356 307 357 308 - David Joyner and William Stein … … 652 603 _init() 653 604 return RDF(meval("airy_bi(%s)"%RDF(x))) 654 605 655 656 def bessel_I(nu,z,algorithm = "pari",prec=53):657 r"""658 Implements the "I-Bessel function", or "modified Bessel function,659 1st kind", with index (or "order") nu and argument z.660 661 INPUT:662 663 664 - ``nu`` - a real (or complex, for pari) number665 666 - ``z`` - a real (positive) algorithm - "pari" or667 "maxima" or "scipy" prec - real precision (for PARI only)668 669 670 DEFINITION::671 672 Maxima:673 inf674 ==== - nu - 2 k nu + 2 k675 \ 2 z676 > -------------------677 / k! Gamma(nu + k + 1)678 ====679 k = 0680 681 PARI:682 683 inf684 ==== - 2 k 2 k685 \ 2 z Gamma(nu + 1)686 > -----------------------687 / k! Gamma(nu + k + 1)688 ====689 k = 0690 691 692 693 Sometimes ``bessel_I(nu,z)`` is denoted694 ``I_nu(z)`` in the literature.695 696 .. warning::697 698 In Maxima (the manual says) i0 is deprecated but699 ``bessel_i(0,*)`` is broken. (Was fixed in recent CVS patch700 though.)701 702 EXAMPLES::703 704 sage: bessel_I(1,1,"pari",500)705 0.565159103992485027207696027609863307328899621621092009480294489479255640964371134092664997766814410064677886055526302676857637684917179812041131208121706 sage: bessel_I(1,1)707 0.565159103992485708 sage: bessel_I(2,1.1,"maxima")709 0.16708949925104...710 sage: bessel_I(0,1.1,"maxima")711 1.32616018371265...712 sage: bessel_I(0,1,"maxima")713 1.2660658777520...714 sage: bessel_I(1,1,"scipy")715 0.565159103992...716 717 Check whether the return value is real whenever the argument is real (#10251)::718 719 sage: bessel_I(5, 1.5, algorithm='scipy') in RR720 True721 722 """723 if algorithm=="pari":724 from sage.libs.pari.all import pari725 try:726 R = RealField(prec)727 nu = R(nu)728 z = R(z)729 except TypeError:730 C = ComplexField(prec)731 nu = C(nu)732 z = C(z)733 K = C734 K = z.parent()735 return K(pari(nu).besseli(z, precision=prec))736 elif algorithm=="scipy":737 if prec != 53:738 raise ValueError, "for the scipy algorithm the precision must be 53"739 import scipy.special740 ans = str(scipy.special.iv(float(nu),complex(real(z),imag(z))))741 ans = ans.replace("(","")742 ans = ans.replace(")","")743 ans = ans.replace("j","*I")744 ans = sage_eval(ans)745 return real(ans) if z in RR else ans # Return real value when arg is real746 elif algorithm == "maxima":747 if prec != 53:748 raise ValueError, "for the maxima algorithm the precision must be 53"749 return sage_eval(maxima.eval("bessel_i(%s,%s)"%(float(nu),float(z))))750 else:751 raise ValueError, "unknown algorithm '%s'"%algorithm752 753 def bessel_J(nu,z,algorithm="pari",prec=53):754 r"""755 Return value of the "J-Bessel function", or "Bessel function, 1st756 kind", with index (or "order") nu and argument z.757 758 ::759 760 Defn:761 Maxima:762 inf763 ==== - nu - 2 k nu + 2 k764 \ (-1)^k 2 z765 > -------------------------766 / k! Gamma(nu + k + 1)767 ====768 k = 0769 770 PARI:771 772 inf773 ==== - 2k 2k774 \ (-1)^k 2 z Gamma(nu + 1)775 > ----------------------------776 / k! Gamma(nu + k + 1)777 ====778 k = 0779 780 781 Sometimes bessel_J(nu,z) is denoted J_nu(z) in the literature.782 783 .. warning::784 785 Inaccurate for small values of z.786 787 EXAMPLES::788 789 sage: bessel_J(2,1.1)790 0.136564153956658791 sage: bessel_J(0,1.1)792 0.719622018527511793 sage: bessel_J(0,1)794 0.765197686557967795 sage: bessel_J(0,0)796 1.00000000000000797 sage: bessel_J(0.1,0.1)798 0.777264368097005799 800 We check consistency of PARI and Maxima::801 802 sage: n(bessel_J(3,10,"maxima"))803 0.0583793793051...804 sage: n(bessel_J(3,10,"pari"))805 0.0583793793051868806 sage: bessel_J(3,10,"scipy")807 0.0583793793052...808 809 Check whether the return value is real whenever the argument is real (#10251)::810 sage: bessel_J(5, 1.5, algorithm='scipy') in RR811 True812 """813 814 if algorithm=="pari":815 from sage.libs.pari.all import pari816 try:817 R = RealField(prec)818 nu = R(nu)819 z = R(z)820 except TypeError:821 C = ComplexField(prec)822 nu = C(nu)823 z = C(z)824 K = C825 if nu == 0:826 nu = ZZ(0)827 K = z.parent()828 return K(pari(nu).besselj(z, precision=prec))829 elif algorithm=="scipy":830 if prec != 53:831 raise ValueError, "for the scipy algorithm the precision must be 53"832 import scipy.special833 ans = str(scipy.special.jv(float(nu),complex(real(z),imag(z))))834 ans = ans.replace("(","")835 ans = ans.replace(")","")836 ans = ans.replace("j","*I")837 ans = sage_eval(ans)838 return real(ans) if z in RR else ans839 elif algorithm == "maxima":840 if prec != 53:841 raise ValueError, "for the maxima algorithm the precision must be 53"842 return maxima_function("bessel_j")(nu, z)843 else:844 raise ValueError, "unknown algorithm '%s'"%algorithm845 846 def bessel_K(nu,z,algorithm="pari",prec=53):847 r"""848 Implements the "K-Bessel function", or "modified Bessel function,849 2nd kind", with index (or "order") nu and argument z. Defn::850 851 pi*(bessel_I(-nu, z) - bessel_I(nu, z))852 ----------------------------------------853 2*sin(pi*nu)854 855 856 if nu is not an integer and by taking a limit otherwise.857 858 Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In859 PARI, nu can be complex and z must be real and positive.860 861 EXAMPLES::862 863 sage: bessel_K(3,2,"scipy")864 0.64738539094...865 sage: bessel_K(3,2)866 0.64738539094...867 sage: bessel_K(1,1)868 0.60190723019...869 sage: bessel_K(1,1,"pari",10)870 0.60871 sage: bessel_K(1,1,"pari",100)872 0.60190723019723457473754000154873 874 TESTS::875 876 sage: bessel_K(2,1.1, algorithm="maxima")877 Traceback (most recent call last):878 ...879 NotImplementedError: The K-Bessel function is only implemented for the pari and scipy algorithms880 881 Check whether the return value is real whenever the argument is real (#10251)::882 883 sage: bessel_K(5, 1.5, algorithm='scipy') in RR884 True885 886 """887 if algorithm=="scipy":888 if prec != 53:889 raise ValueError, "for the scipy algorithm the precision must be 53"890 import scipy.special891 ans = str(scipy.special.kv(float(nu),float(z)))892 ans = ans.replace("(","")893 ans = ans.replace(")","")894 ans = ans.replace("j","*I")895 ans = sage_eval(ans)896 return real(ans) if z in RR else ans897 elif algorithm == 'pari':898 from sage.libs.pari.all import pari899 try:900 R = RealField(prec)901 nu = R(nu)902 z = R(z)903 except TypeError:904 C = ComplexField(prec)905 nu = C(nu)906 z = C(z)907 K = C908 K = z.parent()909 return K(pari(nu).besselk(z, precision=prec))910 elif algorithm == 'maxima':911 raise NotImplementedError, "The K-Bessel function is only implemented for the pari and scipy algorithms"912 else:913 raise ValueError, "unknown algorithm '%s'"%algorithm914 915 916 def bessel_Y(nu,z,algorithm="maxima", prec=53):917 r"""918 Implements the "Y-Bessel function", or "Bessel function of the 2nd919 kind", with index (or "order") nu and argument z.920 921 .. note::922 923 Currently only prec=53 is supported.924 925 Defn::926 927 cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)928 -------------------------------------------------929 sin(nu*pi)930 931 if nu is not an integer and by taking a limit otherwise.932 933 Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.934 935 This is computed using Maxima by default.936 937 EXAMPLES::938 939 sage: bessel_Y(2,1.1,"scipy")940 -1.4314714939...941 sage: bessel_Y(2,1.1)942 -1.4314714939590...943 sage: bessel_Y(3.001,2.1)944 -1.0299574976424...945 946 TESTS::947 948 sage: bessel_Y(2,1.1, algorithm="pari")949 Traceback (most recent call last):950 ...951 NotImplementedError: The Y-Bessel function is only implemented for the maxima and scipy algorithms952 """953 if algorithm=="scipy":954 if prec != 53:955 raise ValueError, "for the scipy algorithm the precision must be 53"956 import scipy.special957 ans = str(scipy.special.yv(float(nu),complex(real(z),imag(z))))958 ans = ans.replace("(","")959 ans = ans.replace(")","")960 ans = ans.replace("j","*I")961 ans = sage_eval(ans)962 return real(ans) if z in RR else ans963 elif algorithm == "maxima":964 if prec != 53:965 raise ValueError, "for the maxima algorithm the precision must be 53"966 return RR(maxima.eval("bessel_y(%s,%s)"%(float(nu),float(z))))967 elif algorithm == "pari":968 raise NotImplementedError, "The Y-Bessel function is only implemented for the maxima and scipy algorithms"969 else:970 raise ValueError, "unknown algorithm '%s'"%algorithm971 972 class Bessel():973 """974 A class implementing the I, J, K, and Y Bessel functions.975 976 EXAMPLES::977 978 sage: g = Bessel(2); g979 J_{2}980 sage: print g981 J-Bessel function of order 2982 sage: g.plot(0,10)983 984 ::985 986 sage: Bessel(2, typ='I')(pi)987 2.61849485263445988 sage: Bessel(2, typ='J')(pi)989 0.485433932631509990 sage: Bessel(2, typ='K')(pi)991 0.0510986902537926992 sage: Bessel(2, typ='Y')(pi)993 -0.0999007139289404994 """995 def __init__(self, nu, typ = "J", algorithm = None, prec = 53):996 """997 Initializes new instance of the Bessel class.998 999 INPUT:1000 1001 - ``typ`` -- (default: J) the type of Bessel function: 'I', 'J', 'K'1002 or 'Y'.1003 1004 - ``algorithm`` -- (default: maxima for type Y, pari for other types)1005 algorithm to use to compute the Bessel function: 'pari', 'maxima' or1006 'scipy'. Note that type K is not implemented in Maxima and type Y1007 is not implemented in PARI.1008 1009 - ``prec`` -- (default: 53) precision in bits of the Bessel function.1010 Only supported for the PARI algorithm.1011 1012 EXAMPLES::1013 1014 sage: g = Bessel(2); g1015 J_{2}1016 sage: Bessel(1,'I')1017 I_{1}1018 sage: Bessel(6, prec=120)(pi)1019 0.0145459669825055605736603696040018041020 sage: Bessel(6, algorithm="pari")(pi)1021 0.01454596698250561022 1023 For the Bessel J-function, Maxima returns a symbolic result. For1024 types I and Y, we always get a numeric result::1025 1026 sage: b = Bessel(6, algorithm="maxima")(pi); b1027 bessel_j(6, pi)1028 sage: b.n(53)1029 0.01454596698250561030 sage: Bessel(6, typ='I', algorithm="maxima")(pi)1031 0.02946198400595681032 sage: Bessel(6, typ='Y', algorithm="maxima")(pi)1033 -4.339328189390381034 1035 SciPy usually gives less precise results::1036 1037 sage: Bessel(6, algorithm="scipy")(pi)1038 0.0145459669825000...1039 1040 TESTS::1041 1042 sage: Bessel(1,'Z')1043 Traceback (most recent call last):1044 ...1045 ValueError: typ must be one of I, J, K, Y1046 """1047 if not (typ in ['I', 'J', 'K', 'Y']):1048 raise ValueError, "typ must be one of I, J, K, Y"1049 1050 # Did the user ask for the default algorithm?1051 if algorithm is None:1052 if typ == 'Y':1053 algorithm = 'maxima'1054 else:1055 algorithm = 'pari'1056 1057 self._system = algorithm1058 self._order = nu1059 self._type = typ1060 prec = int(prec)1061 if prec < 0:1062 raise ValueError, "prec must be a positive integer"1063 self._prec = int(prec)1064 1065 def __str__(self):1066 """1067 Returns a string representation of this Bessel object.1068 1069 TEST::1070 1071 sage: a = Bessel(1,'I')1072 sage: str(a)1073 'I-Bessel function of order 1'1074 """1075 return self.type()+"-Bessel function of order "+str(self.order())1076 1077 def __repr__(self):1078 """1079 Returns a string representation of this Bessel object.1080 1081 TESTS::1082 1083 sage: Bessel(1,'I')1084 I_{1}1085 """1086 return self.type()+"_{"+str(self.order())+"}"1087 1088 def type(self):1089 """1090 Returns the type of this Bessel object.1091 1092 TEST::1093 1094 sage: a = Bessel(3,'K')1095 sage: a.type()1096 'K'1097 """1098 return self._type1099 1100 def prec(self):1101 """1102 Returns the precision (in number of bits) used to represent this1103 Bessel function.1104 1105 TESTS::1106 1107 sage: a = Bessel(3,'K')1108 sage: a.prec()1109 531110 sage: B = Bessel(20,prec=100); B1111 J_{20}1112 sage: B.prec()1113 1001114 """1115 return self._prec1116 1117 def order(self):1118 """1119 Returns the order of this Bessel function.1120 1121 TEST::1122 1123 sage: a = Bessel(3,'K')1124 sage: a.order()1125 31126 """1127 return self._order1128 1129 def system(self):1130 """1131 Returns the package used, e.g. Maxima, PARI, or SciPy, to compute with1132 this Bessel function.1133 1134 TESTS::1135 1136 sage: Bessel(20,algorithm='maxima').system()1137 'maxima'1138 sage: Bessel(20,prec=100).system()1139 'pari'1140 """1141 return self._system1142 1143 def __call__(self,z):1144 """1145 Implements evaluation of all the Bessel functions directly1146 from the Bessel class. This essentially allows a Bessel object to1147 behave like a function that can be invoked.1148 1149 TESTS::1150 1151 sage: Bessel(3,'K')(5.0)1152 0.008291768415230931153 sage: Bessel(20,algorithm='maxima')(5.0)1154 2.77033005213e-111155 sage: Bessel(20,prec=100)(5.0101010101010101)1156 2.8809188227195382093062257967e-111157 sage: B = Bessel(2,'Y',algorithm='scipy',prec=50)1158 sage: B(2.0)1159 Traceback (most recent call last):1160 ...1161 ValueError: for the scipy algorithm the precision must be 531162 """1163 nu = self.order()1164 t = self.type()1165 s = self.system()1166 p = self.prec()1167 if t == "I":1168 return bessel_I(nu,z,algorithm=s,prec=p)1169 if t == "J":1170 return bessel_J(nu,z,algorithm=s,prec=p)1171 if t == "K":1172 return bessel_K(nu,z,algorithm=s,prec=p)1173 if t == "Y":1174 return bessel_Y(nu,z,algorithm=s,prec=p)1175 1176 def plot(self,a,b):1177 """1178 Enables easy plotting of all the Bessel functions directly1179 from the Bessel class.1180 1181 TESTS::1182 1183 sage: plot(Bessel(2),3,4)1184 sage: Bessel(2).plot(3,4)1185 sage: P = Bessel(2,'I').plot(1,5)1186 sage: P += Bessel(2,'J').plot(1,5)1187 sage: P += Bessel(2,'K').plot(1,5)1188 sage: P += Bessel(2,'Y').plot(1,5)1189 sage: show(P)1190 """1191 nu = self.order()1192 s = self.system()1193 t = self.type()1194 if t == "I":1195 f = lambda z: bessel_I(nu,z,s)1196 P = plot(f,a,b)1197 if t == "J":1198 f = lambda z: bessel_J(nu,z,s)1199 P = plot(f,a,b)1200 if t == "K":1201 f = lambda z: bessel_K(nu,z,s)1202 P = plot(f,a,b)1203 if t == "Y":1204 f = lambda z: bessel_Y(nu,z,s)1205 P = plot(f,a,b)1206 return P1207 1208 606 def hypergeometric_U(alpha,beta,x,algorithm="pari",prec=53): 1209 607 r""" 1210 608 Default is a wrap of PARI's hyperu(alpha,beta,x) function.