Hypergeometric functions, Bessel/Airy, orthogonal polynomials

Bessel functions (besselj(), bessely())

mpmath.functions.besselj(n, x)

besselj(n,x) computes the Bessel function of the first kind J_n(x). Bessel functions of the first kind are defined as solutions of the differential equation

x^2 y'' + x y' + (x^2 - n^2) y = 0

which appears, among other things, when solving the radial part of Laplace’s equation in cylindrical coordinates. This equation has two solutions for given n, where the J_n-function is the solution that is nonsingular at x = 0. For positive integer n, J_n(x) behaves roughly like a sine (odd n) or cosine (even n) multiplied by a magnitude factor that decays slowly as x \to \pm\infty.

Generally, J_n is a special case of the hypergeometric function \,_0F_1:

J_n(x) = \frac{x^n}{2^n \Gamma(n+1)}
         \,_0F_1\left(n+1,-\frac{x^2}{4}\right)

Examples

Evaluation is supported for arbitrary arguments, and at arbitrary precision:

>>> from mpmath import *
>>> mp.dps = 15
>>> print besselj(2, 1000)
-0.024777229528606
>>> print besselj(4, 0.75)
0.000801070086542314
>>> print besselj(2, 1000j)
(-2.48071721019185e+432 + 0.0j)
>>> mp.dps = 25
>>> print besselj(0.75j, 3+4j)
(-2.778118364828153309919653 - 1.5863603889018621585533j)
>>> mp.dps = 50
>>> print besselj(1, pi)
0.28461534317975275734531059968613140570981118184947

The Bessel functions of the first kind satisfy simple symmetries around x = 0:

>>> mp.dps = 15
>>> nprint([besselj(n,0) for n in range(5)])
[1.0, 0.0, 0.0, 0.0, 0.0]
>>> nprint([besselj(n,pi) for n in range(5)])
[-0.304242, 0.284615, 0.485434, 0.333458, 0.151425]
>>> nprint([besselj(n,-pi) for n in range(5)])
[-0.304242, -0.284615, 0.485434, -0.333458, 0.151425]

Roots of Bessel functions are often used:

>>> nprint([findroot(j0, k) for k in [2, 5, 8, 11, 14]])
[2.40483, 5.52008, 8.65373, 11.7915, 14.9309]
>>> nprint([findroot(j1, k) for k in [3, 7, 10, 13, 16]])
[3.83171, 7.01559, 10.1735, 13.3237, 16.4706]

The roots are not periodic, but the distance between successive roots asymptotically approaches 2 \pi. Bessel functions of the first kind have the following normalization:

>>> print quadosc(j0, [0, inf], period=2*pi)
1.0
>>> print quadosc(j1, [0, inf], period=2*pi)
1.0

For n = 1/2 or n = -1/2, the Bessel function reduces to a trigonometric function:

>>> x = 10
>>> print besselj(0.5, x), sqrt(2/(pi*x))*sin(x)
-0.13726373575505 -0.13726373575505
>>> print besselj(-0.5, x), sqrt(2/(pi*x))*cos(x)
-0.211708866331398 -0.211708866331398
mpmath.functions.j0(x)
Computes the Bessel function J_0(x). See besselj().
mpmath.functions.j1(x)
Computes the Bessel function J_1(x). See besselj().
mpmath.functions.bessely(n, x)

bessely(n,x) computes the Bessel function of the second kind,

Y_n(x) = \frac{J_n(x) \cos(\pi n) - J_{-n}(x)}{\sin(\pi n)}.

For n an integer, this formula should be understood as a limit.

Examples

Some values of Y_n(x):

>>> from mpmath import *
>>> mp.dps = 25
>>> print bessely(0,0), bessely(1,0), bessely(2,0)
-inf -inf -inf
>>> print bessely(1, pi)
0.3588729167767189594679827
>>> print bessely(0.5, 3+4j)
(9.242861436961450520325216 - 3.085042824915332562522402j)

Modified Bessel functions (besseli(), besselk())

mpmath.functions.besseli(n, x)

besseli(n,x) computes the modified Bessel function of the first kind,

I_n(x) = i^{-n} J_n(ix)

Examples

Some values of I_n(x):

>>> from mpmath import *
>>> mp.dps = 25
>>> print besseli(0,0)
1.0
>>> print besseli(1,0)
0.0
>>> print besseli(0,1)
1.266065877752008335598245
>>> print besseli(3.5, 2+3j)
(-0.2904369752642538144289025 - 0.4469098397654815837307006j)

