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

Glossary

Preface

Although Johann Bernoulli (1667--1748) came in 1694 to an equation that is a particular case of what we call now the Bessel equation, it was his son who actually introduced it. Daniel Bernoulli (1700--1782) is generally credited with being the first to introduce the concept of Bessels functions in 1732 (Saint Petersburg). He used the function of zero order J0(x) as a solution to the problem of an oscillating chain suspended at one end. He also discovered that the equation J0(x) = 0 has infinuitely many zeroes. In 1732 (Saint Petersburg), Leonhard Euler (1707--1783) employed Bessel functions of both zero and integral orders in an analysis of vibrations of a stretched membrane. He also found Maclaurin series for Jn(x) with integer values of n. Also, Leonhard showed that Bessel's function of half integer order can be expressed through elementary functions.

Bessel functions were used by Lagrange in 1770, in the theory of planetary motion, by Fourier in his theory of heat flow (1822), by Poisson in the theory of heat flow in spherical bodies (1823), and by Bessel, who studied these functions in detail around 1824. It was Lord Rayleigh who demonstrated that Bessel's functions are particular cases of Laplace's functions in 1878. Watson (1966) provided the most comprehensive study of the Bessel functions. See Dutka's historical notes.

This section is about Bessel's equation and its solutions, known as Bessel functions or cylinder functions. The origin of the term cylinder is due to the fact that these functions are encountered in studying the boundary–value problems of potential theory for cylindrical coordi-nates. The Bessel differential equation is the linear second-order ordinary differential equation given by

\begin{equation} \label{EqBessel.1} 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 in self-adjoint form}\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 , \end{equation}
where ν is a real constant, called the order of the Bessel equation. Eq.\eqref{EqBessel.1} has a regular singularity at x = 0. 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 . \]

The systematic analysis of solutions to equation \eqref{EqBessel.1} was conducted around 1817 by the German astronomer Friedrich Wilhelm Bessel (1784--1846) during an investigation of solutions of one of Kepler’s equations of planetary motion. The Bessel function was the result of Bessel's 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.

Bessel Equation

Bessel's equation \eqref{EqBessel.1} 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. A generalized series for bounded solution of Bessel's equation was found in section of tutorial I.

It is a custom to separate bounded solutions from unbounded solutions of Bessel's equation. A bounded solution is usually denoted by Jν(x) and called the Bessel function of the first kind or just Bessel's function. All other bounded solutions of Eq.\eqref{EqBessel.1} are constant multiples of this function. On the other hand, an unbounded solution is refferred to as the Bessel function of the second kind. Their complex combination is called the Bessel function of the third kind or "Hunkel's function. Bessel functions of complex argument are referred to as the Kelvin functions.

 

Bessel function of first kind

The Bessel function of first kind of order ν of real argument x is usually denoted by Jν(x). Its generalized series expansion was found previously (see section of tutorial I):
\begin{equation} \label{EqBessel.2} J_{\nu} (x) = \sum_{k\ge 0} \frac{(-1)^k}{k!\,\Gamma (k+\nu +1)} \left( \frac{x}{2} \right)^{2k+\nu} , \end{equation}

The notation Jz,n was first used by a Danish-born German astronomer Peter Hansen (1795--1874) in 1843 and subsequently by Oskar Xavier Schlömilch in 1857 and later modified to Jn(2z) by Watson (1922). Subsequent studies of Bessel functions included the works of Mathews in 1895, “A treatise on Bessel functions and their applications to physics” written in collaboration with Andrew Gray. It was the first major treatise on Bessel functions in English and covered topics such as applications of Bessel functions to electricity, hydrodynamics and diffraction. In 1922, Watson first published his comprehensive examination of Bessel functions “A Treatise on the Theory of Bessel Functions”.

Since Bessel's equation \eqref{EqBessel.1} has a singular point at the origin, initial conditions are not specified for it because they do not identify a unique solution of the equation and cannot be chosen arbitrary. Nevertheless, they make sense for a bounded solution although the initial conditions for the function y(x) and its derivative at x = 0 are coupled (meaning that one of them forced another one to have a specific value). Then to find a bounded solution of Eq.\eqref{EqBessel.1} we apply the Laplace transformation:

