import numpy as np
from scipy.special import jv, jn_zeros
import matplotlib.pyplot as plt

Applications of Bessel functions

Documentation of each function

jv?

Signature:       jv(*args, **kwargs)
Type:            ufunc
String form:     <ufunc 'jv'>
File:            ~/.venv/default/lib/python3.12/site-packages/numpy/__init__.py
Docstring:      
jv(x1, x2, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])


jv(v, z, out=None)

Bessel function of the first kind of real order and complex argument.

Parameters
----------
v : array_like
    Order (float).
z : array_like
    Argument (float or complex).
out : ndarray, optional
    Optional output array for the function values

Returns
-------
J : scalar or ndarray
    Value of the Bessel function, :math:`J_v(z)`.

See Also
--------
jve : :math:`J_v` with leading exponential behavior stripped off.
spherical_jn : spherical Bessel functions.
j0 : faster version of this function for order 0.
j1 : faster version of this function for order 1.

Notes
-----
For positive `v` values, the computation is carried out using the AMOS
[1]_ `zbesj` routine, which exploits the connection to the modified
Bessel function :math:`I_v`,

.. math::
    J_v(z) = \exp(v\pi\imath/2) I_v(-\imath z)\qquad (\Im z > 0)

    J_v(z) = \exp(-v\pi\imath/2) I_v(\imath z)\qquad (\Im z < 0)

For negative `v` values the formula,

.. math:: J_{-v}(z) = J_v(z) \cos(\pi v) - Y_v(z) \sin(\pi v)

is used, where :math:`Y_v(z)` is the Bessel function of the second
kind, computed using the AMOS routine `zbesy`.  Note that the second
term is exactly zero for integer `v`; to improve accuracy the second
term is explicitly omitted for `v` values such that `v = floor(v)`.

Not to be confused with the spherical Bessel functions (see `spherical_jn`).

References
----------
.. [1] Donald E. Amos, "AMOS, A Portable Package for Bessel Functions
       of a Complex Argument and Nonnegative Order",
       http://netlib.org/amos/

Examples
--------
Evaluate the function of order 0 at one point.

>>> from scipy.special import jv
>>> jv(0, 1.)
0.7651976865579666

Evaluate the function at one point for different orders.

>>> jv(0, 1.), jv(1, 1.), jv(1.5, 1.)
(0.7651976865579666, 0.44005058574493355, 0.24029783912342725)

The evaluation for different orders can be carried out in one call by
providing a list or NumPy array as argument for the `v` parameter:

>>> jv([0, 1, 1.5], 1.)
array([0.76519769, 0.44005059, 0.24029784])

Evaluate the function at several points for order 0 by providing an
array for `z`.

>>> import numpy as np
>>> points = np.array([-2., 0., 3.])
>>> jv(0, points)
array([ 0.22389078,  1.        , -0.26005195])

If `z` is an array, the order parameter `v` must be broadcastable to
the correct shape if different orders shall be computed in one call.
To calculate the orders 0 and 1 for an 1D array:

>>> orders = np.array([[0], [1]])
>>> orders.shape
(2, 1)

>>> jv(orders, points)
array([[ 0.22389078,  1.        , -0.26005195],
       [-0.57672481,  0.        ,  0.33905896]])

Plot the functions of order 0 to 3 from -10 to 10.

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots()
>>> x = np.linspace(-10., 10., 1000)
>>> for i in range(4):
...     ax.plot(x, jv(i, x), label=f'$J_{i!r}$')
>>> ax.legend()
>>> plt.show()
Class docstring:
Functions that operate element by element on whole arrays.

To see the documentation for a specific ufunc, use `info`.  For
example, ``np.info(np.sin)``.  Because ufuncs are written in C
(for speed) and linked into Python with NumPy's ufunc facility,
Python's help() function finds this page whenever help() is called
on a ufunc.

A detailed explanation of ufuncs can be found in the docs for :ref:`ufuncs`.

**Calling ufuncs:** ``op(*x[, out], where=True, **kwargs)``

Apply `op` to the arguments `*x` elementwise, broadcasting the arguments.

The broadcasting rules are:

* Dimensions of length 1 may be prepended to either array.
* Arrays may be repeated along dimensions of length 1.

Parameters
----------
*x : array_like
    Input arrays.
out : ndarray, None, or tuple of ndarray and None, optional
    Alternate array object(s) in which to put the result; if provided, it
    must have a shape that the inputs broadcast to. A tuple of arrays
    (possible only as a keyword argument) must have length equal to the
    number of outputs; use None for uninitialized outputs to be
    allocated by the ufunc.
where : array_like, optional
    This condition is broadcast over the input. At locations where the
    condition is True, the `out` array will be set to the ufunc result.
    Elsewhere, the `out` array will retain its original value.
    Note that if an uninitialized `out` array is created via the default
    ``out=None``, locations within it where the condition is False will
    remain uninitialized.
**kwargs
    For other keyword-only arguments, see the :ref:`ufunc docs <ufuncs.kwargs>`.

Returns
-------
r : ndarray or tuple of ndarray
    `r` will have the shape that the arrays in `x` broadcast to; if `out` is
    provided, it will be returned. If not, `r` will be allocated and
    may contain uninitialized values. If the function has more than one
    output, then the result will be a tuple of arrays.

jn_zeros?

Signature: jn_zeros(n, nt)
Docstring:
Compute zeros of integer-order Bessel functions Jn.

Compute `nt` zeros of the Bessel functions :math:`J_n(x)` on the
interval :math:`(0, \infty)`. The zeros are returned in ascending
order. Note that this interval excludes the zero at :math:`x = 0`
that exists for :math:`n > 0`.

Parameters
----------
n : int
    Order of Bessel function
