10.5 Special Functions

Module: sage.functions.special

Special Functions

Author Log:

This module provides easy access to many of Maxima and PARI's special functions.

Maxima's special functions package (which includes spherical harmonic functions, spherical Bessel functions (of the 1st and 2nd kind), and spherical Hankel functions (of the 1st and 2nd kind)) was written by Barton Willis of the University of Nebraska at Kearney. It is released under the terms of the General Public License (GPL).

Support for elliptic functions and integrals was written by Raymond Toy. It is placed under the terms of the General Public License (GPL) that governs the distribution of Maxima.

The (usual) Bessel functions and Airy functions are part of the standard Maxima package. Some Bessel functions also are implemented in Pari. (Caution: The Pari versions are sometimes different than the Maxima version.) For example, the K-Bessel function $ K_\nu (z)$ can be computed using either Maxima or Pari, depending on an optional variable you pass to bessel_K.

Next, we summarize some of the properties of the functions implemented here.

Methods implemented:
    * Bessel functions and Airy functions
    * spherical harmonic functions
    * spherical Bessel functions (of the 1st and 2nd kind)
    * spherical Hankel functions (of the 1st and 2nd kind)
    * Jacobi elliptic functions
    * complete/incomplete elliptic integrals
    * hyperbolic trig functions (for completeness, since
      they are special cases of elliptic functions)
    * Kummer confluent $U$-hypergeometric functions.
    
REFERENCE:
    * Abramowitz and Stegun: Handbook of Mathematical Functions,
      http://www.math.sfu.ca/~cbm/aands/
    * http://en.wikipedia.org/wiki/Bessel_function
    * http://en.wikipedia.org/wiki/Airy_function
    * http://en.wikipedia.org/wiki/Spherical_harmonics
    * http://en.wikipedia.org/wiki/Helmholtz_equation
    * http://en.wikipedia.org/wiki/Jacobi's_elliptic_functions
    * A. Khare, U. Sukhatme, "Cyclic Identities Involving
      Jacobi Elliptic Functions", Math ArXiv, math-ph/0201004
    * Online Encyclopedia of Special Function
      http://algo.inria.fr/esf/index.html

TODO: Resolve weird bug in commented out code in hypergeometric_U below.

Author: David Joyner and William Stein

Added 16-02-2008 (wdj): optional calls to scipy and replace all "#random" by "..." (both at the request of William Stein)

WARNING: SciPy's versions are poorly documented and seem less accurate than the Maxima and Pari versions.

Module-level Functions

bessel_I( nu, z, [algorithm=pari], [prec=53])

Implements the "I-Bessel function", or "modified Bessel function, 1st kind", with index (or "order") nu and argument z.

Input:

nu
- a real (or complex, for pari) number
z
- a real (positive) algorithm - "pari" or "maxima" or "scipy" prec - real precision (for Pari only)

DEFINITION:

    Maxima:
                     inf
                    ====   - nu - 2 k  nu + 2 k
                    \     2          z
                     >    -------------------
                    /     k! Gamma(nu + k + 1)
                    ====
                    k = 0

    Pari:
    
                     inf
                    ====   - 2 k  2 k
                    \     2      z    Gamma(nu + 1)
                     >    -----------------------
                    /       k! Gamma(nu + k + 1)
                    ====
                    k = 0
Sometimes bessel_I(nu,z) is denoted I_nu(z) in the literature.

WARNING: In Maxima (the manual says) i0 is deprecated but bessel_i(0,*) is broken. (Was fixed in recent CVS patch though.)

sage: bessel_I(1,1,"pari",500)
0.5651591039924850272076960276098633073288996216210920094802944894792556409
643711340926649977668144100646778860555263026768576376849171798120411312081
21
sage: bessel_I(1,1)
0.565159103992485
sage: bessel_I(2,1.1,"maxima")  
0.16708949925104...
sage: bessel_I(0,1.1,"maxima") 
1.32616018371265...
sage: bessel_I(0,1,"maxima")   
1.26606587775200...
sage: bessel_I(1,1,"scipy")
0.565159103992...

bessel_J( nu, z, [algorithm=pari], [prec=53])

Return value of the "J-Bessel function", or "Bessel function, 1st kind", with index (or "order") nu and argument z.

    Defn:
    Maxima:
                     inf
                    ====          - nu - 2 k  nu + 2 k
                    \     (-1)^k 2           z
                     >    -------------------------
                    /        k! Gamma(nu + k + 1)
                    ====
                    k = 0

    Pari:
    
                     inf
                    ====          - 2k    2k
                    \     (-1)^k 2      z    Gamma(nu + 1)
                     >    ----------------------------
                    /         k! Gamma(nu + k + 1)
                    ====
                    k = 0

Sometimes bessel_J(nu,z) is denoted J_nu(z) in the literature.

WARNING: Inaccurate for small values of z.

sage: bessel_J(2,1.1)
0.13656415395665...
sage: bessel_J(0,1.1)
0.71962201852751...
sage: bessel_J(0,1) 
0.76519768655796...

