Navigation

  • index
  • next |
  • previous |
  • mpmath v0.13 documentation »
  • Mathematical functions »

Bessel functions and related functions¶

Bessel functions¶

besselj() / j0() / j1()¶

mpmath.functions.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.

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
mpmath.functions.j0(x)¶
Computes the Bessel function J_0(x). See besselj().
mpmath.functions.j1(x)¶
Computes the Bessel function J_1(x). See besselj().

bessely()¶

mpmath.functions.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.

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.functions.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.

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.functions.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.

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

Hankel functions¶

hankel1()¶

mpmath.functions.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).

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.functions.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).

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.functions.ber(ctx, n, z)¶

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().

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()¶

mpmath.functions.bei(ctx, n, z)¶
Computes the Kelvin function bei, which for real arguments gives the imaginary part of the Bessel J function of a rotated argument. See ber().

ker()¶

mpmath.functions.ker(ctx, n, z)¶

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().

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)

kei()¶

mpmath.functions.kei(ctx, n, z)¶
Computes the Kelvin function kei, which for real arguments gives the imaginary part of the (rescaled) Bessel K function of a rotated argument. See ker().

Struve functions¶

struveh()¶

mpmath.functions.struveh(ctx, n, z)¶

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.functions.struvel(ctx, n, z)¶

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

Airy functions¶

airyai()¶

mpmath.functions.airyai(x)¶

Computes the Airy function \mathrm{Ai}(x), which is a solution of the Airy differential equation y''-xy=0. The Ai-function behaves roughly like a slowly decaying sine wave for x < 0 and like a decreasing exponential for x > 0.

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 x = 0 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 x < 0, 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

airybi()¶

mpmath.functions.airybi(x)¶

Computes the Airy function \mathrm{Bi}(x), which is a solution of the Airy differential equation y''-xy=0. The Bi-function behaves roughly like a slowly decaying sine wave for x < 0 and like an increasing exponential for x > 0.

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 x = 0 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 x < 0, 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

Coulomb wave functions¶

coulombf()¶

mpmath.functions.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).

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 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

  1. I.J. Thompson & A.R. Barnett, “Coulomb and Bessel Functions of Complex Arguments and Order”, J. Comp. Phys., vol 64, no. 2, June 1986.
  2. N. Michel, “Precise Coulomb wave functions for a wide range of complex l, \eta and z“, http://arxiv.org/abs/physics/0702051v1

coulombg()¶

mpmath.functions.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.

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 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.functions.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.

Table Of Contents

  • Bessel functions and related functions
    • Bessel functions
      • besselj() / j0() / j1()
      • bessely()
      • besseli()
      • besselk()
    • Hankel functions
      • hankel1()
      • hankel2()
    • Kelvin functions
      • ber()
      • bei()
      • ker()
      • kei()
    • Struve functions
      • struveh()
      • struvel()
    • Airy functions
      • airyai()
      • airybi()
    • Coulomb wave functions
      • coulombf()
      • coulombg()
      • coulombc()

Previous topic

Exponential integrals and error functions

Next topic

Orthogonal polynomials

This Page

  • Show Source

Quick search

Enter search terms or a module, class or function name.

Navigation

  • index
  • next |
  • previous |
  • mpmath v0.13 documentation »
  • Mathematical functions »
© Copyright 2009, Fredrik Johansson. Last updated on Sep 18, 2009. Created using Sphinx 0.6.3.