For integers n, the following integral representation holds:

>>> mp.dps = 15
>>> n = 3
>>> x = 2.3
>>> print quad(lambda t: exp(x*cos(t))*cos(n*t), [0,pi])/pi
0.349223221159309
>>> print besseli(n,x)
0.349223221159309
mpmath.functions.besselk(n, x)

besselk(n,x) computes the modified Bessel function of the second kind,

K_n(x) = \frac{\pi}{2} \frac{I_{-n}(x)-I_{n}(x)}{\sin(\pi n)}

For n an integer, this formula should be understood as a limit.

Examples

Some values and limits of K_n(x):

>>> from mpmath import *
>>> mp.dps = 25
>>> print besselk(0,0)
+inf
>>> print besselk(1,0)
+inf
>>> print besselk(0,1)
0.4210244382407083333356274
>>> print besselk(3.5, 2+3j)
(-0.02090732889633760668464128 + 0.2464022641351420167819697j)

Hankel functions (hankel1(), hankel2())

mpmath.functions.hankel1(n, x)

hankel1(n,x) computes the Hankel function of the first kind, which is the complex combination of Bessel functions given by

H_n^{(1)}(x) = J_n(x) + i Y_n(x).

Examples

The Hankel function is generally complex-valued:

>>> from mpmath import *
>>> mp.dps = 25
>>> print hankel1(2, pi)
(0.4854339326315091097054957 - 0.0999007139290278787734903j)
>>> print hankel1(3.5, pi)
(0.2340002029630507922628888 - 0.6419643823412927142424049j)
mpmath.functions.hankel2(n, x)

hankel2(n,x) computes the Hankel function of the second kind, which is the complex combination of Bessel functions given by

H_n^{(2)}(x) = J_n(x) - i Y_n(x).

Examples

The Hankel function is generally complex-valued:

>>> from mpmath import *
>>> mp.dps = 25
>>> print hankel2(2, pi)
(0.4854339326315091097054957 + 0.0999007139290278787734903j)
>>> print hankel2(3.5, pi)
(0.2340002029630507922628888 + 0.6419643823412927142424049j)

Airy functions (airyai(), airybi())

mpmath.functions.airyai(x)

Computes the Airy function \mathrm{Ai}(x), which is a solution of the Airy differential equation y''-xy=0. The Ai-function behaves roughly like a slowly decaying sine wave for x < 0 and like a decreasing exponential for x > 0.

Limits and values include:

>>> from mpmath import *
>>> mp.dps = 15
>>> print airyai(0), 1/(3**(2/3.)*gamma(2/3.))
0.355028053887817 0.355028053887817
>>> print airyai(1)
0.135292416312881
>>> print airyai(-1)
0.535560883292352
>>> print airyai(inf)
0.0
>>> print airyai(-inf)
0.0

airyai() uses a series expansion around x = 0, so it is slow for extremely large arguments. Here are some evaluations for moderately large arguments:

>>> print airyai(-100)
0.176753393239553
>>> print airyai(100)
2.63448215208818e-291
>>> print airyai(50+50j)
(-5.31790195707456e-68 - 1.16358800377071e-67j)
>>> print airyai(-50+50j)
(1.04124253736317e+158 + 3.3475255449236e+157j)

The first negative root is:

>>> print findroot(airyai, -2)
-2.33810741045977

We can verify the differential equation:

>>> for x in [-3.4, 0, 2.5, 1+2j]:
...     print abs(diff(airyai, x, 2) - x*airyai(x)) < eps
...
True
True
True
True

The Taylor series expansion around x = 0 starts with the following coefficients (note that every third term is zero):

>>> nprint(chop(taylor(airyai, 0, 5)))
[0.355028, -0.258819, 0.0, 5.91713e-2, -2.15683e-2, 0.0]

The Airy functions are a special case of Bessel functions. For x < 0, we have:

>>> x = 3
>>> print airyai(-x)
-0.378814293677658
>>> p = 2*(x**1.5)/3
>>> print sqrt(x)*(besselj(1/3.,p) + besselj(-1/3.,p))/3
-0.378814293677658
mpmath.functions.airybi(x)

Computes the Airy function \mathrm{Bi}(x), which is a solution of the Airy differential equation y''-xy=0. The Bi-function behaves roughly like a slowly decaying sine wave for x < 0 and like an increasing exponential for x > 0.