We check consistency of PARI and Maxima:

sage: n(bessel_J(3,10,"maxima"))
0.0583793793051...
sage: n(bessel_J(3,10,"pari"))  
0.0583793793051...
sage: bessel_J(3,10,"scipy")
0.0583793793052...

bessel_K( nu, z, [algorithm=pari], [prec=53])

Implements the "K-Bessel function", or "modified Bessel function, 2nd kind", with index (or "order") nu and argument z. Defn:

            pi*(bessel_I(-nu, z) - bessel_I(nu, z))
           ----------------------------------------
                        2*sin(pi*nu)

if nu is not an integer and by taking a limit otherwise.

Sometimes bessel_K(nu,z) is denoted K_nu(z) in the literature. In Pari, nu can be complex and z must be real and positive.

sage: bessel_K(3,2,"scipy")
0.64738539094...
sage: bessel_K(3,2)
0.64738539094...
sage: bessel_K(1,1)
0.60190723019...
sage: bessel_K(1,1,"pari",10)
0.60
sage: bessel_K(1,1,"pari",100)
0.60190723019723457473754000154

bessel_Y( nu, z, [algorithm=maxima], [prec=53])

Implements the "Y-Bessel function", or "Bessel function of the 2nd kind", with index (or "order") nu and argument z.

NOTE: Currently only prec=53 is supported.

Defn:

            cos(pi n)*bessel_J(nu, z) - bessel_J(-nu, z)
           -------------------------------------------------
                             sin(nu*pi)
if nu is not an integer and by taking a limit otherwise.

Sometimes bessel_Y(n,z) is denoted Y_n(z) in the literature.

This is computed using Maxima by default.

sage: bessel_Y(2,1.1,"scipy")
-1.4314714939...
sage: bessel_Y(2,1.1)   
-1.4314714939590...
sage: bessel_Y(3.001,2.1) 
-1.0299574976424...

NOTE: Adding "0"+ inside sage_eval as a temporary bug work-around.

elliptic_e( phi, m)

This returns the value of the incomplete elliptic integral of the second kind, $ \int_0^\phi \sqrt(1 - m\sin(x)^2)\, dx$ , ie, integrate(sqrt(1 - m*sin(x)2), x, 0, phi). Taking $ \phi =\pi/2$ gives elliptic_ec.

sage: z = var("z")
sage: elliptic_e (z, 1)
sin(z)
sage: elliptic_e (z, 0)
z
sage: elliptic_e (0.5, 0.1)
0.498011394498832

elliptic_ec( m)

This returns the value of the complete elliptic integral of the second kind, $ \int_0^{\pi/2} \sqrt(1 - m\sin(x)^2)\, dx$ .

sage: elliptic_ec (0.1)
1.530757636897763
sage: elliptic_ec (x).diff()
(elliptic_ec(x) - elliptic_kc(x))/(2*x)

elliptic_eu( u, m)

This returns the value of the incomplete elliptic integral of the second kind defined by $ \int_0^u jacobi_dn(x,m)^2)\, dx$ .

sage: elliptic_eu (0.5, 0.1)
0.496054551286597

elliptic_f( phi, m)

This returns the value of the "incomplete elliptic integral of the first kind", $ \int_0^\phi \frac{dx}{\sqrt{1 - m\sin(x)^2}}$ , ie, integrate(1/sqrt(1 - m*sin(x)2), x, 0, phi). Taking $ \phi =\pi/2$ gives elliptic_kc.

sage: z = var("z")
sage: elliptic_f (z, 0)
z
sage: elliptic_f (z, 1)
log(tan(z/2 + pi/4))
sage: elliptic_f (0.2, 0.1)
0.200132506747543

elliptic_kc( m)

This returns the value of the "complete elliptic integral of the first kind", $ \int_0^{\pi/2} \frac{dx}{\sqrt{1 - m\sin(x)^2}}$ .

sage: elliptic_kc (0.5)
1.854074677301372
sage: elliptic_f (RR(pi/2), 0.5)
1.854074677301378

elliptic_pi( n, phi, m)

This returns the value of the "incomplete elliptic integral of the third kind",

$\displaystyle \int_0^\phi \frac{dx}{\sqrt{(1 - m\sin(x)^2)(1 - n\sin(x)^2)}}$

.

sage: elliptic_pi(0.1, 0.2, 0.3)
0.200665068220979

error_fcn( t)

The complementary error function $ \frac{2}{\sqrt{\pi}}\int_t^\infty e^{-x^2} dx$ (t belongs to RR).

exp_int( t)

The exponential integral $ \int_x^\infty e^{-x}/x dx$ (t belongs to RR).

hypergeometric_U( alpha, beta, x, [algorithm=pari], [prec=53])

Default is a wrap of Pari's hyperu(alpha,beta,x) function. Optionally, algorithm = "scipy" can be used.

