Preface


This section concerns about Bessel equation and its solutions, known as Bessel functions. Daniel Bernoulli is generally credited with being the first to introduce the concept of Bessels functions in 1732. He used the function of zero order as a solution to the problem of an oscillating chain suspended at one end. In 1764 Leonhard Euler employed Bessel functions of both zero and integral orders in an analysis of vibrations of a stretched membrane, an investigation which was further developed by Lord Rayleigh in 1878, where he demonstrated that Bessels functions are particular cases of Laplaces functions.

Bessel, while receiving named credit for these functions, did not incorporate them into his work as an astronomer until 1817. The Bessel function was the result of Bessels study of a problem of Kepler for determining the motion of three bodies moving under mutual gravita- tion. In 1824, he incorporated Bessel functions in a study of planetary perturbations where the Bessel functions appear as coefficients in a series expansion of the indirect perturbation of a planet, that is the motion of the Sun caused by the perturbing body. It was likely Lagrange’s work on elliptical orbits that first suggested to Bessel to work on the Bessel functions.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VII of the course APMA0340
Introduction to Linear Algebra with Mathematica

Bessel Functions


The Bessel differential equation is the linear second-order ordinary differential equation given by
\[ x^2 \frac{{\text d}^2 y}{{\text d} x^2} + x\,\frac{{\text d} y}{{\text d} x} + \left( x^2 - \nu^2 \right) y(x) = 0 \qquad \mbox{or} \qquad \frac{\text d}{{\text d}x} \left( x \,\frac{{\text d} y}{{\text d} x} \right) + \left( x - \frac{\nu^2}{x} \right) y(x) = 0 , \]
where ν is a real constant, called the order of the Bessel equation. Upon substitution \( y = u\, x^{-1/2} , \) it is reduced to the equation without first derivative:
\[ u'' + \left( 1 + \frac{1 - 4\nu^2}{4x^2} \right) u(x) = 0 . \]
A Bessel equation is a special case of a confluent hypergeometric equation. Since x = 0 is a regular singular point for the Bessel equation, one of its solution can be bounded at this point but another linearly independent solution should be unbounded.

 

Transformed version of the Bessel equation


The self-adjoint differential equation
\[ \frac{\text d}{{\text d}x} \left( x^a \,\frac{{\text d} y}{{\text d} x} \right) + b\, x^c y(x) = 0 , \]
where a, b, and c are real numbers, can be transformed to a Bessel equation by transforming both independent and dependent variables. Because of the powers of x in the differential equation, it seems promissing to change variables from x, y(x) to t, u(t) according to
\[ t = A\, x^B , \qquad u = x^C y , \]
and to try to find A, B, C so that the new equation becomes a Bessel equation for u(t). It turns out that the plan works and one finds that under the change of variables
\[ t = \alpha \sqrt{b} \,x^{1/\alpha} \qquad \mbox{and} \qquad u = x^{-\nu /\alpha} y \]
the given self-adjoint equation is transformed into the Bessel equation of order ν,
\[ t^2 \frac{{\text d}^2 u}{{\text d} t^2} + t\,\frac{{\text d} u}{{\text d} t} + \left( t^2 - \nu^2 \right) u(t) = 0 , \]
if we choose
\[ \alpha = \frac{2}{c-a+2} \qquad\mbox{and} \qquad \nu = \frac{1-a}{c-a+2} . \]
The latter is meaningless when c - a + 2 = 0, but in that case the original equation
\[ \frac{\text d}{{\text d}x} \left( x^a \,\frac{{\text d} y}{{\text d} x} \right) + b\, x^{a-2} y(x) = 0 \]
becomes an Euler equation
\[ x^{a-2} \left[ x^2 \frac{{\text d}^2 y}{{\text d}x^2} + ax \, \frac{{\text d}y}{{\text d}x} + b\,y(x) \right] = 0 . \]
Thus, if
\[ Z_{\nu} (x) = C_1 J_{\nu} (x) + C_2 Y_{\nu} (x) \]
denotes the general solution of the Bessel equation with arbitrary constants C1 and C2, then the general solution of the given self-adjoint equation becomes
\[ y(x) = x^{\nu /\alpha} Z_{\nu} \left( \alpha \sqrt{|b|} \, x^{1/\alpha} \right) . \]
If b < 0, then we have to choose
\[ Z_{\nu} (x) = C_1 I_{\nu} (x) + C_2 K_{\nu} (x) \]
as the general solution to the modified Bessel equation.

