When you solve an initial boundary value problem, the boundary condition is used to find eigenvalues while initial conditions are used to determined Fourier coefficients of the series solution. In the previous discussion, we find the solution $R(r)=AJ_n(\lambda r)$ of Bessel’s equation. The boundary condition $R(1)=0$ results the equation $$J_n(\lambda)=0.$$ This is an eigenvalue equation and its solutions are called Bessel zeros. We can easily find Bessel zeros using Maple. First open Maple and type

Bess:=(n,m)->evalf(BesselJZeros(n,m));

and press enter. This defines a function Bess that returns $m$-th zero $\lambda_{nm}$ of $J_n(\lambda)$ for given positive integer $m$ as shown in the screenshot.

To calculate, for example, 1st zero of $J_0(\lambda)$ i.e. $\lambda_{01}$ simply type

Bess(0,1);

and press enter. It will return the value $2.40482558$ as shown in the screenshot.

In Maxima, unfortunately such a command (as BesselJZeros in Maple) does not appear to be available. But we can still find Bessel zeros using an elementary technique from calculus, Newton’s Method. If the initial approximation $x_0$ to a zero of a differentiable function $f(x)$ is given, then the $(n+1)$-th approximation is given by $$x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.$$ In Maxima, Newton’s Method can be performed as follows. As an example we calculate the 1st Bessel zero of $J_0(\lambda)$.

First we define $J_0(x)$ as a function $f(x)$:

(%i1) f(x):=bessel_j(0,x);

Newton’s method applied to $f(x)$ is performed by

(%i2) newton(t):=subst(t,x,x-f(x)/diff(f(x),x));

To use Newton’s method we need the initial approximation. We can make a guess from its graph. As you can see the first zero is near 2, so we choose $x_0=2$. If you just run

(%i3) newton(2);

the output will not be a numerical value. In order to have a numerical value as output, run

(%i4) ev(newton(2),numer);

The output will be

(%o4)                          2.388210765567796

In order to calculate next approximation $x_1$, simply run

(%i5) ev(newton(%),numer);

and its output is

(%o5)                          2.404769548213657

For the next approximation $x_2$, again run

(%i6) ev(newton(%),numer);

and its output is

(%o6)                          2.404825557043583

This value is accurate to nine decimal places. If you think this is good enough, you can stop. The whole process is seen in the following screenshot.

Let us denote the eigenvalues by $\lambda_{nm}$. Then the radial eigenfunctions are given by $$R_{nm}(r)=J_n(\lambda_{nm}r).$$
The solution $U(r,\theta)$ to the Helmholtz boundary value problem is
$$U_{nm}(r,\theta)=J_n(\lambda_{nm}r)[A\cos(n\theta)+B\sin(n\theta)].$$
Using a simple trigonometric identity, one can write
$$A\cos(n\theta)+B\sin(n\theta)=\sqrt{A^2+B^2}\cos(n(\theta-\psi))$$ for some $\psi$. Thus $U_{nm}(r,\theta)$ can be written as
$$U_{nm}(r,\theta)=AJ_n(\lambda_{nm}r)\cos(n\theta).$$
Oscillatory facotrs satisfy the equation $$\ddot{T}_{nm}+\lambda_{nm}^2c^2T_{nm}=0$$ and they are
$$T_{nm}(t)=C\cos(\lambda_{nm}ct)+D\sin(\lambda ct).$$
Now we are ready to put pieces together to write down the solution $u(r,\theta,t)$ as the following Fourier series:
$$u(r,\theta,t)=\sum_{n=0}^\infty\sum_{m=1}^\infty J_n(\lambda_{nm}r)\cos(n\theta)[A_{nm}\cos(\lambda_{nm} ct)+B_{nm}\sin(\lambda_{nm}ct)].$$

The only things remained for us to do are to determine the unknown Fourier coefficients $A_{nm}$ and $B_{nm}$. We can determine the Fourier coefficients using the initial conditions $u(r,\theta,0)$ and $u_t(r,\theta,0)$. This will be discussed in the next lecture.

A vibrating drumhead can be modeled by wave equation in polar coordinates and appropriate boundary and initial conditions. The modeling problem can be set up as follows. We assume that our drumhead is a round disk with unit radius ($r=1$).

PDE: $u_{tt}=c^2\left(u_{rr}+\frac{1}{r}u_r+\frac{1}{r^2}u_{\theta\theta}\right),\ 0<r<1$

BC: $u=0$ when $r=1$, $0<t<\infty$

