pyerrors.special

beta(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

beta(a, b, out=None)

Beta function.

This function is defined in ¹ as

$$B(a, b) = \int_0^1 t^{a-1}(1-t)^{b-1}dt = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)},$$

where \( \Gamma \) is the gamma function.

Parameters

• a, b (array_like): Real-valued arguments
• out (ndarray, optional): Optional output array for the function result

Returns

• scalar or ndarray: Value of the beta function

See Also

gamma: the gamma function
betainc: the regularized incomplete beta function
betaln: the natural logarithm of the absolute value of the beta function

References

Examples

>>> import scipy.special as sc

The beta function relates to the gamma function by the definition given above:

>>> sc.beta(2, 3)
0.08333333333333333
>>> sc.gamma(2)*sc.gamma(3)/sc.gamma(2 + 3)
0.08333333333333333

As this relationship demonstrates, the beta function is symmetric:

>>> sc.beta(1.7, 2.4)
0.16567527689031739
>>> sc.beta(2.4, 1.7)
0.16567527689031739

This function satisfies \( B(1, b) = 1/b \):

>>> sc.beta(1, 4)
0.25

betainc(x1, x2, x3, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

betainc(a, b, x, out=None)

Regularized incomplete beta function.

Computes the regularized incomplete beta function, defined as ¹:

$$I_x(a, b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \int_0^x t^{a-1}(1-t)^{b-1}dt,$$

for \( 0 \leq x \leq 1 \).

This function is the cumulative distribution function for the beta distribution; its range is [0, 1].

Parameters

• a, b (array_like): Positive, real-valued parameters
• x (array_like): Real-valued such that \( 0 \leq x \leq 1 \), the upper limit of integration
• out (ndarray, optional): Optional output array for the function values

Returns

• scalar or ndarray: Value of the regularized incomplete beta function

See Also

beta: beta function
betaincinv: inverse of the regularized incomplete beta function
betaincc: complement of the regularized incomplete beta function
scipy.stats.beta: beta distribution

Notes

The term regularized in the name of this function refers to the scaling of the function by the gamma function terms shown in the formula. When not qualified as regularized, the name incomplete beta function often refers to just the integral expression, without the gamma terms. One can use the function beta from scipy.special to get this "nonregularized" incomplete beta function by multiplying the result of betainc(a, b, x) by beta(a, b).

References

Examples

Let \( B(a, b) \) be the beta function.

>>> import scipy.special as sc

The coefficient in terms of gamma is equal to \( 1/B(a, b) \). Also, when \( x=1 \) the integral is equal to \( B(a, b) \). Therefore, \( I_{x=1}(a, b) = 1 \) for any \( a, b \).

>>> sc.betainc(0.2, 3.5, 1.0)
1.0

It satisfies \( I_x(a, b) = x^a F(a, 1-b, a+1, x)/ (aB(a, b)) \), where \( F \) is the hypergeometric function hyp2f1:

>>> a, b, x = 1.4, 3.1, 0.5
>>> x**a * sc.hyp2f1(a, 1 - b, a + 1, x)/(a * sc.beta(a, b))
0.8148904036225295
>>> sc.betainc(a, b, x)
0.8148904036225296

This functions satisfies the relationship \( I_x(a, b) = 1 - I_{1-x}(b, a) \):

>>> sc.betainc(2.2, 3.1, 0.4)
0.49339638807619446
>>> 1 - sc.betainc(3.1, 2.2, 1 - 0.4)
0.49339638807619446

betaln(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

betaln(a, b, out=None)

Natural logarithm of absolute value of beta function.

Computes ln(abs(beta(a, b))).

Parameters

• a, b (array_like): Positive, real-valued parameters
• out (ndarray, optional): Optional output array for function values

Returns

• scalar or ndarray: Value of the betaln function

See Also

gamma: the gamma function
betainc: the regularized incomplete beta function
beta: the beta function

Examples

>>> import numpy as np
>>> from scipy.special import betaln, beta

Verify that, for moderate values of a and b, betaln(a, b) is the same as log(beta(a, b)):

>>> betaln(3, 4)
-4.0943445622221

>>> np.log(beta(3, 4))
-4.0943445622221

In the following beta(a, b) underflows to 0, so we can't compute the logarithm of the actual value.

>>> a = 400
>>> b = 900
>>> beta(a, b)
0.0

We can compute the logarithm of beta(a, b) by using betaln:

>>> betaln(a, b)
-804.3069951764146

Polygamma functions.

Defined as \( \psi^{(n)}(x) \) where \( \psi \) is the digamma function. See [dlmf]_ for details.

Parameters

• n (array_like): The order of the derivative of the digamma function; must be integral
• x (array_like): Real valued input

Returns

• ndarray: Function results

See Also

digamma

References

.. [dlmf] NIST, Digital Library of Mathematical Functions, https://dlmf.nist.gov/5.15

Examples

>>> from scipy import special
>>> x = [2, 3, 25.5]
>>> special.polygamma(1, x)
array([ 0.64493407,  0.39493407,  0.03999467])
>>> special.polygamma(0, x) == special.psi(x)
array([ True,  True,  True], dtype=bool)

psi(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

psi(z, out=None)

The digamma function.

The logarithmic derivative of the gamma function evaluated at z.

Parameters

• z (array_like): Real or complex argument.
• out (ndarray, optional): Array for the computed values of psi.

Returns

• digamma (scalar or ndarray): Computed values of psi.

Notes

For large values not close to the negative real axis, psi is computed using the asymptotic series (5.11.2) from ¹. For small arguments not close to the negative real axis, the recurrence relation (5.5.2) from ² is used until the argument is large enough to use the asymptotic series. For values close to the negative real axis, the reflection formula (5.5.4) from ³ is used first. Note that psi has a family of zeros on the negative real axis which occur between the poles at nonpositive integers. Around the zeros the reflection formula suffers from cancellation and the implementation loses precision. The sole positive zero and the first negative zero, however, are handled separately by precomputing series expansions using ⁴, so the function should maintain full accuracy around the origin.

References

Examples

>>> from scipy.special import psi
>>> z = 3 + 4j
>>> psi(z)
(1.55035981733341+1.0105022091860445j)

Verify psi(z) = psi(z + 1) - 1/z:

>>> psi(z + 1) - 1/z
(1.55035981733341+1.0105022091860445j)

psi(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

psi(z, out=None)

The digamma function.

The logarithmic derivative of the gamma function evaluated at z.

Parameters

• z (array_like): Real or complex argument.
• out (ndarray, optional): Array for the computed values of psi.

Returns

• digamma (scalar or ndarray): Computed values of psi.

Notes

For large values not close to the negative real axis, psi is computed using the asymptotic series (5.11.2) from ¹. For small arguments not close to the negative real axis, the recurrence relation (5.5.2) from ² is used until the argument is large enough to use the asymptotic series. For values close to the negative real axis, the reflection formula (5.5.4) from ³ is used first. Note that psi has a family of zeros on the negative real axis which occur between the poles at nonpositive integers. Around the zeros the reflection formula suffers from cancellation and the implementation loses precision. The sole positive zero and the first negative zero, however, are handled separately by precomputing series expansions using ⁴, so the function should maintain full accuracy around the origin.

References

Examples

>>> from scipy.special import psi
>>> z = 3 + 4j
>>> psi(z)
(1.55035981733341+1.0105022091860445j)

Verify psi(z) = psi(z + 1) - 1/z:

>>> psi(z + 1) - 1/z
(1.55035981733341+1.0105022091860445j)

gamma(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

gamma(z, out=None)

gamma function.

The gamma function is defined as

$$\Gamma(z) = \int_0^\infty t^{z-1} e^{-t} dt$$

for \( \Re(z) > 0 \) and is extended to the rest of the complex plane by analytic continuation. See [dlmf]_ for more details.

Parameters

• z (array_like): Real or complex valued argument
• out (ndarray, optional): Optional output array for the function values

Returns

• scalar or ndarray: Values of the gamma function

Notes

The gamma function is often referred to as the generalized factorial since \( \Gamma(n + 1) = n! \) for natural numbers \( n \). More generally it satisfies the recurrence relation \( \Gamma(z + 1) = z \cdot \Gamma(z) \) for complex \( z \), which, combined with the fact that \( \Gamma(1) = 1 \), implies the above identity for \( z = n \).

References

.. [dlmf] NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5.2#E1

Examples

>>> import numpy as np
>>> from scipy.special import gamma, factorial

>>> gamma([0, 0.5, 1, 5])
array([         inf,   1.77245385,   1.        ,  24.        ])

>>> z = 2.5 + 1j
>>> gamma(z)
(0.77476210455108352+0.70763120437959293j)
>>> gamma(z+1), z*gamma(z)  # Recurrence property
((1.2292740569981171+2.5438401155000685j),
 (1.2292740569981158+2.5438401155000658j))

>>> gamma(0.5)**2  # gamma(0.5) = sqrt(pi)
3.1415926535897927

Plot gamma(x) for real x

>>> x = np.linspace(-3.5, 5.5, 2251)
>>> y = gamma(x)

>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'b', alpha=0.6, label='gamma(x)')
>>> k = np.arange(1, 7)
>>> plt.plot(k, factorial(k-1), 'k*', alpha=0.6,
...          label='(x-1)!, x = 1, 2, ...')
>>> plt.xlim(-3.5, 5.5)
>>> plt.ylim(-10, 25)
>>> plt.grid()
>>> plt.xlabel('x')
>>> plt.legend(loc='lower right')
>>> plt.show()

gammaln(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

gammaln(x, out=None)

Logarithm of the absolute value of the gamma function.

Defined as

$$\ln(\lvert\Gamma(x)\rvert)$$

where \( \Gamma \) is the gamma function. For more details on the gamma function, see [dlmf]_.

Parameters

• x (array_like): Real argument
• out (ndarray, optional): Optional output array for the function results

Returns

• scalar or ndarray: Values of the log of the absolute value of gamma

See Also

gammasgn: sign of the gamma function
loggamma: principal branch of the logarithm of the gamma function

Notes

It is the same function as the Python standard library function math.lgamma().

When used in conjunction with gammasgn, this function is useful for working in logspace on the real axis without having to deal with complex numbers via the relation exp(gammaln(x)) = gammasgn(x) * gamma(x).

For complex-valued log-gamma, use loggamma instead of gammaln.

References

.. [dlmf] NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5

Examples

>>> import numpy as np
>>> import scipy.special as sc

It has two positive zeros.

>>> sc.gammaln([1, 2])
array([0., 0.])

It has poles at nonpositive integers.

>>> sc.gammaln([0, -1, -2, -3, -4])
array([inf, inf, inf, inf, inf])

It asymptotically approaches x * log(x) (Stirling's formula).

>>> x = np.array([1e10, 1e20, 1e40, 1e80])
>>> sc.gammaln(x)
array([2.20258509e+11, 4.50517019e+21, 9.11034037e+41, 1.83206807e+82])
>>> x * np.log(x)
array([2.30258509e+11, 4.60517019e+21, 9.21034037e+41, 1.84206807e+82])

gammainc(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

gammainc(a, x, out=None)

Regularized lower incomplete gamma function.

It is defined as

$$P(a, x) = \frac{1}{\Gamma(a)} \int_0^x t^{a - 1}e^{-t} dt$$

for \( a > 0 \) and \( x \geq 0 \). See [dlmf]_ for details.

Parameters

• a (array_like): Positive parameter
• x (array_like): Nonnegative argument
• out (ndarray, optional): Optional output array for the function values

Returns

• scalar or ndarray: Values of the lower incomplete gamma function

See Also

gammaincc: regularized upper incomplete gamma function
gammaincinv: inverse of the regularized lower incomplete gamma function
gammainccinv: inverse of the regularized upper incomplete gamma function

Notes

The function satisfies the relation gammainc(a, x) + gammaincc(a, x) = 1 where gammaincc is the regularized upper incomplete gamma function.

The implementation largely follows that of [boost]_.

References

.. [dlmf] NIST Digital Library of Mathematical functions https://dlmf.nist.gov/8.2#E4 .. [boost] Maddock et. al., "Incomplete Gamma Functions", https://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html

Examples

>>> import scipy.special as sc

It is the CDF of the gamma distribution, so it starts at 0 and monotonically increases to 1.

>>> sc.gammainc(0.5, [0, 1, 10, 100])
array([0.        , 0.84270079, 0.99999226, 1.        ])

It is equal to one minus the upper incomplete gamma function.

>>> a, x = 0.5, 0.4
>>> sc.gammainc(a, x)
0.6289066304773024
>>> 1 - sc.gammaincc(a, x)
0.6289066304773024

gammaincc(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

gammaincc(a, x, out=None)

Regularized upper incomplete gamma function.

It is defined as

$$Q(a, x) = \frac{1}{\Gamma(a)} \int_x^\infty t^{a - 1}e^{-t} dt$$

for \( a > 0 \) and \( x \geq 0 \). See [dlmf]_ for details.

Parameters

• a (array_like): Positive parameter
• x (array_like): Nonnegative argument
• out (ndarray, optional): Optional output array for the function values

Returns

• scalar or ndarray: Values of the upper incomplete gamma function

See Also

gammainc: regularized lower incomplete gamma function
gammaincinv: inverse of the regularized lower incomplete gamma function
gammainccinv: inverse of the regularized upper incomplete gamma function

Notes

The function satisfies the relation gammainc(a, x) + gammaincc(a, x) = 1 where gammainc is the regularized lower incomplete gamma function.

The implementation largely follows that of [boost]_.

References

.. [dlmf] NIST Digital Library of Mathematical functions https://dlmf.nist.gov/8.2#E4 .. [boost] Maddock et. al., "Incomplete Gamma Functions", https://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html

Examples

>>> import scipy.special as sc

It is the survival function of the gamma distribution, so it starts at 1 and monotonically decreases to 0.

>>> sc.gammaincc(0.5, [0, 1, 10, 100, 1000])
array([1.00000000e+00, 1.57299207e-01, 7.74421643e-06, 2.08848758e-45,
       0.00000000e+00])

It is equal to one minus the lower incomplete gamma function.

>>> a, x = 0.5, 0.4
>>> sc.gammaincc(a, x)
0.37109336952269756
>>> 1 - sc.gammainc(a, x)
0.37109336952269756

gammasgn(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

gammasgn(x, out=None)

Sign of the gamma function.

It is defined as

$$\text{gammasgn}(x) = \begin{cases} +1 & \Gamma(x) > 0 \ -1 & \Gamma(x) < 0 \end{cases}$$

where \( \Gamma \) is the gamma function; see gamma. This definition is complete since the gamma function is never zero; see the discussion after [dlmf]_.

Parameters

• x (array_like): Real argument
• out (ndarray, optional): Optional output array for the function values

Returns

• scalar or ndarray: Sign of the gamma function

See Also

gamma: the gamma function
gammaln: log of the absolute value of the gamma function
loggamma: analytic continuation of the log of the gamma function

Notes

The gamma function can be computed as gammasgn(x) * np.exp(gammaln(x)).

References

.. [dlmf] NIST Digital Library of Mathematical Functions https://dlmf.nist.gov/5.2#E1

Examples

>>> import numpy as np
>>> import scipy.special as sc

It is 1 for x > 0.

>>> sc.gammasgn([1, 2, 3, 4])
array([1., 1., 1., 1.])

It alternates between -1 and 1 for negative integers.

>>> sc.gammasgn([-0.5, -1.5, -2.5, -3.5])
array([-1.,  1., -1.,  1.])

It can be used to compute the gamma function.

>>> x = [1.5, 0.5, -0.5, -1.5]
>>> sc.gammasgn(x) * np.exp(sc.gammaln(x))
array([ 0.88622693,  1.77245385, -3.5449077 ,  2.3632718 ])
>>> sc.gamma(x)
array([ 0.88622693,  1.77245385, -3.5449077 ,  2.3632718 ])

rgamma(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

rgamma(z, out=None)

Reciprocal of the gamma function.

Defined as \( 1 / \Gamma(z) \), where \( \Gamma \) is the gamma function. For more on the gamma function see gamma.

Parameters

• z (array_like): Real or complex valued input
• out (ndarray, optional): Optional output array for the function results

Returns

• scalar or ndarray: Function results

See Also

gamma,, gammaln,, loggamma

Notes

The gamma function has no zeros and has simple poles at nonpositive integers, so rgamma is an entire function with zeros at the nonpositive integers. See the discussion in [dlmf]_ for more details.

References

.. [dlmf] Nist, Digital Library of Mathematical functions, https://dlmf.nist.gov/5.2#i

Examples

>>> import scipy.special as sc

It is the reciprocal of the gamma function.

>>> sc.rgamma([1, 2, 3, 4])
array([1.        , 1.        , 0.5       , 0.16666667])
>>> 1 / sc.gamma([1, 2, 3, 4])
array([1.        , 1.        , 0.5       , 0.16666667])

It is zero at nonpositive integers.

>>> sc.rgamma([0, -1, -2, -3])
array([0., 0., 0., 0.])

It rapidly underflows to zero along the positive real axis.

>>> sc.rgamma([10, 100, 179])
array([2.75573192e-006, 1.07151029e-156, 0.00000000e+000])

Returns the log of multivariate gamma, also sometimes called the generalized gamma.

Parameters

• a (ndarray): The multivariate gamma is computed for each item of a.
• d (int): The dimension of the space of integration.

Returns

• res (ndarray): The values of the log multivariate gamma at the given points a.

Notes

The formal definition of the multivariate gamma of dimension d for a real a is

$$\Gamma_d(a) = \int_{A>0} e^{-tr(A)} |A|^{a - (d+1)/2} dA$$

with the condition \( a > (d-1)/2 \), and \( A > 0 \) being the set of all the positive definite matrices of dimension d. Note that a is a scalar: the integrand only is multivariate, the argument is not (the function is defined over a subset of the real set).

This can be proven to be equal to the much friendlier equation

$$\Gamma_d(a) = \pi^{d(d-1)/4} \prod_{i=1}^{d} \Gamma(a - (i-1)/2).$$

References

R. J. Muirhead, Aspects of multivariate statistical theory (Wiley Series in probability and mathematical statistics).

Examples

>>> import numpy as np
>>> from scipy.special import multigammaln, gammaln
>>> a = 23.5
>>> d = 10
>>> multigammaln(a, d)
454.1488605074416

Verify that the result agrees with the logarithm of the equation shown above:

>>> d*(d-1)/4*np.log(np.pi) + gammaln(a - 0.5*np.arange(0, d)).sum()
454.1488605074416

Modified Bessel function of the second kind of integer order n

j0(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

j0(x, out=None)

Bessel function of the first kind of order 0.

Parameters

• x (array_like): Argument (float).
• out (ndarray, optional): Optional output array for the function values

Returns

• J (scalar or ndarray): Value of the Bessel function of the first kind of order 0 at x.

See Also

jv: Bessel function of real order and complex argument.
spherical_jn: spherical Bessel functions.

Notes

The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval the following rational approximation is used:

$$J_0(x) \approx (w - r_1^2)(w - r_2^2) \frac{P_3(w)}{Q_8(w)},$$

where \( w = x^2 \) and \( r_1 \), \( r_2 \) are the zeros of \( J_0 \), and \( P_3 \) and \( Q_8 \) are polynomials of degrees 3 and 8, respectively.

In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7.

This function is a wrapper for the Cephes ¹ routine j0. It should not be confused with the spherical Bessel functions (see spherical_jn).

References

Examples

Calculate the function at one point:

>>> from scipy.special import j0
>>> j0(1.)
0.7651976865579665

Calculate the function at several points:

>>> import numpy as np
>>> j0(np.array([-2., 0., 4.]))
array([ 0.22389078,  1.        , -0.39714981])

Plot the function from -20 to 20.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-20., 20., 1000)
>>> y = j0(x)
>>> ax.plot(x, y)
>>> plt.show()

y0(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

y0(x, out=None)

Bessel function of the second kind of order 0.

Parameters

• x (array_like): Argument (float).
• out (ndarray, optional): Optional output array for the function results

Returns

• Y (scalar or ndarray): Value of the Bessel function of the second kind of order 0 at x.

See Also

j0: Bessel function of the first kind of order 0
yv: Bessel function of the first kind

Notes

The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval a rational approximation \( R(x) \) is employed to compute,

$$Y_0(x) = R(x) + \frac{2 \log(x) J_0(x)}{\pi},$$

where \( J_0 \) is the Bessel function of the first kind of order 0.

In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7.

This function is a wrapper for the Cephes ¹ routine y0.

References

Examples

Calculate the function at one point:

>>> from scipy.special import y0
>>> y0(1.)
0.08825696421567697

Calculate at several points:

>>> import numpy as np
>>> y0(np.array([0.5, 2., 3.]))
array([-0.44451873,  0.51037567,  0.37685001])

Plot the function from 0 to 10.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(0., 10., 1000)
>>> y = y0(x)
>>> ax.plot(x, y)
>>> plt.show()

j1(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

j1(x, out=None)

Bessel function of the first kind of order 1.

Parameters

• x (array_like): Argument (float).
• out (ndarray, optional): Optional output array for the function values

Returns

• J (scalar or ndarray): Value of the Bessel function of the first kind of order 1 at x.

See Also

jv: Bessel function of the first kind
spherical_jn: spherical Bessel functions.

Notes

The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 24 term Chebyshev expansion is used. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5.

This function is a wrapper for the Cephes ¹ routine j1. It should not be confused with the spherical Bessel functions (see spherical_jn).

References

Examples

Calculate the function at one point:

>>> from scipy.special import j1
>>> j1(1.)
0.44005058574493355

Calculate the function at several points:

>>> import numpy as np
>>> j1(np.array([-2., 0., 4.]))
array([-0.57672481,  0.        , -0.06604333])

Plot the function from -20 to 20.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-20., 20., 1000)
>>> y = j1(x)
>>> ax.plot(x, y)
>>> plt.show()

y1(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

y1(x, out=None)

Bessel function of the second kind of order 1.

Parameters

• x (array_like): Argument (float).
• out (ndarray, optional): Optional output array for the function results

Returns

• Y (scalar or ndarray): Value of the Bessel function of the second kind of order 1 at x.

See Also

j1: Bessel function of the first kind of order 1
yn: Bessel function of the second kind
yv: Bessel function of the second kind

Notes

The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 25 term Chebyshev expansion is used, and computing \( J_1 \) (the Bessel function of the first kind) is required. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5.

This function is a wrapper for the Cephes ¹ routine y1.

References

Examples

Calculate the function at one point:

>>> from scipy.special import y1
>>> y1(1.)
-0.7812128213002888

Calculate at several points:

>>> import numpy as np
>>> y1(np.array([0.5, 2., 3.]))
array([-1.47147239, -0.10703243,  0.32467442])

Plot the function from 0 to 10.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(0., 10., 1000)
>>> y = y1(x)
>>> ax.plot(x, y)
>>> plt.show()

jv(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

jv(v, z, out=None)

Bessel function of the first kind of real order and complex argument.

Parameters

• v (array_like): Order (float).
• z (array_like): Argument (float or complex).
• out (ndarray, optional): Optional output array for the function values

Returns

• J (scalar or ndarray): Value of the Bessel function, \( J_v(z) \).

See Also

jve: \( J_v \) with leading exponential behavior stripped off.
spherical_jn: spherical Bessel functions.
j0: faster version of this function for order 0.
j1: faster version of this function for order 1.

Notes

For positive v values, the computation is carried out using the AMOS ¹ zbesj routine, which exploits the connection to the modified Bessel function \( I_v \),

$$J_v(z) = \exp(v\pi\imath/2) I_v(-\imath z)\qquad (\Im z > 0)

J_v(z) = \exp(-v\pi\imath/2) I_v(\imath z)\qquad (\Im z < 0)$$

For negative v values the formula,

$$J_{-v}(z) = J_v(z) \cos(\pi v) - Y_v(z) \sin(\pi v)$$

is used, where \( Y_v(z) \) is the Bessel function of the second kind, computed using the AMOS routine zbesy. Note that the second term is exactly zero for integer v; to improve accuracy the second term is explicitly omitted for v values such that v = floor(v).

Not to be confused with the spherical Bessel functions (see spherical_jn).

References

Examples

Evaluate the function of order 0 at one point.

>>> from scipy.special import jv
>>> jv(0, 1.)
0.7651976865579666

Evaluate the function at one point for different orders.

>>> jv(0, 1.), jv(1, 1.), jv(1.5, 1.)
(0.7651976865579666, 0.44005058574493355, 0.24029783912342725)

The evaluation for different orders can be carried out in one call by providing a list or NumPy array as argument for the v parameter:

>>> jv([0, 1, 1.5], 1.)
array([0.76519769, 0.44005059, 0.24029784])

Evaluate the function at several points for order 0 by providing an array for z.

>>> import numpy as np
>>> points = np.array([-2., 0., 3.])
>>> jv(0, points)
array([ 0.22389078,  1.        , -0.26005195])

If z is an array, the order parameter v must be broadcastable to the correct shape if different orders shall be computed in one call. To calculate the orders 0 and 1 for an 1D array:

>>> orders = np.array([[0], [1]])
>>> orders.shape
(2, 1)

>>> jv(orders, points)
array([[ 0.22389078,  1.        , -0.26005195],
       [-0.57672481,  0.        ,  0.33905896]])

Plot the functions of order 0 to 3 from -10 to 10.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-10., 10., 1000)
>>> for i in range(4):
...     ax.plot(x, jv(i, x), label=f'$J_{i!r}$')
>>> ax.legend()
>>> plt.show()

yn(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

yn(n, x, out=None)

Bessel function of the second kind of integer order and real argument.

Parameters

• n (array_like): Order (integer).
• x (array_like): Argument (float).
• out (ndarray, optional): Optional output array for the function results

Returns

• Y (scalar or ndarray): Value of the Bessel function, \( Y_n(x) \).

See Also

yv: For real order and real or complex argument.
y0: faster implementation of this function for order 0
y1: faster implementation of this function for order 1

Notes

Wrapper for the Cephes ¹ routine yn.

The function is evaluated by forward recurrence on n, starting with values computed by the Cephes routines y0 and y1. If n = 0 or 1, the routine for y0 or y1 is called directly.

References

Examples

Evaluate the function of order 0 at one point.

>>> from scipy.special import yn
>>> yn(0, 1.)
0.08825696421567697

Evaluate the function at one point for different orders.

>>> yn(0, 1.), yn(1, 1.), yn(2, 1.)
(0.08825696421567697, -0.7812128213002888, -1.6506826068162546)

The evaluation for different orders can be carried out in one call by providing a list or NumPy array as argument for the v parameter:

>>> yn([0, 1, 2], 1.)
array([ 0.08825696, -0.78121282, -1.65068261])

Evaluate the function at several points for order 0 by providing an array for z.

>>> import numpy as np
>>> points = np.array([0.5, 3., 8.])
>>> yn(0, points)
array([-0.44451873,  0.37685001,  0.22352149])

If z is an array, the order parameter v must be broadcastable to the correct shape if different orders shall be computed in one call. To calculate the orders 0 and 1 for an 1D array:

>>> orders = np.array([[0], [1]])
>>> orders.shape
(2, 1)

>>> yn(orders, points)
array([[-0.44451873,  0.37685001,  0.22352149],
       [-1.47147239,  0.32467442, -0.15806046]])

Plot the functions of order 0 to 3 from 0 to 10.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(0., 10., 1000)
>>> for i in range(4):
...     ax.plot(x, yn(i, x), label=f'$Y_{i!r}$')
>>> ax.set_ylim(-3, 1)
>>> ax.legend()
>>> plt.show()

i0(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

i0(x, out=None)

Modified Bessel function of order 0.

Defined as,

$$I_0(x) = \sum_{k=0}^\infty \frac{(x^2/4)^k}{(k!)^2} = J_0(\imath x),$$

where \( J_0 \) is the Bessel function of the first kind of order 0.

Parameters

• x (array_like): Argument (float)
• out (ndarray, optional): Optional output array for the function values

Returns

• I (scalar or ndarray): Value of the modified Bessel function of order 0 at x.

See Also

iv: Modified Bessel function of any order
i0e: Exponentially scaled modified Bessel function of order 0

Notes

The range is partitioned into the two intervals [0, 8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval.

This function is a wrapper for the Cephes ¹ routine i0.

References

Examples

Calculate the function at one point:

>>> from scipy.special import i0
>>> i0(1.)
1.2660658777520082

Calculate at several points:

>>> import numpy as np
>>> i0(np.array([-2., 0., 3.5]))
array([2.2795853 , 1.        , 7.37820343])

Plot the function from -10 to 10.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-10., 10., 1000)
>>> y = i0(x)
>>> ax.plot(x, y)
>>> plt.show()

i1(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

i1(x, out=None)

Modified Bessel function of order 1.

Defined as,

$$I_1(x) = \frac{1}{2}x \sum_{k=0}^\infty \frac{(x^2/4)^k}{k! (k + 1)!} = -\imath J_1(\imath x),$$

where \( J_1 \) is the Bessel function of the first kind of order 1.

Parameters

• x (array_like): Argument (float)
• out (ndarray, optional): Optional output array for the function values

Returns

• I (scalar or ndarray): Value of the modified Bessel function of order 1 at x.

See Also

iv: Modified Bessel function of the first kind
i1e: Exponentially scaled modified Bessel function of order 1

Notes

The range is partitioned into the two intervals [0, 8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval.

This function is a wrapper for the Cephes ¹ routine i1.

References

Examples

Calculate the function at one point:

>>> from scipy.special import i1
>>> i1(1.)
0.5651591039924851

Calculate the function at several points:

>>> import numpy as np
>>> i1(np.array([-2., 0., 6.]))
array([-1.59063685,  0.        , 61.34193678])

Plot the function between -10 and 10.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-10., 10., 1000)
>>> y = i1(x)
>>> ax.plot(x, y)
>>> plt.show()

iv(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature, extobj])

iv(v, z, out=None)

Modified Bessel function of the first kind of real order.

Parameters

• v (array_like): Order. If z is of real type and negative, v must be integer valued.
• z (array_like of float or complex): Argument.
• out (ndarray, optional): Optional output array for the function values

Returns

• scalar or ndarray: Values of the modified Bessel function.

See Also

ive: This function with leading exponential behavior stripped off.
i0: Faster version of this function for order 0.
i1: Faster version of this function for order 1.

Notes

For real z and \( v \in [-50, 50] \), the evaluation is carried out using Temme's method ¹. For larger orders, uniform asymptotic expansions are applied.

For complex z and positive v, the AMOS ² zbesi routine is called. It uses a power series for small z, the asymptotic expansion for large abs(z), the Miller algorithm normalized by the Wronskian and a Neumann series for intermediate magnitudes, and the uniform asymptotic expansions for \( I_v(z) \) and \( J_v(z) \) for large orders. Backward recurrence is used to generate sequences or reduce orders when necessary.

The calculations above are done in the right half plane and continued into the left half plane by the formula,

$$I_v(z \exp(\pm\imath\pi)) = \exp(\pm\pi v) I_v(z)$$

(valid when the real part of z is positive). For negative v, the formula

$$I_{-v}(z) = I_v(z) + \frac{2}{\pi} \sin(\pi v) K_v(z)$$

is used, where \( K_v(z) \) is the modified Bessel function of the second kind, evaluated using the AMOS routine zbesk.

References

Examples

Evaluate the function of order 0 at one point.

>>> from scipy.special import iv
>>> iv(0, 1.)
1.2660658777520084

Evaluate the function at one point for different orders.

>>> iv(0, 1.), iv(1, 1.), iv(1.5, 1.)
(1.2660658777520084, 0.565159103992485, 0.2935253263474798)

The evaluation for different orders can be carried out in one call by providing a list or NumPy array as argument for the v parameter:

>>> iv([0, 1, 1.5], 1.)
array([1.26606588, 0.5651591 , 0.29352533])

Evaluate the function at several points for order 0 by providing an array for z.

>>> import numpy as np
>>> points = np.array([-2., 0., 3.])
>>> iv(0, points)
array([2.2795853 , 1.        , 4.88079259])

If z is an array, the order parameter v must be broadcastable to the correct shape if different orders shall be computed in one call. To calculate the orders 0 and 1 for an 1D array:

>>> orders = np.array([[0], [1]])
>>> orders.shape
(2, 1)

>>> iv(orders, points)
array([[ 2.2795853 ,  1.        ,  4.88079259],
       [-1.59063685,  0.        ,  3.95337022]])

Plot the functions of order 0 to 3 from -5 to 5.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-5., 5., 1000)
>>> for i in range(4):
...     ax.plot(x, iv(i, x), label=f'$I_{i!r}$')
>>> ax.legend()
>>> plt.show()