Bessel Functions¶
This module provides symbolic Bessel Functions. These functions use the mpmath library for numerical evaluation and Maxima, GiNaC, Pynac for symbolics.
The main objects which are exported from this module are:
bessel_J
– The Bessel J functionbessel_Y
– The Bessel Y functionbessel_I
– The Bessel I functionbessel_K
– The Bessel K functionBessel
– A factory function for producing Bessel functions of various kinds and orders
Bessel functions, first defined by the Swiss mathematician Daniel Bernoulli and named after Friedrich Bessel, are canonical solutions y(x) of Bessel’s differential equation:
for an arbitrary complex number
(the order).
In this module,
denotes the unique solution of Bessel’s equation which is non-singular at
. This function is known as the Bessel Function of the First Kind. This function also arises as a special case of the hypergeometric function
:
The second linearly independent solution to Bessel’s equation (which is singular at
) is denoted by
and is called the Bessel Function of the Second Kind:
There are also two commonly used combinations of the Bessel J and Y Functions. The Bessel I Function, or the Modified Bessel Function of the First Kind, is defined by:
The Bessel K Function, or the Modified Bessel Function of the Second Kind, is defined by:
We should note here that the above formulas for Bessel Y and K functions should be understood as limits when
is an integer.
It follows from Bessel’s differential equation that the derivative of
with respect to
is:
Another important formulation of the two linearly independent solutions to Bessel’s equation are the Hankel functions
and
, defined by:
where
is the imaginary unit (and
and
are the usual J- and Y-Bessel functions). These linear combinations are also known as Bessel functions of the third kind; they are also two linearly independent solutions of Bessel’s differential equation. They are named for Hermann Hankel.
EXAMPLES:
Evaluate the Bessel J function symbolically and numerically:
sage: bessel_J(0, x) bessel_J(0, x) sage: bessel_J(0, 0) bessel_J(0, 0) sage: bessel_J(0, x).diff(x) -1/2*bessel_J(1, x) + 1/2*bessel_J(-1, x) sage: N(bessel_J(0, 0), digits = 20) 1.0000000000000000000 sage: find_root(bessel_J(0,x), 0, 5) 2.404825557695773Plot the Bessel J function:
sage: f(x) = Bessel(0)(x); f x |--> bessel_J(0, x) sage: plot(f, (x, 1, 10)) Graphics object consisting of 1 graphics primitiveVisualize the Bessel Y function on the complex plane (set plot_points to a higher value to get more detail):
sage: complex_plot(bessel_Y(0, x), (-5, 5), (-5, 5), plot_points=20) Graphics object consisting of 1 graphics primitiveEvaluate a combination of Bessel functions:
sage: f(x) = bessel_J(1, x) - bessel_Y(0, x) sage: f(pi) bessel_J(1, pi) - bessel_Y(0, pi) sage: f(pi).n() -0.0437509653365599 sage: f(pi).n(digits=50) -0.043750965336559909054985168023342675387737118378169Symbolically solve a second order differential equation with initial conditions
and
in terms of Bessel functions:
sage: y = function('y', x) sage: a, b = var('a, b') sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0 sage: f = desolve(diffeq, y, [1, a, b]); f (a*bessel_Y(1, 1) + b*bessel_Y(0, 1))*bessel_J(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) - (a*bessel_J(1, 1) + b*bessel_J(0, 1))*bessel_Y(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1))For more examples, see the docstring for
Bessel()
.
AUTHORS:
- Benjamin Jones (2012-12-27): initial version
- Some of the documentation here has been adapted from David Joyner’s original documentation of Sage’s special functions module (2006).
REFERENCES:
- Abramowitz and Stegun: Handbook of Mathematical Functions, http://www.math.sfu.ca/~cbm/aands/
- http://en.wikipedia.org/wiki/Bessel_function
- mpmath Library Bessel Functions
-
sage.functions.bessel.
Bessel
(*args, **kwds)¶ A function factory that produces symbolic I, J, K, and Y Bessel functions. There are several ways to call this function:
Bessel(order, type)
Bessel(order)
– type defaults to ‘J’Bessel(order, typ=T)
Bessel(typ=T)
– order is unspecified, this is a 2-parameter functionBessel()
– order is unspecified, type is ‘J’
where
order
can be any integer and T must be one of the strings ‘I’, ‘J’, ‘K’, or ‘Y’.See the EXAMPLES below.
EXAMPLES:
Construction of Bessel functions with various orders and types:
sage: Bessel() bessel_J sage: Bessel(1)(x) bessel_J(1, x) sage: Bessel(1, 'Y')(x) bessel_Y(1, x) sage: Bessel(-2, 'Y')(x) bessel_Y(-2, x) sage: Bessel(typ='K') bessel_K sage: Bessel(0, typ='I')(x) bessel_I(0, x)
Evaluation:
sage: f = Bessel(1) sage: f(3.0) 0.339058958525936 sage: f(3) bessel_J(1, 3) sage: f(3).n(digits=50) 0.33905895852593645892551459720647889697308041819801 sage: g = Bessel(typ='J') sage: g(1,3) bessel_J(1, 3) sage: g(2, 3+I).n() 0.634160370148554 + 0.0253384000032695*I sage: abs(numerical_integral(1/pi*cos(3*sin(x)), 0.0, pi)[0] - Bessel(0, 'J')(3.0)) < 1e-15 True
Symbolic calculus:
sage: f(x) = Bessel(0, 'J')(x) sage: derivative(f, x) x |--> -1/2*bessel_J(1, x) + 1/2*bessel_J(-1, x) sage: derivative(f, x, x) x |--> 1/4*bessel_J(2, x) - 1/2*bessel_J(0, x) + 1/4*bessel_J(-2, x)
Verify that
satisfies Bessel’s differential equation numerically using the
test_relation()
method:sage: y = bessel_J(0, x) sage: diffeq = x^2*derivative(y,x,x) + x*derivative(y,x) + x^2*y == 0 sage: diffeq.test_relation(proof=False) True
Conversion to other systems:
sage: x,y = var('x,y') sage: f = maxima(Bessel(typ='K')(x,y)) sage: f.derivative('_SAGE_VAR_x') %pi*csc(%pi*_SAGE_VAR_x)*('diff(bessel_i(-_SAGE_VAR_x,_SAGE_VAR_y),_SAGE_VAR_x,1)-'diff(bessel_i(_SAGE_VAR_x,_SAGE_VAR_y),_SAGE_VAR_x,1))/2-%pi*bessel_k(_SAGE_VAR_x,_SAGE_VAR_y)*cot(%pi*_SAGE_VAR_x) sage: f.derivative('_SAGE_VAR_y') -(bessel_k(_SAGE_VAR_x+1,_SAGE_VAR_y)+bessel_k(_SAGE_VAR_x-1,_SAGE_VAR_y))/2
Compute the particular solution to Bessel’s Differential Equation that satisfies
and
, then verify the initial conditions and plot it:
sage: y = function('y', x) sage: diffeq = x^2*diff(y,x,x) + x*diff(y,x) + x^2*y == 0 sage: f = desolve(diffeq, y, [1, 1, 1]); f (bessel_Y(1, 1) + bessel_Y(0, 1))*bessel_J(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) - (bessel_J(1, 1) + bessel_J(0, 1))*bessel_Y(0, x)/(bessel_J(0, 1)*bessel_Y(1, 1) - bessel_J(1, 1)*bessel_Y(0, 1)) sage: f.subs(x=1).n() # numerical verification 1.00000000000000 sage: fp = f.diff(x) sage: fp.subs(x=1).n() 1.00000000000000 sage: f.subs(x=1).simplify_full() # symbolic verification 1 sage: fp = f.diff(x) sage: fp.subs(x=1).simplify_full() 1 sage: plot(f, (x,0,5)) Graphics object consisting of 1 graphics primitive
Plotting:
sage: f(x) = Bessel(0)(x); f x |--> bessel_J(0, x) sage: plot(f, (x, 1, 10)) Graphics object consisting of 1 graphics primitive sage: plot([ Bessel(i, 'J') for i in range(5) ], 2, 10) Graphics object consisting of 5 graphics primitives sage: G = Graphics() sage: G += sum([ plot(Bessel(i), 0, 4*pi, rgbcolor=hue(sin(pi*i/10))) for i in range(5) ]) sage: show(G)
A recreation of Abramowitz and Stegun Figure 9.1:
sage: G = plot(Bessel(0, 'J'), 0, 15, color='black') sage: G += plot(Bessel(0, 'Y'), 0, 15, color='black') sage: G += plot(Bessel(1, 'J'), 0, 15, color='black', linestyle='dotted') sage: G += plot(Bessel(1, 'Y'), 0, 15, color='black', linestyle='dotted') sage: show(G, ymin=-1, ymax=1)
-
class
sage.functions.bessel.
Function_Bessel_I
¶ Bases:
sage.symbolic.function.BuiltinFunction
The Bessel I function, or the Modified Bessel Function of the First Kind.
DEFINITION:
EXAMPLES:
sage: bessel_I(1, x) bessel_I(1, x) sage: bessel_I(1.0, 1.0) 0.565159103992485 sage: n = var('n') sage: bessel_I(n, x) bessel_I(n, x) sage: bessel_I(2, I).n() -0.114903484931900
Examples of symbolic manipulation:
sage: a = bessel_I(pi, bessel_I(1, I)) sage: N(a, digits=20) 0.00026073272117205890528 - 0.0011528954889080572266*I sage: f = bessel_I(2, x) sage: f.diff(x) 1/2*bessel_I(3, x) + 1/2*bessel_I(1, x)
Special identities that bessel_I satisfies:
sage: bessel_I(1/2, x) sqrt(2)*sqrt(1/(pi*x))*sinh(x) sage: eq = bessel_I(1/2, x) == bessel_I(0.5, x) sage: eq.test_relation() True sage: bessel_I(-1/2, x) sqrt(2)*sqrt(1/(pi*x))*cosh(x) sage: eq = bessel_I(-1/2, x) == bessel_I(-0.5, x) sage: eq.test_relation() True
Examples of asymptotic behavior:
sage: limit(bessel_I(0, x), x=oo) +Infinity sage: limit(bessel_I(0, x), x=0) 1
High precision and complex valued inputs:
sage: bessel_I(0, 1).n(128) 1.2660658777520083355982446252147175376 sage: bessel_I(0, RealField(200)(1)) 1.2660658777520083355982446252147175376076703113549622068081 sage: bessel_I(0, ComplexField(200)(0.5+I)) 0.80644357583493619472428518415019222845373366024179916785502 + 0.22686958987911161141397453401487525043310874687430711021434*I
Visualization (set plot_points to a higher value to get more detail):
sage: plot(bessel_I(1,x), (x,0,5), color='blue') Graphics object consisting of 1 graphics primitive sage: complex_plot(bessel_I(1, x), (-5, 5), (-5, 5), plot_points=20) Graphics object consisting of 1 graphics primitive
ALGORITHM:
Numerical evaluation is handled by the mpmath library. Symbolics are handled by a combination of Maxima and Sage (Ginac/Pynac).TESTS:
sage: N(bessel_I(1,1),500) 0.565159103992485027207696027609863307328899621621092009480294489479255640964371134092664997766814410064677886055526302676857637684917179812041131208121
Check whether the return value is real whenever the argument is real (trac ticket #10251):
sage: bessel_I(5, 1.5) in RR True
-
class
sage.functions.bessel.
Function_Bessel_J
¶ Bases:
sage.symbolic.function.BuiltinFunction
The Bessel J Function, denoted by bessel_J(
, x) or
. As a Taylor series about
it is equal to:
The parameter
is called the order and may be any real or complex number; however, integer and half-integer values are most common. It is defined for all complex numbers
when
is an integer or greater than zero and it diverges as
for negative non-integer values of
.
For integer orders
there is an integral representation:
This function also arises as a special case of the hypergeometric function
:
EXAMPLES:
sage: bessel_J(1.0, 1.0) 0.440050585744933 sage: bessel_J(2, I).n(digits=30) -0.135747669767038281182852569995 sage: bessel_J(1, x) bessel_J(1, x) sage: n = var('n') sage: bessel_J(n, x) bessel_J(n, x)
Examples of symbolic manipulation:
sage: a = bessel_J(pi, bessel_J(1, I)); a bessel_J(pi, bessel_J(1, I)) sage: N(a, digits=20) 0.00059023706363796717363 - 0.0026098820470081958110*I sage: f = bessel_J(2, x) sage: f.diff(x) -1/2*bessel_J(3, x) + 1/2*bessel_J(1, x)
Comparison to a well-known integral representation of
:
sage: A = numerical_integral(1/pi*cos(x - sin(x)), 0, pi) sage: A[0] # abs tol 1e-14 0.44005058574493355 sage: bessel_J(1.0, 1.0) - A[0] < 1e-15 True
Integration is supported directly and through Maxima:
sage: f = bessel_J(2, x) sage: f.integrate(x) 1/24*x^3*hypergeometric((3/2,), (5/2, 3), -1/4*x^2) sage: m = maxima(bessel_J(2, x)) sage: m.integrate(x) hypergeometric([3/2],[5/2,3],-_SAGE_VAR_x^2/4)*_SAGE_VAR_x^3/24
Visualization (set plot_points to a higher value to get more detail):
sage: plot(bessel_J(1,x), (x,0,5), color='blue') Graphics object consisting of 1 graphics primitive sage: complex_plot(bessel_J(1, x), (-5, 5), (-5, 5), plot_points=20) Graphics object consisting of 1 graphics primitive
ALGORITHM:
Numerical evaluation is handled by the mpmath library. Symbolics are handled by a combination of Maxima and Sage (Ginac/Pynac).Check whether the return value is real whenever the argument is real (trac ticket #10251):
sage: bessel_J(5, 1.5) in RR True
-
class
sage.functions.bessel.
Function_Bessel_K
¶ Bases:
sage.symbolic.function.BuiltinFunction
The Bessel K function, or the modified Bessel function of the second kind.
DEFINITION:
EXAMPLES:
sage: bessel_K(1, x) bessel_K(1, x) sage: bessel_K(1.0, 1.0) 0.601907230197235 sage: n = var('n') sage: bessel_K(n, x) bessel_K(n, x) sage: bessel_K(2, I).n() -2.59288617549120 + 0.180489972066962*I
Examples of symbolic manipulation:
sage: a = bessel_K(pi, bessel_K(1, I)); a bessel_K(pi, bessel_K(1, I)) sage: N(a, digits=20) 3.8507583115005220157 + 0.068528298579883425792*I sage: f = bessel_K(2, x) sage: f.diff(x) -1/2*bessel_K(3, x) - 1/2*bessel_K(1, x) sage: bessel_K(1/2, x) bessel_K(1/2, x) sage: bessel_K(1/2, -1) bessel_K(1/2, -1) sage: bessel_K(1/2, 1) sqrt(1/2)*sqrt(pi)*e^(-1)
Examples of asymptotic behavior:
sage: bessel_K(0, 0.0) +infinity sage: limit(bessel_K(0, x), x=0) +Infinity sage: limit(bessel_K(0, x), x=oo) 0
High precision and complex valued inputs:
sage: bessel_K(0, 1).n(128) 0.42102443824070833333562737921260903614 sage: bessel_K(0, RealField(200)(1)) 0.42102443824070833333562737921260903613621974822666047229897 sage: bessel_K(0, ComplexField(200)(0.5+I)) 0.058365979093103864080375311643360048144715516692187818271179 - 0.67645499731334483535184142196073004335768129348518210260256*I
Visualization (set plot_points to a higher value to get more detail):
sage: plot(bessel_K(1,x), (x,0,5), color='blue') Graphics object consisting of 1 graphics primitive sage: complex_plot(bessel_K(1, x), (-5, 5), (-5, 5), plot_points=20) Graphics object consisting of 1 graphics primitive
ALGORITHM:
Numerical evaluation is handled by the mpmath library. Symbolics are handled by a combination of Maxima and Sage (Ginac/Pynac).TESTS:
Verify that trac ticket #3426 is fixed:
The Bessel K function can be evaluated numerically at complex orders:
sage: bessel_K(10 * I, 10).n() 9.82415743819925e-8
For a fixed imaginary order and increasing, real, second component the value of Bessel K is exponentially decaying:
sage: for x in [10, 20, 50, 100, 200]: print bessel_K(5*I, x).n() 5.27812176514912e-6 3.11005908421801e-10 2.66182488515423e-23 - 8.59622057747552e-58*I 4.11189776828337e-45 - 1.01494840019482e-80*I 1.15159692553603e-88 - 6.75787862113718e-125*I
Check whether the return value is real whenever the argument is real (trac ticket #10251):
sage: bessel_K(5, 1.5) in RR True
-
class
sage.functions.bessel.
Function_Bessel_Y
¶ Bases:
sage.symbolic.function.BuiltinFunction
The Bessel Y functions, also known as the Bessel functions of the second kind, Weber functions, or Neumann functions.
is a holomorphic function of
on the complex plane, cut along the negative real axis. It is singular at
. When
is fixed,
is an entire function of the order
.
DEFINITION:
Its derivative with respect to
is:
EXAMPLES:
sage: bessel_Y(1, x) bessel_Y(1, x) sage: bessel_Y(1.0, 1.0) -0.781212821300289 sage: n = var('n') sage: bessel_Y(n, x) bessel_Y(n, x) sage: bessel_Y(2, I).n() 1.03440456978312 - 0.135747669767038*I sage: bessel_Y(0, 0).n() -infinity sage: bessel_Y(0, 1).n(128) 0.088256964215676957982926766023515162828
Examples of symbolic manipulation:
sage: a = bessel_Y(pi, bessel_Y(1, I)); a bessel_Y(pi, bessel_Y(1, I)) sage: N(a, digits=20) 4.2059146571791095708 + 21.307914215321993526*I sage: f = bessel_Y(2, x) sage: f.diff(x) -1/2*bessel_Y(3, x) + 1/2*bessel_Y(1, x)
High precision and complex valued inputs (see trac ticket #4230):
sage: bessel_Y(0, 1).n(128) 0.088256964215676957982926766023515162828 sage: bessel_Y(0, RealField(200)(1)) 0.088256964215676957982926766023515162827817523090675546711044 sage: bessel_Y(0, ComplexField(200)(0.5+I)) 0.077763160184438051408593468823822434235010300228009867784073 + 1.0142336049916069152644677682828326441579314239591288411739*I
Visualization (set plot_points to a higher value to get more detail):
sage: plot(bessel_Y(1,x), (x,0,5), color='blue') Graphics object consisting of 1 graphics primitive sage: complex_plot(bessel_Y(1, x), (-5, 5), (-5, 5), plot_points=20) Graphics object consisting of 1 graphics primitive
ALGORITHM:
Numerical evaluation is handled by the mpmath library. Symbolics are handled by a combination of Maxima and Sage (Ginac/Pynac).TESTS:
Check whether the return value is real whenever the argument is real (trac ticket #10251):
sage: bessel_Y(5, 1.5) in RR True
Coercion works correctly (see trac ticket #17130):
sage: r = bessel_Y(RealField(200)(1), 1.0); r -0.781212821300289 sage: parent(r) Real Field with 53 bits of precision sage: r = bessel_Y(RealField(200)(1), 1); r -0.78121282130028871654715000004796482054990639071644460784383 sage: parent(r) Real Field with 200 bits of precision