Skip to main content

Section 2.6 Bessel equations

Subsection 2.6.1 Bessel's differential equation

One of the most important differential equations in applied mathematics and physical science is the Bessel differential equation, which has the form

\begin{equation} x^2 y''+ x y' + (x^2 - \alpha^2) y = 0\label{eq-bessel}\tag{2.6.1} \end{equation}

for some complex constant \(\alpha\text{,}\) which is called the order of the solutions arising from the associated equation. Bessel's equation and the associated solutions are a core technique in the analysis of the propagation of waves, for example.

In general, we have to write the solutions to this equation in terms of integrals; that is, there are not closed form solutions. But since (2.6.1) has a regular singular point at \(0\text{,}\) in the case where \(\alpha\) is a nonnegative real constant, we can use Frobenius techniques to write series solutions.

Since \(p(0) = 1\) and \(q(0) = -\alpha^2\text{,}\) we get the indicial equation

\begin{equation*} r(r-1) + r - \alpha^2, \end{equation*}

which has solutions \(\pm \alpha\text{.}\) As long as \(2\alpha\) is not an integer, we can find two Frobenius series solutions of the form

\begin{equation*} y(x) = x^r \sum_{n=0}^\infty a_n x^n, \end{equation*}

and again it is useful to recall the derivative formulas (2.4.4) and (2.4.5).

We proceed with our standard technique and recover a series solution from a recurrence relation. Substituting into the (2.6.1), we get

\begin{equation*} \sum_{n=0}^\infty [(r + n)^2 - \alpha^2]a_n x^{n+r} + \sum_{n=2}^\infty a_{n-2} x^{n+r} = 0. \end{equation*}

We can peel off the first two terms of the left sum. For \(n=0\text{,}\) we get the indicial equation. For \(n = 1\text{,}\) we get

\begin{equation*} [(r + 1)^2 - \alpha^2]a_1 = 0. \end{equation*}

For \(n \geq 2\text{,}\) we get the recurrence relation

\begin{equation*} [(n+r)^2 - \alpha^2]a_n = a_{n -2}. \end{equation*}

To construct the first solution, consider the root \(r = \alpha\text{.}\) The recurrence relation becomes

\begin{equation*} a_n = -\frac{1}{n(2\alpha + n)} a_{n-2}. \end{equation*}

Since the gap between the terms in the relation is two, there will be two families of coefficients. Since \(a_1 = 0\text{,}\) it must be that \(a_3, a_5, \ldots = 0\text{.}\) On the other hand, with a little work we get that the even terms are characterized by

\begin{equation*} a_{2k} = \frac{(-1)^k}{2 \cdot 4 \ldots \cdot (2k) \cdot (2\alpha + 2) \ldots \cdot (2 \alpha + 2k)} a_0, \end{equation*}

which is a bit of a bear. Next, we'll try to reduce this into a more palatable form. First, we can pull out a whole pile of 2s to get something a bit more condensed, and a form that might help guide a useful choice for \(a_0\text{.}\)

\begin{equation*} a_{2k} = \frac{(-1)^k}{2^{2k} k! (\alpha + 1) \cdots (\alpha + k)} a_0, k = 1, 2, \ldots. \end{equation*}

Then we have a Frobenius series solution to (2.6.1) of the form

\begin{equation} y_1(x) = a_0 x^\alpha \left[1 + \frac{(-1)^k}{2^{2k} k! (\alpha + 1) \cdots (\alpha + k)} x^{2k}\right].\label{eq-bessel-series}\tag{2.6.2} \end{equation}

Subsection 2.6.2 Bessel functions of the first kind

When \(\alpha = N\) is an integer, the formula collapses even further if we make a clever choice for the arbitrary constant \(a_0\text{.}\) Let

\begin{equation*} a_0 = \frac{!}{N!2^N}. \end{equation*}

Then the solution is denoted \(J_N(x)\) and is called the Bessel function of the first kind of order \(N\) and has the form

\begin{equation*} J_N(x) = \sum_{k=0}^\infty \frac{(-1)^k}{k!(N + k)!} \left(\frac{x}{2}\right)^{2k + N}. \end{equation*}

The following code uses prebuilt libraries to construct and plot the first few examples of \(J_N\)

The graphs above should suggest why these functions are useful at modeling wave propagation - each function oscillates with decreasing amplitude and has an infinite number of zeros on \((0,\infty)\text{.}\)

To deal with the case where \(\alpha \in (0,\infty)\) is not an integer, we need to introduce a generalization of the factorial function called the gamma function, one of the most famous special functions in mathematics. The basic idea of the gamma function is to “fill in the spaces” between the factorials with a nice, smooth function that retains the special properties of the factorial function (that is, a nice function \(y\) for which \(y(x) = x y(x-1)\)). It turns out that the best way to do that is with a function defined in terms of an integral.

Definition 2.6.1.

The gamma function \(\Gamma: (0, \infty) \to \R\) is defined by the integral formula

\begin{equation} \Gamma(x) = \int_0^\infty t^{x-1} e^{-t} \, dt.\label{eq-gamma-prop}\tag{2.6.3} \end{equation}

The gamma function is actually defined on all complex numbers outside of the negative integers, but we do not need that fact here.

This is an immediate consequence of integration by parts.

\begin{align*} \Gamma(x + 1) \amp= \int_0^\infty t^{x} e^{-t} \, dt\\ \amp= \left[-t^{x} e^{-t}\right]_0^\infty + x \int_0^\infty t^{x - 1} e^{-t} \, dt\\ \amp= x \Gamma(x) \end{align*}

If we had a starting point, we could show that the gamma function generates the factorials. Fortunately, it is very easy to show that

\begin{equation*} \Gamma(1) = \int_0^\infty e^{-t} \, dt = 1. \end{equation*}

Then induction on (2.6.3) gets us

\begin{equation*} \Gamma(2) = 2\Gamma(1) = 2!, \,\,\, \Gamma(3) = 3\Gamma(2) = 3!, \end{equation*}

and in general that for a positive integer \(N\) that

\begin{equation*} \Gamma(N) = N!. \end{equation*}

At the same time, we can get the factorial-like formula

\begin{equation} \Gamma(x + k + 1) = [(x+1)(x+2)\cdots(x+k)] \Gamma(x).\label{eq-gamma-factorial-like}\tag{2.6.4} \end{equation}

Now we can extend the gamma function to negative numbers that aren't integers by using the relation

\begin{equation} \Gamma(x) = \frac{\Gamma(x+1)}{x}\label{eq-gamma-extension}\tag{2.6.5} \end{equation}

If \(x \in (-1,0)\text{,}\) then \(x + 1 \in (0,1)\text{,}\) and so for \(x \in (-1,0)\) we use (2.6.5) as the definition of \(\Gamma(x)\text{.}\) We can continue in this way inductively: if \(x \in (-2,-1)\text{,}\) then \(x+ 1 \in (-1,0)\text{,}\) which is now defined, and again we use formula (2.6.5). The resulting function is continuous off the negative integers and has the graph below.

The gamma function allows us to define the Bessel functions for non-integer orders. First, choose

\begin{equation*} a_0 = \frac{1}{2^\alpha \Gamma(\alpha + 1)} \end{equation*}

and plug into (2.6.2) to get

\begin{equation} J_\alpha(x) = \sum_{k=0}^\infty \frac{(-1)^k}{\Gamma(k+1) \, \Gamma(\alpha + k + 1)} \left(\frac{x}{2}\right)^{2k + \alpha},\label{eq-bessel-first}\tag{2.6.6} \end{equation}

which by the definition of \(\Gamma\) reduces to \(J_N(x)\) when \(N\) is an integer.

Subsection 2.6.3 Bessel functions of the second kind

Because the differential equation (2.6.1) has lots of solutions, there are many ways to write down a pair of linearly independent functions that solve the equation. The function that we will consider here, we choose for convenience by convention of physicists and engineers, but there are several equivalent ways of formulating the general solution. First, note that as long as \(2\alpha\) isn't an integer, the second root to the indicial equation, \(r = -\alpha\) will also have an associated and linearly independent solution. It is straightforward that we can just replace \(\alpha\) with \(-\alpha\) in (2.6.6) to get

\begin{equation*} J_{-\alpha} = \sum_{k=0}^\infty \frac{(-1)^k}{\Gamma(k + 1) \, \Gamma(k - \alpha + 1)} \left(\frac{x}{2}\right)^{2k - \alpha} \end{equation*}

and a general solution

\begin{equation*} y(x) = c_1 J_\alpha(x) + c_2 J_{-\alpha}(x). \end{equation*}

In practice, we don't use \(J_{-\alpha}\) as the second solution, but rather a linear combination of \(J_{\alpha}\) and \(J_{-\alpha}\text{.}\) Thus we define the Bessel function of the second kind of order \(\alpha\) by

\begin{equation*} Y_\alpha(x) = \frac{J_\alpha(x) \cos \alpha \pi - J_{-\alpha}(x)}{\sin \alpha \pi}. \end{equation*}

The function is construced this way because it extends naturally to the case where \(\alpha\) is an integer, avoiding the complications of the reduction of order technique needed to derive the second solution from the first in the \(\alpha\) is an integer case. Without doing the computations, it is enough to know that if \(\alpha = n\text{,}\) then we can define

\begin{equation*} Y_n(x) = \lim_{\alpha \to n} \left[\frac{J_\alpha(x) \cos \alpha \pi - J_{-\alpha}(x)}{\sin \alpha \pi}\right] \end{equation*}

and that this sum can be computed.

Thus the general solution to Bessel's equation is

\begin{equation*} y(x) = c_1 J_\alpha(x) + c_2 Y_\alpha(x) \end{equation*}

for \(\alpha \in (0,\infty)\text{.}\)

Notice that the \(Y_\alpha\) are unbounded at \(x = 0\text{,}\) which is typically not a desired characteristic in a physical problem.

Subsection 2.6.4 Bessel-Fourier series

We've already seen that Legendre polynomials from an orthgonal basis for the \(C^1\) functions on a finite interval with respect to the inner product

\begin{equation*} \ip{f}{g} = \int f(x) g(x) \, dx. \end{equation*}

It turns out that Bessel functions can be used to form a different orthogonal system of functions that serve as a basis for the \(C^1\) functions (one of of an infinite family), but with a slightly different inner product. In fact, just one function \(J_\alpha\) can be used to generate the basis (where the choice of \(\alpha\) lets us describe functions in terms of a solution to a particular physical scenario).

The graphs of \(J_\alpha\) should be evidence that for a given \(\alpha\text{,}\) \(J_\alpha\) has infinite zeros (as it wiggles back and forth across the \(x\)-axis). For a fixed \(\alpha\text{,}\) let \(\la_{\alpha, n}\) denote the sequence of positive zeros of \(J_\alpha\text{.}\)

To construct an orthonormal family, we'll need to use a modification of the standard inner product by introducing a weight function. One way of thinking about the standard inner product \(x \cdot y = \sum x_i y_i\) or \(\ip{f}{g} = \int fg\) is that we give all terms equal importance. But there's no reason that we have to do that. A weight function \(w(x) \geq 0\) essentially turns an inner product into a weighted average instead of a standard average. On \(\R^n\text{,}\) a weighted inner product has the form \(x \cdot y = \sum x_i y_i w(i)\text{.}\) On function spaces, we have the weighted inner product with respect to \(w\) which is given by

\begin{equation*} \ip{f}{g} = \int_0^1 f(x)g(x)w(x)\, dx. \end{equation*}

We do this because we change the geomerty of the space when we use a new inner product. New vectors will be orthogonal. We will choose the weight function \(w(x) = x\text{,}\) and we get the following fact about the family \(J_\alpha(\la_{\alpha,n}x)\text{.}\)

What we're really saying here is that for any \(n\neq m\) that

\begin{equation*} \ip{J_\alpha(\la_{\alpha,n}x)}{J_\alpha(\la_{\alpha,m}x)} = \int_0^1 x J_\alpha(\la_{\alpha,n}) J_\alpha(\la_{\alpha,m}x) \, dx = 0. \end{equation*}

We can also show (with quite a bit of elbow grease) that

\begin{equation*} \norm{J_\alpha(\la_{\alpha,n}x)} = \frac{1}{2} \left[ J_{\alpha + 1}(\la_{\alpha,n})\right]^2. \end{equation*}

Even better, the family \(J_\alpha(\la_{\alpha,n}x)\) turns out to be a basis for \(C^1\text{,}\) and so we can express any \(C^1\) function in terms of a series of Bessel functions.

Again, the proof of convergence of these types of series is beyond the scope of the course and more suited to a course in real analysis or partial differential equations. We keep emphasizing the vector notation not just to impress with our ability to assemble huge blocks of mathematics, but to point out that all we're really doing here is linear algebra with fancier vectors.