# HG changeset patch
# User Alexandru Ghitza <aghitza@alum.mit.edu>
# Date 1222843216 36000
# Node ID ed46740a68d3602b822352abf6378c70281b2b50
# Parent 4805d6f2365704ea92de52d00a97375299950add
trac 3426: fix Bessel functions in special.py
diff r 4805d6f23657 r ed46740a68d3 sage/functions/special.py
a

b


380  380  index (or "order") nu and argument z. 
381  381  
382  382  INPUT: 
383   nu  a real (or complex, for pari) number 
384   z  a real (positive) 
385   algorithm  "pari" or "maxima" or "scipy" 
386   prec  real precision (for Pari only) 
 383  nu  a real or complex number 
 384  z  a real or complex number 
 385  algorithm  "pari" (default) or "maxima" or "scipy"; note 
 386  that "maxima" and "scipy" can only deal with real nu and 
 387  positive real z 
 388  prec  real precision (for Pari only) 
387  389  
388  390  DEFINITION: 
389  391  \begin{verbatim} 
… 
… 

428  430  sage: bessel_I(1,1,"scipy") 
429  431  0.565159103992... 
430  432  
 433  A corner case: 
 434  sage: bessel_I(0,0) 
 435  1.00000000000000 
 436  
 437  If both nu and z are real, then bessel_I is known to be real if z 
 438  is positive or if nu is an integer: 
 439  sage: bessel_I(2.1, 5.3) 
 440  22.5840990831073 
 441  sage: bessel_I(2.1, 5.3) 
 442  21.4787545976446  6.97887041932783*I 
 443  sage: bessel_I(2, 5.3) 
 444  23.5424857046048 
 445  sage: bessel_I(10*I, 10) 
 446  1.19390297337705e6  688460.429969423*I 
431  447  """ 
432  448  if algorithm=="pari": 
433  449  from sage.libs.pari.all import pari 
 450  R = RealField(prec) 
 451  # this is a workaround for a bug in Pari 2.3.3 
 452  if nu.is_zero() and z.is_zero(): 
 453  return R(1) 
 454  C = ComplexField(prec) 
434  455  try: 
435   R = RealField(prec) 
 456  nu = ZZ(nu) 
 457  # pari doesn't like negative integers here, 
 458  # but luckily I(nu, z) = I(nu, z) if nu is integral 
 459  nu = nu.abs() 
 460  K = R 
 461  except TypeError: 
 462  K = C 
 463  try: 
436  464  nu = R(nu) 
437  465  z = R(z) 
 466  if (z > 0): 
 467  K = R 
438  468  except TypeError: 
439   C = ComplexField(prec) 
440  469  nu = C(nu) 
441  470  z = C(z) 
442  471  K = C 
443   K = z.parent() 
444   return K(pari(nu).besseli(z, precision=prec)) 
 472  pari_bes = pari(nu).besseli(z, precision=prec) 
 473  if K is R: 
 474  return K(pari_bes.real()) 
 475  else: 
 476  return K(pari_bes) 
445  477  elif algorithm=="scipy": 
446  478  if prec != 53: 
447  479  raise ValueError, "for the scipy algorithm the precision must be 53" 
… 
… 

508  540  0.0583793793051868 
509  541  sage: bessel_J(3,10,"scipy") 
510  542  0.0583793793052... 
 543  
 544  A corner case: 
 545  sage: bessel_I(0,0) 
 546  1.00000000000000 
 547  
 548  If both nu and z are real, then bessel_J is known to be real if z 
 549  is positive or if nu is an integer: 
 550  sage: bessel_J(2.1, 5.3) 
 551  0.123120872446240 
 552  sage: bessel_J(2.1, 5.3) 
 553  0.117094908031941 + 0.0380464419481583*I 
 554  sage: bessel_J(2, 5.3) 
 555  0.0547481464525993 
 556  sage: bessel_J(10*I, 10) 
 557  119728.409042333  693848.987632644*I 
511  558  """ 
512  559  if algorithm=="pari": 
513  560  from sage.libs.pari.all import pari 
 561  R = RealField(prec) 
 562  # this is a workaround for a bug in Pari 2.3.3 
 563  if nu.is_zero() and z.is_zero(): 
 564  return R(1) 
 565  C = ComplexField(prec) 
 566  fudge = 1 