Limits and values include:

>>> from mpmath import *
>>> mp.dps = 15
>>> print airybi(0), 1/(3**(1/6.)*gamma(2/3.))
0.614926627446001 0.614926627446001
>>> print airybi(1)
1.20742359495287
>>> print airybi(-1)
0.103997389496945
>>> print airybi(inf)
+inf
>>> print airybi(-inf)
0.0

airyai() uses a series expansion around x = 0, so it is slow for extremely large arguments. Here are some evaluations for moderately large arguments:

>>> print airybi(-100)
0.0242738876801601
>>> print airybi(100)
6.0412239966702e+288
>>> print airybi(50+50j)
(-5.32207626732144e+63 + 1.47845029116524e+65j)
>>> print airybi(-50+50j)
(-3.3475255449236e+157 + 1.04124253736317e+158j)

The first negative root is:

>>> print findroot(airybi, -1)
-1.17371322270913

We can verify the differential equation:

>>> for x in [-3.4, 0, 2.5, 1+2j]:
...     print abs(diff(airybi, x, 2) - x*airybi(x)) < eps
...
True
True
True
True

The Taylor series expansion around x = 0 starts with the following coefficients (note that every third term is zero):

>>> nprint(chop(taylor(airybi, 0, 5)))
[0.614927, 0.448288, 0.0, 0.102488, 3.73574e-2, 0.0]

The Airy functions are a special case of Bessel functions. For x < 0, we have:

>>> x = 3
>>> print airybi(-x)
-0.198289626374927
>>> p = 2*(x**1.5)/3
>>> print sqrt(x/3)*(besselj(-1/3.,p) - besselj(1/3.,p))
-0.198289626374926

Orthogonal polynomials (jacobi(), legendre(), chebyt(), chebyu())

mpmath.functions.jacobi(n, a, b, x)

jacobi(n, a, b, x) evaluates the Jacobi polynomial P_n^{(a,b)}(x). The Jacobi polynomials are a special case of the hypergeometric function \,_2F_1 given by:

P_n^{(a,b)}(x) = {n+a \choose n}
  \,_2F_1\left(-n,1+a+b+n,a+1,\frac{1-x}{2}\right).

Note that this definition generalizes to nonintegral values of n. When n is an integer, the hypergeometric series terminates after a finite number of terms, giving a polynomial in x.

Evaluation of Jacobi polynomials

A special evaluation is P_n^{(a,b)}(1) = {n+a \choose n}:

>>> from mpmath import *
>>> mp.dps = 15
>>> print jacobi(4, 0.5, 0.25, 1)
2.4609375
>>> print binomial(4+0.5, 4)
2.4609375

A Jacobi polynomial of degree n is equal to its Taylor polynomial of degree n. The explicit coefficients of Jacobi polynomials can therefore be recovered easily using taylor():

>>> for n in range(5):
...     nprint(taylor(lambda x: jacobi(n,1,2,x), 0, n))
...
[1.0]
[-0.5, 2.5]
[-0.75, -1.5, 5.25]
[0.5, -3.5, -3.5, 10.5]
[0.625, 2.5, -11.25, -7.5, 20.625]

For nonintegral n, the Jacobi “polynomial” is no longer a polynomial:

>>> nprint(taylor(lambda x: jacobi(0.5,1,2,x), 0, 4))
[0.309983, 1.84119, -1.26933, 1.26699, -1.34808]

Orthogonality

The Jacobi polynomials are orthogonal on the interval [-1, 1] with respect to the weight function w(x) = (1-x)^a (1+x)^b. That is, w(x) P_n^{(a,b)}(x) P_m^{(a,b)}(x) integrates to zero if m \ne n and to a nonzero number if m = n.

The orthogonality is easy to verify using numerical quadrature:

>>> P = jacobi
>>> f = lambda x: (1-x)**a * (1+x)**b * P(m,a,b,x) * P(n,a,b,x)
>>> a = 2
>>> b = 3
>>> m, n = 3, 4
>>> print chop(quad(f, [-1, 1]), 1)
0.0
>>> m, n = 4, 4
>>> print quad(f, [-1, 1])
1.9047619047619

Differential equation

The Jacobi polynomials are solutions of the differential equation

(1-x^2) y'' + (b-a-(a+b+2)x) y' + n (n+a+b+1) y = 0.

