Bessel functions and related functions¶
The functions in this section arise as solutions to various differential equations in physics, typically describing wavelike oscillatory behavior or a combination of oscillation and exponential decay or growth. Mathematically, they are special cases of the confluent hypergeometric functions \(\,_0F_1\), \(\,_1F_1\) and \(\,_1F_2\) (see Hypergeometric functions).
Bessel functions¶
besselj()¶
- mpmath.besselj(n, x, derivative=0)¶
besselj(n, x, derivative=0)gives the Bessel function of the first kind \(J_n(x)\). Bessel functions of the first kind are defined as solutions of the differential equation\[x^2 y'' + x y' + (x^2 - n^2) y = 0\]which appears, among other things, when solving the radial part of Laplace’s equation in cylindrical coordinates. This equation has two solutions for given \(n\), where the \(J_n\)-function is the solution that is nonsingular at \(x = 0\). For positive integer \(n\), \(J_n(x)\) behaves roughly like a sine (odd \(n\)) or cosine (even \(n\)) multiplied by a magnitude factor that decays slowly as \(x \to \pm\infty\).
Generally, \(J_n\) is a special case of the hypergeometric function \(\,_0F_1\):
\[J_n(x) = \frac{x^n}{2^n \Gamma(n+1)} \,_0F_1\left(n+1,-\frac{x^2}{4}\right)\]With derivative = \(m \ne 0\), the \(m\)-th derivative
\[\frac{d^m}{dx^m} J_n(x)\]is computed.
Plots
# Bessel function J_n(x) on the real line for n=0,1,2,3 j0 = lambda x: besselj(0,x) j1 = lambda x: besselj(1,x) j2 = lambda x: besselj(2,x) j3 = lambda x: besselj(3,x) plot([j0,j1,j2,j3],[0,14])
# Bessel function J_n(z) in the complex plane cplot(lambda z: besselj(1,z), [-8,8], [-8,8], points=50000)
Examples
Evaluation is supported for arbitrary arguments, and at arbitrary precision:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> besselj(2, 1000) -0.024777229528606 >>> besselj(4, 0.75) 0.000801070086542314 >>> besselj(2, 1000j) (-2.48071721019185e+432 + 6.41567059811949e-437j) >>> mp.dps = 25 >>> besselj(0.75j, 3+4j) (-2.778118364828153309919653 - 1.5863603889018621585533j) >>> mp.dps = 50 >>> besselj(1, pi) 0.28461534317975275734531059968613140570981118184947
Arguments may be large:
>>> mp.dps = 25 >>> besselj(0, 10000) -0.007096160353388801477265164 >>> besselj(0, 10**10) 0.000002175591750246891726859055 >>> besselj(2, 10**100) 7.337048736538615712436929e-51 >>> besselj(2, 10**5*j) (-3.540725411970948860173735e+43426 + 4.4949812409615803110051e-43433j)
The Bessel functions of the first kind satisfy simple symmetries around \(x = 0\):
>>> mp.dps = 15 >>> nprint([besselj(n,0) for n in range(5)]) [1.0, 0.0, 0.0, 0.0, 0.0] >>> nprint([besselj(n,pi) for n in range(5)]) [-0.304242, 0.284615, 0.485434, 0.333458, 0.151425] >>> nprint([besselj(n,-pi) for n in range(5)]) [-0.304242, -0.284615, 0.485434, -0.333458, 0.151425]
Roots of Bessel functions are often used:
>>> nprint([findroot(j0, k) for k in [2, 5, 8, 11, 14]]) [2.40483, 5.52008, 8.65373, 11.7915, 14.9309] >>> nprint([findroot(j1, k) for k in [3, 7, 10, 13, 16]]) [3.83171, 7.01559, 10.1735, 13.3237, 16.4706]
The roots are not periodic, but the distance between successive roots asymptotically approaches \(2 \pi\). Bessel functions of the first kind have the following normalization:
>>> quadosc(j0, [0, inf], period=2*pi) 1.0 >>> quadosc(j1, [0, inf], period=2*pi) 1.0
For \(n = 1/2\) or \(n = -1/2\), the Bessel function reduces to a trigonometric function:
>>> x = 10 >>> besselj(0.5, x), sqrt(2/(pi*x))*sin(x) (-0.13726373575505, -0.13726373575505) >>> besselj(-0.5, x), sqrt(2/(pi*x))*cos(x) (-0.211708866331398, -0.211708866331398)
Derivatives of any order can be computed (negative orders correspond to integration):
>>> mp.dps = 25 >>> besselj(0, 7.5, 1) -0.1352484275797055051822405 >>> diff(lambda x: besselj(0,x), 7.5) -0.1352484275797055051822405 >>> besselj(0, 7.5, 10) -0.1377811164763244890135677 >>> diff(lambda x: besselj(0,x), 7.5, 10) -0.1377811164763244890135677 >>> besselj(0,7.5,-1) - besselj(0,3.5,-1) -0.1241343240399987693521378 >>> quad(j0, [3.5, 7.5]) -0.1241343240399987693521378
Differentiation with a noninteger order gives the fractional derivative in the sense of the Riemann-Liouville differintegral, as computed by
differint():>>> mp.dps = 15 >>> besselj(1, 3.5, 0.75) -0.385977722939384 >>> differint(lambda x: besselj(1, x), 3.5, 0.75) -0.385977722939384
bessely()¶
- mpmath.bessely(n, x, derivative=0)¶
bessely(n, x, derivative=0)gives the Bessel function of the second kind,\[Y_n(x) = \frac{J_n(x) \cos(\pi n) - J_{-n}(x)}{\sin(\pi n)}.\]For \(n\) an integer, this formula should be understood as a limit. With derivative = \(m \ne 0\), the \(m\)-th derivative
\[\frac{d^m}{dx^m} Y_n(x)\]is computed.
Plots
# Bessel function of 2nd kind Y_n(x) on the real line for n=0,1,2,3 y0 = lambda x: bessely(0,x) y1 = lambda x: bessely(1,x) y2 = lambda x: bessely(2,x) y3 = lambda x: bessely(3,x) plot([y0,y1,y2,y3],[0,10],[-4,1])
# Bessel function of 2nd kind Y_n(z) in the complex plane cplot(lambda z: bessely(1,z), [-8,8], [-8,8], points=50000)
Examples
Some values of \(Y_n(x)\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> bessely(0,0), bessely(1,0), bessely(2,0) (-inf, -inf, -inf) >>> bessely(1, pi) 0.3588729167767189594679827 >>> bessely(0.5, 3+4j) (9.242861436961450520325216 - 3.085042824915332562522402j)
Arguments may be large:
>>> bessely(0, 10000) 0.00364780555898660588668872 >>> bessely(2.5, 10**50) -4.8952500412050989295774e-26 >>> bessely(2.5, -10**50) (0.0 + 4.8952500412050989295774e-26j)
Derivatives and antiderivatives of any order can be computed:
>>> bessely(2, 3.5, 1) 0.3842618820422660066089231 >>> diff(lambda x: bessely(2, x), 3.5) 0.3842618820422660066089231 >>> bessely(0.5, 3.5, 1) -0.2066598304156764337900417 >>> diff(lambda x: bessely(0.5, x), 3.5) -0.2066598304156764337900417 >>> diff(lambda x: bessely(2, x), 0.5, 10) -208173867409.5547350101511 >>> bessely(2, 0.5, 10) -208173867409.5547350101511 >>> bessely(2, 100.5, 100) 0.02668487547301372334849043 >>> quad(lambda x: bessely(2,x), [1,3]) -1.377046859093181969213262 >>> bessely(2,3,-1) - bessely(2,1,-1) -1.377046859093181969213262
besseli()¶
- mpmath.besseli(n, x, derivative=0)¶
besseli(n, x, derivative=0)gives the modified Bessel function of the first kind,\[I_n(x) = i^{-n} J_n(ix).\]With derivative = \(m \ne 0\), the \(m\)-th derivative
\[\frac{d^m}{dx^m} I_n(x)\]is computed.
Plots
# Modified Bessel function I_n(x) on the real line for n=0,1,2,3 i0 = lambda x: besseli(0,x) i1 = lambda x: besseli(1,x) i2 = lambda x: besseli(2,x) i3 = lambda x: besseli(3,x) plot([i0,i1,i2,i3],[0,5],[0,5])
# Modified Bessel function I_n(z) in the complex plane cplot(lambda z: besseli(1,z), [-8,8], [-8,8], points=50000)
Examples
Some values of \(I_n(x)\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> besseli(0,0) 1.0 >>> besseli(1,0) 0.0 >>> besseli(0,1) 1.266065877752008335598245 >>> besseli(3.5, 2+3j) (-0.2904369752642538144289025 - 0.4469098397654815837307006j)
Arguments may be large:
>>> besseli(2, 1000) 2.480717210191852440616782e+432 >>> besseli(2, 10**10) 4.299602851624027900335391e+4342944813 >>> besseli(2, 6000+10000j) (-2.114650753239580827144204e+2603 + 4.385040221241629041351886e+2602j)
For integers \(n\), the following integral representation holds:
>>> mp.dps = 15 >>> n = 3 >>> x = 2.3 >>> quad(lambda t: exp(x*cos(t))*cos(n*t), [0,pi])/pi 0.349223221159309 >>> besseli(n,x) 0.349223221159309
Derivatives and antiderivatives of any order can be computed:
>>> mp.dps = 25 >>> besseli(2, 7.5, 1) 195.8229038931399062565883 >>> diff(lambda x: besseli(2,x), 7.5) 195.8229038931399062565883 >>> besseli(2, 7.5, 10) 153.3296508971734525525176 >>> diff(lambda x: besseli(2,x), 7.5, 10) 153.3296508971734525525176 >>> besseli(2,7.5,-1) - besseli(2,3.5,-1) 202.5043900051930141956876 >>> quad(lambda x: besseli(2,x), [3.5, 7.5]) 202.5043900051930141956876
besselk()¶
- mpmath.besselk(n, x)¶
besselk(n, x)gives the modified Bessel function of the second kind,\[K_n(x) = \frac{\pi}{2} \frac{I_{-n}(x)-I_{n}(x)}{\sin(\pi n)}\]For \(n\) an integer, this formula should be understood as a limit.
Plots
# Modified Bessel function of 2nd kind K_n(x) on the real line for n=0,1,2,3 k0 = lambda x: besselk(0,x) k1 = lambda x: besselk(1,x) k2 = lambda x: besselk(2,x) k3 = lambda x: besselk(3,x) plot([k0,k1,k2,k3],[0,8],[0,5])
# Modified Bessel function of 2nd kind K_n(z) in the complex plane cplot(lambda z: besselk(1,z), [-8,8], [-8,8], points=50000)
Examples
Evaluation is supported for arbitrary complex arguments:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> besselk(0,1) 0.4210244382407083333356274 >>> besselk(0, -1) (0.4210244382407083333356274 - 3.97746326050642263725661j) >>> besselk(3.5, 2+3j) (-0.02090732889633760668464128 + 0.2464022641351420167819697j) >>> besselk(2+3j, 0.5) (0.9615816021726349402626083 + 0.1918250181801757416908224j)
Arguments may be large:
>>> besselk(0, 100) 4.656628229175902018939005e-45 >>> besselk(1, 10**6) 4.131967049321725588398296e-434298 >>> besselk(1, 10**6*j) (0.001140348428252385844876706 - 0.0005200017201681152909000961j) >>> besselk(4.5, fmul(10**50, j, exact=True)) (1.561034538142413947789221e-26 + 1.243554598118700063281496e-25j)
The point \(x = 0\) is a singularity (logarithmic if \(n = 0\)):
>>> besselk(0,0) +inf >>> besselk(1,0) +inf >>> for n in range(-4, 5): ... print(besselk(n, '1e-1000')) ... 4.8e+4001 8.0e+3000 2.0e+2000 1.0e+1000 2302.701024509704096466802 1.0e+1000 2.0e+2000 8.0e+3000 4.8e+4001
Bessel function zeros¶
besseljzero()¶
- mpmath.besseljzero(v, m, derivative=0)¶
For a real order \(\nu \ge 0\) and a positive integer \(m\), returns \(j_{\nu,m}\), the \(m\)-th positive zero of the Bessel function of the first kind \(J_{\nu}(z)\) (see
besselj()). Alternatively, with derivative=1, gives the first nonnegative simple zero \(j'_{\nu,m}\) of \(J'_{\nu}(z)\).The indexing convention is that used by Abramowitz & Stegun and the DLMF. Note the special case \(j'_{0,1} = 0\), while all other zeros are positive. In effect, only simple zeros are counted (all zeros of Bessel functions are simple except possibly \(z = 0\)) and \(j_{\nu,m}\) becomes a monotonic function of both \(\nu\) and \(m\).
The zeros are interlaced according to the inequalities
\[ \begin{align}\begin{aligned}j'_{\nu,k} < j_{\nu,k} < j'_{\nu,k+1}\\j_{\nu,1} < j_{\nu+1,2} < j_{\nu,2} < j_{\nu+1,2} < j_{\nu,3} < \cdots\end{aligned}\end{align} \]Examples
Initial zeros of the Bessel functions \(J_0(z), J_1(z), J_2(z)\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> besseljzero(0,1); besseljzero(0,2); besseljzero(0,3) 2.404825557695772768621632 5.520078110286310649596604 8.653727912911012216954199 >>> besseljzero(1,1); besseljzero(1,2); besseljzero(1,3) 3.831705970207512315614436 7.01558666981561875353705 10.17346813506272207718571 >>> besseljzero(2,1); besseljzero(2,2); besseljzero(2,3) 5.135622301840682556301402 8.417244140399864857783614 11.61984117214905942709415
Initial zeros of \(J'_0(z), J'_1(z), J'_2(z)\):
0.0 3.831705970207512315614436 7.01558666981561875353705 >>> besseljzero(1,1,1); besseljzero(1,2,1); besseljzero(1,3,1) 1.84118378134065930264363 5.331442773525032636884016 8.536316366346285834358961 >>> besseljzero(2,1,1); besseljzero(2,2,1); besseljzero(2,3,1) 3.054236928227140322755932 6.706133194158459146634394 9.969467823087595793179143
Zeros with large index:
>>> besseljzero(0,100); besseljzero(0,1000); besseljzero(0,10000) 313.3742660775278447196902 3140.807295225078628895545 31415.14114171350798533666 >>> besseljzero(5,100); besseljzero(5,1000); besseljzero(5,10000) 321.1893195676003157339222 3148.657306813047523500494 31422.9947255486291798943 >>> besseljzero(0,100,1); besseljzero(0,1000,1); besseljzero(0,10000,1) 311.8018681873704508125112 3139.236339643802482833973 31413.57032947022399485808
Zeros of functions with large order:
>>> besseljzero(50,1) 57.11689916011917411936228 >>> besseljzero(50,2) 62.80769876483536093435393 >>> besseljzero(50,100) 388.6936600656058834640981 >>> besseljzero(50,1,1) 52.99764038731665010944037 >>> besseljzero(50,2,1) 60.02631933279942589882363 >>> besseljzero(50,100,1) 387.1083151608726181086283
Zeros of functions with fractional order:
>>> besseljzero(0.5,1); besseljzero(1.5,1); besseljzero(2.25,4) 3.141592653589793238462643 4.493409457909064175307881 15.15657692957458622921634
Both \(J_{\nu}(z)\) and \(J'_{\nu}(z)\) can be expressed as infinite products over their zeros:
>>> v,z = 2, mpf(1) >>> (z/2)**v/gamma(v+1) * \ ... nprod(lambda k: 1-(z/besseljzero(v,k))**2, [1,inf]) ... 0.1149034849319004804696469 >>> besselj(v,z) 0.1149034849319004804696469 >>> (z/2)**(v-1)/2/gamma(v) * \ ... nprod(lambda k: 1-(z/besseljzero(v,k,1))**2, [1,inf]) ... 0.2102436158811325550203884 >>> besselj(v,z,1) 0.2102436158811325550203884
besselyzero()¶
- mpmath.besselyzero(v, m, derivative=0)¶
For a real order \(\nu \ge 0\) and a positive integer \(m\), returns \(y_{\nu,m}\), the \(m\)-th positive zero of the Bessel function of the second kind \(Y_{\nu}(z)\) (see
bessely()). Alternatively, with derivative=1, gives the first positive zero \(y'_{\nu,m}\) of \(Y'_{\nu}(z)\).The zeros are interlaced according to the inequalities
\[ \begin{align}\begin{aligned}y_{\nu,k} < y'_{\nu,k} < y_{\nu,k+1}\\y_{\nu,1} < y_{\nu+1,2} < y_{\nu,2} < y_{\nu+1,2} < y_{\nu,3} < \cdots\end{aligned}\end{align} \]Examples
Initial zeros of the Bessel functions \(Y_0(z), Y_1(z), Y_2(z)\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> besselyzero(0,1); besselyzero(0,2); besselyzero(0,3) 0.8935769662791675215848871 3.957678419314857868375677 7.086051060301772697623625 >>> besselyzero(1,1); besselyzero(1,2); besselyzero(1,3) 2.197141326031017035149034 5.429681040794135132772005 8.596005868331168926429606 >>> besselyzero(2,1); besselyzero(2,2); besselyzero(2,3) 3.384241767149593472701426 6.793807513268267538291167 10.02347797936003797850539
Initial zeros of \(Y'_0(z), Y'_1(z), Y'_2(z)\):
>>> besselyzero(0,1,1); besselyzero(0,2,1); besselyzero(0,3,1) 2.197141326031017035149034 5.429681040794135132772005 8.596005868331168926429606 >>> besselyzero(1,1,1); besselyzero(1,2,1); besselyzero(1,3,1) 3.683022856585177699898967 6.941499953654175655751944 10.12340465543661307978775 >>> besselyzero(2,1,1); besselyzero(2,2,1); besselyzero(2,3,1) 5.002582931446063945200176 8.350724701413079526349714 11.57419546521764654624265
Zeros with large index:
>>> besselyzero(0,100); besselyzero(0,1000); besselyzero(0,10000) 311.8034717601871549333419 3139.236498918198006794026 31413.57034538691205229188 >>> besselyzero(5,100); besselyzero(5,1000); besselyzero(5,10000) 319.6183338562782156235062 3147.086508524556404473186 31421.42392920214673402828 >>> besselyzero(0,100,1); besselyzero(0,1000,1); besselyzero(0,10000,1) 313.3726705426359345050449 3140.807136030340213610065 31415.14112579761578220175
Zeros of functions with large order:
>>> besselyzero(50,1) 53.50285882040036394680237 >>> besselyzero(50,2) 60.11244442774058114686022 >>> besselyzero(50,100) 387.1096509824943957706835 >>> besselyzero(50,1,1) 56.96290427516751320063605 >>> besselyzero(50,2,1) 62.74888166945933944036623 >>> besselyzero(50,100,1) 388.6923300548309258355475
Zeros of functions with fractional order:
>>> besselyzero(0.5,1); besselyzero(1.5,1); besselyzero(2.25,4) 1.570796326794896619231322 2.798386045783887136720249 13.56721208770735123376018
Hankel functions¶
hankel1()¶
- mpmath.hankel1(n, x)¶
hankel1(n,x)computes the Hankel function of the first kind, which is the complex combination of Bessel functions given by\[H_n^{(1)}(x) = J_n(x) + i Y_n(x).\]Plots
# Hankel function H1_n(x) on the real line for n=0,1,2,3 h0 = lambda x: hankel1(0,x) h1 = lambda x: hankel1(1,x) h2 = lambda x: hankel1(2,x) h3 = lambda x: hankel1(3,x) plot([h0,h1,h2,h3],[0,6],[-2,1])
# Hankel function H1_n(z) in the complex plane cplot(lambda z: hankel1(1,z), [-8,8], [-8,8], points=50000)
Examples
The Hankel function is generally complex-valued:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> hankel1(2, pi) (0.4854339326315091097054957 - 0.0999007139290278787734903j) >>> hankel1(3.5, pi) (0.2340002029630507922628888 - 0.6419643823412927142424049j)
hankel2()¶
- mpmath.hankel2(n, x)¶
hankel2(n,x)computes the Hankel function of the second kind, which is the complex combination of Bessel functions given by\[H_n^{(2)}(x) = J_n(x) - i Y_n(x).\]Plots
# Hankel function H2_n(x) on the real line for n=0,1,2,3 h0 = lambda x: hankel2(0,x) h1 = lambda x: hankel2(1,x) h2 = lambda x: hankel2(2,x) h3 = lambda x: hankel2(3,x) plot([h0,h1,h2,h3],[0,6],[-1,2])
# Hankel function H2_n(z) in the complex plane cplot(lambda z: hankel2(1,z), [-8,8], [-8,8], points=50000)
Examples
The Hankel function is generally complex-valued:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> hankel2(2, pi) (0.4854339326315091097054957 + 0.0999007139290278787734903j) >>> hankel2(3.5, pi) (0.2340002029630507922628888 + 0.6419643823412927142424049j)
Kelvin functions¶
ber()¶
- mpmath.ber(n, z, **kwargs)¶
Computes the Kelvin function ber, which for real arguments gives the real part of the Bessel J function of a rotated argument
\[J_n\left(x e^{3\pi i/4}\right) = \mathrm{ber}_n(x) + i \mathrm{bei}_n(x).\]The imaginary part is given by
bei().Plots
# Kelvin functions ber_n(x) and bei_n(x) on the real line for n=0,2 f0 = lambda x: ber(0,x) f1 = lambda x: bei(0,x) f2 = lambda x: ber(2,x) f3 = lambda x: bei(2,x) plot([f0,f1,f2,f3],[0,10],[-10,10])
Examples
Verifying the defining relation:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> n, x = 2, 3.5 >>> ber(n,x) 1.442338852571888752631129 >>> bei(n,x) -0.948359035324558320217678 >>> besselj(n, x*root(1,8,3)) (1.442338852571888752631129 - 0.948359035324558320217678j)
The ber and bei functions are also defined by analytic continuation for complex arguments:
>>> ber(1+j, 2+3j) (4.675445984756614424069563 - 15.84901771719130765656316j) >>> bei(1+j, 2+3j) (15.83886679193707699364398 + 4.684053288183046528703611j)
bei()¶
ker()¶
- mpmath.ker(n, z, **kwargs)¶
Computes the Kelvin function ker, which for real arguments gives the real part of the (rescaled) Bessel K function of a rotated argument
\[e^{-\pi i/2} K_n\left(x e^{3\pi i/4}\right) = \mathrm{ker}_n(x) + i \mathrm{kei}_n(x).\]The imaginary part is given by
kei().Plots
# Kelvin functions ker_n(x) and kei_n(x) on the real line for n=0,2 f0 = lambda x: ker(0,x) f1 = lambda x: kei(0,x) f2 = lambda x: ker(2,x) f3 = lambda x: kei(2,x) plot([f0,f1,f2,f3],[0,5],[-1,4])
Examples
Verifying the defining relation:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> n, x = 2, 4.5 >>> ker(n,x) 0.02542895201906369640249801 >>> kei(n,x) -0.02074960467222823237055351 >>> exp(-n*pi*j/2) * besselk(n, x*root(1,8,1)) (0.02542895201906369640249801 - 0.02074960467222823237055351j)
The ker and kei functions are also defined by analytic continuation for complex arguments:
>>> ker(1+j, 3+4j) (1.586084268115490421090533 - 2.939717517906339193598719j) >>> kei(1+j, 3+4j) (-2.940403256319453402690132 - 1.585621643835618941044855j)
Struve functions¶
struveh()¶
- mpmath.struveh(n, z, **kwargs)¶
Gives the Struve function
\[\,\mathbf{H}_n(z) = \sum_{k=0}^\infty \frac{(-1)^k}{\Gamma(k+\frac{3}{2}) \Gamma(k+n+\frac{3}{2})} {\left({\frac{z}{2}}\right)}^{2k+n+1}\]which is a solution to the Struve differential equation
\[z^2 f''(z) + z f'(z) + (z^2-n^2) f(z) = \frac{2 z^{n+1}}{\pi (2n-1)!!}.\]Examples
Evaluation for arbitrary real and complex arguments:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> struveh(0, 3.5) 0.3608207733778295024977797 >>> struveh(-1, 10) -0.255212719726956768034732 >>> struveh(1, -100.5) 0.5819566816797362287502246 >>> struveh(2.5, 10000000000000) 3153915652525200060.308937 >>> struveh(2.5, -10000000000000) (0.0 - 3153915652525200060.308937j) >>> struveh(1+j, 1000000+4000000j) (-3.066421087689197632388731e+1737173 - 1.596619701076529803290973e+1737173j)
A Struve function of half-integer order is elementary; for example:
>>> z = 3 >>> struveh(0.5, 3) 0.9167076867564138178671595 >>> sqrt(2/(pi*z))*(1-cos(z)) 0.9167076867564138178671595
Numerically verifying the differential equation:
>>> z = mpf(4.5) >>> n = 3 >>> f = lambda z: struveh(n,z) >>> lhs = z**2*diff(f,z,2) + z*diff(f,z) + (z**2-n**2)*f(z) >>> rhs = 2*z**(n+1)/fac2(2*n-1)/pi >>> lhs 17.40359302709875496632744 >>> rhs 17.40359302709875496632744
struvel()¶
- mpmath.struvel(n, z, **kwargs)¶
Gives the modified Struve function
\[\,\mathbf{L}_n(z) = -i e^{-n\pi i/2} \mathbf{H}_n(i z)\]which solves to the modified Struve differential equation
\[z^2 f''(z) + z f'(z) - (z^2+n^2) f(z) = \frac{2 z^{n+1}}{\pi (2n-1)!!}.\]Examples
Evaluation for arbitrary real and complex arguments:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> struvel(0, 3.5) 7.180846515103737996249972 >>> struvel(-1, 10) 2670.994904980850550721511 >>> struvel(1, -100.5) 1.757089288053346261497686e+42 >>> struvel(2.5, 10000000000000) 4.160893281017115450519948e+4342944819025 >>> struvel(2.5, -10000000000000) (0.0 - 4.160893281017115450519948e+4342944819025j) >>> struvel(1+j, 700j) (-0.1721150049480079451246076 + 0.1240770953126831093464055j) >>> struvel(1+j, 1000000+4000000j) (-2.973341637511505389128708e+434290 - 5.164633059729968297147448e+434290j)
Numerically verifying the differential equation:
>>> z = mpf(3.5) >>> n = 3 >>> f = lambda z: struvel(n,z) >>> lhs = z**2*diff(f,z,2) + z*diff(f,z) - (z**2+n**2)*f(z) >>> rhs = 2*z**(n+1)/fac2(2*n-1)/pi >>> lhs 6.368850306060678353018165 >>> rhs 6.368850306060678353018165
Anger-Weber functions¶
angerj()¶
- mpmath.angerj(v, z, **kwargs)¶
Gives the Anger function
\[\mathbf{J}_{\nu}(z) = \frac{1}{\pi} \int_0^{\pi} \cos(\nu t - z \sin t) dt\]which is an entire function of both the parameter \(\nu\) and the argument \(z\). It solves the inhomogeneous Bessel differential equation
\[f''(z) + \frac{1}{z}f'(z) + \left(1-\frac{\nu^2}{z^2}\right) f(z) = \frac{(z-\nu)}{\pi z^2} \sin(\pi \nu).\]Examples
Evaluation for real and complex parameter and argument:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> angerj(2,3) 0.4860912605858910769078311 >>> angerj(-3+4j, 2+5j) (-5033.358320403384472395612 + 585.8011892476145118551756j) >>> angerj(3.25, 1e6j) (4.630743639715893346570743e+434290 - 1.117960409887505906848456e+434291j) >>> angerj(-1.5, 1e6) 0.0002795719747073879393087011
The Anger function coincides with the Bessel J-function when \(\nu\) is an integer:
>>> angerj(1,3); besselj(1,3) 0.3390589585259364589255146 0.3390589585259364589255146 >>> angerj(1.5,3); besselj(1.5,3) 0.4088969848691080859328847 0.4777182150870917715515015
Verifying the differential equation:
>>> v,z = mpf(2.25), 0.75 >>> f = lambda z: angerj(v,z) >>> diff(f,z,2) + diff(f,z)/z + (1-(v/z)**2)*f(z) -0.6002108774380707130367995 >>> (z-v)/(pi*z**2) * sinpi(v) -0.6002108774380707130367995
Verifying the integral representation:
>>> angerj(v,z) 0.1145380759919333180900501 >>> quad(lambda t: cos(v*t-z*sin(t))/pi, [0,pi]) 0.1145380759919333180900501
References
[DLMF] section 11.10: Anger-Weber Functions
webere()¶
- mpmath.webere(v, z, **kwargs)¶
Gives the Weber function
\[\mathbf{E}_{\nu}(z) = \frac{1}{\pi} \int_0^{\pi} \sin(\nu t - z \sin t) dt\]which is an entire function of both the parameter \(\nu\) and the argument \(z\). It solves the inhomogeneous Bessel differential equation
\[f''(z) + \frac{1}{z}f'(z) + \left(1-\frac{\nu^2}{z^2}\right) f(z) = -\frac{1}{\pi z^2} (z+\nu+(z-\nu)\cos(\pi \nu)).\]Examples
Evaluation for real and complex parameter and argument:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> webere(2,3) -0.1057668973099018425662646 >>> webere(-3+4j, 2+5j) (-585.8081418209852019290498 - 5033.314488899926921597203j) >>> webere(3.25, 1e6j) (-1.117960409887505906848456e+434291 - 4.630743639715893346570743e+434290j) >>> webere(3.25, 1e6) -0.00002812518265894315604914453
Up to addition of a rational function of \(z\), the Weber function coincides with the Struve H-function when \(\nu\) is an integer:
>>> webere(1,3); 2/pi-struveh(1,3) -0.3834897968188690177372881 -0.3834897968188690177372881 >>> webere(5,3); 26/(35*pi)-struveh(5,3) 0.2009680659308154011878075 0.2009680659308154011878075
Verifying the differential equation:
>>> v,z = mpf(2.25), 0.75 >>> f = lambda z: webere(v,z) >>> diff(f,z,2) + diff(f,z)/z + (1-(v/z)**2)*f(z) -1.097441848875479535164627 >>> -(z+v+(z-v)*cospi(v))/(pi*z**2) -1.097441848875479535164627
Verifying the integral representation:
>>> webere(v,z) 0.1486507351534283744485421 >>> quad(lambda t: sin(v*t-z*sin(t))/pi, [0,pi]) 0.1486507351534283744485421
References
[DLMF] section 11.10: Anger-Weber Functions
Lommel functions¶
lommels1()¶
- mpmath.lommels1(u, v, z, **kwargs)¶
Gives the Lommel function \(s_{\mu,\nu}\) or \(s^{(1)}_{\mu,\nu}\)
\[s_{\mu,\nu}(z) = \frac{z^{\mu+1}}{(\mu-\nu+1)(\mu+\nu+1)} \,_1F_2\left(1; \frac{\mu-\nu+3}{2}, \frac{\mu+\nu+3}{2}; -\frac{z^2}{4} \right)\]which solves the inhomogeneous Bessel equation
\[z^2 f''(z) + z f'(z) + (z^2-\nu^2) f(z) = z^{\mu+1}.\]A second solution is given by
lommels2().Plots
# Lommel function s_(u,v)(x) on the real line for a few different u,v f1 = lambda x: lommels1(-1,2.5,x) f2 = lambda x: lommels1(0,0.5,x) f3 = lambda x: lommels1(0,6,x) f4 = lambda x: lommels1(0.5,3,x) plot([f1,f2,f3,f4], [0,20])
Examples
An integral representation:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> u,v,z = 0.25, 0.125, mpf(0.75) >>> lommels1(u,v,z) 0.4276243877565150372999126 >>> (bessely(v,z)*quad(lambda t: t**u*besselj(v,t), [0,z]) - \ ... besselj(v,z)*quad(lambda t: t**u*bessely(v,t), [0,z]))*(pi/2) 0.4276243877565150372999126
A special value:
>>> lommels1(v,v,z) 0.5461221367746048054932553 >>> gamma(v+0.5)*sqrt(pi)*power(2,v-1)*struveh(v,z) 0.5461221367746048054932553
Verifying the differential equation:
>>> f = lambda z: lommels1(u,v,z) >>> z**2*diff(f,z,2) + z*diff(f,z) + (z**2-v**2)*f(z) 0.6979536443265746992059141 >>> z**(u+1) 0.6979536443265746992059141
References
lommels2()¶
- mpmath.lommels2(u, v, z, **kwargs)¶
Gives the second Lommel function \(S_{\mu,\nu}\) or \(s^{(2)}_{\mu,\nu}\)
\[ \begin{align}\begin{aligned}S_{\mu,\nu}(z) = s_{\mu,\nu}(z) + 2^{\mu-1} \Gamma\left(\tfrac{1}{2}(\mu-\nu+1)\right) \Gamma\left(\tfrac{1}{2}(\mu+\nu+1)\right) \times\\ \left[\sin(\tfrac{1}{2}(\mu-\nu)\pi) J_{\nu}(z) - \cos(\tfrac{1}{2}(\mu-\nu)\pi) Y_{\nu}(z) \right]\end{aligned}\end{align} \]which solves the same differential equation as
lommels1().Plots
# Lommel function S_(u,v)(x) on the real line for a few different u,v f1 = lambda x: lommels2(-1,2.5,x) f2 = lambda x: lommels2(1.5,2,x) f3 = lambda x: lommels2(2.5,1,x) f4 = lambda x: lommels2(3.5,-0.5,x) plot([f1,f2,f3,f4], [0,8], [-8,8])
Examples
For large \(|z|\), \(S_{\mu,\nu} \sim z^{\mu-1}\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> lommels2(10,2,30000) 1.968299831601008419949804e+40 >>> power(30000,9) 1.9683e+40
A special value:
>>> u,v,z = 0.5, 0.125, mpf(0.75) >>> lommels2(v,v,z) 0.9589683199624672099969765 >>> (struveh(v,z)-bessely(v,z))*power(2,v-1)*sqrt(pi)*gamma(v+0.5) 0.9589683199624672099969765
Verifying the differential equation:
>>> f = lambda z: lommels2(u,v,z) >>> z**2*diff(f,z,2) + z*diff(f,z) + (z**2-v**2)*f(z) 0.6495190528383289850727924 >>> z**(u+1) 0.6495190528383289850727924
References
Airy and Scorer functions¶
airyai()¶
- mpmath.airyai(z, derivative=0, **kwargs)¶
Computes the Airy function \(\operatorname{Ai}(z)\), which is the solution of the Airy differential equation \(f''(z) - z f(z) = 0\) with initial conditions
\[ \begin{align}\begin{aligned}\operatorname{Ai}(0) = \frac{1}{3^{2/3}\Gamma\left(\frac{2}{3}\right)}\\\operatorname{Ai}'(0) = -\frac{1}{3^{1/3}\Gamma\left(\frac{1}{3}\right)}.\end{aligned}\end{align} \]Other common ways of defining the Ai-function include integrals such as
\[ \begin{align}\begin{aligned}\operatorname{Ai}(x) = \frac{1}{\pi} \int_0^{\infty} \cos\left(\frac{1}{3}t^3+xt\right) dt \qquad x \in \mathbb{R}\\\operatorname{Ai}(z) = \frac{\sqrt{3}}{2\pi} \int_0^{\infty} \exp\left(-\frac{t^3}{3}-\frac{z^3}{3t^3}\right) dt.\end{aligned}\end{align} \]The Ai-function is an entire function with a turning point, behaving roughly like a slowly decaying sine wave for \(z < 0\) and like a rapidly decreasing exponential for \(z > 0\). A second solution of the Airy differential equation is given by \(\operatorname{Bi}(z)\) (see
airybi()).Optionally, with derivative=alpha,
airyai()can compute the \(\alpha\)-th order fractional derivative with respect to \(z\). For \(\alpha = n = 1,2,3,\ldots\) this gives the derivative \(\operatorname{Ai}^{(n)}(z)\), and for \(\alpha = -n = -1,-2,-3,\ldots\) this gives the \(n\)-fold iterated integral\[ \begin{align}\begin{aligned}f_0(z) = \operatorname{Ai}(z)\\f_n(z) = \int_0^z f_{n-1}(t) dt.\end{aligned}\end{align} \]The Ai-function has infinitely many zeros, all located along the negative half of the real axis. They can be computed with
airyaizero().Plots
# Airy function Ai(x), Ai'(x) and int_0^x Ai(t) dt on the real line f = airyai f_diff = lambda z: airyai(z, derivative=1) f_int = lambda z: airyai(z, derivative=-1) plot([f, f_diff, f_int], [-10,5])
# Airy function Ai(z) in the complex plane cplot(airyai, [-8,8], [-8,8], points=50000)
Basic examples
Limits and values include:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> airyai(0); 1/(power(3,'2/3')*gamma('2/3')) 0.3550280538878172392600632 0.3550280538878172392600632 >>> airyai(1) 0.1352924163128814155241474 >>> airyai(-1) 0.5355608832923521187995166 >>> airyai(inf); airyai(-inf) 0.0 0.0
Evaluation is supported for large magnitudes of the argument:
>>> airyai(-100) 0.1767533932395528780908311 >>> airyai(100) 2.634482152088184489550553e-291 >>> airyai(50+50j) (-5.31790195707456404099817e-68 - 1.163588003770709748720107e-67j) >>> airyai(-50+50j) (1.041242537363167632587245e+158 + 3.347525544923600321838281e+157j)
Huge arguments are also fine:
>>> airyai(10**10) 1.162235978298741779953693e-289529654602171 >>> airyai(-10**10) 0.0001736206448152818510510181 >>> w = airyai(10**10*(1+j)) >>> w.real 5.711508683721355528322567e-186339621747698 >>> w.imag 1.867245506962312577848166e-186339621747697
The first root of the Ai-function is:
>>> findroot(airyai, -2) -2.338107410459767038489197 >>> airyaizero(1) -2.338107410459767038489197
Properties and relations
Verifying the Airy differential equation:
>>> for z in [-3.4, 0, 2.5, 1+2j]: ... chop(airyai(z,2) - z*airyai(z)) ... 0.0 0.0 0.0 0.0
The first few terms of the Taylor series expansion around \(z = 0\) (every third term is zero):
>>> nprint(taylor(airyai, 0, 5)) [0.355028, -0.258819, 0.0, 0.0591713, -0.0215683, 0.0]
The Airy functions satisfy the Wronskian relation \(\operatorname{Ai}(z) \operatorname{Bi}'(z) - \operatorname{Ai}'(z) \operatorname{Bi}(z) = 1/\pi\):
>>> z = -0.5 >>> airyai(z)*airybi(z,1) - airyai(z,1)*airybi(z) 0.3183098861837906715377675 >>> 1/pi 0.3183098861837906715377675
The Airy functions can be expressed in terms of Bessel functions of order \(\pm 1/3\). For \(\Re[z] \le 0\), we have:
>>> z = -3 >>> airyai(z) -0.3788142936776580743472439 >>> y = 2*power(-z,'3/2')/3 >>> (sqrt(-z) * (besselj('1/3',y) + besselj('-1/3',y)))/3 -0.3788142936776580743472439
Derivatives and integrals
Derivatives of the Ai-function (directly and using
diff()):>>> airyai(-3,1); diff(airyai,-3) 0.3145837692165988136507873 0.3145837692165988136507873 >>> airyai(-3,2); diff(airyai,-3,2) 1.136442881032974223041732 1.136442881032974223041732 >>> airyai(1000,1); diff(airyai,1000) -2.943133917910336090459748e-9156 -2.943133917910336090459748e-9156
Several derivatives at \(z = 0\):
>>> airyai(0,0); airyai(0,1); airyai(0,2) 0.3550280538878172392600632 -0.2588194037928067984051836 0.0 >>> airyai(0,3); airyai(0,4); airyai(0,5) 0.3550280538878172392600632 -0.5176388075856135968103671 0.0 >>> airyai(0,15); airyai(0,16); airyai(0,17) 1292.30211615165475090663 -3188.655054727379756351861 0.0
The integral of the Ai-function:
>>> airyai(3,-1); quad(airyai, [0,3]) 0.3299203760070217725002701 0.3299203760070217725002701 >>> airyai(-10,-1); quad(airyai, [0,-10]) -0.765698403134212917425148 -0.765698403134212917425148
Integrals of high or fractional order:
>>> airyai(-2,0.5); differint(airyai,-2,0.5,0) (0.0 + 0.2453596101351438273844725j) (0.0 + 0.2453596101351438273844725j) >>> airyai(-2,-4); differint(airyai,-2,-4,0) 0.2939176441636809580339365 0.2939176441636809580339365 >>> airyai(0,-1); airyai(0,-2); airyai(0,-3) 0.0 0.0 0.0
Integrals of the Ai-function can be evaluated at limit points:
>>> airyai(-1000000,-1); airyai(-inf,-1) -0.6666843728311539978751512 -0.6666666666666666666666667 >>> airyai(10,-1); airyai(+inf,-1) 0.3333333332991690159427932 0.3333333333333333333333333 >>> airyai(+inf,-2); airyai(+inf,-3) +inf +inf >>> airyai(-1000000,-2); airyai(-inf,-2) 666666.4078472650651209742 +inf >>> airyai(-1000000,-3); airyai(-inf,-3) -333333074513.7520264995733 -inf
References
[DLMF] Chapter 9: Airy and Related Functions
[WolframFunctions] section: Bessel-Type Functions
airybi()¶
- mpmath.airybi(z, derivative=0, **kwargs)¶
Computes the Airy function \(\operatorname{Bi}(z)\), which is the solution of the Airy differential equation \(f''(z) - z f(z) = 0\) with initial conditions
\[ \begin{align}\begin{aligned}\operatorname{Bi}(0) = \frac{1}{3^{1/6}\Gamma\left(\frac{2}{3}\right)}\\\operatorname{Bi}'(0) = \frac{3^{1/6}}{\Gamma\left(\frac{1}{3}\right)}.\end{aligned}\end{align} \]Like the Ai-function (see
airyai()), the Bi-function is oscillatory for \(z < 0\), but it grows rather than decreases for \(z > 0\).Optionally, as for
airyai(), derivatives, integrals and fractional derivatives can be computed with the derivative parameter.The Bi-function has infinitely many zeros along the negative half-axis, as well as complex zeros, which can all be computed with
airybizero().Plots
# Airy function Bi(x), Bi'(x) and int_0^x Bi(t) dt on the real line f = airybi f_diff = lambda z: airybi(z, derivative=1) f_int = lambda z: airybi(z, derivative=-1) plot([f, f_diff, f_int], [-10,2], [-1,2])
# Airy function Bi(z) in the complex plane cplot(airybi, [-8,8], [-8,8], points=50000)
Basic examples
Limits and values include:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> airybi(0); 1/(power(3,'1/6')*gamma('2/3')) 0.6149266274460007351509224 0.6149266274460007351509224 >>> airybi(1) 1.207423594952871259436379 >>> airybi(-1) 0.10399738949694461188869 >>> airybi(inf); airybi(-inf) +inf 0.0
Evaluation is supported for large magnitudes of the argument:
>>> airybi(-100) 0.02427388768016013160566747 >>> airybi(100) 6.041223996670201399005265e+288 >>> airybi(50+50j) (-5.322076267321435669290334e+63 + 1.478450291165243789749427e+65j) >>> airybi(-50+50j) (-3.347525544923600321838281e+157 + 1.041242537363167632587245e+158j)
Huge arguments:
>>> airybi(10**10) 1.369385787943539818688433e+289529654602165 >>> airybi(-10**10) 0.001775656141692932747610973 >>> w = airybi(10**10*(1+j)) >>> w.real -6.559955931096196875845858e+186339621747689 >>> w.imag -6.822462726981357180929024e+186339621747690
The first real root of the Bi-function is:
>>> findroot(airybi, -1); airybizero(1) -1.17371322270912792491998 -1.17371322270912792491998
Properties and relations
Verifying the Airy differential equation:
>>> for z in [-3.4, 0, 2.5, 1+2j]: ... chop(airybi(z,2) - z*airybi(z)) ... 0.0 0.0 0.0 0.0
The first few terms of the Taylor series expansion around \(z = 0\) (every third term is zero):
>>> nprint(taylor(airybi, 0, 5)) [0.614927, 0.448288, 0.0, 0.102488, 0.0373574, 0.0]
The Airy functions can be expressed in terms of Bessel functions of order \(\pm 1/3\). For \(\Re[z] \le 0\), we have:
>>> z = -3 >>> airybi(z) -0.1982896263749265432206449 >>> p = 2*power(-z,'3/2')/3 >>> sqrt(-mpf(z)/3)*(besselj('-1/3',p) - besselj('1/3',p)) -0.1982896263749265432206449
Derivatives and integrals
Derivatives of the Bi-function (directly and using
diff()):>>> airybi(-3,1); diff(airybi,-3) -0.675611222685258537668032 -0.675611222685258537668032 >>> airybi(-3,2); diff(airybi,-3,2) 0.5948688791247796296619346 0.5948688791247796296619346 >>> airybi(1000,1); diff(airybi,1000) 1.710055114624614989262335e+9156 1.710055114624614989262335e+9156
Several derivatives at \(z = 0\):
>>> airybi(0,0); airybi(0,1); airybi(0,2) 0.6149266274460007351509224 0.4482883573538263579148237 0.0 >>> airybi(0,3); airybi(0,4); airybi(0,5) 0.6149266274460007351509224 0.8965767147076527158296474 0.0 >>> airybi(0,15); airybi(0,16); airybi(0,17) 2238.332923903442675949357 5522.912562599140729510628 0.0
The integral of the Bi-function:
>>> airybi(3,-1); quad(airybi, [0,3]) 10.06200303130620056316655 10.06200303130620056316655 >>> airybi(-10,-1); quad(airybi, [0,-10]) -0.01504042480614002045135483 -0.01504042480614002045135483
Integrals of high or fractional order:
>>> airybi(-2,0.5); differint(airybi, -2, 0.5, 0) (0.0 + 0.5019859055341699223453257j) (0.0 + 0.5019859055341699223453257j) >>> airybi(-2,-4); differint(airybi,-2,-4,0) 0.2809314599922447252139092 0.2809314599922447252139092 >>> airybi(0,-1); airybi(0,-2); airybi(0,-3) 0.0 0.0 0.0
Integrals of the Bi-function can be evaluated at limit points:
>>> airybi(-1000000,-1); airybi(-inf,-1) 0.000002191261128063434047966873 0.0 >>> airybi(10,-1); airybi(+inf,-1) 147809803.1074067161675853 +inf >>> airybi(+inf,-2); airybi(+inf,-3) +inf +inf >>> airybi(-1000000,-2); airybi(-inf,-2) 0.4482883750599908479851085 0.4482883573538263579148237 >>> gamma('2/3')*power(3,'2/3')/(2*pi) 0.4482883573538263579148237 >>> airybi(-100000,-3); airybi(-inf,-3) -44828.52827206932872493133 -inf >>> airybi(-100000,-4); airybi(-inf,-4) 2241411040.437759489540248 +inf
airyaizero()¶
- mpmath.airyaizero(k, derivative=0)¶
Gives the \(k\)-th zero of the Airy Ai-function, i.e. the \(k\)-th number \(a_k\) ordered by magnitude for which \(\operatorname{Ai}(a_k) = 0\).
Optionally, with derivative=1, the corresponding zero \(a'_k\) of the derivative function, i.e. \(\operatorname{Ai}'(a'_k) = 0\), is computed.
Examples
Some values of \(a_k\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> airyaizero(1) -2.338107410459767038489197 >>> airyaizero(2) -4.087949444130970616636989 >>> airyaizero(3) -5.520559828095551059129856 >>> airyaizero(1000) -281.0315196125215528353364
Some values of \(a'_k\):
>>> airyaizero(1,1) -1.018792971647471089017325 >>> airyaizero(2,1) -3.248197582179836537875424 >>> airyaizero(3,1) -4.820099211178735639400616 >>> airyaizero(1000,1) -280.9378080358935070607097
Verification:
>>> chop(airyai(airyaizero(1))) 0.0 >>> chop(airyai(airyaizero(1,1),1)) 0.0
airybizero()¶
- mpmath.airybizero(k, derivative=0, complex=0)¶
With complex=False, gives the \(k\)-th real zero of the Airy Bi-function, i.e. the \(k\)-th number \(b_k\) ordered by magnitude for which \(\operatorname{Bi}(b_k) = 0\).
With complex=True, gives the \(k\)-th complex zero in the upper half plane \(\beta_k\). Also the conjugate \(\overline{\beta_k}\) is a zero.
Optionally, with derivative=1, the corresponding zero \(b'_k\) or \(\beta'_k\) of the derivative function, i.e. \(\operatorname{Bi}'(b'_k) = 0\) or \(\operatorname{Bi}'(\beta'_k) = 0\), is computed.
Examples
Some values of \(b_k\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> airybizero(1) -1.17371322270912792491998 >>> airybizero(2) -3.271093302836352715680228 >>> airybizero(3) -4.830737841662015932667709 >>> airybizero(1000) -280.9378112034152401578834
Some values of \(b_k\):
>>> airybizero(1,1) -2.294439682614123246622459 >>> airybizero(2,1) -4.073155089071828215552369 >>> airybizero(3,1) -5.512395729663599496259593 >>> airybizero(1000,1) -281.0315164471118527161362
Some values of \(\beta_k\):
>>> airybizero(1,complex=True) (0.9775448867316206859469927 + 2.141290706038744575749139j) >>> airybizero(2,complex=True) (1.896775013895336346627217 + 3.627291764358919410440499j) >>> airybizero(3,complex=True) (2.633157739354946595708019 + 4.855468179979844983174628j) >>> airybizero(1000,complex=True) (140.4978560578493018899793 + 243.3907724215792121244867j)
Some values of \(\beta'_k\):
>>> airybizero(1,1,complex=True) (0.2149470745374305676088329 + 1.100600143302797880647194j) >>> airybizero(2,1,complex=True) (1.458168309223507392028211 + 2.912249367458445419235083j) >>> airybizero(3,1,complex=True) (2.273760763013482299792362 + 4.254528549217097862167015j) >>> airybizero(1000,1,complex=True) (140.4509972835270559730423 + 243.3096175398562811896208j)
Verification:
>>> chop(airybi(airybizero(1))) 0.0 >>> chop(airybi(airybizero(1,1),1)) 0.0 >>> u = airybizero(1,complex=True) >>> chop(airybi(u)) 0.0 >>> chop(airybi(conj(u))) 0.0
The complex zeros (in the upper and lower half-planes respectively) asymptotically approach the rays \(z = R \exp(\pm i \pi /3)\):
>>> arg(airybizero(1,complex=True)) 1.142532510286334022305364 >>> arg(airybizero(1000,complex=True)) 1.047271114786212061583917 >>> arg(airybizero(1000000,complex=True)) 1.047197624741816183341355 >>> pi/3 1.047197551196597746154214
scorergi()¶
- mpmath.scorergi(z, **kwargs)¶
Evaluates the Scorer function
\[\operatorname{Gi}(z) = \operatorname{Ai}(z) \int_0^z \operatorname{Bi}(t) dt + \operatorname{Bi}(z) \int_z^{\infty} \operatorname{Ai}(t) dt\]which gives a particular solution to the inhomogeneous Airy differential equation \(f''(z) - z f(z) = 1/\pi\). Another particular solution is given by the Scorer Hi-function (
scorerhi()). The two functions are related as \(\operatorname{Gi}(z) + \operatorname{Hi}(z) = \operatorname{Bi}(z)\).Plots
# Scorer function Gi(x) and Gi'(x) on the real line plot([scorergi, diffun(scorergi)], [-10,10])
# Scorer function Gi(z) in the complex plane cplot(scorergi, [-8,8], [-8,8], points=50000)
Examples
Some values and limits:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> scorergi(0); 1/(power(3,'7/6')*gamma('2/3')) 0.2049755424820002450503075 0.2049755424820002450503075 >>> diff(scorergi, 0); 1/(power(3,'5/6')*gamma('1/3')) 0.1494294524512754526382746 0.1494294524512754526382746 >>> scorergi(+inf); scorergi(-inf) 0.0 0.0 >>> scorergi(1) 0.2352184398104379375986902 >>> scorergi(-1) -0.1166722172960152826494198
Evaluation for large arguments:
>>> scorergi(10) 0.03189600510067958798062034 >>> scorergi(100) 0.003183105228162961476590531 >>> scorergi(1000000) 0.0000003183098861837906721743873 >>> 1/(pi*1000000) 0.0000003183098861837906715377675 >>> scorergi(-1000) -0.08358288400262780392338014 >>> scorergi(-100000) 0.02886866118619660226809581 >>> scorergi(50+10j) (0.0061214102799778578790984 - 0.001224335676457532180747917j) >>> scorergi(-50-10j) (5.236047850352252236372551e+29 - 3.08254224233701381482228e+29j) >>> scorergi(100000j) (-8.806659285336231052679025e+6474077 + 8.684731303500835514850962e+6474077j)
Verifying the connection between Gi and Hi:
>>> z = 0.25 >>> scorergi(z) + scorerhi(z) 0.7287469039362150078694543 >>> airybi(z) 0.7287469039362150078694543
Verifying the differential equation:
>>> for z in [-3.4, 0, 2.5, 1+2j]: ... chop(diff(scorergi,z,2) - z*scorergi(z)) ... -0.3183098861837906715377675 -0.3183098861837906715377675 -0.3183098861837906715377675 -0.3183098861837906715377675
Verifying the integral representation:
>>> z = 0.5 >>> scorergi(z) 0.2447210432765581976910539 >>> Ai,Bi = airyai,airybi >>> Bi(z)*(Ai(inf,-1)-Ai(z,-1)) + Ai(z)*(Bi(z,-1)-Bi(0,-1)) 0.2447210432765581976910539
References
[DLMF] section 9.12: Scorer Functions
scorerhi()¶
- mpmath.scorerhi(z, **kwargs)¶
Evaluates the second Scorer function
\[\operatorname{Hi}(z) = \operatorname{Bi}(z) \int_{-\infty}^z \operatorname{Ai}(t) dt - \operatorname{Ai}(z) \int_{-\infty}^z \operatorname{Bi}(t) dt\]which gives a particular solution to the inhomogeneous Airy differential equation \(f''(z) - z f(z) = 1/\pi\). See also
scorergi().Plots
# Scorer function Hi(x) and Hi'(x) on the real line plot([scorerhi, diffun(scorerhi)], [-10,2], [0,2])
# Scorer function Hi(z) in the complex plane cplot(scorerhi, [-8,8], [-8,8], points=50000)
Examples
Some values and limits:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> scorerhi(0); 2/(power(3,'7/6')*gamma('2/3')) 0.4099510849640004901006149 0.4099510849640004901006149 >>> diff(scorerhi,0); 2/(power(3,'5/6')*gamma('1/3')) 0.2988589049025509052765491 0.2988589049025509052765491 >>> scorerhi(+inf); scorerhi(-inf) +inf 0.0 >>> scorerhi(1) 0.9722051551424333218376886 >>> scorerhi(-1) 0.2206696067929598945381098
Evaluation for large arguments:
>>> scorerhi(10) 455641153.5163291358991077 >>> scorerhi(100) 6.041223996670201399005265e+288 >>> scorerhi(1000000) 7.138269638197858094311122e+289529652 >>> scorerhi(-10) 0.0317685352825022727415011 >>> scorerhi(-100) 0.003183092495767499864680483 >>> scorerhi(100j) (-6.366197716545672122983857e-9 + 0.003183098861710582761688475j) >>> scorerhi(50+50j) (-5.322076267321435669290334e+63 + 1.478450291165243789749427e+65j) >>> scorerhi(-1000-1000j) (0.0001591549432510502796565538 - 0.000159154943091895334973109j)
Verifying the differential equation:
>>> for z in [-3.4, 0, 2, 1+2j]: ... chop(diff(scorerhi,z,2) - z*scorerhi(z)) ... 0.3183098861837906715377675 0.3183098861837906715377675 0.3183098861837906715377675 0.3183098861837906715377675
Verifying the integral representation:
>>> z = 0.5 >>> scorerhi(z) 0.6095559998265972956089949 >>> Ai,Bi = airyai,airybi >>> Bi(z)*(Ai(z,-1)-Ai(-inf,-1)) - Ai(z)*(Bi(z,-1)-Bi(-inf,-1)) 0.6095559998265972956089949
Coulomb wave functions¶
coulombf()¶
- mpmath.coulombf(l, eta, z)¶
Calculates the regular Coulomb wave function
\[F_l(\eta,z) = C_l(\eta) z^{l+1} e^{-iz} \,_1F_1(l+1-i\eta, 2l+2, 2iz)\]where the normalization constant \(C_l(\eta)\) is as calculated by
coulombc(). This function solves the differential equation\[f''(z) + \left(1-\frac{2\eta}{z}-\frac{l(l+1)}{z^2}\right) f(z) = 0.\]A second linearly independent solution is given by the irregular Coulomb wave function \(G_l(\eta,z)\) (see
coulombg()) and thus the general solution is \(f(z) = C_1 F_l(\eta,z) + C_2 G_l(\eta,z)\) for arbitrary constants \(C_1\), \(C_2\). Physically, the Coulomb wave functions give the radial solution to the Schrodinger equation for a point particle in a \(1/z\) potential; \(z\) is then the radius and \(l\), \(\eta\) are quantum numbers.The Coulomb wave functions with real parameters are defined in Abramowitz & Stegun, section 14. However, all parameters are permitted to be complex in this implementation (see references).
Plots
# Regular Coulomb wave functions -- equivalent to figure 14.3 in A&S F1 = lambda x: coulombf(0,0,x) F2 = lambda x: coulombf(0,1,x) F3 = lambda x: coulombf(0,5,x) F4 = lambda x: coulombf(0,10,x) F5 = lambda x: coulombf(0,x/2,x) plot([F1,F2,F3,F4,F5], [0,25], [-1.2,1.6])
# Regular Coulomb wave function in the complex plane cplot(lambda z: coulombf(1,1,z), points=50000)
Examples
Evaluation is supported for arbitrary magnitudes of \(z\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> coulombf(2, 1.5, 3.5) 0.4080998961088761187426445 >>> coulombf(-2, 1.5, 3.5) 0.7103040849492536747533465 >>> coulombf(2, 1.5, '1e-10') 4.143324917492256448770769e-33 >>> coulombf(2, 1.5, 1000) 0.4482623140325567050716179 >>> coulombf(2, 1.5, 10**10) -0.066804196437694360046619
Verifying the differential equation:
>>> l, eta, z = 2, 3, mpf(2.75) >>> A, B = 1, 2 >>> f = lambda z: A*coulombf(l,eta,z) + B*coulombg(l,eta,z) >>> chop(diff(f,z,2) + (1-2*eta/z - l*(l+1)/z**2)*f(z)) 0.0
A Wronskian relation satisfied by the Coulomb wave functions:
>>> l = 2 >>> eta = 1.5 >>> F = lambda z: coulombf(l,eta,z) >>> G = lambda z: coulombg(l,eta,z) >>> for z in [3.5, -1, 2+3j]: ... chop(diff(F,z)*G(z) - F(z)*diff(G,z)) ... 1.0 1.0 1.0
Another Wronskian relation:
>>> F = coulombf >>> G = coulombg >>> for z in [3.5, -1, 2+3j]: ... chop(F(l-1,eta,z)*G(l,eta,z)-F(l,eta,z)*G(l-1,eta,z) - l/sqrt(l**2+eta**2)) ... 0.0 0.0 0.0
An integral identity connecting the regular and irregular wave functions:
>>> l, eta, z = 4+j, 2-j, 5+2j >>> coulombf(l,eta,z) + j*coulombg(l,eta,z) (0.7997977752284033239714479 + 0.9294486669502295512503127j) >>> g = lambda t: exp(-t)*t**(l-j*eta)*(t+2*j*z)**(l+j*eta) >>> j*exp(-j*z)*z**(-l)/fac(2*l+1)/coulombc(l,eta)*quad(g, [0,inf]) (0.7997977752284033239714479 + 0.9294486669502295512503127j)
Some test case with complex parameters, taken from Michel [2]:
>>> mp.dps = 15 >>> coulombf(1+0.1j, 50+50j, 100.156) (-1.02107292320897e+15 - 2.83675545731519e+15j) >>> coulombg(1+0.1j, 50+50j, 100.156) (2.83675545731519e+15 - 1.02107292320897e+15j) >>> coulombf(1e-5j, 10+1e-5j, 0.1+1e-6j) (4.30566371247811e-14 - 9.03347835361657e-19j) >>> coulombg(1e-5j, 10+1e-5j, 0.1+1e-6j) (778709182061.134 + 18418936.2660553j)
The following reproduces a table in Abramowitz & Stegun, at twice the precision:
>>> mp.dps = 10 >>> eta = 2; z = 5 >>> for l in [5, 4, 3, 2, 1, 0]: ... print("%s %s %s" % (l, coulombf(l,eta,z), ... diff(lambda z: coulombf(l,eta,z), z))) ... 5 0.09079533488 0.1042553261 4 0.2148205331 0.2029591779 3 0.4313159311 0.320534053 2 0.7212774133 0.3952408216 1 0.9935056752 0.3708676452 0 1.143337392 0.2937960375
References
I.J. Thompson & A.R. Barnett, “Coulomb and Bessel Functions of Complex Arguments and Order”, J. Comp. Phys., vol 64, no. 2, June 1986.
N. Michel, “Precise Coulomb wave functions for a wide range of complex \(l\), \(\eta\) and \(z\)”, http://arxiv.org/abs/physics/0702051v1
coulombg()¶
- mpmath.coulombg(l, eta, z)¶
Calculates the irregular Coulomb wave function
\[G_l(\eta,z) = \frac{F_l(\eta,z) \cos(\chi) - F_{-l-1}(\eta,z)}{\sin(\chi)}\]where \(\chi = \sigma_l - \sigma_{-l-1} - (l+1/2) \pi\) and \(\sigma_l(\eta) = (\ln \Gamma(1+l+i\eta)-\ln \Gamma(1+l-i\eta))/(2i)\).
See
coulombf()for additional information.Plots
# Irregular Coulomb wave functions -- equivalent to figure 14.5 in A&S F1 = lambda x: coulombg(0,0,x) F2 = lambda x: coulombg(0,1,x) F3 = lambda x: coulombg(0,5,x) F4 = lambda x: coulombg(0,10,x) F5 = lambda x: coulombg(0,x/2,x) plot([F1,F2,F3,F4,F5], [0,30], [-2,2])
# Irregular Coulomb wave function in the complex plane cplot(lambda z: coulombg(1,1,z), points=50000)
Examples
Evaluation is supported for arbitrary magnitudes of \(z\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> coulombg(-2, 1.5, 3.5) 1.380011900612186346255524 >>> coulombg(2, 1.5, 3.5) 1.919153700722748795245926 >>> coulombg(-2, 1.5, '1e-10') 201126715824.7329115106793 >>> coulombg(-2, 1.5, 1000) 0.1802071520691149410425512 >>> coulombg(-2, 1.5, 10**10) 0.652103020061678070929794
The following reproduces a table in Abramowitz & Stegun, at twice the precision:
>>> mp.dps = 10 >>> eta = 2; z = 5 >>> for l in [1, 2, 3, 4, 5]: ... print("%s %s %s" % (l, coulombg(l,eta,z), ... -diff(lambda z: coulombg(l,eta,z), z))) ... 1 1.08148276 0.6028279961 2 1.496877075 0.5661803178 3 2.048694714 0.7959909551 4 3.09408669 1.731802374 5 5.629840456 4.549343289
Evaluation close to the singularity at \(z = 0\):
>>> mp.dps = 15 >>> coulombg(0,10,1) 3088184933.67358 >>> coulombg(0,10,'1e-10') 5554866000719.8 >>> coulombg(0,10,'1e-100') 5554866221524.1
Evaluation with a half-integer value for \(l\):
>>> coulombg(1.5, 1, 10) 0.852320038297334
coulombc()¶
- mpmath.coulombc(l, eta)¶
Gives the normalizing Gamow constant for Coulomb wave functions,
\[C_l(\eta) = 2^l \exp\left(-\pi \eta/2 + [\ln \Gamma(1+l+i\eta) + \ln \Gamma(1+l-i\eta)]/2 - \ln \Gamma(2l+2)\right),\]where the log gamma function with continuous imaginary part away from the negative half axis (see
loggamma()) is implied.This function is used internally for the calculation of Coulomb wave functions, and automatically cached to make multiple evaluations with fixed \(l\), \(\eta\) fast.
Confluent U and Whittaker functions¶
hyperu()¶
- mpmath.hyperu(a, b, z)¶
Gives the Tricomi confluent hypergeometric function \(U\), also known as the Kummer or confluent hypergeometric function of the second kind. This function gives a second linearly independent solution to the confluent hypergeometric differential equation (the first is provided by \(\,_1F_1\) – see
hyp1f1()).Examples
Evaluation for arbitrary complex arguments:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> hyperu(2,3,4) 0.0625 >>> hyperu(0.25, 5, 1000) 0.1779949416140579573763523 >>> hyperu(0.25, 5, -1000) (0.1256256609322773150118907 - 0.1256256609322773150118907j)
The \(U\) function may be singular at \(z = 0\):
>>> hyperu(1.5, 2, 0) +inf >>> hyperu(1.5, -2, 0) 0.1719434921288400112603671
Verifying the differential equation:
>>> a, b = 1.5, 2 >>> f = lambda z: hyperu(a,b,z) >>> for z in [-10, 3, 3+4j]: ... chop(z*diff(f,z,2) + (b-z)*diff(f,z) - a*f(z)) ... 0.0 0.0 0.0
An integral representation:
>>> a,b,z = 2, 3.5, 4.25 >>> hyperu(a,b,z) 0.06674960718150520648014567 >>> quad(lambda t: exp(-z*t)*t**(a-1)*(1+t)**(b-a-1),[0,inf]) / gamma(a) 0.06674960718150520648014567
whitm()¶
- mpmath.whitm(k, m, z)¶
Evaluates the Whittaker function \(M(k,m,z)\), which gives a solution to the Whittaker differential equation
\[\frac{d^2f}{dz^2} + \left(-\frac{1}{4}+\frac{k}{z}+ \frac{(\frac{1}{4}-m^2)}{z^2}\right) f = 0.\]A second solution is given by
whitw().The Whittaker functions are defined in Abramowitz & Stegun, section 13.1. They are alternate forms of the confluent hypergeometric functions \(\,_1F_1\) and \(U\):
\[ \begin{align}\begin{aligned}M(k,m,z) = e^{-\frac{1}{2}z} z^{\frac{1}{2}+m} \,_1F_1(\tfrac{1}{2}+m-k, 1+2m, z)\\W(k,m,z) = e^{-\frac{1}{2}z} z^{\frac{1}{2}+m} U(\tfrac{1}{2}+m-k, 1+2m, z).\end{aligned}\end{align} \]Examples
Evaluation for arbitrary real and complex arguments is supported:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> whitm(1, 1, 1) 0.7302596799460411820509668 >>> whitm(1, 1, -1) (0.0 - 1.417977827655098025684246j) >>> whitm(j, j/2, 2+3j) (3.245477713363581112736478 - 0.822879187542699127327782j) >>> whitm(2, 3, 100000) 4.303985255686378497193063e+21707
Evaluation at zero:
>>> whitm(1,-1,0); whitm(1,-0.5,0); whitm(1,0,0) +inf nan 0.0
We can verify that
whitm()numerically satisfies the differential equation for arbitrarily chosen values:>>> k = mpf(0.25) >>> m = mpf(1.5) >>> f = lambda z: whitm(k,m,z) >>> for z in [-1, 2.5, 3, 1+2j]: ... chop(diff(f,z,2) + (-0.25 + k/z + (0.25-m**2)/z**2)*f(z)) ... 0.0 0.0 0.0 0.0
An integral involving both
whitm()andwhitw(), verifying evaluation along the real axis:>>> quad(lambda x: exp(-x)*whitm(3,2,x)*whitw(1,-2,x), [0,inf]) 3.438869842576800225207341 >>> 128/(21*sqrt(pi)) 3.438869842576800225207341
whitw()¶
- mpmath.whitw(k, m, z)¶
Evaluates the Whittaker function \(W(k,m,z)\), which gives a second solution to the Whittaker differential equation. (See
whitm().)Examples
Evaluation for arbitrary real and complex arguments is supported:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> whitw(1, 1, 1) 1.19532063107581155661012 >>> whitw(1, 1, -1) (-0.9424875979222187313924639 - 0.2607738054097702293308689j) >>> whitw(j, j/2, 2+3j) (0.1782899315111033879430369 - 0.01609578360403649340169406j) >>> whitw(2, 3, 100000) 1.887705114889527446891274e-21705 >>> whitw(-1, -1, 100) 1.905250692824046162462058e-24
Evaluation at zero:
>>> for m in [-1, -0.5, 0, 0.5, 1]: ... whitw(1, m, 0) ... +inf nan 0.0 nan +inf
We can verify that
whitw()numerically satisfies the differential equation for arbitrarily chosen values:>>> k = mpf(0.25) >>> m = mpf(1.5) >>> f = lambda z: whitw(k,m,z) >>> for z in [-1, 2.5, 3, 1+2j]: ... chop(diff(f,z,2) + (-0.25 + k/z + (0.25-m**2)/z**2)*f(z)) ... 0.0 0.0 0.0 0.0
Parabolic cylinder functions¶
pcfd()¶
- mpmath.pcfd(n, z, **kwargs)¶
Gives the parabolic cylinder function in Whittaker’s notation \(D_n(z) = U(-n-1/2, z)\) (see
pcfu()). It solves the differential equation\[y'' + \left(n + \frac{1}{2} - \frac{1}{4} z^2\right) y = 0.\]and can be represented in terms of Hermite polynomials (see
hermite()) as\[D_n(z) = 2^{-n/2} e^{-z^2/4} H_n\left(\frac{z}{\sqrt{2}}\right).\]Plots
# Parabolic cylinder function D_n(x) on the real line for n=0,1,2,3,4 d0 = lambda x: pcfd(0,x) d1 = lambda x: pcfd(1,x) d2 = lambda x: pcfd(2,x) d3 = lambda x: pcfd(3,x) d4 = lambda x: pcfd(4,x) plot([d0,d1,d2,d3,d4],[-7,7])
Examples
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> pcfd(0,0); pcfd(1,0); pcfd(2,0); pcfd(3,0) 1.0 0.0 -1.0 0.0 >>> pcfd(4,0); pcfd(-3,0) 3.0 0.6266570686577501256039413 >>> pcfd('1/2', 2+3j) (-5.363331161232920734849056 - 3.858877821790010714163487j) >>> pcfd(2, -10) 1.374906442631438038871515e-9
Verifying the differential equation:
>>> n = mpf(2.5) >>> y = lambda z: pcfd(n,z) >>> z = 1.75 >>> chop(diff(y,z,2) + (n+0.5-0.25*z**2)*y(z)) 0.0
Rational Taylor series expansion when \(n\) is an integer:
>>> taylor(lambda z: pcfd(5,z), 0, 7) [0.0, 15.0, 0.0, -13.75, 0.0, 3.96875, 0.0, -0.6015625]
pcfu()¶
- mpmath.pcfu(a, z, **kwargs)¶
Gives the parabolic cylinder function \(U(a,z)\), which may be defined for \(\Re(z) > 0\) in terms of the confluent U-function (see
hyperu()) by\[U(a,z) = 2^{-\frac{1}{4}-\frac{a}{2}} e^{-\frac{1}{4} z^2} U\left(\frac{a}{2}+\frac{1}{4}, \frac{1}{2}, \frac{1}{2}z^2\right)\]or, for arbitrary \(z\),
\[e^{-\frac{1}{4}z^2} U(a,z) = U(a,0) \,_1F_1\left(-\tfrac{a}{2}+\tfrac{1}{4}; \tfrac{1}{2}; -\tfrac{1}{2}z^2\right) + U'(a,0) z \,_1F_1\left(-\tfrac{a}{2}+\tfrac{3}{4}; \tfrac{3}{2}; -\tfrac{1}{2}z^2\right).\]Examples
Connection to other functions:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> z = mpf(3) >>> pcfu(0.5,z) 0.03210358129311151450551963 >>> sqrt(pi/2)*exp(z**2/4)*erfc(z/sqrt(2)) 0.03210358129311151450551963 >>> pcfu(0.5,-z) 23.75012332835297233711255 >>> sqrt(pi/2)*exp(z**2/4)*erfc(-z/sqrt(2)) 23.75012332835297233711255 >>> pcfu(0.5,-z) 23.75012332835297233711255 >>> sqrt(pi/2)*exp(z**2/4)*erfc(-z/sqrt(2)) 23.75012332835297233711255
pcfv()¶
- mpmath.pcfv(a, z, **kwargs)¶
Gives the parabolic cylinder function \(V(a,z)\), which can be represented in terms of
pcfu()as\[V(a,z) = \frac{\Gamma(a+\tfrac{1}{2}) (U(a,-z)-\sin(\pi a) U(a,z)}{\pi}.\]Examples
Wronskian relation between \(U\) and \(V\):
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> a, z = 2, 3 >>> pcfu(a,z)*diff(pcfv,(a,z),(0,1))-diff(pcfu,(a,z),(0,1))*pcfv(a,z) 0.7978845608028653558798921 >>> sqrt(2/pi) 0.7978845608028653558798921 >>> a, z = 2.5, 3 >>> pcfu(a,z)*diff(pcfv,(a,z),(0,1))-diff(pcfu,(a,z),(0,1))*pcfv(a,z) 0.7978845608028653558798921 >>> a, z = 0.25, -1 >>> pcfu(a,z)*diff(pcfv,(a,z),(0,1))-diff(pcfu,(a,z),(0,1))*pcfv(a,z) 0.7978845608028653558798921 >>> a, z = 2+1j, 2+3j >>> chop(pcfu(a,z)*diff(pcfv,(a,z),(0,1))-diff(pcfu,(a,z),(0,1))*pcfv(a,z)) 0.7978845608028653558798921
pcfw()¶
- mpmath.pcfw(a, z, **kwargs)¶
Gives the parabolic cylinder function \(W(a,z)\) defined in (DLMF 12.14).
Examples
Value at the origin:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> a = mpf(0.25) >>> pcfw(a,0) 0.9722833245718180765617104 >>> power(2,-0.75)*sqrt(abs(gamma(0.25+0.5j*a)/gamma(0.75+0.5j*a))) 0.9722833245718180765617104 >>> diff(pcfw,(a,0),(0,1)) -0.5142533944210078966003624 >>> -power(2,-0.25)*sqrt(abs(gamma(0.75+0.5j*a)/gamma(0.25+0.5j*a))) -0.5142533944210078966003624