514  567  try: 
515   R = RealField(prec) 
 568  nu = ZZ(nu) 
 569  # pari doesn't like negative integers here, 
 570  # but luckily J(nu, z) = (1)**nu J(nu, z) if nu is integral 
 571  if nu < 0: 
 572  nu = nu 
 573  fudge = (1)**nu 
 574  K = R 
 575  except TypeError: 
 576  K = C 
 577  try: 
516  578  nu = R(nu) 
517  579  z = R(z) 
 580  if (z > 0): 
 581  K = R 
518  582  except TypeError: 
519   C = ComplexField(prec) 
520  583  nu = C(nu) 
521  584  z = C(z) 
522  585  K = C 
523   K = z.parent() 
524   return K(pari(nu).besselj(z, precision=prec)) 
 586  pari_bes = pari(nu).besselj(z, precision=prec) 
 587  if K is R: 
 588  return fudge*K(pari_bes.real()) 
 589  else: 
 590  return fudge*K(pari_bes) 
525  591  elif algorithm=="scipy": 
526  592  if prec != 53: 
527  593  raise ValueError, "for the scipy algorithm the precision must be 53" 
… 
… 

553  619  otherwise. 
554  620  
555  621  Sometimes bessel_K(nu,z) is denoted K_nu(z) in the 
556   literature. In Pari, nu can be complex and 
557   z must be real and positive. 
 622  literature. In Pari, both nu and z can be complex; in SciPy they 
 623  must be real and z mush be positive. 
558  624  
559  625  EXAMPLES: 
560  626  sage: bessel_K(3,2,"scipy") 
… 
… 

567  633  0.60 
568  634  sage: bessel_K(1,1,"pari",100) 
569  635  0.60190723019723457473754000154 
 636  
 637  sage: bessel_K(0,0) 
 638  +Infinity 
 639  
 640  If both nu and z are real, then bessel_K is known to be real if 
 641  z is positive: 
 642  sage: bessel_K(2.1, 5.3) 
 643  0.00388720082599840 
 644  sage: bessel_K(2.1, 5.3) 
 645  0.00369694767571369  70.9488385563182*I 
 646  sage: bessel_K(2, 5.3) 
 647  0.00375340283798982  73.9609001368291*I 
 648  sage: bessel_K(10*I, 10) 
 649  9.82415743819925e8 
 650  sage: bessel_K(2, 1.121) 
 651  1.23414145962983  0.547231658266306*I 
570  652  """ 
571  653  if algorithm=="scipy": 
572  654  if prec != 53: 
… 
… 

579  661  return sage_eval(ans) 
580  662  elif algorithm == 'pari': 
581  663  from sage.libs.pari.all import pari 
 664  R = RealField(prec) 
 665  C = ComplexField(prec) 
 666  if z.is_zero() and nu.is_zero(): 
 667  from sage.rings.infinity import Infinity 
 668  return Infinity 
 669  #TODO: sort out cases z = 0, nu != 0 
582  670  try: 
583   R = RealField(prec) 
584  671  nu = R(nu) 
 672  # K(nu, z) = K(nu, z) for all nu 
 673  nu = nu.abs() 
585  674  z = R(z) 
 675  if (z > 0): 
 676  K = R 
 677  else: 
 678  K = C 
586  679  except TypeError: 
587   C = ComplexField(prec) 
588  680  nu = C(nu) 
589  681  z = C(z) 
590  682  K = C 
591   K = z.parent() 
592   return K(pari(nu).besselk(z, precision=prec)) 
 683  pari_bes = pari(nu).besselk(z, precision=prec) 
 684  if K is R: 
 685  return K(pari_bes.real()) 
 686  else: 
 687  return K(pari_bes) 
593  688  else: 
594  689  raise ValueError, "unknown algorithm '%s'"%algorithm 
595  690  