In the case of the general Bessel equation given by Bowman (1958)

\[ x^2 \frac{{\text d}^2 y}{{\text d}x^2} + \left( 2p +1 \right) \frac{{\text d}y}{{\text d}x} + \left( \alpha^2 x^{2r} + \beta^2 \right) y =0 , \]
its solution is
\[ y(x) = x^{-p} \left[ C_1 J_{q/r} \left( \frac{\alpha}{r}\,x^r \right) + C_2 Y_{q/r} \left( \frac{\alpha}{r}\,x^r \right) \right] , \]
where
\[ q = \sqrt{p^2 - \beta^2} , \]
Jn(x) and Yn(x) are the Bessel functions of the first and second kind, and C1 and C2 are arbitrary constants. Another form is given by letting \( y= x^{\alpha} J_n \left( \beta\, x^{\gamma} \right) , \quad \eta = y\,x^{\gamma} , \ \mbox{ and } \ \xi = \beta\,x^{\gamma} \) (see Bowman, 1958, p. 117), then
\[ \frac{{\text d}^2 y}{{\text d}x^2} - \frac{2\alpha -1}{x} \, \frac{{\text d}y}{{\text d}x} + \left( \beta^2 \gamma^2 x^{2\gamma -2} + \frac{\alpha^2 - n^2 \gamma^2}{x^2} \right) y(x) = 0 . \]
The solution is
\[ y(x) = \begin{cases} x^{\alpha} \left[ C_1 J_n \left( \beta\, x^{\gamma} \right) + C_2 Y_n \left( \beta\, x^{\gamma} \right) \right] \quad \mbox{for integer } n , \\ x^{\alpha} \left[ C_1 J_n \left( \beta\, x^{\gamma} \right) + C_2 J_{-n} \left( \beta\, x^{\gamma} \right) \right] \quad \mbox{for noninteger } n , \end{cases} \]

 

Bessel function of first kind


The Bessel function of first kind of order ν is
\[ J_{\nu} (x) = \sum_{k\ge 0} \frac{(-1)^k}{k!\,\Gamma (k+\nu +1)} \left( \frac{x}{2} \right)^{2k+\nu} , \]
Bessel functions of the first kind Jν(x) are bounded at the origin for any ν. When ν is not an integer, Bessel functions Jν(x) and J(x) are linearly independent because their Wronskian is
\[ W\left[ J_{\nu} (x), J_{-\nu} (x) \right] = J_{\nu +1} (x)\,J_{-\nu} + J_{\nu } (x) (x)\,J_{-\nu -1} (x) = -\frac{2\,\sin \left( \nu \pi \right)}{\pi\, x} . \]
When ν = n an integer,
J_{-n} (x) = (-1)^n J_n (x) . \]
plot = Plot[{BesselJ[0, x], BesselJ[1, x], BesselJ[2, x]}, {x, 0.01, 14}, AxesLabel -> {x, Bessel}, PlotStyle -> {Thick}, PlotLegends -> "Expressions"];
left = Graphics[ Text[Style["\!\(\*SubscriptBox[\(J\), \(0\)]\)(x)"], {1.5, 0.95}]];
center = Graphics[ Text[Style["\!\(\*SubscriptBox[\(J\), \(1\)]\)(x)"], {2.2, 0.65}]];
right = Graphics[ Text[Style["\!\(\*SubscriptBox[\(J\), \(2\)]\)(x)"], {4.0, 0.55}]];
Show[plot, left, center, right]
Graphs of Bessel functions of first kind.

Integral Representations

