Miscellaneous special functions#
This module provides easy access to many of Maxima and PARI’s special functions.
Maxima’s special functions package (which includes spherical harmonic functions, spherical Bessel functions (of the 1st and 2nd kind), and spherical Hankel functions (of the 1st and 2nd kind)) was written by Barton Willis of the University of Nebraska at Kearney. It is released under the terms of the General Public License (GPL).
Support for elliptic functions and integrals was written by Raymond Toy. It is placed under the terms of the General Public License (GPL) that governs the distribution of Maxima.
Next, we summarize some of the properties of the functions implemented here.
Spherical harmonics: Laplace’s equation in spherical coordinates is:
\[\frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial f}{\partial r} \right) + \frac{1}{r^2\sin\theta} \frac{\partial}{\partial \theta} \left( \sin\theta \frac{\partial f}{\partial \theta} \right) + \frac{1}{r^2\sin^2\theta} \frac{\partial^2 f}{\partial \varphi^2} = 0.\]Note that the spherical coordinates \(\theta\) and \(\varphi\) are defined here as follows: \(\theta\) is the colatitude or polar angle, ranging from \(0\leq\theta\leq\pi\) and \(\varphi\) the azimuth or longitude, ranging from \(0\leq\varphi<2\pi\).
The general solution which remains finite towards infinity is a linear combination of functions of the form
\[r^{-1-\ell} \cos (m \varphi) P_\ell^m (\cos{\theta} )\]and
\[r^{-1-\ell} \sin (m \varphi) P_\ell^m (\cos{\theta} )\]where \(P_\ell^m\) are the associated Legendre polynomials (cf.
Func_assoc_legendre_P), and with integer parameters \(\ell \ge 0\) and \(m\) from \(0\) to \(\ell\). Put in another way, the solutions with integer parameters \(\ell \ge 0\) and \(- \ell\leq m\leq \ell\), can be written as linear combinations of:\[U_{\ell,m}(r,\theta , \varphi ) = r^{-1-\ell} Y_\ell^m( \theta , \varphi )\]where the functions \(Y\) are the spherical harmonic functions with parameters \(\ell\), \(m\), which can be written as:
\[Y_\ell^m( \theta , \varphi ) = \sqrt{ \frac{(2\ell+1)}{4\pi} \frac{(\ell-m)!}{(\ell+m)!} } \, e^{i m \varphi } \, P_\ell^m ( \cos{\theta} ) .\]The spherical harmonics obey the normalisation condition
\[\int_{\theta=0}^\pi\int_{\varphi=0}^{2\pi} Y_\ell^mY_{\ell'}^{m'*}\,d\Omega = \delta_{\ell\ell'}\delta_{mm'}\quad\quad d\Omega = \sin\theta\,d\varphi\,d\theta .\]The incomplete elliptic integrals (of the first kind, etc.) are:
\[\begin{split}\begin{array}{c} \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2}}\, dx,\\ \displaystyle\int_0^\phi \sqrt{1 - m\sin(x)^2}\, dx,\\ \displaystyle\int_0^\phi \frac{\sqrt{1-mt^2}}{\sqrt(1 - t^2)}\, dx,\\ \displaystyle\int_0^\phi \frac{1}{\sqrt{1 - m\sin(x)^2\sqrt{1 - n\sin(x)^2}}}\, dx, \end{array}\end{split}\]and the complete ones are obtained by taking \(\phi =\pi/2\).
Warning
SciPy’s versions are poorly documented and seem less accurate than the Maxima and PARI versions. Typically they are limited by hardware floats precision.
REFERENCES:
Abramowitz and Stegun: Handbook of Mathematical Functions [AS1964]
AUTHORS:
David Joyner (2006-13-06): initial version
David Joyner (2006-30-10): bug fixes to pari wrappers of Bessel functions, hypergeometric_U
William Stein (2008-02): Impose some sanity checks.
David Joyner (2008-02-16): optional calls to scipy and replace all
#randomby...David Joyner (2008-04-23): addition of elliptic integrals
Eviatar Bach (2013): making elliptic integrals symbolic
Eric Gourgoulhon (2022): add Condon-Shortley phase to spherical harmonics
- class sage.functions.special.EllipticE[source]#
Bases:
BuiltinFunctionReturn the incomplete elliptic integral of the second kind:
\[E(\varphi\,|\,m)=\int_0^\varphi \sqrt{1 - m\sin(x)^2}\, dx.\]EXAMPLES:
sage: z = var("z") # needs sage.symbolic sage: elliptic_e(z, 1) # needs sage.symbolic elliptic_e(z, 1) sage: elliptic_e(z, 1).simplify() # not tested # needs sage.symbolic 2*round(z/pi) - sin(pi*round(z/pi) - z) sage: elliptic_e(z, 0) # needs sage.symbolic z sage: elliptic_e(0.5, 0.1) # abs tol 2e-15 # needs mpmath 0.498011394498832 sage: elliptic_e(1/2, 1/10).n(200) # needs sage.symbolic 0.4980113944988315331154610406...
>>> from sage.all import * >>> z = var("z") # needs sage.symbolic >>> elliptic_e(z, Integer(1)) # needs sage.symbolic elliptic_e(z, 1) >>> elliptic_e(z, Integer(1)).simplify() # not tested # needs sage.symbolic 2*round(z/pi) - sin(pi*round(z/pi) - z) >>> elliptic_e(z, Integer(0)) # needs sage.symbolic z >>> elliptic_e(RealNumber('0.5'), RealNumber('0.1')) # abs tol 2e-15 # needs mpmath 0.498011394498832 >>> elliptic_e(Integer(1)/Integer(2), Integer(1)/Integer(10)).n(Integer(200)) # needs sage.symbolic 0.4980113944988315331154610406...
See also
Taking \(\varphi = \pi/2\) gives
elliptic_ec().Taking \(\varphi = \operatorname{arc\,sin}(\operatorname{sn}(u,m))\) gives
elliptic_eu().
REFERENCES:
- class sage.functions.special.EllipticEC[source]#
Bases:
BuiltinFunctionReturn the complete elliptic integral of the second kind:
\[E(m)=\int_0^{\pi/2} \sqrt{1 - m\sin(x)^2}\, dx.\]EXAMPLES:
sage: elliptic_ec(0.1) # needs mpmath 1.53075763689776 sage: elliptic_ec(x).diff() # needs sage.symbolic 1/2*(elliptic_ec(x) - elliptic_kc(x))/x
>>> from sage.all import * >>> elliptic_ec(RealNumber('0.1')) # needs mpmath 1.53075763689776 >>> elliptic_ec(x).diff() # needs sage.symbolic 1/2*(elliptic_ec(x) - elliptic_kc(x))/x
See also
REFERENCES:
- class sage.functions.special.EllipticEU[source]#
Bases:
BuiltinFunctionReturn Jacobi’s form of the incomplete elliptic integral of the second kind:
\[E(u,m)= \int_0^u \mathrm{dn}(x,m)^2\, dx = \int_0^\tau \frac{\sqrt{1-m x^2}}{\sqrt{1-x^2}}\, dx.\]where \(\tau = \mathrm{sn}(u, m)\).
Also,
elliptic_eu(u, m) = elliptic_e(asin(sn(u,m)),m).EXAMPLES:
sage: elliptic_eu(0.5, 0.1) # needs mpmath 0.496054551286597
>>> from sage.all import * >>> elliptic_eu(RealNumber('0.5'), RealNumber('0.1')) # needs mpmath 0.496054551286597
See also
REFERENCES:
- class sage.functions.special.EllipticF[source]#
Bases:
BuiltinFunctionReturn the incomplete elliptic integral of the first kind.
\[F(\varphi\,|\,m)=\int_0^\varphi \frac{dx}{\sqrt{1 - m\sin(x)^2}},\]Taking \(\varphi = \pi/2\) gives
elliptic_kc().EXAMPLES:
sage: z = var("z") # needs sage.symbolic sage: elliptic_f(z, 0) # needs sage.symbolic z sage: elliptic_f(z, 1).simplify() # needs sage.symbolic log(tan(1/4*pi + 1/2*z)) sage: elliptic_f(0.2, 0.1) # needs mpmath 0.200132506747543
>>> from sage.all import * >>> z = var("z") # needs sage.symbolic >>> elliptic_f(z, Integer(0)) # needs sage.symbolic z >>> elliptic_f(z, Integer(1)).simplify() # needs sage.symbolic log(tan(1/4*pi + 1/2*z)) >>> elliptic_f(RealNumber('0.2'), RealNumber('0.1')) # needs mpmath 0.200132506747543
See also
REFERENCES:
- class sage.functions.special.EllipticKC[source]#
Bases:
BuiltinFunctionReturn the complete elliptic integral of the first kind:
\[K(m)=\int_0^{\pi/2} \frac{dx}{\sqrt{1 - m\sin(x)^2}}.\]EXAMPLES:
sage: elliptic_kc(0.5) # needs mpmath 1.85407467730137
>>> from sage.all import * >>> elliptic_kc(RealNumber('0.5')) # needs mpmath 1.85407467730137
See also
REFERENCES:
- class sage.functions.special.EllipticPi[source]#
Bases:
BuiltinFunctionReturn the incomplete elliptic integral of the third kind:
\[\Pi(n, t, m) = \int_0^t \frac{dx}{(1 - n \sin(x)^2)\sqrt{1 - m \sin(x)^2}}.\]INPUT:
n– a real number, called the “characteristic”t– a real number, called the “amplitude”m– a real number, called the “parameter”
EXAMPLES:
sage: N(elliptic_pi(1, pi/4, 1)) # needs sage.symbolic 1.14779357469632
>>> from sage.all import * >>> N(elliptic_pi(Integer(1), pi/Integer(4), Integer(1))) # needs sage.symbolic 1.14779357469632
Compare the value computed by Maxima to the definition as a definite integral (using GSL):
sage: elliptic_pi(0.1, 0.2, 0.3) # needs mpmath 0.200665068220979 sage: numerical_integral(1/(1-0.1*sin(x)^2)/sqrt(1-0.3*sin(x)^2), 0.0, 0.2) # needs sage.symbolic (0.2006650682209791, 2.227829789769088e-15)
>>> from sage.all import * >>> elliptic_pi(RealNumber('0.1'), RealNumber('0.2'), RealNumber('0.3')) # needs mpmath 0.200665068220979 >>> numerical_integral(Integer(1)/(Integer(1)-RealNumber('0.1')*sin(x)**Integer(2))/sqrt(Integer(1)-RealNumber('0.3')*sin(x)**Integer(2)), RealNumber('0.0'), RealNumber('0.2')) # needs sage.symbolic (0.2006650682209791, 2.227829789769088e-15)
REFERENCES:
- class sage.functions.special.SphericalHarmonic[source]#
Bases:
BuiltinFunctionReturns the spherical harmonic function \(Y_n^m(\theta, \varphi)\).
For integers \(n > -1\), \(|m| \leq n\), simplification is done automatically. Numeric evaluation is supported for complex \(n\) and \(m\).
EXAMPLES:
sage: # needs sage.symbolic sage: x, y = var('x, y') sage: spherical_harmonic(3, 2, x, y) 1/8*sqrt(30)*sqrt(7)*cos(x)*e^(2*I*y)*sin(x)^2/sqrt(pi) sage: spherical_harmonic(3, 2, 1, 2) 1/8*sqrt(30)*sqrt(7)*cos(1)*e^(4*I)*sin(1)^2/sqrt(pi) sage: spherical_harmonic(3 + I, 2., 1, 2) -0.351154337307488 - 0.415562233975369*I sage: latex(spherical_harmonic(3, 2, x, y, hold=True)) Y_{3}^{2}\left(x, y\right) sage: spherical_harmonic(1, 2, x, y) 0
>>> from sage.all import * >>> # needs sage.symbolic >>> x, y = var('x, y') >>> spherical_harmonic(Integer(3), Integer(2), x, y) 1/8*sqrt(30)*sqrt(7)*cos(x)*e^(2*I*y)*sin(x)^2/sqrt(pi) >>> spherical_harmonic(Integer(3), Integer(2), Integer(1), Integer(2)) 1/8*sqrt(30)*sqrt(7)*cos(1)*e^(4*I)*sin(1)^2/sqrt(pi) >>> spherical_harmonic(Integer(3) + I, RealNumber('2.'), Integer(1), Integer(2)) -0.351154337307488 - 0.415562233975369*I >>> latex(spherical_harmonic(Integer(3), Integer(2), x, y, hold=True)) Y_{3}^{2}\left(x, y\right) >>> spherical_harmonic(Integer(1), Integer(2), x, y) 0
The degree \(n\) and the order \(m\) can be symbolic:
sage: # needs sage.symbolic sage: n, m = var('n m') sage: spherical_harmonic(n, m, x, y) spherical_harmonic(n, m, x, y) sage: latex(spherical_harmonic(n, m, x, y)) Y_{n}^{m}\left(x, y\right) sage: diff(spherical_harmonic(n, m, x, y), x) m*cot(x)*spherical_harmonic(n, m, x, y) + sqrt(-(m + n + 1)*(m - n))*e^(-I*y)*spherical_harmonic(n, m + 1, x, y) sage: diff(spherical_harmonic(n, m, x, y), y) I*m*spherical_harmonic(n, m, x, y)
>>> from sage.all import * >>> # needs sage.symbolic >>> n, m = var('n m') >>> spherical_harmonic(n, m, x, y) spherical_harmonic(n, m, x, y) >>> latex(spherical_harmonic(n, m, x, y)) Y_{n}^{m}\left(x, y\right) >>> diff(spherical_harmonic(n, m, x, y), x) m*cot(x)*spherical_harmonic(n, m, x, y) + sqrt(-(m + n + 1)*(m - n))*e^(-I*y)*spherical_harmonic(n, m + 1, x, y) >>> diff(spherical_harmonic(n, m, x, y), y) I*m*spherical_harmonic(n, m, x, y)
The convention regarding the Condon-Shortley phase \((-1)^m\) is the same as for SymPy’s spherical harmonics and Wikipedia article Spherical_harmonics:
sage: # needs sage.symbolic sage: spherical_harmonic(1, 1, x, y) -1/4*sqrt(3)*sqrt(2)*e^(I*y)*sin(x)/sqrt(pi) sage: from sympy import Ynm # needs sympy sage: Ynm(1, 1, x, y).expand(func=True) # needs sympy -sqrt(6)*exp(I*y)*sin(x)/(4*sqrt(pi)) sage: spherical_harmonic(1, 1, x, y) - Ynm(1, 1, x, y) # needs sympy 0
>>> from sage.all import * >>> # needs sage.symbolic >>> spherical_harmonic(Integer(1), Integer(1), x, y) -1/4*sqrt(3)*sqrt(2)*e^(I*y)*sin(x)/sqrt(pi) >>> from sympy import Ynm # needs sympy >>> Ynm(Integer(1), Integer(1), x, y).expand(func=True) # needs sympy -sqrt(6)*exp(I*y)*sin(x)/(4*sqrt(pi)) >>> spherical_harmonic(Integer(1), Integer(1), x, y) - Ynm(Integer(1), Integer(1), x, y) # needs sympy 0
It also agrees with SciPy’s spherical harmonics:
sage: spherical_harmonic(1, 1, pi/2, pi).n() # abs tol 1e-14 # needs sage.symbolic 0.345494149471335 sage: from scipy.special import sph_harm # NB: arguments x and y are swapped # needs scipy sage: sph_harm(1, 1, pi.n(), (pi/2).n()) # abs tol 1e-14 # needs scipy sage.symbolic (0.3454941494713355-4.231083042742082e-17j)
>>> from sage.all import * >>> spherical_harmonic(Integer(1), Integer(1), pi/Integer(2), pi).n() # abs tol 1e-14 # needs sage.symbolic 0.345494149471335 >>> from scipy.special import sph_harm # NB: arguments x and y are swapped # needs scipy >>> sph_harm(Integer(1), Integer(1), pi.n(), (pi/Integer(2)).n()) # abs tol 1e-14 # needs scipy sage.symbolic (0.3454941494713355-4.231083042742082e-17j)
Note that this convention differs from the one in Maxima, as revealed by the sign difference for odd values of \(m\):
sage: maxima.spherical_harmonic(1, 1, x, y).sage() # needs sage.symbolic 1/2*sqrt(3/2)*e^(I*y)*sin(x)/sqrt(pi)
>>> from sage.all import * >>> maxima.spherical_harmonic(Integer(1), Integer(1), x, y).sage() # needs sage.symbolic 1/2*sqrt(3/2)*e^(I*y)*sin(x)/sqrt(pi)
It follows that, contrary to Maxima, SageMath uses the same sign convention for spherical harmonics as SymPy, SciPy, Mathematica and Wikipedia article Table_of_spherical_harmonics.
REFERENCES:
- sage.functions.special.elliptic_eu_f(u, m)[source]#
Internal function for numeric evaluation of
elliptic_eu, defined as \(E\left(\operatorname{am}(u, m)|m\right)\), where \(E\) is the incomplete elliptic integral of the second kind and \(\operatorname{am}\) is the Jacobi amplitude function.EXAMPLES:
sage: from sage.functions.special import elliptic_eu_f sage: elliptic_eu_f(0.5, 0.1) # needs mpmath mpf('0.49605455128659691')
>>> from sage.all import * >>> from sage.functions.special import elliptic_eu_f >>> elliptic_eu_f(RealNumber('0.5'), RealNumber('0.1')) # needs mpmath mpf('0.49605455128659691')
- sage.functions.special.elliptic_j(z, prec=53)[source]#
Returns the elliptic modular \(j\)-function evaluated at \(z\).
INPUT:
z(complex) – a complex number with positive imaginary part.prec(default: 53) – precision in bits for the complex field.
OUTPUT:
(complex) The value of \(j(z)\).
ALGORITHM:
Calls the
parifunctionellj().AUTHOR:
John Cremona
EXAMPLES:
sage: elliptic_j(CC(i)) # needs sage.rings.real_mpfr 1728.00000000000 sage: elliptic_j(sqrt(-2.0)) # needs sage.rings.complex_double 8000.00000000000 sage: z = ComplexField(100)(1, sqrt(11))/2 # needs sage.rings.real_mpfr sage.symbolic sage: elliptic_j(z) # needs sage.rings.real_mpfr sage.symbolic -32768.000... sage: elliptic_j(z).real().round() # needs sage.rings.real_mpfr sage.symbolic -32768
>>> from sage.all import * >>> elliptic_j(CC(i)) # needs sage.rings.real_mpfr 1728.00000000000 >>> elliptic_j(sqrt(-RealNumber('2.0'))) # needs sage.rings.complex_double 8000.00000000000 >>> z = ComplexField(Integer(100))(Integer(1), sqrt(Integer(11)))/Integer(2) # needs sage.rings.real_mpfr sage.symbolic >>> elliptic_j(z) # needs sage.rings.real_mpfr sage.symbolic -32768.000... >>> elliptic_j(z).real().round() # needs sage.rings.real_mpfr sage.symbolic -32768
sage: tau = (1 + sqrt(-163))/2 # needs sage.symbolic sage: (-elliptic_j(tau.n(100)).real().round())^(1/3) # needs sage.symbolic 640320
>>> from sage.all import * >>> tau = (Integer(1) + sqrt(-Integer(163)))/Integer(2) # needs sage.symbolic >>> (-elliptic_j(tau.n(Integer(100))).real().round())**(Integer(1)/Integer(3)) # needs sage.symbolic 640320
This example shows the need for higher precision than the default one of the \(ComplexField\), see Issue #28355:
sage: # needs sage.symbolic sage: -elliptic_j(tau) # rel tol 1e-2 2.62537412640767e17 - 732.558854258998*I sage: -elliptic_j(tau, 75) # rel tol 1e-2 2.625374126407680000000e17 - 0.0001309913593909879441262*I sage: -elliptic_j(tau, 100) # rel tol 1e-2 2.6253741264076799999999999999e17 - 1.3012822400356887122945119790e-12*I sage: (-elliptic_j(tau, 100).real().round())^(1/3) 640320
>>> from sage.all import * >>> # needs sage.symbolic >>> -elliptic_j(tau) # rel tol 1e-2 2.62537412640767e17 - 732.558854258998*I >>> -elliptic_j(tau, Integer(75)) # rel tol 1e-2 2.625374126407680000000e17 - 0.0001309913593909879441262*I >>> -elliptic_j(tau, Integer(100)) # rel tol 1e-2 2.6253741264076799999999999999e17 - 1.3012822400356887122945119790e-12*I >>> (-elliptic_j(tau, Integer(100)).real().round())**(Integer(1)/Integer(3)) 640320