Special

DiracDelta

class sympy.functions.special.delta_functions.DiracDelta

The DiracDelta function and its derivatives.

DiracDelta function has the following properties:

  1. diff(Heaviside(x),x) = DiracDelta(x)
  2. integrate(DiracDelta(x-a)*f(x),(x,-oo,oo)) = f(a) and integrate(DiracDelta(x-a)*f(x),(x,a-e,a+e)) = f(a)
  3. DiracDelta(x) = 0 for all x != 0
  4. DiracDelta(g(x)) = Sum_i(DiracDelta(x-x_i)/abs(g'(x_i))) Where x_i-s are the roots of g

Derivatives of k-th order of DiracDelta have the following property:

  1. DiracDelta(x,k) = 0, for all x != 0

See also

Heaviside, simplify, is_simple, sympy.functions.special.tensor_functions.KroneckerDelta

References

[R67]http://mathworld.wolfram.com/DeltaFunction.html
is_simple(self, x)

Tells whether the argument(args[0]) of DiracDelta is a linear expression in x.

x can be:

  • a symbol

See also

simplify, Directdelta

Examples

>>> from sympy import DiracDelta, cos
>>> from sympy.abc import x, y
>>> DiracDelta(x*y).is_simple(x)
True
>>> DiracDelta(x*y).is_simple(y)
True
>>> DiracDelta(x**2+x-2).is_simple(x)
False
>>> DiracDelta(cos(x)).is_simple(x)
False
simplify(self, x)

Compute a simplified representation of the function using property number 4.

x can be:

  • a symbol

See also

is_simple, Directdelta

Examples

>>> from sympy import DiracDelta
>>> from sympy.abc import x, y
>>> DiracDelta(x*y).simplify(x)
DiracDelta(x)/Abs(y)
>>> DiracDelta(x*y).simplify(y)
DiracDelta(y)/Abs(x)
>>> DiracDelta(x**2 + x - 2).simplify(x)
DiracDelta(x - 1)/3 + DiracDelta(x + 2)/3

Heaviside

class sympy.functions.special.delta_functions.Heaviside

Heaviside Piecewise function

Heaviside function has the following properties [*]:

  1. diff(Heaviside(x),x) = DiracDelta(x)

    ( 0, if x < 0

  2. Heaviside(x) = < ( 1/2 if x==0 [*]

    ( 1, if x > 0

[*]Regarding to the value at 0, Mathematica defines H(0) = 1, but Maple uses H(0) = undefined

I think is better to have H(0) = 1/2, due to the following:

integrate(DiracDelta(x), x) = Heaviside(x)
integrate(DiracDelta(x), (x, -oo, oo)) = 1

and since DiracDelta is a symmetric function, integrate(DiracDelta(x), (x, 0, oo)) should be 1/2 (which is what Maple returns).

If we take Heaviside(0) = 1/2, we would have integrate(DiracDelta(x), (x, 0, oo)) = `` ``Heaviside(oo) - Heaviside(0) = 1 - 1/2 = 1/2 and integrate(DiracDelta(x), (x, -oo, 0)) = `` ``Heaviside(0) - Heaviside(-oo) = 1/2 - 0 = 1/2

If we consider, instead Heaviside(0) = 1, we would have integrate(DiracDelta(x), (x, 0, oo)) = Heaviside(oo) - Heaviside(0) = 0 and integrate(DiracDelta(x), (x, -oo, 0)) = Heaviside(0) - Heaviside(-oo) = 1

See also

DiracDelta

References

[R68]http://mathworld.wolfram.com/HeavisideStepFunction.html

beta

sympy.functions.special.gamma_functions.beta(x, y)

Euler Beta function

beta(x, y) == gamma(x)*gamma(y) / gamma(x+y)

Special Cases of the Incomplete Gamma Functions

class sympy.functions.special.error_functions.Ei

The classical exponential integral.

For use in SymPy, this function is defined as

\[\operatorname{Ei}(x) = \sum_{n=1}^\infty \frac{x^n}{n\, n!} + \log(x) + \gamma,\]

where \(\gamma\) is the Euler-Mascheroni constant.

If \(x\) is a polar number, this defines an analytic function on the riemann surface of the logarithm. Otherwise this defines an analytic function in the cut plane \(\mathbb{C} \setminus (-\infty, 0]\).

Background

The name exponential integral comes from the following statement:

\[\operatorname{Ei}(x) = \int_{-\infty}^x \frac{e^t}{t} \mathrm{d}t\]

If the integral is interpreted as a Cauchy principal value, this statement holds for \(x > 0\) and \(\operatorname{Ei}(x)\) as defined above.

Note that we carefully avoided defining \(\operatorname{Ei}(x)\) for negative real \(x\). This is because above integral formula does not hold for any polar lift of such \(x\), indeed all branches of \(\operatorname{Ei}(x)\) above the negative reals are imaginary.

However, the following statement holds for all \(x \in \mathbb{R}^*\):

\[\int_{-\infty}^x \frac{e^t}{t} \mathrm{d}t = \frac{\operatorname{Ei}\left(|x|e^{i \arg(x)}\right) + \operatorname{Ei}\left(|x|e^{- i \arg(x)}\right)}{2},\]

where the integral is again understood to be a principal value if \(x > 0\), and \(|x|e^{i \arg(x)}\), \(|x|e^{- i \arg(x)}\) denote two conjugate polar lifts of \(x\).

See also

expint
Generalised exponential integral.
E1
Special case of the generalised exponential integral.
li
Logarithmic integral.
Li
Offset logarithmic integral.
Si
Sine integral.
Ci
Cosine integral.
Shi
Hyperbolic sine integral.
Chi
Hyperbolic cosine integral.

sympy.functions.special.gamma_functions.uppergamma

References

[R79]http://dlmf.nist.gov/6.6
[R80]http://en.wikipedia.org/wiki/Exponential_integral
[R81]Abramowitz & Stegun, section 5: http://www.math.sfu.ca/~cbm/aands/page_228.htm

Examples

>>> from sympy import Ei, polar_lift, exp_polar, I, pi
>>> from sympy.abc import x

The exponential integral in SymPy is strictly undefined for negative values of the argument. For convenience, exponential integrals with negative arguments are immediately converted into an expression that agrees with the classical integral definition:

>>> Ei(-1)
-I*pi + Ei(exp_polar(I*pi))

This yields a real value:

>>> Ei(-1).n(chop=True)
-0.219383934395520

On the other hand the analytic continuation is not real:

>>> Ei(polar_lift(-1)).n(chop=True)
-0.21938393439552 + 3.14159265358979*I

The exponential integral has a logarithmic branch point at the origin:

>>> Ei(x*exp_polar(2*I*pi))
Ei(x) + 2*I*pi

Differentiation is supported:

>>> Ei(x).diff(x)
exp(x)/x

The exponential integral is related to many other special functions. For example:

>>> from sympy import uppergamma, expint, Shi
>>> Ei(x).rewrite(expint)
-expint(1, x*exp_polar(I*pi)) - I*pi
>>> Ei(x).rewrite(Shi)
Chi(x) + Shi(x)
class sympy.functions.special.error_functions.expint

Generalized exponential integral.

This function is defined as

\[\operatorname{E}_\nu(z) = z^{\nu - 1} \Gamma(1 - \nu, z),\]

where \(\Gamma(1 - \nu, z)\) is the upper incomplete gamma function (uppergamma).

Hence for \(z\) with positive real part we have

\[\operatorname{E}_\nu(z) = \int_1^\infty \frac{e^{-zt}}{z^\nu} \mathrm{d}t,\]

which explains the name.

The representation as an incomplete gamma function provides an analytic continuation for \(\operatorname{E}_\nu(z)\). If \(\nu\) is a non-positive integer the exponential integral is thus an unbranched function of \(z\), otherwise there is a branch point at the origin. Refer to the incomplete gamma function documentation for details of the branching behavior.

See also

Ei
Another related function called exponential integral.
E1
The classical case, returns expint(1, z).
li
Logarithmic integral.
Li
Offset logarithmic integral.
Si
Sine integral.
Ci
Cosine integral.
Shi
Hyperbolic sine integral.
Chi
Hyperbolic cosine integral.

sympy.functions.special.gamma_functions.uppergamma

References

[R82]http://dlmf.nist.gov/8.19
[R83]http://functions.wolfram.com/GammaBetaErf/ExpIntegralE/
[R84]http://en.wikipedia.org/wiki/Exponential_integral

Examples

>>> from sympy import expint, S
>>> from sympy.abc import nu, z

Differentiation is supported. Differentiation with respect to z explains further the name: for integral orders, the exponential integral is an iterated integral of the exponential function.

>>> expint(nu, z).diff(z)
-expint(nu - 1, z)

Differentiation with respect to nu has no classical expression:

>>> expint(nu, z).diff(nu)
-z**(nu - 1)*meijerg(((), (1, 1)), ((0, 0, -nu + 1), ()), z)

At non-postive integer orders, the exponential integral reduces to the exponential function:

>>> expint(0, z)
exp(-z)/z
>>> expint(-1, z)
exp(-z)/z + exp(-z)/z**2

At half-integers it reduces to error functions:

>>> expint(S(1)/2, z)
-sqrt(pi)*erf(sqrt(z))/sqrt(z) + sqrt(pi)/sqrt(z)

At positive integer orders it can be rewritten in terms of exponentials and expint(1, z). Use expand_func() to do this:

>>> from sympy import expand_func
>>> expand_func(expint(5, z))
z**4*expint(1, z)/24 + (-z**3 + z**2 - 2*z + 6)*exp(-z)/24

The generalised exponential integral is essentially equivalent to the incomplete gamma function:

>>> from sympy import uppergamma
>>> expint(nu, z).rewrite(uppergamma)
z**(nu - 1)*uppergamma(-nu + 1, z)

As such it is branched at the origin:

>>> from sympy import exp_polar, pi, I
>>> expint(4, z*exp_polar(2*pi*I))
I*pi*z**3/3 + expint(4, z)
>>> expint(nu, z*exp_polar(2*pi*I))
z**(nu - 1)*(exp(2*I*pi*nu) - 1)*gamma(-nu + 1) + expint(nu, z)
sympy.functions.special.error_functions.E1(z)

Classical case of the generalized exponential integral.

This is equivalent to expint(1, z).

See also

Ei
Exponential integral.
expint
Generalised exponential integral.
li
Logarithmic integral.
Li
Offset logarithmic integral.
Si
Sine integral.
Ci
Cosine integral.
Shi
Hyperbolic sine integral.
Chi
Hyperbolic cosine integral.
class sympy.functions.special.error_functions.li

The classical logarithmic integral.

For the use in SymPy, this function is defined as

\[\operatorname{li}(x) = \int_0^x \frac{1}{\log(t)} \mathrm{d}t \,.\]

See also

Li
Offset logarithmic integral.
Ei
Exponential integral.
expint
Generalised exponential integral.
E1
Special case of the generalised exponential integral.
Si
Sine integral.
Ci
Cosine integral.
Shi
Hyperbolic sine integral.
Chi
Hyperbolic cosine integral.

References

[R85]http://en.wikipedia.org/wiki/Logarithmic_integral
[R86]http://mathworld.wolfram.com/LogarithmicIntegral.html
[R87]http://dlmf.nist.gov/6
[R88]http://mathworld.wolfram.com/SoldnersConstant.html

Examples

>>> from sympy import I, oo, li
>>> from sympy.abc import z

Several special values are known:

>>> li(0)
0
>>> li(1)
-oo
>>> li(oo)
oo

Differentiation with respect to z is supported:

>>> from sympy import diff
>>> diff(li(z), z)
1/log(z)

Defining the \(li\) function via an integral:

The logarithmic integral can also be defined in terms of Ei:

>>> from sympy import Ei
>>> li(z).rewrite(Ei)
Ei(log(z))
>>> diff(li(z).rewrite(Ei), z)
1/log(z)

We can numerically evaluate the logarithmic integral to arbitrary precision on the whole complex plane (except the singular points):

>>> li(2).evalf(30)
1.04516378011749278484458888919
>>> li(2*I).evalf(30)
1.0652795784357498247001125598 + 3.08346052231061726610939702133*I

We can even compute Soldner’s constant by the help of mpmath:

>>> from sympy.mpmath import findroot
>>> findroot(li, 2)
1.45136923488338

Further transformations include rewriting \(li\) in terms of the trigonometric integrals \(Si\), \(Ci\), \(Shi\) and \(Chi\):

>>> from sympy import Si, Ci, Shi, Chi
>>> li(z).rewrite(Si)
-log(I*log(z)) - log(1/log(z))/2 + log(log(z))/2 + Ci(I*log(z)) + Shi(log(z))
>>> li(z).rewrite(Ci)
-log(I*log(z)) - log(1/log(z))/2 + log(log(z))/2 + Ci(I*log(z)) + Shi(log(z))
>>> li(z).rewrite(Shi)
-log(1/log(z))/2 + log(log(z))/2 + Chi(log(z)) - Shi(log(z))
>>> li(z).rewrite(Chi)
-log(1/log(z))/2 + log(log(z))/2 + Chi(log(z)) - Shi(log(z))
class sympy.functions.special.error_functions.Li

The offset logarithmic integral.

For the use in SymPy, this function is defined as

\[\operatorname{Li}(x) = \operatorname{li}(x) - \operatorname{li}(2)\]

See also

li
Logarithmic integral.
Ei
Exponential integral.
expint
Generalised exponential integral.
E1
Special case of the generalised exponential integral.
Si
Sine integral.
Ci
Cosine integral.
Shi
Hyperbolic sine integral.
Chi
Hyperbolic cosine integral.

References

[R89]http://en.wikipedia.org/wiki/Logarithmic_integral
[R90]http://mathworld.wolfram.com/LogarithmicIntegral.html
[R91]http://dlmf.nist.gov/6

Examples

>>> from sympy import I, oo, Li
>>> from sympy.abc import z

The following special value is known:

>>> Li(2)
0

Differentiation with respect to z is supported:

>>> from sympy import diff
>>> diff(Li(z), z)
1/log(z)

The shifted logarithmic integral can be written in terms of \(li(z)\):

>>> from sympy import li
>>> Li(z).rewrite(li)
li(z) - li(2)

We can numerically evaluate the logarithmic integral to arbitrary precision on the whole complex plane (except the singular points):

>>> Li(2).evalf(30)
0
>>> Li(4).evalf(30)
1.92242131492155809316615998938
class sympy.functions.special.error_functions.Si

Sine integral.

This function is defined by

\[\operatorname{Si}(z) = \int_0^z \frac{\sin{t}}{t} \mathrm{d}t.\]

It is an entire function.

See also

Ci
Cosine integral.
Shi
Hyperbolic sine integral.
Chi
Hyperbolic cosine integral.
Ei
Exponential integral.
expint
Generalised exponential integral.
E1
Special case of the generalised exponential integral.
li
Logarithmic integral.
Li
Offset logarithmic integral.

References

[R92]http://en.wikipedia.org/wiki/Trigonometric_integral

Examples

>>> from sympy import Si
>>> from sympy.abc import z

The sine integral is an antiderivative of sin(z)/z:

>>> Si(z).diff(z)
sin(z)/z

It is unbranched:

>>> from sympy import exp_polar, I, pi
>>> Si(z*exp_polar(2*I*pi))
Si(z)

Sine integral behaves much like ordinary sine under multiplication by I:

>>> Si(I*z)
I*Shi(z)
>>> Si(-z)
-Si(z)

It can also be expressed in terms of exponential integrals, but beware that the latter is branched:

>>> from sympy import expint
>>> Si(z).rewrite(expint)
-I*(-expint(1, z*exp_polar(-I*pi/2))/2 +
     expint(1, z*exp_polar(I*pi/2))/2) + pi/2
class sympy.functions.special.error_functions.Ci

Cosine integral.

This function is defined for positive \(x\) by

\[\operatorname{Ci}(x) = \gamma + \log{x} + \int_0^x \frac{\cos{t} - 1}{t} \mathrm{d}t = -\int_x^\infty \frac{\cos{t}}{t} \mathrm{d}t,\]

where \(\gamma\) is the Euler-Mascheroni constant.

We have

\[\operatorname{Ci}(z) = -\frac{\operatorname{E}_1\left(e^{i\pi/2} z\right) + \operatorname{E}_1\left(e^{-i \pi/2} z\right)}{2}\]

which holds for all polar \(z\) and thus provides an analytic continuation to the Riemann surface of the logarithm.

The formula also holds as stated for \(z \in \mathbb{C}\) with \(\Re(z) > 0\). By lifting to the principal branch we obtain an analytic function on the cut complex plane.

See also

Si
Sine integral.
Shi
Hyperbolic sine integral.
Chi
Hyperbolic cosine integral.
Ei
Exponential integral.
expint
Generalised exponential integral.
E1
Special case of the generalised exponential integral.
li
Logarithmic integral.
Li
Offset logarithmic integral.

References

[R93]http://en.wikipedia.org/wiki/Trigonometric_integral

Examples

>>> from sympy import Ci
>>> from sympy.abc import z

The cosine integral is a primitive of \(\cos(z)/z\):

>>> Ci(z).diff(z)
cos(z)/z

It has a logarithmic branch point at the origin:

>>> from sympy import exp_polar, I, pi
>>> Ci(z*exp_polar(2*I*pi))
Ci(z) + 2*I*pi

The cosine integral behaves somewhat like ordinary \(\cos\) under multiplication by \(i\):

>>> from sympy import polar_lift
>>> Ci(polar_lift(I)*z)
Chi(z) + I*pi/2
>>> Ci(polar_lift(-1)*z)
Ci(z) + I*pi

It can also be expressed in terms of exponential integrals:

>>> from sympy import expint
>>> Ci(z).rewrite(expint)
-expint(1, z*exp_polar(-I*pi/2))/2 - expint(1, z*exp_polar(I*pi/2))/2
class sympy.functions.special.error_functions.Shi

Sinh integral.

This function is defined by

\[\operatorname{Shi}(z) = \int_0^z \frac{\sinh{t}}{t} \mathrm{d}t.\]

It is an entire function.

See also

Si
Sine integral.
Ci
Cosine integral.
Chi
Hyperbolic cosine integral.
Ei
Exponential integral.
expint
Generalised exponential integral.
E1
Special case of the generalised exponential integral.
li
Logarithmic integral.
Li
Offset logarithmic integral.

References

[R94]http://en.wikipedia.org/wiki/Trigonometric_integral

Examples

>>> from sympy import Shi
>>> from sympy.abc import z

The Sinh integral is a primitive of \(\sinh(z)/z\):

>>> Shi(z).diff(z)
sinh(z)/z

It is unbranched:

>>> from sympy import exp_polar, I, pi
>>> Shi(z*exp_polar(2*I*pi))
Shi(z)

The \(\sinh\) integral behaves much like ordinary \(\sinh\) under multiplication by \(i\):

>>> Shi(I*z)
I*Si(z)
>>> Shi(-z)
-Shi(z)

It can also be expressed in terms of exponential integrals, but beware that the latter is branched:

>>> from sympy import expint
>>> Shi(z).rewrite(expint)
expint(1, z)/2 - expint(1, z*exp_polar(I*pi))/2 - I*pi/2
class sympy.functions.special.error_functions.Chi

Cosh integral.

This function is defined for positive \(x\) by

\[\operatorname{Chi}(x) = \gamma + \log{x} + \int_0^x \frac{\cosh{t} - 1}{t} \mathrm{d}t,\]

where \(\gamma\) is the Euler-Mascheroni constant.

We have

\[\operatorname{Chi}(z) = \operatorname{Ci}\left(e^{i \pi/2}z\right) - i\frac{\pi}{2},\]

which holds for all polar \(z\) and thus provides an analytic continuation to the Riemann surface of the logarithm. By lifting to the principal branch we obtain an analytic function on the cut complex plane.

See also

Si
Sine integral.
Ci
Cosine integral.
Shi
Hyperbolic sine integral.
Ei
Exponential integral.
expint
Generalised exponential integral.
E1
Special case of the generalised exponential integral.
li
Logarithmic integral.
Li
Offset logarithmic integral.

References

[R95]http://en.wikipedia.org/wiki/Trigonometric_integral

Examples

>>> from sympy import Chi
>>> from sympy.abc import z

The \(\cosh\) integral is a primitive of \(\cosh(z)/z\):

>>> Chi(z).diff(z)
cosh(z)/z

It has a logarithmic branch point at the origin:

>>> from sympy import exp_polar, I, pi
>>> Chi(z*exp_polar(2*I*pi))
Chi(z) + 2*I*pi

The \(\cosh\) integral behaves somewhat like ordinary \(\cosh\) under multiplication by \(i\):

>>> from sympy import polar_lift
>>> Chi(polar_lift(I)*z)
Ci(z) + I*pi/2
>>> Chi(polar_lift(-1)*z)
Chi(z) + I*pi

It can also be expressed in terms of exponential integrals:

>>> from sympy import expint
>>> Chi(z).rewrite(expint)
-expint(1, z)/2 - expint(1, z*exp_polar(I*pi))/2 - I*pi/2
class sympy.functions.special.error_functions.FresnelIntegral

Base class for the Fresnel integrals.

class sympy.functions.special.error_functions.fresnels

Fresnel integral S.

This function is defined by

\[\operatorname{S}(z) = \int_0^z \sin{\frac{\pi}{2} t^2} \mathrm{d}t.\]

It is an entire function.

See also

fresnelc

References

[R96]http://en.wikipedia.org/wiki/Fresnel_integral
[R97]http://dlmf.nist.gov/7
[R98]http://mathworld.wolfram.com/FresnelIntegrals.html
[R99]http://functions.wolfram.com/GammaBetaErf/FresnelS

Examples

>>> from sympy import I, oo, fresnels
>>> from sympy.abc import z

Several special values are known:

>>> fresnels(0)
0
>>> fresnels(oo)
1/2
>>> fresnels(-oo)
-1/2
>>> fresnels(I*oo)
-I/2
>>> fresnels(-I*oo)
I/2

In general one can pull out factors of -1 and \(i\) from the argument:

>>> fresnels(-z)
-fresnels(z)
>>> fresnels(I*z)
-I*fresnels(z)

The Fresnel S integral obeys the mirror symmetry \(\overline{S(z)} = S(\bar{z})\):

>>> from sympy import conjugate
>>> conjugate(fresnels(z))
fresnels(conjugate(z))

Differentiation with respect to \(z\) is supported:

>>> from sympy import diff
>>> diff(fresnels(z), z)
sin(pi*z**2/2)

Defining the Fresnel functions via an integral

>>> from sympy import integrate, pi, sin, gamma, expand_func
>>> integrate(sin(pi*z**2/2), z)
3*fresnels(z)*gamma(3/4)/(4*gamma(7/4))
>>> expand_func(integrate(sin(pi*z**2/2), z))
fresnels(z)

We can numerically evaluate the Fresnel integral to arbitrary precision on the whole complex plane:

>>> fresnels(2).evalf(30)
0.343415678363698242195300815958
>>> fresnels(-2*I).evalf(30)
0.343415678363698242195300815958*I
class sympy.functions.special.error_functions.fresnelc

Fresnel integral C.

This function is defined by

\[\operatorname{C}(z) = \int_0^z \cos{\frac{\pi}{2} t^2} \mathrm{d}t.\]

It is an entire function.

See also

fresnels

References

[R100]http://en.wikipedia.org/wiki/Fresnel_integral
[R101]http://dlmf.nist.gov/7
[R102]http://mathworld.wolfram.com/FresnelIntegrals.html
[R103]http://functions.wolfram.com/GammaBetaErf/FresnelC

Examples

>>> from sympy import I, oo, fresnelc
>>> from sympy.abc import z

Several special values are known:

>>> fresnelc(0)
0
>>> fresnelc(oo)
1/2
>>> fresnelc(-oo)
-1/2
>>> fresnelc(I*oo)
I/2
>>> fresnelc(-I*oo)
-I/2

In general one can pull out factors of -1 and \(i\) from the argument:

>>> fresnelc(-z)
-fresnelc(z)
>>> fresnelc(I*z)
I*fresnelc(z)

The Fresnel C integral obeys the mirror symmetry \(\overline{C(z)} = C(\bar{z})\):

>>> from sympy import conjugate
>>> conjugate(fresnelc(z))
fresnelc(conjugate(z))

Differentiation with respect to \(z\) is supported:

>>> from sympy import diff
>>> diff(fresnelc(z), z)
cos(pi*z**2/2)

Defining the Fresnel functions via an integral

>>> from sympy import integrate, pi, cos, gamma, expand_func
>>> integrate(cos(pi*z**2/2), z)
fresnelc(z)*gamma(1/4)/(4*gamma(5/4))
>>> expand_func(integrate(cos(pi*z**2/2), z))
fresnelc(z)

We can numerically evaluate the Fresnel integral to arbitrary precision on the whole complex plane:

>>> fresnelc(2).evalf(30)
0.488253406075340754500223503357
>>> fresnelc(-2*I).evalf(30)
-0.488253406075340754500223503357*I

Error Functions

class sympy.functions.special.error_functions.erf

The Gauss error function. This function is defined as:

\[\mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{-t^2} \mathrm{d}t.\]

See also

erfc
Complementary error function.
erfi
Imaginary error function.
erf2
Two-argument error function.
erfinv
Inverse error function.
erfcinv
Inverse Complementary error function.
erf2inv
Inverse two-argument error function.

References

[R104]http://en.wikipedia.org/wiki/Error_function
[R105]http://dlmf.nist.gov/7
[R106]http://mathworld.wolfram.com/Erf.html
[R107]http://functions.wolfram.com/GammaBetaErf/Erf

Examples

>>> from sympy import I, oo, erf
>>> from sympy.abc import z

Several special values are known:

>>> erf(0)
0
>>> erf(oo)
1
>>> erf(-oo)
-1
>>> erf(I*oo)
oo*I
>>> erf(-I*oo)
-oo*I

In general one can pull out factors of -1 and I from the argument:

>>> erf(-z)
-erf(z)

The error function obeys the mirror symmetry:

>>> from sympy import conjugate
>>> conjugate(erf(z))
erf(conjugate(z))

Differentiation with respect to z is supported:

>>> from sympy import diff
>>> diff(erf(z), z)
2*exp(-z**2)/sqrt(pi)

We can numerically evaluate the error function to arbitrary precision on the whole complex plane:

>>> erf(4).evalf(30)
0.999999984582742099719981147840
>>> erf(-4*I).evalf(30)
-1296959.73071763923152794095062*I
class sympy.functions.special.error_functions.erfc

Complementary Error Function. The function is defined as:

\[\mathrm{erfc}(x) = \frac{2}{\sqrt{\pi}} \int_x^\infty e^{-t^2} \mathrm{d}t\]

See also

erf
Gaussian error function.
erfi
Imaginary error function.
erf2
Two-argument error function.
erfinv
Inverse error function.
erfcinv
Inverse Complementary error function.
erf2inv
Inverse two-argument error function.

References

[R108]http://en.wikipedia.org/wiki/Error_function
[R109]http://dlmf.nist.gov/7
[R110]http://mathworld.wolfram.com/Erfc.html
[R111]http://functions.wolfram.com/GammaBetaErf/Erfc

Examples

>>> from sympy import I, oo, erfc
>>> from sympy.abc import z

Several special values are known:

>>> erfc(0)
1
>>> erfc(oo)
0
>>> erfc(-oo)
2
>>> erfc(I*oo)
-oo*I
>>> erfc(-I*oo)
oo*I

The error function obeys the mirror symmetry:

>>> from sympy import conjugate
>>> conjugate(erfc(z))
erfc(conjugate(z))

Differentiation with respect to z is supported:

>>> from sympy import diff
>>> diff(erfc(z), z)
-2*exp(-z**2)/sqrt(pi)

It also follows

>>> erfc(-z)
-erfc(z) + 2

We can numerically evaluate the complementary error function to arbitrary precision on the whole complex plane:

>>> erfc(4).evalf(30)
0.0000000154172579002800188521596734869
>>> erfc(4*I).evalf(30)
1.0 - 1296959.73071763923152794095062*I
class sympy.functions.special.error_functions.erfi

Imaginary error function. The function erfi is defined as:

\[\mathrm{erfi}(x) = \frac{2}{\sqrt{\pi}} \int_0^x e^{t^2} \mathrm{d}t\]

See also

erf
Gaussian error function.
erfc
Complementary error function.
erf2
Two-argument error function.
erfinv
Inverse error function.
erfcinv
Inverse Complementary error function.
erf2inv
Inverse two-argument error function.

References

[R112]http://en.wikipedia.org/wiki/Error_function
[R113]http://mathworld.wolfram.com/Erfi.html
[R114]http://functions.wolfram.com/GammaBetaErf/Erfi

Examples

>>> from sympy import I, oo, erfi
>>> from sympy.abc import z

Several special values are known:

>>> erfi(0)
0
>>> erfi(oo)
oo
>>> erfi(-oo)
-oo
>>> erfi(I*oo)
I
>>> erfi(-I*oo)
-I

In general one can pull out factors of -1 and I from the argument:

>>> erfi(-z)
-erfi(z)
>>> from sympy import conjugate
>>> conjugate(erfi(z))
erfi(conjugate(z))

Differentiation with respect to z is supported:

>>> from sympy import diff
>>> diff(erfi(z), z)
2*exp(z**2)/sqrt(pi)

We can numerically evaluate the imaginary error function to arbitrary precision on the whole complex plane:

>>> erfi(2).evalf(30)
18.5648024145755525987042919132
>>> erfi(-2*I).evalf(30)
-0.995322265018952734162069256367*I
class sympy.functions.special.error_functions.erf2

Two-argument error function. This function is defined as:

\[\mathrm{erf2}(x, y) = \frac{2}{\sqrt{\pi}} \int_x^y e^{-t^2} \mathrm{d}t\]

See also

erf
Gaussian error function.
erfc
Complementary error function.
erfi
Imaginary error function.
erfinv
Inverse error function.
erfcinv
Inverse Complementary error function.
erf2inv
Inverse two-argument error function.

References

[R115]http://functions.wolfram.com/GammaBetaErf/Erf2/

Examples

>>> from sympy import I, oo, erf2
>>> from sympy.abc import x, y

Several special values are known:

>>> erf2(0, 0)
0
>>> erf2(x, x)
0
>>> erf2(x, oo)
-erf(x) + 1
>>> erf2(x, -oo)
-erf(x) - 1
>>> erf2(oo, y)
erf(y) - 1
>>> erf2(-oo, y)
erf(y) + 1

In general one can pull out factors of -1:

>>> erf2(-x, -y)
-erf2(x, y)

The error function obeys the mirror symmetry:

>>> from sympy import conjugate
>>> conjugate(erf2(x, y))
erf2(conjugate(x), conjugate(y))

Differentiation with respect to x, y is supported:

>>> from sympy import diff
>>> diff(erf2(x, y), x)
-2*exp(-x**2)/sqrt(pi)
>>> diff(erf2(x, y), y)
2*exp(-y**2)/sqrt(pi)
class sympy.functions.special.error_functions.erfinv

Inverse Error Function. The erfinv function is defined as:

\[\mathrm{erf}(x) = y \quad \Rightarrow \quad \mathrm{erfinv}(y) = x\]

See also

erf
Gaussian error function.
erfc
Complementary error function.
erfi
Imaginary error function.
erf2
Two-argument error function.
erfcinv
Inverse Complementary error function.
erf2inv
Inverse two-argument error function.

References

[R116]http://en.wikipedia.org/wiki/Error_function#Inverse_functions
[R117]http://functions.wolfram.com/GammaBetaErf/InverseErf/

Examples

>>> from sympy import I, oo, erfinv
>>> from sympy.abc import x

Several special values are known:

>>> erfinv(0)
0
>>> erfinv(1)
oo

Differentiation with respect to x is supported:

>>> from sympy import diff
>>> diff(erfinv(x), x)
sqrt(pi)*exp(erfinv(x)**2)/2

We can numerically evaluate the inverse error function to arbitrary precision on [-1, 1]:

>>> erfinv(0.2).evalf(30)
0.179143454621291692285822705344
class sympy.functions.special.error_functions.erfcinv

Inverse Complementary Error Function. The erfcinv function is defined as:

\[\mathrm{erfc}(x) = y \quad \Rightarrow \quad \mathrm{erfcinv}(y) = x\]

See also

erf
Gaussian error function.
erfc
Complementary error function.
erfi
Imaginary error function.
erf2
Two-argument error function.
erfinv
Inverse error function.
erf2inv
Inverse two-argument error function.

References

[R118]http://en.wikipedia.org/wiki/Error_function#Inverse_functions
[R119]http://functions.wolfram.com/GammaBetaErf/InverseErfc/

Examples

>>> from sympy import I, oo, erfcinv
>>> from sympy.abc import x

Several special values are known:

>>> erfcinv(1)
0
>>> erfcinv(0)
oo

Differentiation with respect to x is supported:

>>> from sympy import diff
>>> diff(erfcinv(x), x)
-sqrt(pi)*exp(erfcinv(x)**2)/2
class sympy.functions.special.error_functions.erf2inv

Two-argument Inverse error function. The erf2inv function is defined as:

\[\mathrm{erf2}(x, w) = y \quad \Rightarrow \quad \mathrm{erf2inv}(x, y) = w\]

See also

erf
Gaussian error function.
erfc
Complementary error function.
erfi
Imaginary error function.
erf2
Two-argument error function.
erfinv
Inverse error function.
erfcinv
Inverse complementary error function.

References

[R120]http://functions.wolfram.com/GammaBetaErf/InverseErf2/

Examples

>>> from sympy import I, oo, erf2inv, erfinv, erfcinv
>>> from sympy.abc import x, y

Several special values are known:

>>> erf2inv(0, 0)
0
>>> erf2inv(1, 0)
1
>>> erf2inv(0, 1)
oo
>>> erf2inv(0, y)
erfinv(y)
>>> erf2inv(oo, y)
erfcinv(-y)

Differentiation with respect to x and y is supported:

>>> from sympy import diff
>>> diff(erf2inv(x, y), x)
exp(-x**2 + erf2inv(x, y)**2)
>>> diff(erf2inv(x, y), y)
sqrt(pi)*exp(erf2inv(x, y)**2)/2

Bessel Type Functions

class sympy.functions.special.bessel.BesselBase

Abstract base class for bessel-type functions.

This class is meant to reduce code duplication. All Bessel type functions can 1) be differentiated, and the derivatives expressed in terms of similar functions and 2) be rewritten in terms of other bessel-type functions.

Here “bessel-type functions” are assumed to have one complex parameter.

To use this base class, define class attributes _a and _b such that 2*F_n' = -_a*F_{n+1} + b*F_{n-1}.

argument

The argument of the bessel-type function.

order

The order of the bessel-type function.

class sympy.functions.special.bessel.besselj

Bessel function of the first kind.

The Bessel \(J\) function of order \(\nu\) is defined to be the function satisfying Bessel’s differential equation

\[z^2 \frac{\mathrm{d}^2 w}{\mathrm{d}z^2} + z \frac{\mathrm{d}w}{\mathrm{d}z} + (z^2 - \nu^2) w = 0,\]

with Laurent expansion

\[J_\nu(z) = z^\nu \left(\frac{1}{\Gamma(\nu + 1) 2^\nu} + O(z^2) \right),\]

if \(\nu\) is not a negative integer. If \(\nu=-n \in \mathbb{Z}_{<0}\) is a negative integer, then the definition is

\[J_{-n}(z) = (-1)^n J_n(z).\]

See also

bessely, besseli, besselk

References

[R121]Abramowitz, Milton; Stegun, Irene A., eds. (1965), “Chapter 9”, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables
[R122]Luke, Y. L. (1969), The Special Functions and Their Approximations, Volume 1
[R123]http://en.wikipedia.org/wiki/Bessel_function
[R124]http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/

Examples

Create a Bessel function object:

>>> from sympy import besselj, jn
>>> from sympy.abc import z, n
>>> b = besselj(n, z)

Differentiate it:

>>> b.diff(z)
besselj(n - 1, z)/2 - besselj(n + 1, z)/2

Rewrite in terms of spherical Bessel functions:

>>> b.rewrite(jn)
sqrt(2)*sqrt(z)*jn(n - 1/2, z)/sqrt(pi)

Access the parameter and argument:

>>> b.order
n
>>> b.argument
z
class sympy.functions.special.bessel.bessely

Bessel function of the second kind.

The Bessel \(Y\) function of order \(\nu\) is defined as

\[Y_\nu(z) = \lim_{\mu \to \nu} \frac{J_\mu(z) \cos(\pi \mu) - J_{-\mu}(z)}{\sin(\pi \mu)},\]

where \(J_\mu(z)\) is the Bessel function of the first kind.

It is a solution to Bessel’s equation, and linearly independent from \(J_\nu\).

See also

besselj, besseli, besselk

References

[R125]http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/

Examples

>>> from sympy import bessely, yn
>>> from sympy.abc import z, n
>>> b = bessely(n, z)
>>> b.diff(z)
bessely(n - 1, z)/2 - bessely(n + 1, z)/2
>>> b.rewrite(yn)
sqrt(2)*sqrt(z)*yn(n - 1/2, z)/sqrt(pi)
class sympy.functions.special.bessel.besseli

Modified Bessel function of the first kind.

The Bessel I function is a solution to the modified Bessel equation

\[z^2 \frac{\mathrm{d}^2 w}{\mathrm{d}z^2} + z \frac{\mathrm{d}w}{\mathrm{d}z} + (z^2 + \nu^2)^2 w = 0.\]

It can be defined as

\[I_\nu(z) = i^{-\nu} J_\nu(iz),\]

where \(J_\nu(z)\) is the Bessel function of the first kind.

See also

besselj, bessely, besselk

References

[R126]http://functions.wolfram.com/Bessel-TypeFunctions/BesselI/

Examples

>>> from sympy import besseli
>>> from sympy.abc import z, n
>>> besseli(n, z).diff(z)
besseli(n - 1, z)/2 + besseli(n + 1, z)/2
class sympy.functions.special.bessel.besselk

Modified Bessel function of the second kind.

The Bessel K function of order \(\nu\) is defined as

\[K_\nu(z) = \lim_{\mu \to \nu} \frac{\pi}{2} \frac{I_{-\mu}(z) -I_\mu(z)}{\sin(\pi \mu)},\]

where \(I_\mu(z)\) is the modified Bessel function of the first kind.

It is a solution of the modified Bessel equation, and linearly independent from \(Y_\nu\).

See also

besselj, besseli, bessely

References

[R127]http://functions.wolfram.com/Bessel-TypeFunctions/BesselK/

Examples

>>> from sympy import besselk
>>> from sympy.abc import z, n
>>> besselk(n, z).diff(z)
-besselk(n - 1, z)/2 - besselk(n + 1, z)/2
class sympy.functions.special.bessel.hankel1

Hankel function of the first kind.

This function is defined as

\[H_\nu^{(1)} = J_\nu(z) + iY_\nu(z),\]

where \(J_\nu(z)\) is the Bessel function of the first kind, and \(Y_\nu(z)\) is the Bessel function of the second kind.

It is a solution to Bessel’s equation.

See also

hankel2, besselj, bessely

References

[R128]http://functions.wolfram.com/Bessel-TypeFunctions/HankelH1/

Examples

>>> from sympy import hankel1
>>> from sympy.abc import z, n
>>> hankel1(n, z).diff(z)
hankel1(n - 1, z)/2 - hankel1(n + 1, z)/2
class sympy.functions.special.bessel.hankel2

Hankel function of the second kind.

This function is defined as

\[H_\nu^{(2)} = J_\nu(z) - iY_\nu(z),\]

where \(J_\nu(z)\) is the Bessel function of the first kind, and \(Y_\nu(z)\) is the Bessel function of the second kind.

It is a solution to Bessel’s equation, and linearly independent from \(H_\nu^{(1)}\).

See also

hankel1, besselj, bessely

References

[R129]http://functions.wolfram.com/Bessel-TypeFunctions/HankelH2/

Examples

>>> from sympy import hankel2
>>> from sympy.abc import z, n
>>> hankel2(n, z).diff(z)
hankel2(n - 1, z)/2 - hankel2(n + 1, z)/2
class sympy.functions.special.bessel.jn

Spherical Bessel function of the first kind.

This function is a solution to the spherical Bessel equation

\[z^2 \frac{\mathrm{d}^2 w}{\mathrm{d}z^2} + 2z \frac{\mathrm{d}w}{\mathrm{d}z} + (z^2 - \nu(\nu + 1)) w = 0.\]

It can be defined as

\[j_\nu(z) = \sqrt{\frac{\pi}{2z}} J_{\nu + \frac{1}{2}}(z),\]

where \(J_\nu(z)\) is the Bessel function of the first kind.

See also

besselj, bessely, besselk, yn

Examples

>>> from sympy import Symbol, jn, sin, cos, expand_func
>>> z = Symbol("z")
>>> print(jn(0, z).expand(func=True))
sin(z)/z
>>> jn(1, z).expand(func=True) == sin(z)/z**2 - cos(z)/z
True
>>> expand_func(jn(3, z))
(-6/z**2 + 15/z**4)*sin(z) + (1/z - 15/z**3)*cos(z)

The spherical Bessel functions of integral order are calculated using the formula:

\[j_n(z) = f_n(z) \sin{z} + (-1)^{n+1} f_{-n-1}(z) \cos{z},\]

where the coefficients \(f_n(z)\) are available as polys.orthopolys.spherical_bessel_fn().

class sympy.functions.special.bessel.yn

Spherical Bessel function of the second kind.

This function is another solution to the spherical Bessel equation, and linearly independent from \(j_n\). It can be defined as

\[j_\nu(z) = \sqrt{\frac{\pi}{2z}} Y_{\nu + \frac{1}{2}}(z),\]

where \(Y_\nu(z)\) is the Bessel function of the second kind.

See also

besselj, bessely, besselk, jn

Examples

>>> from sympy import Symbol, yn, sin, cos, expand_func
>>> z = Symbol("z")
>>> print(expand_func(yn(0, z)))
-cos(z)/z
>>> expand_func(yn(1, z)) == -cos(z)/z**2-sin(z)/z
True

For integral orders \(n\), \(y_n\) is calculated using the formula:

\[y_n(z) = (-1)^{n+1} j_{-n-1}(z)\]
sympy.functions.special.bessel.jn_zeros(n, k, method='sympy', dps=15)

Zeros of the spherical Bessel function of the first kind.

This returns an array of zeros of jn up to the k-th zero.

  • method = “sympy”: uses mpmath.besseljzero()
  • method = “scipy”: uses the SciPy’s sph_jn and newton to find all roots, which is faster than computing the zeros using a general numerical solver, but it requires SciPy and only works with low precision floating point numbers. [The function used with method=”sympy” is a recent addition to mpmath, before that a general solver was used.]

See also

jn, yn, besselj, besselk, bessely

Examples

>>> from sympy import jn_zeros
>>> jn_zeros(2, 4, dps=5)
[5.7635, 9.095, 12.323, 15.515]

B-Splines

sympy.functions.special.bsplines.bspline_basis(d, knots, n, x, close=True)

The \(n\)-th B-spline at \(x\) of degree \(d\) with knots.

B-Splines are piecewise polynomials of degree \(d\) [R130]. They are defined on a set of knots, which is a sequence of integers or floats.

The 0th degree splines have a value of one on a single interval:

>>> from sympy import bspline_basis
>>> from sympy.abc import x
>>> d = 0
>>> knots = range(5)
>>> bspline_basis(d, knots, 0, x)
Piecewise((1, And(x <= 1, x >= 0)), (0, True))

For a given (d, knots) there are len(knots)-d-1 B-splines defined, that are indexed by n (starting at 0).

Here is an example of a cubic B-spline:

>>> bspline_basis(3, range(5), 0, x)
Piecewise((x**3/6, And(x < 1, x >= 0)),
          (-x**3/2 + 2*x**2 - 2*x + 2/3, And(x < 2, x >= 1)),
          (x**3/2 - 4*x**2 + 10*x - 22/3, And(x < 3, x >= 2)),
          (-x**3/6 + 2*x**2 - 8*x + 32/3, And(x <= 4, x >= 3)),
          (0, True))

By repeating knot points, you can introduce discontinuities in the B-splines and their derivatives:

>>> d = 1
>>> knots = [0,0,2,3,4]
>>> bspline_basis(d, knots, 0, x)
Piecewise((-x/2 + 1, And(x <= 2, x >= 0)), (0, True))

It is quite time consuming to construct and evaluate B-splines. If you need to evaluate a B-splines many times, it is best to lambdify them first:

>>> from sympy import lambdify
>>> d = 3
>>> knots = range(10)
>>> b0 = bspline_basis(d, knots, 0, x)
>>> f = lambdify(x, b0)
>>> y = f(0.5)

See also

bsplines_basis_set

References

[R130](1, 2) http://en.wikipedia.org/wiki/B-spline
sympy.functions.special.bsplines.bspline_basis_set(d, knots, x)

Return the len(knots)-d-1 B-splines at x of degree d with knots.

This function returns a list of Piecewise polynomials that are the len(knots)-d-1 B-splines of degree d for the given knots. This function calls bspline_basis(d, knots, n, x) for different values of n.

See also

bsplines_basis

Examples

>>> from sympy import bspline_basis_set
>>> from sympy.abc import x
>>> d = 2
>>> knots = range(5)
>>> splines = bspline_basis_set(d, knots, x)
>>> splines
[Piecewise((x**2/2, And(x < 1, x >= 0)),
           (-x**2 + 3*x - 3/2, And(x < 2, x >= 1)),
           (x**2/2 - 3*x + 9/2, And(x <= 3, x >= 2)),
           (0, True)),
 Piecewise((x**2/2 - x + 1/2, And(x < 2, x >= 1)),
           (-x**2 + 5*x - 11/2, And(x < 3, x >= 2)),
           (x**2/2 - 4*x + 8, And(x <= 4, x >= 3)),
           (0, True))]

Riemann Zeta and Related Functions

class sympy.functions.special.zeta_functions.zeta

Hurwitz zeta function (or Riemann zeta function).

For \(\operatorname{Re}(a) > 0\) and \(\operatorname{Re}(s) > 1\), this function is defined as

\[\zeta(s, a) = \sum_{n=0}^\infty \frac{1}{(n + a)^s},\]

where the standard choice of argument for \(n + a\) is used. For fixed \(a\) with \(\operatorname{Re}(a) > 0\) the Hurwitz zeta function admits a meromorphic continuation to all of \(\mathbb{C}\), it is an unbranched function with a simple pole at \(s = 1\).

Analytic continuation to other \(a\) is possible under some circumstances, but this is not typically done.

The Hurwitz zeta function is a special case of the Lerch transcendent:

\[\zeta(s, a) = \Phi(1, s, a).\]

This formula defines an analytic continuation for all possible values of \(s\) and \(a\) (also \(\operatorname{Re}(a) < 0\)), see the documentation of lerchphi for a description of the branching behavior.

If no value is passed for \(a\), by this function assumes a default value of \(a = 1\), yielding the Riemann zeta function.

References

[R131]http://dlmf.nist.gov/25.11
[R132]http://en.wikipedia.org/wiki/Hurwitz_zeta_function

Examples

For \(a = 1\) the Hurwitz zeta function reduces to the famous Riemann zeta function:

\[\zeta(s, 1) = \zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s}.\]
>>> from sympy import zeta
>>> from sympy.abc import s
>>> zeta(s, 1)
zeta(s)
>>> zeta(s)
zeta(s)

The Riemann zeta function can also be expressed using the Dirichlet eta function:

>>> from sympy import dirichlet_eta
>>> zeta(s).rewrite(dirichlet_eta)
dirichlet_eta(s)/(-2**(-s + 1) + 1)

The Riemann zeta function at positive even integer and negative odd integer values is related to the Bernoulli numbers:

>>> zeta(2)
pi**2/6
>>> zeta(4)
pi**4/90
>>> zeta(-1)
-1/12

The specific formulae are:

\[\zeta(2n) = (-1)^{n+1} \frac{B_{2n} (2\pi)^{2n}}{2(2n)!}\]
\[\zeta(-n) = -\frac{B_{n+1}}{n+1}\]

At negative even integers the Riemann zeta function is zero:

>>> zeta(-4)
0

No closed-form expressions are known at positive odd integers, but numerical evaluation is possible:

>>> zeta(3).n()
1.20205690315959

The derivative of \(\zeta(s, a)\) with respect to \(a\) is easily computed:

>>> from sympy.abc import a
>>> zeta(s, a).diff(a)
-s*zeta(s + 1, a)

However the derivative with respect to \(s\) has no useful closed form expression:

>>> zeta(s, a).diff(s)
Derivative(zeta(s, a), s)

The Hurwitz zeta function can be expressed in terms of the Lerch transcendent, sympy.functions.special.lerchphi:

>>> from sympy import lerchphi
>>> zeta(s, a).rewrite(lerchphi)
lerchphi(1, s, a)
class sympy.functions.special.zeta_functions.dirichlet_eta

Dirichlet eta function.

For \(\operatorname{Re}(s) > 0\), this function is defined as

\[\eta(s) = \sum_{n=1}^\infty \frac{(-1)^n}{n^s}.\]

It admits a unique analytic continuation to all of \(\mathbb{C}\). It is an entire, unbranched function.

See also

zeta

References

[R133]http://en.wikipedia.org/wiki/Dirichlet_eta_function

Examples

The Dirichlet eta function is closely related to the Riemann zeta function:

>>> from sympy import dirichlet_eta, zeta
>>> from sympy.abc import s
>>> dirichlet_eta(s).rewrite(zeta)
(-2**(-s + 1) + 1)*zeta(s)
class sympy.functions.special.zeta_functions.polylog

Polylogarithm function.

For \(|z| < 1\) and \(s \in \mathbb{C}\), the polylogarithm is defined by

\[\operatorname{Li}_s(z) = \sum_{n=1}^\infty \frac{z^n}{n^s},\]

where the standard branch of the argument is used for \(n\). It admits an analytic continuation which is branched at \(z=1\) (notably not on the sheet of initial definition), \(z=0\) and \(z=\infty\).

The name polylogarithm comes from the fact that for \(s=1\), the polylogarithm is related to the ordinary logarithm (see examples), and that

\[\operatorname{Li}_{s+1}(z) = \int_0^z \frac{\operatorname{Li}_s(t)}{t} \mathrm{d}t.\]

The polylogarithm is a special case of the Lerch transcendent:

\[\operatorname{Li}_{s}(z) = z \Phi(z, s, 1)\]

See also

zeta, lerchphi

Examples

For \(z \in \{0, 1, -1\}\), the polylogarithm is automatically expressed using other functions:

>>> from sympy import polylog
>>> from sympy.abc import s
>>> polylog(s, 0)
0
>>> polylog(s, 1)
zeta(s)
>>> polylog(s, -1)
dirichlet_eta(s)

If \(s\) is a negative integer, \(0\) or \(1\), the polylogarithm can be expressed using elementary functions. This can be done using expand_func():

>>> from sympy import expand_func
>>> from sympy.abc import z
>>> expand_func(polylog(1, z))
-log(z*exp_polar(-I*pi) + 1)
>>> expand_func(polylog(0, z))
z/(-z + 1)

The derivative with respect to \(z\) can be computed in closed form:

>>> polylog(s, z).diff(z)
polylog(s - 1, z)/z

The polylogarithm can be expressed in terms of the lerch transcendent:

>>> from sympy import lerchphi
>>> polylog(s, z).rewrite(lerchphi)
z*lerchphi(z, s, 1)
class sympy.functions.special.zeta_functions.lerchphi

Lerch transcendent (Lerch phi function).

For \(\operatorname{Re}(a) > 0\), \(|z| < 1\) and \(s \in \mathbb{C}\), the Lerch transcendent is defined as

\[\Phi(z, s, a) = \sum_{n=0}^\infty \frac{z^n}{(n + a)^s},\]

where the standard branch of the argument is used for \(n + a\), and by analytic continuation for other values of the parameters.

A commonly used related function is the Lerch zeta function, defined by

\[L(q, s, a) = \Phi(e^{2\pi i q}, s, a).\]

Analytic Continuation and Branching Behavior

It can be shown that

\[\Phi(z, s, a) = z\Phi(z, s, a+1) + a^{-s}.\]

This provides the analytic continuation to \(\operatorname{Re}(a) \le 0\).

Assume now \(\operatorname{Re}(a) > 0\). The integral representation

\[\Phi_0(z, s, a) = \int_0^\infty \frac{t^{s-1} e^{-at}}{1 - ze^{-t}} \frac{\mathrm{d}t}{\Gamma(s)}\]

provides an analytic continuation to \(\mathbb{C} - [1, \infty)\). Finally, for \(x \in (1, \infty)\) we find

\[\lim_{\epsilon \to 0^+} \Phi_0(x + i\epsilon, s, a) -\lim_{\epsilon \to 0^+} \Phi_0(x - i\epsilon, s, a) = \frac{2\pi i \log^{s-1}{x}}{x^a \Gamma(s)},\]

using the standard branch for both \(\log{x}\) and \(\log{\log{x}}\) (a branch of \(\log{\log{x}}\) is needed to evaluate \(\log{x}^{s-1}\)). This concludes the analytic continuation. The Lerch transcendent is thus branched at \(z \in \{0, 1, \infty\}\) and \(a \in \mathbb{Z}_{\le 0}\). For fixed \(z, a\) outside these branch points, it is an entire function of \(s\).

See also

polylog, zeta

References

[R134]Bateman, H.; Erdelyi, A. (1953), Higher Transcendental Functions, Vol. I, New York: McGraw-Hill. Section 1.11.
[R135]http://dlmf.nist.gov/25.14
[R136]http://en.wikipedia.org/wiki/Lerch_transcendent

Examples

The Lerch transcendent is a fairly general function, for this reason it does not automatically evaluate to simpler functions. Use expand_func() to achieve this.

If \(z=1\), the Lerch transcendent reduces to the Hurwitz zeta function:

>>> from sympy import lerchphi, expand_func
>>> from sympy.abc import z, s, a
>>> expand_func(lerchphi(1, s, a))
zeta(s, a)

More generally, if \(z\) is a root of unity, the Lerch transcendent reduces to a sum of Hurwitz zeta functions:

>>> expand_func(lerchphi(-1, s, a))
2**(-s)*zeta(s, a/2) - 2**(-s)*zeta(s, a/2 + 1/2)

If \(a=1\), the Lerch transcendent reduces to the polylogarithm:

>>> expand_func(lerchphi(z, s, 1))
polylog(s, z)/z

More generally, if \(a\) is rational, the Lerch transcendent reduces to a sum of polylogarithms:

>>> from sympy import S
>>> expand_func(lerchphi(z, s, S(1)/2))
2**(s - 1)*(polylog(s, sqrt(z))/sqrt(z) -
            polylog(s, sqrt(z)*exp_polar(I*pi))/sqrt(z))
>>> expand_func(lerchphi(z, s, S(3)/2))
-2**s/z + 2**(s - 1)*(polylog(s, sqrt(z))/sqrt(z) -
                      polylog(s, sqrt(z)*exp_polar(I*pi))/sqrt(z))/z

The derivatives with respect to \(z\) and \(a\) can be computed in closed form:

>>> lerchphi(z, s, a).diff(z)
(-a*lerchphi(z, s, a) + lerchphi(z, s - 1, a))/z
>>> lerchphi(z, s, a).diff(a)
-s*lerchphi(z, s + 1, a)

Hypergeometric Functions

class sympy.functions.special.hyper.hyper

The (generalized) hypergeometric function is defined by a series where the ratios of successive terms are a rational function of the summation index. When convergent, it is continued analytically to the largest possible domain.

The hypergeometric function depends on two vectors of parameters, called the numerator parameters \(a_p\), and the denominator parameters \(b_q\). It also has an argument \(z\). The series definition is

\[\begin{split}{}_pF_q\left(\begin{matrix} a_1, \dots, a_p \\ b_1, \dots, b_q \end{matrix} \middle| z \right) = \sum_{n=0}^\infty \frac{(a_1)_n \dots (a_p)_n}{(b_1)_n \dots (b_q)_n} \frac{z^n}{n!},\end{split}\]

where \((a)_n = (a)(a+1)\dots(a+n-1)\) denotes the rising factorial.

If one of the \(b_q\) is a non-positive integer then the series is undefined unless one of the \(a_p\) is a larger (i.e. smaller in magnitude) non-positive integer. If none of the \(b_q\) is a non-positive integer and one of the \(a_p\) is a non-positive integer, then the series reduces to a polynomial. To simplify the following discussion, we assume that none of the \(a_p\) or \(b_q\) is a non-positive integer. For more details, see the references.

The series converges for all \(z\) if \(p \le q\), and thus defines an entire single-valued function in this case. If \(p = q+1\) the series converges for \(|z| < 1\), and can be continued analytically into a half-plane. If \(p > q+1\) the series is divergent for all \(z\).

Note: The hypergeometric function constructor currently does not check if the parameters actually yield a well-defined function.

References

[R137]Luke, Y. L. (1969), The Special Functions and Their Approximations, Volume 1
[R138]http://en.wikipedia.org/wiki/Generalized_hypergeometric_function

Examples

The parameters \(a_p\) and \(b_q\) can be passed as arbitrary iterables, for example:

>>> from sympy.functions import hyper
>>> from sympy.abc import x, n, a
>>> hyper((1, 2, 3), [3, 4], x)
hyper((1, 2, 3), (3, 4), x)

There is also pretty printing (it looks better using unicode):

>>> from sympy import pprint
>>> pprint(hyper((1, 2, 3), [3, 4], x), use_unicode=False)
  _
 |_  /1, 2, 3 |  \
 |   |        | x|
3  2 \  3, 4  |  /

The parameters must always be iterables, even if they are vectors of length one or zero:

>>> hyper((1, ), [], x)
hyper((1,), (), x)

But of course they may be variables (but if they depend on x then you should not expect much implemented functionality):

>>> hyper((n, a), (n**2,), x)
hyper((n, a), (n**2,), x)

The hypergeometric function generalizes many named special functions. The function hyperexpand() tries to express a hypergeometric function using named special functions. For example:

>>> from sympy import hyperexpand
>>> hyperexpand(hyper([], [], x))
exp(x)

You can also use expand_func:

>>> from sympy import expand_func
>>> expand_func(x*hyper([1, 1], [2], -x))
log(x + 1)

More examples:

>>> from sympy import S
>>> hyperexpand(hyper([], [S(1)/2], -x**2/4))
cos(x)
>>> hyperexpand(x*hyper([S(1)/2, S(1)/2], [S(3)/2], x**2))
asin(x)

We can also sometimes hyperexpand parametric functions:

>>> from sympy.abc import a
>>> hyperexpand(hyper([-a], [], x))
(-x + 1)**a
ap

Numerator parameters of the hypergeometric function.

argument

Argument of the hypergeometric function.

bq

Denominator parameters of the hypergeometric function.

convergence_statement

Return a condition on z under which the series converges.

eta

A quantity related to the convergence of the series.

radius_of_convergence

Compute the radius of convergence of the defining series.

Note that even if this is not oo, the function may still be evaluated outside of the radius of convergence by analytic continuation. But if this is zero, then the function is not actually defined anywhere else.

>>> from sympy.functions import hyper
>>> from sympy.abc import z
>>> hyper((1, 2), [3], z).radius_of_convergence
1
>>> hyper((1, 2, 3), [4], z).radius_of_convergence
0
>>> hyper((1, 2), (3, 4), z).radius_of_convergence
oo
class sympy.functions.special.hyper.meijerg

The Meijer G-function is defined by a Mellin-Barnes type integral that resembles an inverse Mellin transform. It generalizes the hypergeometric functions.

The Meijer G-function depends on four sets of parameters. There are “numerator parameters\(a_1, \dots, a_n\) and \(a_{n+1}, \dots, a_p\), and there are “denominator parameters\(b_1, \dots, b_m\) and \(b_{m+1}, \dots, b_q\). Confusingly, it is traditionally denoted as follows (note the position of \(m\), \(n\), \(p\), \(q\), and how they relate to the lengths of the four parameter vectors):

\[\begin{split}G_{p,q}^{m,n} \left(\begin{matrix}a_1, \dots, a_n & a_{n+1}, \dots, a_p \\ b_1, \dots, b_m & b_{m+1}, \dots, b_q \end{matrix} \middle| z \right).\end{split}\]

However, in sympy the four parameter vectors are always available separately (see examples), so that there is no need to keep track of the decorating sub- and super-scripts on the G symbol.

The G function is defined as the following integral:

\[\frac{1}{2 \pi i} \int_L \frac{\prod_{j=1}^m \Gamma(b_j - s) \prod_{j=1}^n \Gamma(1 - a_j + s)}{\prod_{j=m+1}^q \Gamma(1- b_j +s) \prod_{j=n+1}^p \Gamma(a_j - s)} z^s \mathrm{d}s,\]

where \(\Gamma(z)\) is the gamma function. There are three possible contours which we will not describe in detail here (see the references). If the integral converges along more than one of them the definitions agree. The contours all separate the poles of \(\Gamma(1-a_j+s)\) from the poles of \(\Gamma(b_k-s)\), so in particular the G function is undefined if \(a_j - b_k \in \mathbb{Z}_{>0}\) for some \(j \le n\) and \(k \le m\).

The conditions under which one of the contours yields a convergent integral are complicated and we do not state them here, see the references.

Note: Currently the Meijer G-function constructor does not check any convergence conditions.

See also

hyper, sympy.simplify.hyperexpand

References

[R139]Luke, Y. L. (1969), The Special Functions and Their Approximations, Volume 1
[R140]http://en.wikipedia.org/wiki/Meijer_G-function

Examples

You can pass the parameters either as four separate vectors:

>>> from sympy.functions import meijerg
>>> from sympy.abc import x, a
>>> from sympy.core.containers import Tuple
>>> from sympy import pprint
>>> pprint(meijerg((1, 2), (a, 4), (5,), [], x), use_unicode=False)
 __1, 2 /1, 2  a, 4 |  \
/__     |           | x|
\_|4, 1 \ 5         |  /

or as two nested vectors:

>>> pprint(meijerg([(1, 2), (3, 4)], ([5], Tuple()), x), use_unicode=False)
 __1, 2 /1, 2  3, 4 |  \
/__     |           | x|
\_|4, 1 \ 5         |  /

As with the hypergeometric function, the parameters may be passed as arbitrary iterables. Vectors of length zero and one also have to be passed as iterables. The parameters need not be constants, but if they depend on the argument then not much implemented functionality should be expected.

All the subvectors of parameters are available:

>>> from sympy import pprint
>>> g = meijerg([1], [2], [3], [4], x)
>>> pprint(g, use_unicode=False)
 __1, 1 /1  2 |  \
/__     |     | x|
\_|2, 2 \3  4 |  /
>>> g.an
(1,)
>>> g.ap
(1, 2)
>>> g.aother
(2,)
>>> g.bm
(3,)
>>> g.bq
(3, 4)
>>> g.bother
(4,)

The Meijer G-function generalizes the hypergeometric functions. In some cases it can be expressed in terms of hypergeometric functions, using Slater’s theorem. For example:

>>> from sympy import hyperexpand
>>> from sympy.abc import a, b, c
>>> hyperexpand(meijerg([a], [], [c], [b], x), allow_hyper=True)
x**c*gamma(-a + c + 1)*hyper((-a + c + 1,),
                             (-b + c + 1,), -x)/gamma(-b + c + 1)

Thus the Meijer G-function also subsumes many named functions as special cases. You can use expand_func or hyperexpand to (try to) rewrite a Meijer G-function in terms of named special functions. For example:

>>> from sympy import expand_func, S
>>> expand_func(meijerg([[],[]], [[0],[]], -x))
exp(x)
>>> hyperexpand(meijerg([[],[]], [[S(1)/2],[0]], (x/2)**2))
sin(x)/sqrt(pi)
an

First set of numerator parameters.

aother

Second set of numerator parameters.

ap

Combined numerator parameters.

argument

Argument of the Meijer G-function.

bm

First set of denominator parameters.

bother

Second set of denominator parameters.

bq

Combined denominator parameters.

delta

A quantity related to the convergence region of the integral, c.f. references.

get_period()

Return a number P such that G(x*exp(I*P)) == G(x).

>>> from sympy.functions.special.hyper import meijerg
>>> from sympy.abc import z
>>> from sympy import pi, S
>>> meijerg([1], [], [], [], z).get_period()
2*pi
>>> meijerg([pi], [], [], [], z).get_period()
oo
>>> meijerg([1, 2], [], [], [], z).get_period()
oo
>>> meijerg([1,1], [2], [1, S(1)/2, S(1)/3], [1], z).get_period()
12*pi
integrand(s)

Get the defining integrand D(s).

nu

A quantity related to the convergence region of the integral, c.f. references.

Elliptic integrals

class sympy.functions.special.elliptic_integrals.elliptic_k

The complete elliptic integral of the first kind, defined by

\[K(z) = F\left(\tfrac{\pi}{2}\middle| z\right)\]

where \(F\left(z\middle| m\right)\) is the Legendre incomplete elliptic integral of the first kind.

The function \(K(z)\) is a single-valued function on the complex plane with branch cut along the interval \((1, \infty)\).

See also

elliptic_f

References

[R141]http://en.wikipedia.org/wiki/Elliptic_integrals
[R142]http://functions.wolfram.com/EllipticIntegrals/EllipticK

Examples

>>> from sympy import elliptic_k, I, pi
>>> from sympy.abc import z
>>> elliptic_k(0)
pi/2
>>> elliptic_k(1.0 + I)
1.50923695405127 + 0.625146415202697*I
>>> elliptic_k(z).series(z, n=3)
pi/2 + pi*z/8 + 9*pi*z**2/128 + O(z**3)
class sympy.functions.special.elliptic_integrals.elliptic_f

The Legendre incomplete elliptic integral of the first kind, defined by

\[F\left(z\middle| m\right) = \int_0^z \frac{dt}{\sqrt{1 - m \sin^2 t}}\]

This function reduces to a complete elliptic integral of the first kind, \(K(m)\), when \(z = \pi/2\).

See also

elliptic_k

References

[R143]http://en.wikipedia.org/wiki/Elliptic_integrals
[R144]http://functions.wolfram.com/EllipticIntegrals/EllipticF

Examples

>>> from sympy import elliptic_f, I, O
>>> from sympy.abc import z, m
>>> elliptic_f(z, m).series(z)
z + z**5*(3*m**2/40 - m/30) + m*z**3/6 + O(z**6)
>>> elliptic_f(3.0 + I/2, 1.0 + I)
2.909449841483 + 1.74720545502474*I
class sympy.functions.special.elliptic_integrals.elliptic_e

Called with two arguments \(z\) and \(m\), evaluates the incomplete elliptic integral of the second kind, defined by

\[E\left(z\middle| m\right) = \int_0^z \sqrt{1 - m \sin^2 t} dt\]

Called with a single argument \(z\), evaluates the Legendre complete elliptic integral of the second kind

\[E(z) = E\left(\tfrac{\pi}{2}\middle| z\right)\]

The function \(E(z)\) is a single-valued function on the complex plane with branch cut along the interval \((1, \infty)\).

References

[R145]http://en.wikipedia.org/wiki/Elliptic_integrals
[R146]http://functions.wolfram.com/EllipticIntegrals/EllipticE2
[R147]http://functions.wolfram.com/EllipticIntegrals/EllipticE

Examples

>>> from sympy import elliptic_e, I, pi, O
>>> from sympy.abc import z, m
>>> elliptic_e(z, m).series(z)
z + z**5*(-m**2/40 + m/30) - m*z**3/6 + O(z**6)
>>> elliptic_e(z).series(z, n=4)
pi/2 - pi*z/8 - 3*pi*z**2/128 - 5*pi*z**3/512 + O(z**4)
>>> elliptic_e(1 + I, 2 - I/2).n()
1.55203744279187 + 0.290764986058437*I
>>> elliptic_e(0)
pi/2
>>> elliptic_e(2.0 - I)
0.991052601328069 + 0.81879421395609*I
class sympy.functions.special.elliptic_integrals.elliptic_pi

Called with three arguments \(n\), \(z\) and \(m\), evaluates the Legendre incomplete elliptic integral of the third kind, defined by

\[\Pi\left(n; z\middle| m\right) = \int_0^z \frac{dt} {\left(1 - n \sin^2 t\right) \sqrt{1 - m \sin^2 t}}\]

Called with two arguments \(n\) and \(m\), evaluates the complete elliptic integral of the third kind:

\[\Pi\left(n\middle| m\right) = \Pi\left(n; \tfrac{\pi}{2}\middle| m\right)\]

References

[R148]http://en.wikipedia.org/wiki/Elliptic_integrals
[R149]http://functions.wolfram.com/EllipticIntegrals/EllipticPi3
[R150]http://functions.wolfram.com/EllipticIntegrals/EllipticPi

Examples

>>> from sympy import elliptic_pi, I, pi, O, S
>>> from sympy.abc import z, n, m
>>> elliptic_pi(n, z, m).series(z, n=4)
z + z**3*(m/6 + n/3) + O(z**4)
>>> elliptic_pi(0.5 + I, 1.0 - I, 1.2)
2.50232379629182 - 0.760939574180767*I
>>> elliptic_pi(0, 0)
pi/2
>>> elliptic_pi(1.0 - I/3, 2.0 + I)
3.29136443417283 + 0.32555634906645*I

Orthogonal Polynomials

This module mainly implements special orthogonal polynomials.

See also functions.combinatorial.numbers which contains some combinatorial polynomials.

Jacobi Polynomials

class sympy.functions.special.polynomials.jacobi

Jacobi polynomial \(P_n^{\left(\alpha, \beta\right)}(x)\)

jacobi(n, alpha, beta, x) gives the nth Jacobi polynomial in x, \(P_n^{\left(\alpha, \beta\right)}(x)\).

The Jacobi polynomials are orthogonal on \([-1, 1]\) with respect to the weight \(\left(1-x\right)^\alpha \left(1+x\right)^\beta\).

References

[R151]http://en.wikipedia.org/wiki/Jacobi_polynomials
[R152]http://mathworld.wolfram.com/JacobiPolynomial.html
[R153]http://functions.wolfram.com/Polynomials/JacobiP/

Examples

>>> from sympy import jacobi, S, conjugate, diff
>>> from sympy.abc import n,a,b,x
>>> jacobi(0, a, b, x)
1
>>> jacobi(1, a, b, x)
a/2 - b/2 + x*(a/2 + b/2 + 1)
>>> jacobi(2, a, b, x)   
(a**2/8 - a*b/4 - a/8 + b**2/8 - b/8 + x**2*(a**2/8 + a*b/4 + 7*a/8 +
b**2/8 + 7*b/8 + 3/2) + x*(a**2/4 + 3*a/4 - b**2/4 - 3*b/4) - 1/2)
>>> jacobi(n, a, b, x)
jacobi(n, a, b, x)
>>> jacobi(n, a, a, x)
RisingFactorial(a + 1, n)*gegenbauer(n,
    a + 1/2, x)/RisingFactorial(2*a + 1, n)
>>> jacobi(n, 0, 0, x)
legendre(n, x)
>>> jacobi(n, S(1)/2, S(1)/2, x)
RisingFactorial(3/2, n)*chebyshevu(n, x)/factorial(n + 1)
>>> jacobi(n, -S(1)/2, -S(1)/2, x)
RisingFactorial(1/2, n)*chebyshevt(n, x)/factorial(n)
>>> jacobi(n, a, b, -x)
(-1)**n*jacobi(n, b, a, x)
>>> jacobi(n, a, b, 0)
2**(-n)*gamma(a + n + 1)*hyper((-b - n, -n), (a + 1,), -1)/(factorial(n)*gamma(a + 1))
>>> jacobi(n, a, b, 1)
RisingFactorial(a + 1, n)/factorial(n)
>>> conjugate(jacobi(n, a, b, x))
jacobi(n, conjugate(a), conjugate(b), conjugate(x))
>>> diff(jacobi(n,a,b,x), x)
(a/2 + b/2 + n/2 + 1/2)*jacobi(n - 1, a + 1, b + 1, x)
sympy.functions.special.polynomials.jacobi_normalized(n, a, b, x)

Jacobi polynomial \(P_n^{\left(\alpha, \beta\right)}(x)\)

jacobi_normalized(n, alpha, beta, x) gives the nth Jacobi polynomial in x, \(P_n^{\left(\alpha, \beta\right)}(x)\).

The Jacobi polynomials are orthogonal on \([-1, 1]\) with respect to the weight \(\left(1-x\right)^\alpha \left(1+x\right)^\beta\).

This functions returns the polynomials normilzed:

\[\int_{-1}^{1} P_m^{\left(\alpha, \beta\right)}(x) P_n^{\left(\alpha, \beta\right)}(x) (1-x)^{\alpha} (1+x)^{\beta} \mathrm{d}x = \delta_{m,n}\]

References

[R154]http://en.wikipedia.org/wiki/Jacobi_polynomials
[R155]http://mathworld.wolfram.com/JacobiPolynomial.html
[R156]http://functions.wolfram.com/Polynomials/JacobiP/

Examples

>>> from sympy import jacobi_normalized
>>> from sympy.abc import n,a,b,x
>>> jacobi_normalized(n, a, b, x)
jacobi(n, a, b, x)/sqrt(2**(a + b + 1)*gamma(a + n + 1)*gamma(b + n + 1)/((a + b + 2*n + 1)*factorial(n)*gamma(a + b + n + 1)))

Gegenbauer Polynomials

class sympy.functions.special.polynomials.gegenbauer

Gegenbauer polynomial \(C_n^{\left(\alpha\right)}(x)\)

gegenbauer(n, alpha, x) gives the nth Gegenbauer polynomial in x, \(C_n^{\left(\alpha\right)}(x)\).

The Gegenbauer polynomials are orthogonal on \([-1, 1]\) with respect to the weight \(\left(1-x^2\right)^{\alpha-\frac{1}{2}}\).

References

[R157]http://en.wikipedia.org/wiki/Gegenbauer_polynomials
[R158]http://mathworld.wolfram.com/GegenbauerPolynomial.html
[R159]http://functions.wolfram.com/Polynomials/GegenbauerC3/

Examples

>>> from sympy import gegenbauer, conjugate, diff
>>> from sympy.abc import n,a,x
>>> gegenbauer(0, a, x)
1
>>> gegenbauer(1, a, x)
2*a*x
>>> gegenbauer(2, a, x)
-a + x**2*(2*a**2 + 2*a)
>>> gegenbauer(3, a, x)
x**3*(4*a**3/3 + 4*a**2 + 8*a/3) + x*(-2*a**2 - 2*a)
>>> gegenbauer(n, a, x)
gegenbauer(n, a, x)
>>> gegenbauer(n, a, -x)
(-1)**n*gegenbauer(n, a, x)
>>> gegenbauer(n, a, 0)
2**n*sqrt(pi)*gamma(a + n/2)/(gamma(a)*gamma(-n/2 + 1/2)*gamma(n + 1))
>>> gegenbauer(n, a, 1)
gamma(2*a + n)/(gamma(2*a)*gamma(n + 1))
>>> conjugate(gegenbauer(n, a, x))
gegenbauer(n, conjugate(a), conjugate(x))
>>> diff(gegenbauer(n, a, x), x)
2*a*gegenbauer(n - 1, a + 1, x)

Chebyshev Polynomials

class sympy.functions.special.polynomials.chebyshevt

Chebyshev polynomial of the first kind, \(T_n(x)\)

chebyshevt(n, x) gives the nth Chebyshev polynomial (of the first kind) in x, \(T_n(x)\).

The Chebyshev polynomials of the first kind are orthogonal on \([-1, 1]\) with respect to the weight \(\frac{1}{\sqrt{1-x^2}}\).

References

[R160]http://en.wikipedia.org/wiki/Chebyshev_polynomial
[R161]http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html
[R162]http://mathworld.wolfram.com/ChebyshevPolynomialoftheSecondKind.html
[R163]http://functions.wolfram.com/Polynomials/ChebyshevT/
[R164]http://functions.wolfram.com/Polynomials/ChebyshevU/

Examples

>>> from sympy import chebyshevt, chebyshevu, diff
>>> from sympy.abc import n,x
>>> chebyshevt(0, x)
1
>>> chebyshevt(1, x)
x
>>> chebyshevt(2, x)
2*x**2 - 1
>>> chebyshevt(n, x)
chebyshevt(n, x)
>>> chebyshevt(n, -x)
(-1)**n*chebyshevt(n, x)
>>> chebyshevt(-n, x)
chebyshevt(n, x)
>>> chebyshevt(n, 0)
cos(pi*n/2)
>>> chebyshevt(n, -1)
(-1)**n
>>> diff(chebyshevt(n, x), x)
n*chebyshevu(n - 1, x)
class sympy.functions.special.polynomials.chebyshevu

Chebyshev polynomial of the second kind, \(U_n(x)\)

chebyshevu(n, x) gives the nth Chebyshev polynomial of the second kind in x, \(U_n(x)\).

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

References

[R165]http://en.wikipedia.org/wiki/Chebyshev_polynomial
[R166]http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html
[R167]http://mathworld.wolfram.com/ChebyshevPolynomialoftheSecondKind.html
[R168]http://functions.wolfram.com/Polynomials/ChebyshevT/
[R169]http://functions.wolfram.com/Polynomials/ChebyshevU/

Examples

>>> from sympy import chebyshevt, chebyshevu, diff
>>> from sympy.abc import n,x
>>> chebyshevu(0, x)
1
>>> chebyshevu(1, x)
2*x
>>> chebyshevu(2, x)
4*x**2 - 1
>>> chebyshevu(n, x)
chebyshevu(n, x)
>>> chebyshevu(n, -x)
(-1)**n*chebyshevu(n, x)
>>> chebyshevu(-n, x)
-chebyshevu(n - 2, x)
>>> chebyshevu(n, 0)
cos(pi*n/2)
>>> chebyshevu(n, 1)
n + 1
>>> diff(chebyshevu(n, x), x)
(-x*chebyshevu(n, x) + (n + 1)*chebyshevt(n + 1, x))/(x**2 - 1)
class sympy.functions.special.polynomials.chebyshevt_root

chebyshev_root(n, k) returns the kth root (indexed from zero) of the nth Chebyshev polynomial of the first kind; that is, if 0 <= k < n, chebyshevt(n, chebyshevt_root(n, k)) == 0.

Examples

>>> from sympy import chebyshevt, chebyshevt_root
>>> chebyshevt_root(3, 2)
-sqrt(3)/2
>>> chebyshevt(3, chebyshevt_root(3, 2))
0
class sympy.functions.special.polynomials.chebyshevu_root

chebyshevu_root(n, k) returns the kth root (indexed from zero) of the nth Chebyshev polynomial of the second kind; that is, if 0 <= k < n, chebyshevu(n, chebyshevu_root(n, k)) == 0.

Examples

>>> from sympy import chebyshevu, chebyshevu_root
>>> chebyshevu_root(3, 2)
-sqrt(2)/2
>>> chebyshevu(3, chebyshevu_root(3, 2))
0

Legendre Polynomials

class sympy.functions.special.polynomials.legendre

legendre(n, x) gives the nth Legendre polynomial of x, \(P_n(x)\)

The Legendre polynomials are orthogonal on [-1, 1] with respect to the constant weight 1. They satisfy \(P_n(1) = 1\) for all n; further, \(P_n\) is odd for odd n and even for even n.

References

[R170]http://en.wikipedia.org/wiki/Legendre_polynomial
[R171]http://mathworld.wolfram.com/LegendrePolynomial.html
[R172]http://functions.wolfram.com/Polynomials/LegendreP/
[R173]http://functions.wolfram.com/Polynomials/LegendreP2/

Examples

>>> from sympy import legendre, diff
>>> from sympy.abc import x, n
>>> legendre(0, x)
1
>>> legendre(1, x)
x
>>> legendre(2, x)
3*x**2/2 - 1/2
>>> legendre(n, x)
legendre(n, x)
>>> diff(legendre(n,x), x)
n*(x*legendre(n, x) - legendre(n - 1, x))/(x**2 - 1)
class sympy.functions.special.polynomials.assoc_legendre

assoc_legendre(n,m, x) gives \(P_n^m(x)\), where n and m are the degree and order or an expression which is related to the nth order Legendre polynomial, \(P_n(x)\) in the following manner:

\[P_n^m(x) = (-1)^m (1 - x^2)^{\frac{m}{2}} \frac{\mathrm{d}^m P_n(x)}{\mathrm{d} x^m}\]

Associated Legendre polynomial are orthogonal on [-1, 1] with:

  • weight = 1 for the same m, and different n.
  • weight = 1/(1-x**2) for the same n, and different m.

References

[R174]http://en.wikipedia.org/wiki/Associated_Legendre_polynomials
[R175]http://mathworld.wolfram.com/LegendrePolynomial.html
[R176]http://functions.wolfram.com/Polynomials/LegendreP/
[R177]http://functions.wolfram.com/Polynomials/LegendreP2/

Examples

>>> from sympy import assoc_legendre
>>> from sympy.abc import x, m, n
>>> assoc_legendre(0,0, x)
1
>>> assoc_legendre(1,0, x)
x
>>> assoc_legendre(1,1, x)
-sqrt(-x**2 + 1)
>>> assoc_legendre(n,m,x)
assoc_legendre(n, m, x)

Hermite Polynomials

class sympy.functions.special.polynomials.hermite

hermite(n, x) gives the nth Hermite polynomial in x, \(H_n(x)\)

The Hermite polynomials are orthogonal on \((-\infty, \infty)\) with respect to the weight \(\exp\left(-\frac{x^2}{2}\right)\).

References

[R178]http://en.wikipedia.org/wiki/Hermite_polynomial
[R179]http://mathworld.wolfram.com/HermitePolynomial.html
[R180]http://functions.wolfram.com/Polynomials/HermiteH/

Examples

>>> from sympy import hermite, diff
>>> from sympy.abc import x, n
>>> hermite(0, x)
1
>>> hermite(1, x)
2*x
>>> hermite(2, x)
4*x**2 - 2
>>> hermite(n, x)
hermite(n, x)
>>> diff(hermite(n,x), x)
2*n*hermite(n - 1, x)
>>> diff(hermite(n,x), x)
2*n*hermite(n - 1, x)
>>> hermite(n, -x)
(-1)**n*hermite(n, x)

Laguerre Polynomials

class sympy.functions.special.polynomials.laguerre

Returns the nth Laguerre polynomial in x, \(L_n(x)\).

Parameters:

n : int

Degree of Laguerre polynomial. Must be n >= 0.

References

[R181]http://en.wikipedia.org/wiki/Laguerre_polynomial
[R182]http://mathworld.wolfram.com/LaguerrePolynomial.html
[R183]http://functions.wolfram.com/Polynomials/LaguerreL/
[R184]http://functions.wolfram.com/Polynomials/LaguerreL3/

Examples

>>> from sympy import laguerre, diff
>>> from sympy.abc import x, n
>>> laguerre(0, x)
1
>>> laguerre(1, x)
-x + 1
>>> laguerre(2, x)
x**2/2 - 2*x + 1
>>> laguerre(3, x)
-x**3/6 + 3*x**2/2 - 3*x + 1
>>> laguerre(n, x)
laguerre(n, x)
>>> diff(laguerre(n, x), x)
-assoc_laguerre(n - 1, 1, x)
class sympy.functions.special.polynomials.assoc_laguerre

Returns the nth generalized Laguerre polynomial in x, \(L_n(x)\).

Parameters:

n : int

Degree of Laguerre polynomial. Must be n >= 0.

alpha : Expr

Arbitrary expression. For alpha=0 regular Laguerre polynomials will be generated.

References

[R185]http://en.wikipedia.org/wiki/Laguerre_polynomial#Assoc_laguerre_polynomials
[R186]http://mathworld.wolfram.com/AssociatedLaguerrePolynomial.html
[R187]http://functions.wolfram.com/Polynomials/LaguerreL/
[R188]http://functions.wolfram.com/Polynomials/LaguerreL3/

Examples

>>> from sympy import laguerre, assoc_laguerre, diff
>>> from sympy.abc import x, n, a
>>> assoc_laguerre(0, a, x)
1
>>> assoc_laguerre(1, a, x)
a - x + 1
>>> assoc_laguerre(2, a, x)
a**2/2 + 3*a/2 + x**2/2 + x*(-a - 2) + 1
>>> assoc_laguerre(3, a, x)
a**3/6 + a**2 + 11*a/6 - x**3/6 + x**2*(a/2 + 3/2) +
    x*(-a**2/2 - 5*a/2 - 3) + 1
>>> assoc_laguerre(n, a, 0)
binomial(a + n, a)
>>> assoc_laguerre(n, a, x)
assoc_laguerre(n, a, x)
>>> assoc_laguerre(n, 0, x)
laguerre(n, x)
>>> diff(assoc_laguerre(n, a, x), x)
-assoc_laguerre(n - 1, a + 1, x)
>>> diff(assoc_laguerre(n, a, x), a)
Sum(assoc_laguerre(_k, a, x)/(-a + n), (_k, 0, n - 1))

Spherical Harmonics

class sympy.functions.special.spherical_harmonics.Ynm

Spherical harmonics defined as

\[Y_n^m(\theta, \varphi) := \sqrt{\frac{(2n+1)(n-m)!}{4\pi(n+m)!}} \exp(i m \varphi) \mathrm{P}_n^m\left(\cos(\theta)\right)\]

Ynm() gives the spherical harmonic function of order \(n\) and \(m\) in \(\theta\) and \(\varphi\), \(Y_n^m(\theta, \varphi)\). The four parameters are as follows: \(n \geq 0\) an integer and \(m\) an integer such that \(-n \leq m \leq n\) holds. The two angles are real-valued with \(\theta \in [0, \pi]\) and \(\varphi \in [0, 2\pi]\).

See also

Ynm_c, Znm

References

[R189]http://en.wikipedia.org/wiki/Spherical_harmonics
[R190]http://mathworld.wolfram.com/SphericalHarmonic.html
[R191]http://functions.wolfram.com/Polynomials/SphericalHarmonicY/
[R192]http://dlmf.nist.gov/14.30

Examples

>>> from sympy import Ynm, Symbol
>>> from sympy.abc import n,m
>>> theta = Symbol("theta")
>>> phi = Symbol("phi")
>>> Ynm(n, m, theta, phi)
Ynm(n, m, theta, phi)

Several symmetries are known, for the order

>>> from sympy import Ynm, Symbol
>>> from sympy.abc import n,m
>>> theta = Symbol("theta")
>>> phi = Symbol("phi")
>>> Ynm(n, -m, theta, phi)
(-1)**m*exp(-2*I*m*phi)*Ynm(n, m, theta, phi)

as well as for the angles

>>> from sympy import Ynm, Symbol, simplify
>>> from sympy.abc import n,m
>>> theta = Symbol("theta")
>>> phi = Symbol("phi")
>>> Ynm(n, m, -theta, phi)
Ynm(n, m, theta, phi)
>>> Ynm(n, m, theta, -phi)
exp(-2*I*m*phi)*Ynm(n, m, theta, phi)

For specific integers n and m we can evalute the harmonics to more useful expressions

>>> simplify(Ynm(0, 0, theta, phi).expand(func=True))
1/(2*sqrt(pi))
>>> simplify(Ynm(1, -1, theta, phi).expand(func=True))
sqrt(6)*exp(-I*phi)*sin(theta)/(4*sqrt(pi))
>>> simplify(Ynm(1, 0, theta, phi).expand(func=True))
sqrt(3)*cos(theta)/(2*sqrt(pi))
>>> simplify(Ynm(1, 1, theta, phi).expand(func=True))
-sqrt(6)*exp(I*phi)*sin(theta)/(4*sqrt(pi))
>>> simplify(Ynm(2, -2, theta, phi).expand(func=True))
sqrt(30)*exp(-2*I*phi)*sin(theta)**2/(8*sqrt(pi))
>>> simplify(Ynm(2, -1, theta, phi).expand(func=True))
sqrt(30)*exp(-I*phi)*sin(2*theta)/(8*sqrt(pi))
>>> simplify(Ynm(2, 0, theta, phi).expand(func=True))
sqrt(5)*(3*cos(theta)**2 - 1)/(4*sqrt(pi))
>>> simplify(Ynm(2, 1, theta, phi).expand(func=True))
-sqrt(30)*exp(I*phi)*sin(2*theta)/(8*sqrt(pi))
>>> simplify(Ynm(2, 2, theta, phi).expand(func=True))
sqrt(30)*exp(2*I*phi)*sin(theta)**2/(8*sqrt(pi))

We can differentiate the functions with respect to both angles

>>> from sympy import Ynm, Symbol, diff
>>> from sympy.abc import n,m
>>> theta = Symbol("theta")
>>> phi = Symbol("phi")
>>> diff(Ynm(n, m, theta, phi), theta)
m*cot(theta)*Ynm(n, m, theta, phi) + sqrt((-m + n)*(m + n + 1))*exp(-I*phi)*Ynm(n, m + 1, theta, phi)
>>> diff(Ynm(n, m, theta, phi), phi)
I*m*Ynm(n, m, theta, phi)

Further we can compute the complex conjugation

>>> from sympy import Ynm, Symbol, conjugate
>>> from sympy.abc import n,m
>>> theta = Symbol("theta")
>>> phi = Symbol("phi")
>>> conjugate(Ynm(n, m, theta, phi))
(-1)**(2*m)*exp(-2*I*m*phi)*Ynm(n, m, theta, phi)

To get back the well known expressions in spherical coordinates we use full expansion

>>> from sympy import Ynm, Symbol, expand_func
>>> from sympy.abc import n,m
>>> theta = Symbol("theta")
>>> phi = Symbol("phi")
>>> expand_func(Ynm(n, m, theta, phi))
sqrt((2*n + 1)*factorial(-m + n)/factorial(m + n))*exp(I*m*phi)*assoc_legendre(n, m, cos(theta))/(2*sqrt(pi))
sympy.functions.special.spherical_harmonics.Ynm_c(n, m, theta, phi)

Conjugate spherical harmonics defined as

\[\overline{Y_n^m(\theta, \varphi)} := (-1)^m Y_n^{-m}(\theta, \varphi)\]

See also

Ynm, Znm

References

[R193]http://en.wikipedia.org/wiki/Spherical_harmonics
[R194]http://mathworld.wolfram.com/SphericalHarmonic.html
[R195]http://functions.wolfram.com/Polynomials/SphericalHarmonicY/
class sympy.functions.special.spherical_harmonics.Znm

Real spherical harmonics defined as

\[\begin{split}Z_n^m(\theta, \varphi) := \begin{cases} \frac{Y_n^m(\theta, \varphi) + \overline{Y_n^m(\theta, \varphi)}}{\sqrt{2}} &\quad m > 0 \\ Y_n^m(\theta, \varphi) &\quad m = 0 \\ \frac{Y_n^m(\theta, \varphi) - \overline{Y_n^m(\theta, \varphi)}}{i \sqrt{2}} &\quad m < 0 \\ \end{cases}\end{split}\]

which gives in simplified form

\[\begin{split}Z_n^m(\theta, \varphi) = \begin{cases} \frac{Y_n^m(\theta, \varphi) + (-1)^m Y_n^{-m}(\theta, \varphi)}{\sqrt{2}} &\quad m > 0 \\ Y_n^m(\theta, \varphi) &\quad m = 0 \\ \frac{Y_n^m(\theta, \varphi) - (-1)^m Y_n^{-m}(\theta, \varphi)}{i \sqrt{2}} &\quad m < 0 \\ \end{cases}\end{split}\]

See also

Ynm, Ynm_c

References

[R196]http://en.wikipedia.org/wiki/Spherical_harmonics
[R197]http://mathworld.wolfram.com/SphericalHarmonic.html
[R198]http://functions.wolfram.com/Polynomials/SphericalHarmonicY/

Tensor Functions

sympy.functions.special.tensor_functions.Eijk(*args, **kwargs)

Represent the Levi-Civita symbol.

This is just compatibility wrapper to LeviCivita().

See also

LeviCivita

sympy.functions.special.tensor_functions.eval_levicivita(*args)

Evaluate Levi-Civita symbol.

class sympy.functions.special.tensor_functions.LeviCivita

Represent the Levi-Civita symbol.

For even permutations of indices it returns 1, for odd permutations -1, and for everything else (a repeated index) it returns 0.

Thus it represents an alternating pseudotensor.

See also

Eijk

Examples

>>> from sympy import LeviCivita
>>> from sympy.abc import i, j, k
>>> LeviCivita(1, 2, 3)
1
>>> LeviCivita(1, 3, 2)
-1
>>> LeviCivita(1, 2, 2)
0
>>> LeviCivita(i, j, k)
LeviCivita(i, j, k)
>>> LeviCivita(i, j, i)
0

Attributes

nargs  
class sympy.functions.special.tensor_functions.KroneckerDelta

The discrete, or Kronecker, delta function.

A function that takes in two integers \(i\) and \(j\). It returns \(0\) if \(i\) and \(j\) are not equal or it returns \(1\) if \(i\) and \(j\) are equal.

Parameters:

i : Number, Symbol

The first index of the delta function.

j : Number, Symbol

The second index of the delta function.

References

[R199]http://en.wikipedia.org/wiki/Kronecker_delta

Examples

A simple example with integer indices:

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> KroneckerDelta(1, 2)
0
>>> KroneckerDelta(3, 3)
1

Symbolic indices:

>>> from sympy.abc import i, j, k
>>> KroneckerDelta(i, j)
KroneckerDelta(i, j)
>>> KroneckerDelta(i, i)
1
>>> KroneckerDelta(i, i + 1)
0
>>> KroneckerDelta(i, i + 1 + k)
KroneckerDelta(i, i + k + 1)
classmethod eval(i, j)

Evaluates the discrete delta function.

Examples

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> from sympy.abc import i, j, k
>>> KroneckerDelta(i, j)
KroneckerDelta(i, j)
>>> KroneckerDelta(i, i)
1
>>> KroneckerDelta(i, i + 1)
0
>>> KroneckerDelta(i, i + 1 + k)
KroneckerDelta(i, i + k + 1)

# indirect doctest

indices_contain_equal_information

Returns True if indices are either both above or below fermi.

Examples

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> from sympy import Symbol
>>> a = Symbol('a', above_fermi=True)
>>> i = Symbol('i', below_fermi=True)
>>> p = Symbol('p')
>>> q = Symbol('q')
>>> KroneckerDelta(p, q).indices_contain_equal_information
True
>>> KroneckerDelta(p, q+1).indices_contain_equal_information
True
>>> KroneckerDelta(i, p).indices_contain_equal_information
False
is_above_fermi

True if Delta can be non-zero above fermi

See also

is_below_fermi, is_only_below_fermi, is_only_above_fermi

Examples

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> from sympy import Symbol
>>> a = Symbol('a', above_fermi=True)
>>> i = Symbol('i', below_fermi=True)
>>> p = Symbol('p')
>>> q = Symbol('q')
>>> KroneckerDelta(p, a).is_above_fermi
True
>>> KroneckerDelta(p, i).is_above_fermi
False
>>> KroneckerDelta(p, q).is_above_fermi
True
is_below_fermi

True if Delta can be non-zero below fermi

See also

is_above_fermi, is_only_above_fermi, is_only_below_fermi

Examples

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> from sympy import Symbol
>>> a = Symbol('a', above_fermi=True)
>>> i = Symbol('i', below_fermi=True)
>>> p = Symbol('p')
>>> q = Symbol('q')
>>> KroneckerDelta(p, a).is_below_fermi
False
>>> KroneckerDelta(p, i).is_below_fermi
True
>>> KroneckerDelta(p, q).is_below_fermi
True
is_only_above_fermi

True if Delta is restricted to above fermi

See also

is_above_fermi, is_below_fermi, is_only_below_fermi

Examples

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> from sympy import Symbol
>>> a = Symbol('a', above_fermi=True)
>>> i = Symbol('i', below_fermi=True)
>>> p = Symbol('p')
>>> q = Symbol('q')
>>> KroneckerDelta(p, a).is_only_above_fermi
True
>>> KroneckerDelta(p, q).is_only_above_fermi
False
>>> KroneckerDelta(p, i).is_only_above_fermi
False
is_only_below_fermi

True if Delta is restricted to below fermi

See also

is_above_fermi, is_below_fermi, is_only_above_fermi

Examples

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> from sympy import Symbol
>>> a = Symbol('a', above_fermi=True)
>>> i = Symbol('i', below_fermi=True)
>>> p = Symbol('p')
>>> q = Symbol('q')
>>> KroneckerDelta(p, i).is_only_below_fermi
True
>>> KroneckerDelta(p, q).is_only_below_fermi
False
>>> KroneckerDelta(p, a).is_only_below_fermi
False
killable_index

Returns the index which is preferred to substitute in the final expression.

The index to substitute is the index with less information regarding fermi level. If indices contain same information, ‘a’ is preferred before ‘b’.

See also

preferred_index

Examples

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> from sympy import Symbol
>>> a = Symbol('a', above_fermi=True)
>>> i = Symbol('i', below_fermi=True)
>>> j = Symbol('j', below_fermi=True)
>>> p = Symbol('p')
>>> KroneckerDelta(p, i).killable_index
p
>>> KroneckerDelta(p, a).killable_index
p
>>> KroneckerDelta(i, j).killable_index
j
preferred_index

Returns the index which is preferred to keep in the final expression.

The preferred index is the index with more information regarding fermi level. If indices contain same information, ‘a’ is preferred before ‘b’.

See also

killable_index

Examples

>>> from sympy.functions.special.tensor_functions import KroneckerDelta
>>> from sympy import Symbol
>>> a = Symbol('a', above_fermi=True)
>>> i = Symbol('i', below_fermi=True)
>>> j = Symbol('j', below_fermi=True)
>>> p = Symbol('p')
>>> KroneckerDelta(p, i).preferred_index
i
>>> KroneckerDelta(p, a).preferred_index
a
>>> KroneckerDelta(i, j).preferred_index
i