besselj(n, x, derivative=0) gives the Bessel function of the first kind
. Bessel functions of the first kind are defined as
solutions of the differential equation
which appears, among other things, when solving the radial
part of Laplace’s equation in cylindrical coordinates. This
equation has two solutions for given , where the
-function is the solution that is nonsingular at
.
For positive integer
,
behaves roughly like a sine
(odd
) or cosine (even
) multiplied by a magnitude factor
that decays slowly as
.
Generally, is a special case of the hypergeometric
function
:
With derivative = , the
-th derivative
is computed.
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 :
>>> 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 . 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 or
, 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(n, x, derivative=0) gives the Bessel function of the second kind,
For an integer, this formula should be understood as a
limit. With derivative =
, the
-th derivative
is computed.
Examples
Some values of :
>>> 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(n, x, derivative=0) gives the modified Bessel function of the first kind,
With derivative = , the
-th derivative
is computed.
Examples
Some values of :
>>> 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 , 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(n, x) gives the modified Bessel function of the second kind,
For an integer, this formula should be understood as a
limit.
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 is a singularity (logarithmic if
):
>>> 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
hankel1(n,x) computes the Hankel function of the first kind, which is the complex combination of Bessel functions given by
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(n,x) computes the Hankel function of the second kind, which is the complex combination of Bessel functions given by
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)
Computes the Kelvin function ber, which for real arguments gives the real part of the Bessel J function of a rotated argument
The imaginary part is given by bei().
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)
Computes the Kelvin function ker, which for real arguments gives the real part of the (rescaled) Bessel K function of a rotated argument
The imaginary part is given by kei().
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)
Gives the Struve function
which is a solution to the Struve differential equation
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
Gives the modified Struve function
which solves to the modified Struve differential equation
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
Computes the Airy function , which is
a solution of the Airy differential equation
.
The Ai-function behaves roughly like a slowly decaying
sine wave for
and like a decreasing exponential for
.
Limits and values include:
>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> airyai(0), 1/(3**(2/3.)*gamma(2/3.))
(0.355028053887817, 0.355028053887817)
>>> airyai(1)
0.135292416312881
>>> airyai(-1)
0.535560883292352
>>> airyai(inf)
0.0
>>> airyai(-inf)
0.0
Evaluation is supported for large arguments:
>>> airyai(-100)
0.176753393239553
>>> airyai(100)
2.63448215208818e-291
>>> airyai(50+50j)
(-5.31790195707456e-68 - 1.16358800377071e-67j)
>>> airyai(-50+50j)
(1.04124253736317e+158 + 3.3475255449236e+157j)
Huge arguments are also fine:
>>> airyai(10**10)
1.16223597829874e-289529654602171
>>> airyai(-10**10)
0.000173620644815282
>>> airyai(10**10*(1+j))
(5.71150868372136e-186339621747698 + 1.86724550696231e-186339621747697j)
The first negative root of the Ai function is:
>>> findroot(airyai, -2)
-2.33810741045977
We can verify the differential equation:
>>> for x in [-3.4, 0, 2.5, 1+2j]:
... print abs(diff(airyai, x, 2) - x*airyai(x)) < eps
...
True
True
True
True
The Taylor series expansion around starts with
the following coefficients (note that every third term
is zero):
>>> nprint(chop(taylor(airyai, 0, 5)))
[0.355028, -0.258819, 0.0, 5.91713e-2, -2.15683e-2, 0.0]
The Airy functions are a special case of Bessel functions.
For , we have:
>>> x = 3
>>> airyai(-x)
-0.378814293677658
>>> p = 2*(x**1.5)/3
>>> sqrt(x)*(besselj(1/3.,p) + besselj(-1/3.,p))/3
-0.378814293677658
Computes the Airy function , which is
a solution of the Airy differential equation
.
The Bi-function behaves roughly like a slowly decaying
sine wave for
and like an increasing exponential
for
.
Limits and values include:
>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> airybi(0), 1/(3**(1/6.)*gamma(2/3.))
(0.614926627446001, 0.614926627446001)
>>> airybi(1)
1.20742359495287
>>> airybi(-1)
0.103997389496945
>>> airybi(inf)
+inf
>>> airybi(-inf)
0.0
Evaluation is supported for large arguments:
>>> airybi(-100)
0.0242738876801601
>>> airybi(100)
6.0412239966702e+288
>>> airybi(50+50j)
(-5.32207626732144e+63 + 1.47845029116524e+65j)
>>> airybi(-50+50j)
(-3.3475255449236e+157 + 1.04124253736317e+158j)
Huge arguments are also fine:
>>> mp.dps = 15
>>> airybi(10**10)
1.36938578794354e+289529654602165
>>> airybi(-10**10)
0.00177565614169293
>>> airybi(10**10*(1+j))
(-6.5599559310962e+186339621747689 - 6.82246272698136e+186339621747690j)
The first negative root of the Bi function is:
>>> findroot(airybi, -1)
-1.17371322270913
We can verify the differential equation:
>>> for x in [-3.4, 0, 2.5, 1+2j]:
... print abs(diff(airybi, x, 2) - x*airybi(x)) < eps
...
True
True
True
True
The Taylor series expansion around starts with
the following coefficients (note that every third term
is zero):
>>> nprint(chop(taylor(airybi, 0, 5)))
[0.614927, 0.448288, 0.0, 0.102488, 3.73574e-2, 0.0]
The Airy functions are a special case of Bessel functions.
For , we have:
>>> x = 3
>>> airybi(-x)
-0.198289626374927
>>> p = 2*(x**1.5)/3
>>> sqrt(x/3)*(besselj(-1/3.,p) - besselj(1/3.,p))
-0.198289626374926
Calculates the regular Coulomb wave function
where the normalization constant is as calculated by
coulombc(). This function solves the differential equation
A second linearly independent solution is given by the irregular
Coulomb wave function (see coulombg())
and thus the general solution is
for arbitrary
constants
,
.
Physically, the Coulomb wave functions give the radial solution
to the Schrodinger equation for a point particle in a
potential;
is
then the radius and
,
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).
Examples
Evaluation is supported for arbitrary magnitudes of :
>>> 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 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
Calculates the irregular Coulomb wave function
where
and
.
See coulombf() for additional information.
Examples
Evaluation is supported for arbitrary magnitudes of :
>>> 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 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 :
>>> 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 :
>>> coulombg(1.5, 1, 10)
0.852320038297334
Gives the normalizing Gamow constant for Coulomb wave functions,
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 ,
fast.
Exponential integrals and error functions
Enter search terms or a module, class or function name.