We can verify that jacobi() approximately satisfies this equation:

>>> from mpmath import *
>>> mp.dps = 15
>>> a = 2.5
>>> b = 4
>>> n = 3
>>> y = lambda x: jacobi(n,a,b,x)
>>> x = pi
>>> A0 = n*(n+a+b+1)*y(x)
>>> A1 = (b-a-(a+b+2)*x)*diff(y,x)
>>> A2 = (1-x**2)*diff(y,x,2)
>>> nprint(A2 + A1 + A0, 1)
4.0e-12

The difference of order 10^{-12} is as close to zero as it could be at 15-digit working precision, since the terms are large:

>>> print A0, A1, A2
26560.2328981879 -21503.7641037294 -5056.46879445852
mpmath.functions.legendre(n, x)

legendre(n, x) evaluates the Legendre polynomial P_n(x). The Legendre polynomials are given by the formula

P_n(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} (x^2 -1)^n.

Alternatively, they can be computed recursively using

P_0(x) = 1

P_1(x) = x

(n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x).

A third definition is in terms of the hypergeometric function \,_2F_1, whereby they can be generalized to arbitrary n:

P_n(x) = \,_2F_1\left(-n, n+1, 1, \frac{1-x}{2}\right)

Basic evaluation

The Legendre polynomials assume fixed values at the points x = -1 and x = 1:

>>> from mpmath import *
>>> mp.dps = 15
>>> nprint([legendre(n, 1) for n in range(6)])
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
>>> nprint([legendre(n, -1) for n in range(6)])
[1.0, -1.0, 1.0, -1.0, 1.0, -1.0]

The coefficients of Legendre polynomials can be recovered using degree-n Taylor expansion:

>>> for n in range(5):
...     nprint(chop(taylor(lambda x: legendre(n, x), 0, n)))
...
[1.0]
[0.0, 1.0]
[-0.5, 0.0, 1.5]
[0.0, -1.5, 0.0, 2.5]
[0.375, 0.0, -3.75, 0.0, 4.375]

The roots of Legendre polynomials are located symmetrically on the interval [-1, 1]:

>>> for n in range(5):
...     nprint(polyroots(taylor(lambda x: legendre(n, x), 0, n)[::-1]))
...
[]
[0.0]
[-0.57735, 0.57735]
[-0.774597, 0.0, 0.774597]
[-0.861136, -0.339981, 0.339981, 0.861136]

An example of an evaluation for arbitrary n:

>>> print legendre(0.75, 2+4j)
(1.94952805264875 + 2.1071073099422j)

Orthogonality

The Legendre polynomials are orthogonal on [-1, 1] with respect to the trivial weight w(x) = 1. That is, P_m(x) P_n(x) integrates to zero if m \ne n and to 2/(2n+1) if m = n:

>>> m, n = 3, 4
>>> print quad(lambda x: legendre(m,x)*legendre(n,x), [-1, 1])
0.0
>>> m, n = 4, 4
>>> print quad(lambda x: legendre(m,x)*legendre(n,x), [-1, 1])
0.222222222222222

Differential equation

The Legendre polynomials satisfy the differential equation