The confluent hypergeometric function $ y = U(a,b,x)$ is defined to be the solution to Kummer's differential equation

$\displaystyle xy'' + (b-x)y' - ay = 0.
$

This satisfies $ U(a,b,x) \sim x^{-a}$ , as $ x\rightarrow \infty$ , and is sometimes denoted x^{-a}2_F_0(a,1+a-b,-1/x). This is not the same as Kummer's $ M$ -hypergeometric function, denoted sometimes as _1F_1(alpha,beta,x), though it satisfies the same DE that $ U$ does.

WARNING: In the literature, both are called "Kummer confluent hypergeometric" functions.

sage: hypergeometric_U(1,1,1,"scipy")
0.596347362323...
sage: hypergeometric_U(1,1,1)  
0.59634736232319...
sage: hypergeometric_U(1,1,1,"pari",70) 
0.59634736232319407434...

inverse_jacobi( sym, x, m)

Here sym = "pq", where p,q in c,d,n,s. This returns the value of the inverse Jacobi function $ pq^{-1}(x,m)$ . There are a total of 12 functions described by this.

sage: jacobi("sn",1/2,1/2)
jacobi_sn(1/2, 1/2)
sage: float(jacobi("sn",1/2,1/2))
0.4707504736556572
sage: float(inverse_jacobi("sn",0.47,1/2))
0.4990982313222197
sage: float(inverse_jacobi("sn",0.4707504,0.5))
0.49999991146655459
sage: P = plot(inverse_jacobi('sn', x, 0.5), 0, 1, plot_points=20)

Now to view this, just type show(P).

jacobi( sym, x, m)

Here sym = "pq", where p,q in c,d,n,s. This returns the value of the Jacobi function pq(x,m), as described in the documentation for SAGE's "special" module. There are a total of 12 functions described by this.

sage: jacobi("sn",1,1)
tanh(1)
sage: jacobi("cd",1,1/2)
jacobi_cd(1, 1/2)
sage: RDF(jacobi("cd",1,1/2))
0.724009721659
sage: RDF(jacobi("cn",1,1/2)); RDF(jacobi("dn",1,1/2)); RDF(jacobi("cn",1,1/2)/jacobi("dn",1,1/2))
0.595976567672
0.823161001632
0.724009721659
sage: jsn = jacobi("sn",x,1)
sage: P = plot(jsn,0,1)

To view this, type P.show().

lngamma( t)

The principal branch of the logarithm of the Gamma function of t.

meval( x)

spherical_bessel_J( n, var, [algorithm=maxima])

Returns the spherical Bessel function of the first kind for integers n > -1.

Reference: A&S 10.1.8 page 437 and A&S 10.1.15 page 439.

sage: spherical_bessel_J(2,x)
(-(1 - 24/(8*x^2))*sin(x) - 3*cos(x)/x)/x

spherical_bessel_Y( n, var, [algorithm=maxima])

Returns the spherical Bessel function of the second kind for integers n > -1.

Reference: A&S 10.1.9 page 437 and A&S 10.1.15 page 439.

sage: x = PolynomialRing(QQ, 'x').gen()
sage: spherical_bessel_Y(2,x)
-(3*sin(x)/x - (1 - 24/(8*x^2))*cos(x))/x

spherical_hankel1( n, var)

Returns the spherical Hankel function of the first kind for integers $ n > -1$ , written as a string. Reference: A&S 10.1.36 page 439.

sage: spherical_hankel1(2,'x')
-3*I*(-x^2/3 - I*x + 1)*e^(I*x)/x^3

spherical_hankel2( n, x)

Returns the spherical Hankel function of the second kind for integers n > -1, written as a string. Reference: A&S 10.1.17 page 439.

sage: spherical_hankel2(2,'x')
'3*I*(-x^2/3+I*x+1)*%e^-(I*x)/x^3'

Here I = sqrt(-1).

spherical_harmonic( m, n, x, y)

Returns the spherical Harmonic function of the second kind for integers $ n > -1$ , $ \vert m\vert\leq n$ . Reference: Merzbacher 9.64.

sage: x,y = var('x,y') 
sage: spherical_harmonic(3,2,x,y)
15*sqrt(7)*cos(x)*sin(x)^2*e^(2*I*y)/(4*sqrt(30)*sqrt(pi))
sage: spherical_harmonic(3,2,1,2)
15*sqrt(7)*e^(4*I)*cos(1)*sin(1)^2/(4*sqrt(30)*sqrt(pi))

Class: Bessel

class Bessel
A class implementing the I, J, K, and Y Bessel functions.

sage: g = Bessel(2); g
J_{2}
sage: print g
J-Bessel function of order 2
sage: g.plot(0,10)
Bessel( self, nu, [typ=J], [algorithm=pari], [prec=53])

Functions: order,$ \,$ plot,$ \,$ prec,$ \,$ system,$ \,$ type

Special Functions: __call__,$ \,$ __init__,$ \,$ __repr__,$ \,$ __str__

See About this document... for information on suggesting changes.