\[ \frac{{\text d}^2}{{\text d}\lambda^2} \left[ y^L - \lambda y(0) - y'(0) \right] + \frac{{\text d}}{{\text d}\lambda} \left[ y^L - y(0) \right] + \frac{{\text d}^2}{{\text d}\lambda^2} \, y^L - \nu^2 y^L =0 , \]
where
\[ y^L (\lambda ) = {\cal L} \left[ y(t) \right] (\lambda ) = \int_0^{\infty} e^{-\lambda t} y(t)\,{\text d}t \]
is the Laplace transform of unknown solution to Bessel's equation. In its derivation we used the property of the Laplace transformation:
\[ {\cal L} \left[ t\, f(t) \right] = \frac{{\text d}}{{\text d}\lambda} \, f^L \qquad \mbox{and} \qquad {\cal L} \left[ t^2 f(t) \right] = \frac{{\text d}^2}{{\text d}\lambda^2} \, f^L . \]
This yields the following differential equation for \( Y(\lambda ) = y^L : \)
\[ \left( \lambda^2 + 1 \right) y'' + 3\lambda\, Y' + \left( 1 - \nu^2 \right) Y = 0 . \]
To solve this equation, we make the change of variables:
\[ \lambda = \sinh q , \qquad Y(\lambda ) = \frac{1}{\cosh q}\, Z (q) . \]
Then we will get
\[ Z'' - \nu^2 Z = 0 \qquad \Longrightarrow \qquad Z = e^{-\nu\,q} . \]
Upon returning to variables Y and λ, we obtain the explicit formula for the Laplace transform:
\[ y^L = Y = \frac{1}{\sqrt{\lambda^2 +1}} \, e^{-\nu\,\mbox{arcsinh} \lambda} = \frac{1}{\sqrt{\lambda^2 +1}} \cdot \frac{1}{\left( \lambda + \sqrt{\lambda^2 +1} \right)^{\nu}} \]
Using notation \( \omega = \lambda + \sqrt{\lambda^2 +1} , \) we apply the inverse Laplace transform to determine a bounded solution to Bessel's equation:
\[ y(t) = J_{\nu} (x) = {\cal L}^{-1} \left[ Y (\lambda ) \right] (t) = \mbox{P.V.}\,\frac{1}{2\pi {\bf j}} \int_{h- {\bf j}\infty}^{h+ {\bf j}\infty} \frac{1}{\omega^{\nu +1}} \, e^{t \left( \omega - \omega^{-1} \right) /2} {\text d}\omega , \]
where P.V. stands for Principle Value regularization and j is the unit vector in positive verstical direction on the complex plane, so j² = −1. This integral is named after Russian mathematician Nikolay Sonin (1849--1915). Upon setting \( \omega = e^{{\bf j}\zeta} , \) the Sonin integral is expressed as
\[ J_{\nu} (x) = \frac{1}{2\pi} \int_{\Pi} e^{{\bf j} x\,\sin \zeta - {\bf j}\nu \zeta} {\text d}\zeta . \]
Using formal expansion of the integrand into Maclaurin series, we transfer the Sonin integral into the following sum:
\[ J_{\nu} (x) = \sum_{n\ge 0} \frac{(-1)^n}{n!} \left( \frac{x}{2} \right)^{2n+\nu} \frac{1}{2\pi{\bf j}} \int_C e^{\zeta} \zeta^{-\nu -1 - n} {\text d}\zeta \]
This leads to series representation \eqref{EqBessel.2} of Bessel's function of the first kind.

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} . \]
Also
\[ J'_{\nu} (x)\, J_{-\nu} (x) - J'_{-\nu} (x)\, J_{\nu} (x) = \frac{2\,\sin \nu \pi}{\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*}

 

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} \]

============================================ to be modified ===============

Bessel functions

In 1873, Hermite used a sequence of polynomials \( \{ y_n (x) \}_{n\ge 0} \) in his proof of the transcendency of e. Also, in 1880, a German physicist Heinrich Rudolf Hertz (1857--1894) essentially introduced a Bessel polynomial of degree n:

\[ y_n (x) = \left( \frac{\pi x}{2} \right)^{-1/2} e^{1/x} \,K_{n + 1/2} \left( \frac{1}{x} \right) = \sum_{k=0}^n \frac{(n+k)!}{k! \,(n-k)!}\left( \frac{x}{2} \right)^k , \]
where \( K_{n+ 1/2} (x) \) is Macdonald's function. The Bessel polynomials satisfy the differential equation
\[ x^2 y'' + 2 \left( x+1 \right) y' -n \left( n+1 \right) y =0 . \]
From the formula
\[ \frac{\text d}{{\text d}x} \left( x^2 e^{-2/x}\left( y'_n y_m - y_m y'_m \right) \right) = \left( n-m \right) \left( n_m+1 \right) e^{-2/x} \,y_n y_m , \]
it follows that there is no ordinary function ψ that satisfy the orthogonality relation
\[ \int_{-\infty}^{\infty} y_n (x)\,y_m (x)\,{\text d} \psi (x) = K_n \delta_{n,m} , \quad k_n \ne 0, \quad n,m = 0,1,2,\ldots , \]
but the generating function (distribution):
\[ \frac{{\text d}\psi}{{\text d}x} = w(x) , \quad \mbox{where} \quad w(x) = - \sum_{n\ge 0} \frac{2^{n+1} \delta^{(n)} (x)}{n! \,(n+1)!} . \]

If we are looking for bounded solutions on the finite interval \( (0, \ell ) , \) we have to disregard the Weber/Neumann function (by setting the corresponding coefficient to zero) because these functions are unbounded at the origin. Therefore, its solution will be \( J_{\nu} (\alpha x) . \)

There are known three typical types of boundary conditions:

  • Dirichlet boundary conditions: \( y(\ell ) =0 \qquad \Longrightarrow \qquad J_{\nu} (\alpha \ell ) =0. \)
  • Neumann boundary conditions: \( y'(\ell ) =0 \qquad \Longrightarrow \qquad J_{\nu}' (\alpha \ell ) =0. \)
  • Boundary conditions of third kind: \( h\,y(\ell ) + y' (\ell )=0 \qquad \Longrightarrow \qquad h\,J_{\nu} (\alpha \ell ) + \alpha \,J_{\nu}' (\alpha \ell ) =0. \)

 

Bessel--Fourier Series

The Fourier--Bessel series of a function f(x) defined in the finite interval \( (0, \ell ) \) is its representation as infinite series

\[ f (x) = \sum_{n \ge 1} c_n J_{\nu} \left( \mu_{\nu , n} \,\frac{x}{\ell} \right) , \]
where \( \mu_{\nu , n} \) (which we will denote by \( \alpha_{n} \) dropping index ν) is a root, numbered n associated with the Bessel function \( J_{\mu} \) and cn are the assigned coefficients associated with one of three types boundary conditions. Let αn (\( n=1,2, \ldots ) \) be a sequence of roots of one of the types boundary conditions (Dirichlet, Neumann, or third kind), then \( \lambda_n = \alpha_n^2 \) becomes the eigenvalue of the corresponding Sturm--Liouville problem, with eigenfunction \( J_{\nu} \left( \alpha_n \,\frac{x}{\ell} \right) . \)
  1. Dirichlet boundary conditions: \( y(\ell ) =0 \qquad \Longrightarrow \qquad J_{\nu} (\alpha \ell ) =0. \)
    \[ c_n = \frac{2}{\ell^2 J^2_{\nu +1} (\alpha_n \ell )} \, \int_0^{\ell} x\,f(x)\, J_{\nu} (\alpha_n x)\,{\text d}x, \qquad n=1,2,\ldots . \]
  2. Neumann boundary conditions: \( y'(\ell ) =0 \qquad \Longrightarrow \qquad J_{\nu}' (\alpha \ell ) =0. \)
    \[ c_n = \frac{2}{\ell^2 J^2_{\nu} (\alpha_n \ell )} \, \int_0^{\ell} x\,f(x)\, J_{\nu} (\alpha_n x)\,{\text d}x, \qquad n=1,2,\ldots . \]
  3. Boundary conditions of third kind: \( h\,y(\ell ) + y' (\ell )=0 \qquad \Longrightarrow \qquad h\,J_{\nu} (\alpha \ell ) + \alpha \,J_{\nu}' (\alpha \ell ) =0. \)
    \[ c_n = \frac{2\alpha_n^2}{\left( \ell^2 \alpha_n^2 -\nu^2 +h \right) J^2_{\nu} (\alpha_n \ell )} \, \int_0^{\ell} x\,f(x)\, J_{\nu} (\alpha_n x)\,{\text d}x, \qquad n=1,2,\ldots . \]
If the Fourier series of f converges at the point x, then the Fourier--Bessel series of f converges to the same value at the same point.

 

  2. Dutka, J., On the early history of Bessel functions, Archive for History of Exact Sciences, volume 49, pages 105–134 (1995). https://doi.org/10.1007/BF00376544
  4. Olver, E.W.J., Maximon, L.C., Chapter 10 Bessel Functions, National Institute of Standards and Tech-nology, 2010.
  8. Watson, G.N., A Treatise on the Theory of Bessel Functions, Cambridge University Press; 2nd edition (August 1, 1995). ISBN-13 ‏ : ‎ 978-0521483919

 

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