((1-x^2) y')' + n(n+1) y' = 0.

We can verify this numerically:

>>> n = 3.6
>>> x = 0.73
>>> P = legendre
>>> A = diff(lambda t: (1-t**2)*diff(lambda u: P(n,u), t), x)
>>> B = n*(n+1)*P(n,x)
>>> nprint(A+B,1)
9.0e-16
mpmath.functions.chebyt(n, x)

chebyt(n, x) evaluates the Chebyshev polynomial of the first kind T_n(x), defined by the identity

T_n(\cos x) = \cos(n x).

The Chebyshev polynomials of the first kind are a special case of the Jacobi polynomials, and by extension of the hypergeometric function \,_2F_1. They can thus also be evaluated for nonintegral n.

Basic evaluation

The coefficients of the n-th polynomial can be recovered using using degree-n Taylor expansion:

>>> from mpmath import *
>>> mp.dps = 15
>>> for n in range(5):
...     nprint(chop(taylor(lambda x: chebyt(n, x), 0, n)))
...
[1.0]
[0.0, 1.0]
[-1.0, 0.0, 2.0]
[0.0, -3.0, 0.0, 4.0]
[1.0, 0.0, -8.0, 0.0, 8.0]

Orthogonality

The Chebyshev polynomials of the first kind are orthogonal on the interval [-1, 1] with respect to the weight function w(x) = 1/\sqrt{1-x^2}:

>>> f = lambda x: chebyt(m,x)*chebyt(n,x)/sqrt(1-x**2)
>>> m, n = 3, 4
>>> nprint(quad(f, [-1, 1]),1)
0.0
>>> m, n = 4, 4
>>> print quad(f, [-1, 1])
1.57079632596448
mpmath.functions.chebyu(n, x)

chebyu(n, x) evaluates the Chebyshev polynomial of the second kind U_n(x), defined by the identity

U_n(\cos x) = \frac{\sin((n+1)x)}{\sin(x)}.

The Chebyshev polynomials of the second kind are a special case of the Jacobi polynomials, and by extension of the hypergeometric function \,_2F_1. They can thus also be evaluated for nonintegral n.

Basic evaluation

The coefficients of the n-th polynomial can be recovered using using degree-n Taylor expansion:

>>> from mpmath import *
>>> mp.dps = 15
>>> for n in range(5):
...     nprint(chop(taylor(lambda x: chebyu(n, x), 0, n)))
...
[1.0]
[0.0, 2.0]
[-1.0, 0.0, 4.0]
[0.0, -4.0, 0.0, 8.0]
[1.0, 0.0, -12.0, 0.0, 16.0]

Orthogonality

The Chebyshev polynomials of the second kind are orthogonal on the interval [-1, 1] with respect to the weight function w(x) = \sqrt{1-x^2}:

>>> f = lambda x: chebyu(m,x)*chebyu(n,x)*sqrt(1-x**2)
>>> m, n = 3, 4
>>> print quad(f, [-1, 1])
0.0
>>> m, n = 4, 4
>>> print quad(f, [-1, 1])
1.5707963267949

Generic hypergeometric series (hyper(), hyp0f1(), hyp1f1(), hyp2f1())

mpmath.functions.hyper(a_s, b_s, z)

Evaluates the generalized hypergeometric function

\,_pF_q(a_1,\ldots,a_p; b_1,\ldots,b_q; z) =
\sum_{n=0}^\infty \frac{(a_1)_n (a_2)_n \ldots (a_p)_n}
   {(b_1)_n(b_2)_n\ldots(b_q)_n} \frac{z^n}{n!}

where (x)_n denotes the rising factorial (see rf()).

The parameters lists a_s and b_s may contain integers, real numbers, complex numbers, as well as exact fractions given in the form of tuples (p, q). hyper() is optimized to handle integers and fractions more efficiently than arbitrary floating-point parameters (since rational parameters are by far the most common).

Examples

We can compare the output of hyper() with nsum():

>>> from mpmath import *
>>> mp.dps = 25
>>> a,b,c,d = 2,3,4,5
>>> x = 0.25
>>> print hyper([a,b],[c,d],x)
1.078903941164934876086237
>>> fn = lambda n: rf(a,n)*rf(b,n)/rf(c,n)/rf(d,n)*x**n/fac(n)
>>> print nsum(fn, [0, inf])
1.078903941164934876086237

The parameters can be any combination of integers, fractions, floats and complex numbers:

>>> a, b, c, d, e = 1, (-1,2), pi, 3+4j, (2,3)
>>> x = 0.2j
>>> print hyper([a,b],[c,d,e],x)
(0.9923571616434024810831887 - 0.005753848733883879742993122j)
>>> b, e = -0.5, mpf(2)/3
>>> fn = lambda n: rf(a,n)*rf(b,n)/rf(c,n)/rf(d,n)/rf(e,n)*x**n/fac(n)
>>> print nsum(fn, [0, inf])
(0.9923571616434024810831887 - 0.005753848733883879742993122j)
mpmath.functions.hyp0f1(a, z)
Hypergeometric function \,_0F_1. hyp0f1(a,z) is equivalent to hyper([],[a],z); see documentation for hyper() for more information.
mpmath.functions.hyp1f1(a, b, z)
Hypergeometric function \,_1F_1. hyp1f1(a,b,z) is equivalent to hyper([a],[b],z); see documentation for hyper() for more information.
mpmath.functions.hyp2f1(a, b, c, z)
Hypergeometric function \,_2F_1. hyp2f1(a,b,c,z) is equivalent to hyper([a,b],[c],z); see documentation for hyper() for more information.