ICs: \begin{align*}u(r,\theta,0)&=f(r,\theta),\\u_t(r,\theta,0)&=g(r,\theta)
\end{align*}
Note that the expression inside () is the Laplacian $\nabla^2 u$ in polar coordinates, i.e. $$\nabla^2u=u_{rr}+\frac{1}{r}u_r+\frac{1}{r^2}u_{\theta\theta}.$$
Let $u(r,\theta,t)=U(r,\theta)T(t)$. $U(r,\theta)$ is the shape of the drumhead while $T(t)$ is the oscillatory factor. Then the PDE can be broken into two differential equtions, one PDE and one ODE:
\begin{align*}
\nabla^2U+\lambda^2U&=0\ (\mbox{Helmholtz equation})\\\ddot{T}+\lambda^2c^2T&=0\ (\mbox{Simple harmonic motion})
\end{align*}
With the boundary condition, we have the Helmholtz eigenvalue problem:
\begin{align*}
\nabla^2U+\lambda^2U&=0\\
U(1,\theta)&=0.
\end{align*}
To solve this Helmholtz eigenvalue problem, we assume that $U(r,\theta)=R(r)\Theta(\theta)$. Then the Helmholtz equation is broken into two ODEs with BCs:
\begin{align*}
r^2R^{\prime\prime}+rR’+(\lambda^2r^2-n^2)R&=0\ (\mbox{Bessel’s equation}),\\
R(1)&=0
\end{align*}
and
\begin{align*}
\Theta^{\prime\prime}+n^2\Theta&=0,\\
\Theta(0)&=\Theta(2\pi).
\end{align*}
The boundary condition for $\Theta(\theta)$ is imposed because we want $\Theta$ to be periodic with period $2\pi$. We also impose an additional condition for $R(r)$:
$$R(0)<\infty\ (\mbox{physical condition}).$$

Let us discuss the solutions of the Bessel’s equation. The Bessel’s equation has two independent solutions
$J_n(\lambda r)$ and $Y_n(\lambda r)$, where $$J_n(x)=\sum_{k=0}^\infty\frac{(-1)^kx^{2k+n}}{2^{2k+n}k!(n+k)!}$$ and $$Y_n(x)=\frac{2}{\pi}J_n(x)\left[\ln\left(\frac{x}{2}\right)+\gamma\right]+\frac{2^n}{\pi x^n}\sum_{k=0}^\infty\frac{\beta_{mk}}{2^{2k}k!}x^{2k}.$$
Here, $\gamma$ is called the Euler-Macheroni constant and is given by
$$\gamma=\lim_{n\to\infty}\left[\left(\sum_{k=1}^n\frac{1}{k}\right)-\ln n\right]$$
and $\beta_{mk}$ is defined by
$$\beta_{mk}=\left\{\begin{array}{ccc}
-(n-1-k)! & \rm{if} & k\leq n-1,\\(-1)^{k-n-1}\frac{h_{k-n}+h_k}{(k-n)!} & \rm{if} & k>m-1,
\end{array}\right.$$ where $h_p=\sum_{i=1}^p\frac{1}{i}$.
$J_n(x)$ is called Bessel function of the first kind and $Y_n(x)$ is Bessel function of the second kind. $Y_n(x)$ is also called Neumann function. The general solution $R(r)$ of the Bessel’s equation is then written as
$$R(r)=AJ_n(\lambda r)+BY_n(\lambda r).$$ Since $Y_n(\lambda r)$ is not defined at $r=0$, due to the boundary condition $R(0)<\infty$, we require that $B=0$. So we consider only Bessel function of the first kind $J_n(x)$ here. Let us plot Bessel function of the first kind using Maxima. In Maxima, Bessel function of the first kind of order $n$ and argument $x$, $J_n(x)$ is denoted by bessel_j(n,x). To plot $J_1(x)$, $0\leq x\leq 10$, type the command

plot2d(bessel_j(1,x),[x,0,10]);

and press enter as you see in the following screenshot.

Then you will see the plot.

If you want to plot $J_2(x)$, type the command

plot2d(bessel_j(2,x),[x,0,10]);

and press enter.

If you want to plot multiple Bessel functions, for example,  $J_1(x)$, $J_2(x)$, $J_3(x)$, $J_4(x)$,  type the command

plot2d([bessel_j(1,x),bessel_j(2,x),bessel_j(3,x),bessel_j(4,x)],[x,0,10]);

and press enter.