\begin{align*} J_n (x) &= \frac{1}{\pi} \int_0^{\pi} \cos \left( n\theta -x\,\sin\theta \right) {\text d}\theta , \\ J_1 (x) &= \frac{1}{\pi} \int_0^{\pi} \cos \left( \theta - x\,\sin\theta \right) {\text d}\theta = \frac{1}{\pi} \int_0^{\pi} \cos \theta \,\sin \left( x\,\cos\theta \right) {\text d}\theta , \\ J_0 (x) &= \frac{1}{\pi} \int_0^{\pi} \cos \left( x\,\sin\theta \right) {\text d}\theta = \frac{1}{\pi} \int_0^{\pi} \cos \left( x\,\cos\theta \right) {\text d}\theta . \end{align*}

Recurrence Relations

\begin{align*} J_{\nu +1} (x) &= \frac{2\nu}{x}\, J_{\nu} (x) - j_{\nu -1} (x) , \end{align*}

Derivatives

\begin{align*} J'_{\nu +1} (x) &= \frac{1}{2} \left[ J_{\nu -1} (x) - J_{\nu +1} (x) \right] , \\ J'_{\nu} (x) &= J_{\nu -1} (x) - \frac{\nu}{x}\, J_{\nu} (x) , \\ J'_{\nu} (x) &= \frac{\nu}{x}\, J_{\nu} (x) - J_{\nu +1} (x) , \\ \frac{\text d}{{\text d}x} \left[ x^{\nu} J_{\nu} (x) \right] &= x^{\nu} J_{\nu -1} (x) , \\ \frac{\text d}{{\text d}x} \left[ x^{-\nu} J_{\nu} (x) \right] &= -x^{-\nu} J_{\nu +1} (x) . \end{align*}

 

Bessel function of second kind


The Bessel function of the second kind of order ν is defined as follows: \[ Y_{\nu} (x) = \frac{\cos \nu \pi \, J_{\nu} (x) - J_{-\nu} (x)}{\sin \nu \pi} . \]
This function is called either after Heinrich Martin Weber (1875--1913), who introduced Yν(x) in 1873, or after Carl Neumann (1832--1925), with notation Nν(x). Since there is no standard notation, we will distinguish them by defining the Neumann function N(x) = πY(x) as π multiple of the Weber function:
\[ N_{\nu} (x) = \pi\,Y_{\nu} (x) . \]
These functions are clearly solutions of the Bessel equation, as they are linear combinations of Bessel functions of the first kind. However, if ν is an integer, then Y(x), as defined, is the indeterminate form 0/0. Therefore, we need to use l’Hospital’s Rule to determine whether the limit as ν approaches an integer n is nonzero, so that we can obtain a meaningful solution.

Applying l’Hospital’s Rule yields