nt : int
    Number of zeros to return

Returns
-------
ndarray
    First `nt` zeros of the Bessel function.

See Also
--------
jv: Real-order Bessel functions of the first kind
jnp_zeros: Zeros of :math:`Jn'`

References
----------
.. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
       Functions", John Wiley and Sons, 1996, chapter 5.
       https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html

Examples
--------
Compute the first four positive roots of :math:`J_3`.

>>> from scipy.special import jn_zeros
>>> jn_zeros(3, 4)
array([ 6.3801619 ,  9.76102313, 13.01520072, 16.22346616])

Plot :math:`J_3` and its first four positive roots. Note
that the root located at 0 is not returned by `jn_zeros`.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.special import jn, jn_zeros
>>> j3_roots = jn_zeros(3, 4)
>>> xmax = 18
>>> xmin = -1
>>> x = np.linspace(xmin, xmax, 500)
>>> fig, ax = plt.subplots()
>>> ax.plot(x, jn(3, x), label=r'$J_3$')
>>> ax.scatter(j3_roots, np.zeros((4, )), s=30, c='r',
...            label=r"$J_3$_Zeros", zorder=5)
>>> ax.scatter(0, 0, s=30, c='k',
...            label=r"Root at 0", zorder=5)
>>> ax.hlines(0, 0, xmax, color='k')
>>> ax.set_xlim(xmin, xmax)
>>> plt.legend()
>>> plt.show()
File:      ~/.venv/default/lib/python3.12/site-packages/scipy/special/_basic.py
Type:      function

Plots

x = np.linspace(0,20,100)
plt.plot(x,jv(0,x),label="J_0(x)")
plt.plot(x,jv(1,x),label="J_1(x)")
plt.plot(x,jv(2,x),label="J_2(x)")
plt.plot(x,jv(1/2,x),label="J$_{1/2}$(x)")
plt.plot(x,jv(3/2,x),label="J$_{3/2}$(x)")
plt.plot(x,jv(5/2,x),label="J$_{5/2}$(x)")
plt.legend()

<matplotlib.legend.Legend at 0x10e6cf560>

Example: Fraunhofer diffraction integral

(from Arfken) We need to solve the integral

$$\Phi = \int_0^a r dr \int_0^{2\pi} e^{ibr\cos\theta} d\theta$$

We know the integral representation $$J_0(x) = \frac{1}{2\pi}\int_0^{2\pi} e^{ix\cos\theta} d\theta$$

and a recurrence formula

$$ \frac{d}{dx} [x^n J_n(x)] = x^n J_{n-1}(x) $$

The integral on $\theta$ on $\Phi$ reduces to $2\pi J_0(br)$. Therefore, we need to integrate $\Phi = 2\pi \int_0^a J_0(br) r dr$. This can be done with the recurrence formula.

wavelength_green = 5.5e-5 #cm
aperture = 0.5 # cm
alpha = np.linspace(0,np.pi/10000,100) # small angles

b = 2*np.pi/wavelength_green * np.sin(alpha)
a = aperture
r = np.linspace(0,a,1000)
dr = np.diff(r)[0]

Analytical integration, numerical evaluation

Phi = 2*np.pi*a/b * jv(1,a*b)
plt.plot(alpha,Phi,label=r"$\Phi$")
plt.plot(alpha,Phi**2,label=r"$\Phi^2$")
plt.legend()

/var/folders/2m/tgbjydr93_9_19gzzmp_v_3w0000gn/T/ipykernel_52857/1956156308.py:1: RuntimeWarning: divide by zero encountered in divide
  Phi = 2*np.pi*a/b * jv(1,a*b)
/var/folders/2m/tgbjydr93_9_19gzzmp_v_3w0000gn/T/ipykernel_52857/1956156308.py:1: RuntimeWarning: invalid value encountered in multiply
  Phi = 2*np.pi*a/b * jv(1,a*b)

<matplotlib.legend.Legend at 0x10e877500>

Example: zeros of the Bessel function

(From Arfken) We need to solve the Bessel differential equation (that comes out in a problem of standing electromagnetic waves in a hollow metallic cylinder). The solution is $J_{m}(n\rho)$, with $n$ being a constant, $m \in \mathbb N_0$ and $\rho$ is the variable. The boundary conditions at $\rho = a$ requires $J_m(na) = 0$. Evaluate numerically the zeros of the Bessel function.

m = 0

# zeros of the Bessel function. Careful! In the scipy documentation, m <-> n
zeros = jn_zeros(m,6)
zeros

array([ 2.40482556,  5.52007811,  8.65372791, 11.79153444, 14.93091771,
       18.07106397])

If $n\in \mathbb R$, both $\rho$ and $a$ are free. However, if $n\in \mathbb Z$, this fixes the values that $a$ can have to satisfy the Bessel equation. For example, let's take n = 2.

n = 2 # effect: scales rho

rho = np.linspace(0,10,1000) # in the electromagnetism problem, this is the cylindrical radial coordinate

Let's plot the situation

# plot of the Bessel function
plt.plot(rho,jv(m,n*rho))

# plot of the zeros
plt.scatter(zeros/n,zeros*0)

<matplotlib.collections.PathCollection at 0x10e83baa0>

So, $a$ can have the values

zeros/n

array([1.20241278, 2.76003906, 4.32686396, 5.89576722, 7.46545885,
       9.03553198])

This is analogous to the problem of a standing wave on a rope. Given the (spatial) harmonic oscillator equation (i.e., the spatial wave equation) and an angular frequency (which is the same role that $m$ plays here in the Bessel equation, but careful, the Bessel functions have no frequency despite having an oscillatory behavior), only ropes of a certain length can accommodate a given standing wave (i.e., only certain lengths satisfy the boundary conditions).