Rewriting Bessel Equation into Sturm-Liouville Form
The Bessel equation can be made into the Sturm-Liouville form but in a rather unorthodox fashion. Recall (1) of the last tutorial, and now make the change of variable \( x \rightarrow \alpha x \), and hence \( y’ \rightarrow y’/\alpha, y^{\prime\prime} \rightarrow y^{\prime\prime}/\alpha^2 \) (or more explicitly, \( \frac{dy}{dx} \rightarrow \frac{dy}{d(\alpha x)} = \frac{1}{\alpha}\frac{dy}{dx}, \frac{d^2y}{dx^2} \rightarrow \frac{1}{\alpha^2}\frac{d^2y}{dx^2} \)), so it becomes
\begin{align}
(\alpha x)^2(\frac{1}{\alpha^2}\frac{d^2y}{dx^2}) + \alpha x (\frac{1}{\alpha}\frac{dy}{dx}) + ((\alpha x)^2 -\nu^2)y &= 0 \\
x^2\frac{d^2y}{dx^2} + x\frac{dy}{dx} + (\alpha^2 x^2 -\nu^2)y &= 0 \\
x\frac{d^2y}{dx^2} + \frac{dy}{dx} + (\alpha^2x -\frac{\nu^2}{x})y &= 0 \\
\frac{d}{dx}(x\frac{dy}{dx}) + (\alpha^2x -\frac{\nu^2}{x})y &= 0 \tag{1}
\end{align}
The corresponding solution will now be \( y(x) = J_\nu(\alpha x) \). Following the notation of this tutorial, \( p(x) = x, q(x) = -\nu^2/x, \lambda = \alpha^2 , \rho = x \). Again, by the general results derived for the Sturm-Liouville equation in this tutorial, we can immediately deduce some facts, albeit in a different way compared to those special polynomials before. Notice that the eigenfunctions that satisfy the suitable boundary conditions (as seen below) will be the Bessel function of the first kind \( J_\nu(\alpha x) \) with \( \alpha \) (not \(\nu\)!) being the eigenvalue.
Orthogonality
The orthogonality property of Bessel functions then follows, which is
\begin{align}
\int_a^b x J_\nu(\alpha x) J_\nu(\beta x) dx &= 0 & & \alpha \neq \beta \tag{2}
\end{align}
with the weighting \(\rho = x \). Again, notice that the difference lies in \( \alpha \) and \( \beta \), not \(\nu\). Since \(p(x) = x\) does not suggest strongly any natural interval, we will need \( J_\nu(\alpha x) \) and \( J_\nu(\beta x) \) themselves to obey some boundary conditions. To show this, we can start by putting (1) into:
\begin{align}
x^2\frac{d^2J_\nu(\alpha x)}{dx^2} + x\frac{dJ_\nu(\alpha x)}{dx} + (\alpha^2 x^2 -\nu^2)J_\nu(\alpha x) &= 0 \\
\alpha^2 x^2\frac{d^2J_\nu(\alpha x)}{d(\alpha x)^2} + \alpha x\frac{dJ_\nu(\alpha x)}{d(\alpha x)} + (\alpha^2 x^2 -\nu^2)J_\nu(\alpha x) &= 0 \\
\alpha^2 x^2 J_\nu(\alpha x)^{\prime\prime} + \alpha x J_\nu(\alpha x)’ + (\alpha^2 x^2 -\nu^2)J_\nu(\alpha x) &= 0 \tag{3}
\end{align}
and similarly
\begin{align}
\beta^2 x^2 J_\nu(\beta x)^{\prime\prime} + \beta x J_\nu(\beta x)’ + (\beta^2 x^2 -\nu^2)J_\nu(\beta x) &= 0 \tag{4}
\end{align}
Multiplying the first one by \( J_\nu(\beta x) \) and the second one by \( J_\nu(\alpha x) \) then subtracting gives
\begin{align}
\begin{aligned} x^2(\alpha^2 J_\nu(\alpha x)^{\prime\prime}J_\nu(\beta x) -\beta^2 J_\nu(\beta x)^{\prime\prime}J_\nu(\alpha x) ) \\ + x(\alpha J_\nu(\alpha x)’J_\nu(\beta x)-\beta J_\nu(\beta x)’J_\nu(\alpha x)) \\
+(\alpha^2 x^2 -\beta^2 x^2)J_\nu(\alpha x)J_\nu(\beta x) \end{aligned} &= 0 \\
\begin{aligned} x(\alpha \frac{d}{dx}(J_\nu(\alpha x)^{\prime})J_\nu(\beta x) -\beta \frac{d}{dx}(J_\nu(\beta x)^{\prime})J_\nu(\alpha x) ) \\ + (\alpha J_\nu(\alpha x)’J_\nu(\beta x)-\beta J_\nu(\beta x)’J_\nu(\alpha x)) \\
+x(\alpha^2 -\beta^2)J_\nu(\alpha x)J_\nu(\beta x) \end{aligned} &= 0 \\
\begin{aligned} \frac{d}{dx}(x(\alpha J_\nu(\alpha x)’J_\nu(\beta x)- \beta J_\nu(\alpha x)J_\nu(\beta x)’)) \\ + x(\alpha^2-\beta^2)J_\nu(\alpha x)J_\nu(\beta x) \end{aligned} &= 0 \tag{5}
\end{align}
Integrating both sides over certain interval \( [a,b] \), we have
\begin{align}
&\quad -[x(\alpha J_\nu(\alpha x)’J_\nu(\beta x)-\beta J_\nu(\alpha x)J_\nu(\beta x)’)]_a^b \\
&=\int_a^bx(\alpha^2-\beta^2)J_\nu(\alpha x)J_\nu(\beta x)dx \tag{6}
\end{align}
For (2) to hold, the boundary condition then clearly requires that
\begin{align}
[\alpha x J_\nu(\alpha x)’J_\nu(\beta x)- \beta x J_\nu(\alpha x)J_\nu(\beta x)’]_a^b = 0 \tag{7}
\end{align}
So, common possibilities include that both \( J_\nu(\alpha x) \) and \( J_\nu(\beta x) \) vanish at \(x = a,b\), or both of their derivatives \( J_\nu(\alpha x)’ \) and \( J_\nu(\beta x)’ \) vanish at these end-points instead. Note that if \( a = 0 \), then the lower portion of (7) automatically equals \( 0 \) due to \( p(x) = x \) and we just need to enforce the boundary condition for \( x = b \). For convenience, from now on we may take \( a = 0, b = 1 \). (For other values, we just need to make a linear change of variable.)
Normalization Factor
The normalization factor for the Bessel functions is
\begin{align}
\int_a^b xJ_\nu(\alpha x)^2 dx = \frac{1}{2}[(x^2 -\frac{\nu^2}{\alpha^2})J_\nu(\alpha x)^2 + x^2[J_\nu(\alpha x)’]^2 ]_a^b \tag{8}
\end{align}
This can be derived from integration by parts:
\begin{align}
&\quad \int_a^b xJ_\nu(\alpha x)^2 dx \\
&= \frac{1}{\alpha^2} \int_a^b \alpha x J_\nu(\alpha x)^2 d(\alpha x) \\
&= \frac{1}{\alpha^2} ([\frac{1}{2}\alpha^2 x^2 J_\nu(\alpha x)^2]_a^b -\int_a^b \frac{1}{2}\alpha^2 x^2 d(J_\nu(\alpha x)^2)) \\
&= \frac{1}{\alpha^2} ([\frac{1}{2}\alpha^2 x^2 J_\nu(\alpha x)^2]_a^b -\int_a^b \alpha^2 x^2 J_\nu(\alpha x)\frac{dJ_\nu(\alpha x)}{dx} dx) \tag{9}
\end{align}
Plugging in (3) into the last integral of (9):
\begin{align}
&\quad \int_a^b [(-x^2\frac{d^2J_\nu(\alpha x)}{dx^2} -x\frac{dJ_\nu(\alpha x)}{dx} + \nu^2 J_\nu(\alpha x)) \frac{dJ_\nu(\alpha x)}{dx}] dx \\
&= \int_a^b [-\frac{1}{2} \frac{d}{dx}(x^2[\frac{dJ_{\nu}(\alpha x)}{dx}]^2) + \frac{1}{2} \nu^2 \frac{d}{dx} [J_{\nu}(\alpha x)]^2] dx \tag{10}
\end{align}
So (9) becomes
\begin{align}
&\quad \frac{1}{\alpha^2} \left([\frac{1}{2}\alpha^2 x^2 J_\nu(\alpha x)^2]_a^b -\int_a^b -\frac{1}{2} \frac{d}{dx}(x^2[\frac{dJ_{\nu}(\alpha x)}{dx}]^2) + \frac{1}{2} \nu^2 \frac{d}{dx} [J_{\nu}(\alpha x)]^2 dx\right) \\
&= [\frac{1}{2}x^2 J_\nu(\alpha x)^2 + \frac{1}{2}x^2[J_{\nu}(\alpha x)’^2] -\frac{1}{2}\frac{\nu^2}{\alpha^2} J_\nu(\alpha x)^2]_a^b \tag{11}
\end{align}
Recurrence Relations
The major recurrence relations for Bessel Functions (of the first kind) are the following two:
\begin{align}
\begin{aligned}
\frac{d}{dx}(x^\nu J_\nu(x)) &= x^\nu J_{\nu-1}(x) \\
\frac{d}{dx}(x^{-\nu} J_\nu(x)) &= -x^{-\nu} J_{\nu+1}(x)
\end{aligned} \tag{12}
\end{align}
It can be shown by a direct evaluation of the derivatives and here we briefly do that for the first identity. Recall the form of Bessel functions, (8) of the last tutorial:
\begin{align}
\frac{d}{dx}(x^\nu J_\nu(x)) &= \frac{d}{dx}(\sum_{m=0}^{\infty} \frac{(-1)^m x^{2\nu+2m}}{2^{\nu+2m}m!\Gamma(m+1+\nu)}) \\
&= \sum_{m=0}^{\infty} \frac{(-1)^m 2(\nu+m) x^{2\nu+2m-1}}{2^{\nu+2m}m!\Gamma(m+1+\nu)} \\
&= \sum_{m=0}^{\infty} \frac{(-1)^m x^{2\nu+2m-1}}{2^{\nu+2m-1}m!\Gamma(m+\nu)} \\
&= x^\nu \sum_{m=0}^{\infty} \frac{(-1)^m x^{(\nu-1)+2m}}{2^{(\nu-1)+2m}m!\Gamma(m+1+(\nu-1))} \\
&= x^\nu J_{\nu-1}(x) \tag{13}
\end{align}
(assuming that the differentiation can be put into the summation) The proof for the second one is very similar. By expanding the derivatives on the left for (12), we further have
\begin{align}
x^\nu J_\nu(x)’ + \nu x^{\nu-1} J_\nu(x) &= x^\nu J_{\nu-1}(x) \\
x J_\nu(x)’ + \nu J_\nu(x) &= x J_{\nu-1}(x) \tag{14}
\end{align}
And similarly, we have
\begin{align}
x J_\nu(x)’ -\nu J_\nu(x) &= -x J_{\nu+1}(x) \tag{15}
\end{align}
Adding up (14) and (15) then gives
\begin{align}
2xJ_\nu(x)’ &= x J_{\nu-1}(x) -x J_{\nu+1}(x) \\
2J_\nu(x)’ &= J_{\nu-1}(x) -J_{\nu+1}(x) \tag{16}
\end{align}
while a subtraction between them leads to
\begin{align}
2\nu J_\nu(x) &= x J_{\nu-1}(x) + x J_{\nu+1}(x) \\
\frac{2\nu}{x} J_\nu(x) &= J_{\nu-1}(x) + J_{\nu+1}(x) \tag{17}
\end{align}
It is also not difficult to show that the Bessel functions of the second kind \( Y_\nu(x) \) share exactly the same recurrence relations (12)-(17). By the way, (12) also reveals something special about any pair of Bessel functions with \( \nu \) differed by \(1\). Without loss of generality, let’s say they are \( J_\nu(x) \) and \( J_{\nu-1}(x) \). Between any two consecutive zeros of \( J_\nu(x) \) and hence \( x^\nu J_\nu(x) \), by Rolle’s Theorem there will be a point where the derivative \( d/dx(x^\nu J_\nu(x)) = x^\nu J_{\nu-1}(x) \) attains the value of \( 0 \), and thus it is a zero of \( J_{\nu-1}(x) \). The same argument shows that between any two consecutive zeros of \( J_{\nu-1}(x) \) there will be a zero of \( J_\nu(x) \), and we can conclude that the zeros of \( J_\nu(x) \) and \( J_{\nu-1}(x) \) appear alternately.
Generating Function
The generating function of Bessel functions of the first kind with integer order \( n \) is
\begin{align}
\exp(\frac{x}{2} (t-\frac{1}{t})) = \sum_{n=-\infty}^{\infty} J_n(x) t^n \tag{18}
\end{align}
This can be verified also by a direct expansion:
\begin{align}
&\quad \exp(\frac{x}{2} (t-\frac{1}{t})) \\
&= \sum_{n=0}^{\infty} \frac{(\frac{x}{2})^n (t-\frac{1}{t})^n}{n!} \\
&= \sum_{n=0}^{\infty} \frac{(\frac{x}{2})^n}{n!} \sum_{m=0}^{n} \frac{n!}{m!(n-m)!} t^m(-\frac{1}{t})^{n-m} \\
&= \sum_{n=0}^{\infty} (\frac{x}{2})^n \sum_{m=0}^{n} \frac{(-1)^{n-m}}{m!(n-m)!} t^{2m-n} \\
&= \sum_{m=0}^{\infty} \sum_{n=m}^{\infty} (\frac{x}{2})^n \frac{(-1)^{n-m}}{m!(n-m)!} t^{2m-n}\\
&= \sum_{m=0}^{\infty} \sum_{n=-m}^{\infty} (\frac{x}{2})^{n+2m} \frac{(-1)^{n+m}}{m!(n+m)!} t^{-n} \\
&= \sum_{n=-\infty}^{\infty} \sum_{m=\max(-n,0)}^{\infty} (\frac{x}{2})^{n+2m} \frac{(-1)^{n+m}}{m!(n+m)!} t^{-n} \\
&= \sum_{n=-\infty}^{\infty} \sum_{m=\max(n,0)}^{\infty} (\frac{x}{2})^{-n+2m} \frac{(-1)^{-n+m}}{m!\Gamma(-n+m+1)} t^{n} \tag{19}
\end{align}
When \(n \leq 0\), the inner summation is
\begin{align}
\sum_{m=0}^{\infty} (\frac{x}{2})^{-n+2m} \frac{(-1)^{-n+m}}{m!\Gamma(-n+m+1)} &= (-1)^n J_{-n}(x) = J_n(x) \tag{20}
\end{align}
by (8) and (16) of the last tutorial. And if \( n > 0 \), it also yields
\begin{align}
&\quad \sum_{m=n}^{\infty} (\frac{x}{2})^{-n+2m} \frac{(-1)^{-n+m}}{m!\Gamma(-n+m+1)} \\
&= \sum_{m=0}^{\infty} (\frac{x}{2})^{n+2m} \frac{(-1)^{m}}{(m+n)!\Gamma(m+1)} \\
&= \sum_{m=0}^{\infty} (\frac{x}{2})^{n+2m} \frac{(-1)^{m}}{\Gamma(m+n+1)m!} \\
&= J_n(x) \tag{21}
\end{align}
Integral Representation
With the generating function above, we can derive several integral representations for Bessel functions of integer order \(n\), one of them being
\begin{align}
J_n(x) = \frac{1}{\pi} \int_0^\pi \cos(n\theta -x \sin \theta) d\theta \tag{22}
\end{align}
To show this, first use the trigonometric identity to expand the integral in (22):
\begin{align}
&\quad \frac{1}{\pi} \int_0^\pi \cos(n\theta -x \sin \theta) d\theta \\
&= \frac{1}{\pi} \int_0^\pi (\cos(x\sin\theta)\cos n\theta +\sin(x\sin\theta)\sin n\theta) d\theta \tag{23}
\end{align}
Let \( t = \exp(i\theta) \) in (18) to get
\begin{align}
&\quad \exp(\frac{x}{2} (\exp(i\theta) -\exp(-i\theta))) \\
&= \exp(ix \sin \theta)= \sum_{m=-\infty}^{\infty} J_m(x) \exp(im\theta) \tag{24}
\end{align}
Then,
\begin{align}
&\quad \exp(ix \sin \theta)= \cos(x\sin\theta) + i \sin(x\sin\theta) \\
&= \sum_{m=-\infty}^{\infty} J_m(x) (\cos m\theta + i\sin m\theta) \tag{25}
\end{align}
where we have used de Moivre’s theorem. Comparing the real and imaginary parts of (25) then leads to
\begin{align}
\begin{aligned}
\cos(x\sin\theta) &= \sum_{m=-\infty}^{\infty} J_m(x) \cos m\theta \\
\sin(x\sin\theta) &= \sum_{m=-\infty}^{\infty} J_m(x) \sin m\theta
\end{aligned} \tag{26}
\end{align}
Plugging them back into (23) gives
\begin{align}
&\quad \frac{1}{\pi} \sum_{m=-\infty}^{\infty} \int_0^\pi ( J_m(x) \cos m\theta\cos n\theta +J_m(x) \sin m\theta\sin n\theta) d\theta \\
&= \frac{1}{\pi} [\frac{\pi}{2}J_n(x) + \frac{\pi}{2}J_n(x)] = J_n(x) \tag{27}
\end{align}
where we recall the orthogonality and normalization properties for trigonometric functions. Other integral representations can be found in DLMF 10.9.
Bessel Series
An observation is that there are infinitely many (countable) zeros of Bessel functions \( J_\nu(x) \) for any \( \nu \). To show this, we rewrite the Bessel equation into the canonical form as stated in this tutorial:
\begin{align}
u(x) &= \exp(-\frac{1}{2} \int \frac{1}{x} dx) = \exp(-\frac{1}{2} \ln x) = \frac{1}{\sqrt{x}} \tag{28} \\
g(x) &= (1-\frac{\nu^2}{x^2}) -\frac{1}{4} \frac{1}{x^2} -\frac{1}{2} (-\frac{1}{x^2}) \\
&= 1-\frac{\nu^2-1/4}{x^2} \tag{29}
\end{align}
Then we recall Theorem 4 of this tutorial, by noting that \( \int_1^\infty g(x) dx = \infty \) here. Now, according to the Spectral Theorem (see this tutorial), we can expand a given function into Bessel functions of the same order \( \nu \) with eigenvalues \( \lambda_1, \lambda_2, \ldots \) as long as these Bessel functions satisfy the boundary condition required to be an eigenfunction. If we choose the interval to be \( [0,1] \), then this means that the \( \lambda_j \) are just the zeros of Bessel functions \( J_\nu(\lambda_j) = 0 \). Then, for a given function \( f(x) \) defined on the same interval, we have
\begin{align}
f(x) &= \sum_{j=1}^{\infty} c_j J_\nu(\lambda_j x) = c_1 J_\nu(\lambda_1 x) + c_2 J_\nu(\lambda_2 x) + \cdots \tag{30}
\end{align}
as its Bessel Series. The coefficients \( c_j \) are then computed as
\begin{align}
c_j &= \frac{\int_0^1 xf(x) J_\nu(\lambda_j x) dx}{\int_0^1 x J_\nu(\lambda_j x)^2 dx} \\
&= \frac{\int_0^1 xf(x) J_\nu(\lambda_j x) dx}{\frac{1}{2}(J_\nu(\lambda_j)’)^2} \\
&= \frac{2\int_0^1 xf(x) J_\nu(\lambda_j x) dx}{J_{\nu+1}(\lambda_j)^2} \tag{31}
\end{align}
where we have used (8) and (15) with \(J_\nu(\lambda_j) = 0 \).
Exercise
Solve the following ODE
\begin{equation}
xy^{\prime\prime} + y = 0 \tag{32}
\end{equation}
by making the change of variable \( w = 2\sqrt{x} \), \( v = y/\sqrt{x} \), identifying it with the Bessel equation.
Answer
Following the hint, we have
\begin{align}
\frac{dy}{dx} &= \frac{d}{dx} (v\sqrt{x}) \\
&= \frac{d}{dw} (\frac{1}{2}vw) \frac{dw}{dx} \\
&= (\frac{1}{2}v + \frac{w}{2}\frac{dv}{dw}) \frac{1}{\sqrt{x}} \\
&= (\frac{1}{2}v + \frac{w}{2}\frac{dv}{dw}) \frac{2}{w} \\
&= \frac{v}{w} + \frac{dv}{dw} \\
\frac{d^2y}{dx^2} &= \frac{d}{dw} (\frac{v}{w} + \frac{dv}{dw}) (\frac{2}{w}) \\
&= \frac{2}{w}(\frac{1}{w}\frac{dv}{dw} -\frac{v}{w^2} + \frac{d^2v}{dw^2})
\end{align}
Substituting these into (32), we have
\begin{align}
\frac{w^2}{4} \frac{2}{w}(\frac{1}{w}\frac{dv}{dw} -\frac{v}{w^2} + \frac{d^2v}{dw^2}) + \frac{vw}{2} &= 0 \\
\frac{w}{2}\frac{d^2v}{dw^2} + \frac{1}{2}\frac{dv}{dw} -\frac{1}{2}\frac{v}{w} + \frac{vw}{2} &= 0 \\
w^2\frac{d^2v}{dw^2} +w\frac{dv}{dw} + (w^2-1) v &= 0
\end{align}
This is the Bessel equation with \( \nu = 1 \), and hence the general solution is
\begin{align}
v &= c_1 J_1(w) + c_2 Y_1(w) \\
y &= c_1 \sqrt{x} J_1(2\sqrt{x}) + c_2 \sqrt{x} Y_1(2\sqrt{x})
\end{align}
Leave a Reply