\[ \begin{split} Y_{\nu} (x) &= \lim_{\nu \to n} \,\frac{\cos \nu \pi \, J_{\nu} (x) - J_{-\nu} (x)}{\sin \nu \pi} \\ &= \lim_{\nu \to n} \,\frac{-\pi\,\sin \nu\pi \,J_{\nu} (x) + \cos \nu \pi \, \frac{\text d}{{\text d}\nu} \, J_{\nu} (x) - \frac{\text d}{{\text d}\nu} \, J_{-\nu} (x)}{\pi\,\cos \nu\pi} \\ &= \lim_{\nu \to n} \,\frac{1}{\pi} \left[ \frac{\text d}{{\text d}\nu} \, J_{\nu} (x) - \frac{\text d}{{\text d}\nu} \, J_{-\nu} (x) \right] . \end{split} \]
Using the series definition of the Bessel function of the first kind,
\[ J_{\nu} (x) = \sum_{k\ge 0} \frac{(-1)^k}{k!\,\Gamma (k+\nu +1)} \left( \frac{x}{2} \right)^{2k+\nu} , \]
as well as the digamma function
\[ \psi (x) = \frac{\text d}{{\text d}x}\, \ln \Gamma (x) = \frac{\Gamma' (x)}{\Gamma (x)} , \]
we obtain
\[ \begin{split} Y_{n} (x) &= \frac{1}{\pi} \left[ \sum_{k\ge 0} \frac{(-1)^k}{k! \, (k+n)!} \left( \frac{x}{2} \right)^{2k+n} \ln \left( \frac{x}{2} \right) + (-1)^n \sum_{k\ge 0}\frac{(-1)^k}{k! \, (k-n)!} \left( \frac{x}{2} \right)^{2k-n} \ln \left( \frac{x}{2} \right) \right] \\ & \quad - \frac{1}{\pi} \left[ \sum_{k\ge 0} \frac{(-1)^k \Gamma' (k+n+1)}{k! \,(k+n)!^2} \left( \frac{x}{2} \right)^{2k+n} + (-1)^n \sum_{k\ge 0} \frac{(-1)^k \Gamma' (k-n+1)}{k! \,(k-n)!^2} \left( \frac{x}{2} \right)^{2k-n} \right] \\ &= \frac{1}{\pi} \left[ J_n (x)\, \ln \left( \frac{x}{2} \right) + (-1)^n J_{-n} (x)\, \ln \left( \frac{x}{2} \right) \right] \\ & \quad - \frac{1}{\pi} \left[ \sum_{k\ge 0} \frac{(-1)^k \psi (k+n+1)}{k!\,(k+n)!} \left( \frac{x}{2} \right)^{2k+n} + (-1)^n \sum_{k\ge 0} \frac{(-1)^k \psi (k+n+1)}{k!\,(k-n)!} \left( \frac{x}{2} \right)^{2k-n} \right] \\ &= \frac{2}{\pi} \, J_n (x)\, \ln \left( \frac{x}{2} \right) \\ & \quad - \frac{1}{\pi} \left[ \sum_{k\ge 0} \frac{(-1)^k \psi (k+n+1)}{k!\,(k+n)!} \left( \frac{x}{2} \right)^{2k+n} + (-1)^n \sum_{k\ge -n} \frac{(-1)^{k+n} \psi (k+1)}{k!\,(k+n)!} \left( \frac{x}{2} \right)^{2k+n} \right] \\ &= \frac{2}{\pi} \, J_n (x)\, \ln \left( \frac{x}{2} \right) - \frac{(-1)^n}{\pi} \,\sum_{k=-n}^{-1} \frac{(-1)^k \psi (k+1)}{k!\,(k+n)!} \left( \frac{x}{2} \right)^{2k+n} \\ & \quad - \frac{1}{\pi} \sum_{k\ge 0} \frac{(-1)^k}{k! \,(k+n)!} \left( \frac{x}{2} \right)^{2k+n} \left[ \psi (k+n+1) + \psi (k+1) \right] \\ &= \frac{2}{\pi} \, J_n (x)\, \ln \left( \frac{x}{2} \right) - \frac{(-1)^n}{\pi} \,\sum_{k=0}^{n-1} \frac{(-1)^k \psi (k-n+1)}{k!\,(k-n)!} \left( \frac{x}{2} \right)^{2k-n} \\ & \quad - \frac{1}{\pi} \sum_{k\ge 0} \frac{(-1)^k}{k! \,(k+n)!} \left( \frac{x}{2} \right)^{2k+n} \left[ \psi (k+n+1) + \psi (k+1) \right] \\ &= \frac{2}{\pi} \, J_n (x)\, \ln \left( \frac{x}{2} \right) - \frac{1}{\pi} \,\sum_{k=0}^{n-1} \frac{(n-k-1)!}{k!} \left( \frac{x}{2} \right)^{2k-n} \\ & \quad - \frac{1}{\pi} \sum_{k\ge 0} \frac{(-1)^k}{k!\, (k+n)!} \left( \frac{x}{2} \right)^{2k+n} \left[ \psi (k+n+1) + \psi (k+1) \right] , \end{split} \]
where we have used the limit
\[ \lim_{z\to -n} \,\frac{\psi (z)}{\Gamma (z)} = (-1)^{n-1} n! \, . \]
We can see from the above expression for Yn(x) that it is indeed linearly independent of Jn(x), so that we have two linearly independent solutions of the Bessel equation for integer order n. Also, unlike Jn(x), the Weber function is singular at x = 0.

Using the identity \( \psi (n) = H_{n-1} - \gamma , \) we have

\begin{align*} Y_0 (x) &= \frac{2}{\pi} \left[ \gamma + \ln \left( \frac{x}{2} \right) \right] J_0 (x) - \frac{2}{\pi} \sum_{k\ge 0} \frac{(-1)^k}{k!\, k!}\left( \frac{x}{2} \right)^{2k} H_k , \\ Y_n (x) &= \frac{2}{\pi} \left[ \gamma + \ln \left( \frac{x}{2} \right) \right] J_n (x) - \frac{1}{\pi} \sum_{k= 0}^{n-1} \frac{(n-k-1)!}{k!}\left( \frac{x}{2} \right)^{2k-n} - \frac{1}{\pi} \sum_{k\ge 0} \frac{(-1)^k }{k!\, (k+n)!}\left( \frac{x}{2} \right)^{2k} \left[ H_k + H_{k+n} \right] , \end{align*}
where γ ≈ 0.5772156649… is the Euler constant (also called the Euler-Mascheroni constant in honor of Lorenzo Mascheroni (1750--1800)), which is implemented in the Wolfram Language as EulerGamma, and Hn is the n-th harmonic number:
\[ H_n = 1 + \frac{1}{2} + \frac{1}{3} + \cdots + \frac{1}{n} = \sum_{k=1}^n \frac{1}{k} . \]
All Weber functions are unbounded at the origin, as it is seen from the plot.
Plot[{BesselY[0, x], BesselY[1, x], BesselY[2, x]}, {x, 0.01, 14}, AxesLabel -> {x, Bessel}, PlotStyle -> {Thick}, PlotLegends -> "Expressions"]
Graphs of Weber functions.
The Wronskian is
\[ W \left[ J_{\nu} (x), Y_{\nu} (x) \right] = J_{\nu +1} (x)\,Y_{\nu} (x) - J_{\nu} (x)\,Y_{\nu +1} (x) = \frac{2}{\pi\,x} . \]

Integral Representations

Bessel functions of the second kind also have integral representations. We have
\begin{align*} Y_0 (x) &= - \frac{2}{\pi} \int_0^{\infty} \cos \left( x\,\cosh t \right) {\text d} t , \\ Y_n (x) &= \frac{1}{\pi} \int_0^{\pi} \sin \left( x\,\sin t - nt \right) {\text d} t - \frac{1}{\pi} \int_0^{\infty} \left[ e^{nt} +(-1)^n e^{-nt} \right] e^{-x\,\sinh t} {\text d} t . \end{align*}

Recurrence Relations

Bessel functions of the second kind, being solutions of the Bessel equation, satisfy the same recurrence relations as the Bessel functions of the first kind. Specifically,
\begin{align*} Y_{\nu -1} (x) - Y_{\nu +1} (x) &= 2\,Y_{\nu} (x) , \\ Y_{\nu -1} (x) + Y_{\nu +1} (x) &= \frac{2\nu}{x}\, Y_{\nu} (x) . \end{align*}
We also have the relation
\[ Y_{-n} (x) = (-1)^n Y_n (x) , \]
when n is an integer.

 

Bessel function of third kind


We can define two new linearly dependent functions, called the Hankel functions
\begin{align*} H_{\nu}^{(1)} (x) &= J_{\nu} (x) + {\bf j}\,Y_{\nu} (x) = \frac{J_{-\nu} (x) - e^{-\nu\pi{\bf j} J_{\nu} (x)} }{{\bf j}\,\sin \nu\pi} , \\ H_{\nu}^{(2)} (x) &= J_{\nu} (x) - {\bf j}\,Y_{\nu} (x) = \frac{J_{-\nu} (x) - e^{\nu\pi{\bf j} J_{\nu} (x)} }{-{\bf j}\,\sin \nu\pi} , \end{align*}
where j is a unit vector in positive vertical direction on the complex plane ℂ so that j² = -1.

 

Modified Bessel equation


 

Kelvin's functions


 

  1. Grigoryan, V, Parial Differential Equations, 2010, University of California, Santa Barbara.
  2. Olver, E.W.J., Maximon, L.C., Chapter 10 Bessel Functions, National Institute of Standards and Tech-nology, 2010.

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions