Modeling Vibrations by Bessel Functions
=======================================

The wave equation models vibrations.

Vibrations
----------

To model vibrations, use *the wave equation*
:math:`\displaystyle \frac{\partial^2 u}{\partial t^2} = \nabla^2 u`.
For example, to model a vibrating string, use:

.. math::

   \displaystyle \frac{\partial^2 u}{\partial t^2}
   = \frac{\partial^2 u}{\partial x^2},
   \quad \mbox{for} \quad x \in [0,1], \quad t \geq 0,

.. index:: wave equation

subject to the initial conditions:

1. :math:`u(x, 0) = f(x)`, profile of a string,

2. :math:`u_t(x, 0) = 0`, no velocity,

and boundary conditions:

3. :math:`u(0, t) = 0`, string is fixed at the left end,

4. :math:`u(1, t) = 0`, string is fixed at the right end.

The method of :index:`separation of variables` provides 
only a weak solution.

The Wave Equation
-----------------

We apply the wave equation to the following situation.
Strike a drum at its center and 
the :index:`vibrations` are radially symmetric.
Assuming the drum has radius one, in polar coordinates,

.. math::

   \quad u(r, t)

gives deviation from the equilibrium state zero,
at distance :math:`r` from the center of the drum,
at time :math:`t`, starting at the strike of the drum.

As the vibrations are radially symmetric,
we have a one dimensional spatial problem.

We solve the wave equation in polar coordinates

.. math::

   \frac{\partial^2 u}{\partial t^2}
   = \nabla^2 u
   % = \frac{\partial^2 u}{\partial r^2},  this is wrong!
   \quad 0 \leq r < 1, \quad t \geq 0,

with initial conditions

1. :math:`u(r, 0) = f(r)`, the profile of the drum
   at :math:`t = 0`, :math:`r < 1`,

2. :math:`u_t(r, 0) = 0`, we strike once, at :math:`t = 0`,

and boundary condition

3. :math:`u(1, t) = 0`, for :math:`t \geq 0`,
   the drum is fixed at end, at :math:`r = 1`.

In the method of the separation of variables,
we propose :math:`U(r, t)` as the form of the solution

.. math::

   U(r, t) = R(r) T(t),

where 

* :math:`R` is a function depending only on :math:`r`, and

* :math:`T` is a function depending only on :math:`t`.

Subjecting :math:`U` to the conditions of the wave equation
:math:`\displaystyle \frac{\partial^2 u}{\partial t^2} = \nabla^2 u`
we compute

.. math::

   \frac{\partial^2 U}{\partial t^2} = R(r) 
   \frac{\partial^2 T}{\partial t^2} 
   \quad \mbox{and} \quad
   \frac{\partial^2 U}{\partial r^2} = \nabla^2 R \cdot T(t).

Then, seperating to the left and the right:

.. math::

   R(r) \frac{\partial^2 T}{\partial t^2} 
   = \nabla^2 R \cdot T(t)
   \quad \mbox{or equivalently} \quad
   \frac{1}{T(t)} \frac{\partial^2 T}{\partial t^2} 
   = \frac{1}{R(r)} \nabla^2 R.

The assumption of :math:`U(r, t) = R(r) T(t)` for solution 
implies the existence of a separation constant :math:`\lambda`,
wedged in between the equation:

.. math::

   \frac{1}{T(t)} \frac{\partial^2 T}{\partial t^2} 
   = \lambda = \frac{1}{R(r)} \nabla^2 R.

We obtain two decoupled ordinary differential equations

.. math::

   \frac{\partial^2 T}{\partial t^2} = \lambda T(t),

and

.. math::

   \nabla^2 R = \lambda R(r).

*What is :math:`\nabla^2 R`?*

For :math:`u = u(x,y)`, the Laplacian
:math:`\displaystyle \nabla^2 u (x,y) = \frac{\partial^2 u}{\partial x^2}
+  \frac{\partial^2 u}{\partial x^2}`
in polar coordinates, 
for :math:`u = u(x = r\cos(\theta), y = r \sin(\theta))`, is

.. math::

   \nabla^2 u (r, \theta) = 
   \frac{\partial^2 u}{\partial r^2}
   + \frac{1}{r} \frac{\partial u}{\partial r}
   + \frac{1}{r^2} \frac{\partial r^2}{\partial \theta^2}.

Recall exercise 3 of L-36.

For :math:`u = R(r)`, we then have
:math:`\displaystyle \nabla^2 R (r) = \frac{\partial^2 R}{\partial r^2}
+ \frac{1}{r} \frac{\partial R}{\partial r}`
and we can write

.. math::

   \nabla^2 R = \lambda R(r)
   \quad \mbox{as} \quad
   \frac{\partial^2 R}{\partial r^2}
   + \frac{1}{r} \frac{\partial R}{\partial r} = \lambda R(r).

As a solution to

.. math::

   \frac{\partial^2 T}{\partial t^2} = - \alpha^2 T(t),

consider a linear combination of a cosine and a sine function:

.. math::

   \begin{array}{rcl}
   T(t) & = & a \cos(\alpha t) + b \sin(\alpha t) \\
   T'(t) & = & -a~\!\alpha \sin(\alpha t) + b~\!\alpha \cos(\alpha t) \\
   T''(t) & = & -a~\!\alpha^2 \cos(\alpha t) - b~\!\alpha^2 \sin(\alpha t) \\
   & = & -\alpha^2 T(t).
   \end{array}

Therefore, for some constants :math:`a` and :math:`b`,

.. math::

   T(t) = a \cos(\alpha t) + b \sin(\alpha t)

is a solution to the second order differential equation.

Applying the second initial condition, subject :math:`T(t)`
to the second initial condition:

2. :math:`u_t(r, 0) = 0`, we strike once, at :math:`t = 0`,

which gives

.. math::

   T'(t) = - a \sin(\alpha t) + b \cos(\alpha t),

evaluate at :math:`t = 0`: :math:`T'(0) = b`, thus :math:`b = 0`, and

.. math::

   T(t) = a \cos(\alpha t).

Bessel Functions
----------------

The application of separation of variables in polar coordinates
to the wave equation is

.. math::

   \frac{\partial^2 R}{\partial r^2}
   + \frac{1}{r} \frac{\partial R}{\partial r} = -\alpha^2 R,

or equivalently, after multiplication with :math:`r^2`:

.. math::

   r^2 \frac{\partial^2 R}{\partial r^2}
   + r \frac{\partial R}{\partial r} + r^2 \alpha^2 R = 0.

We eliminate :math:`\alpha` by the substitution :math:`r = z/\alpha`:

.. math::

   \left(\frac{z^2}{\alpha^2} \right)
   \alpha^2 \frac{\partial^2 R}{\partial z^2}
   + \left( \frac{z}{\alpha} \right) \alpha \frac{\partial R}{\partial z} 
   + \left( \frac{z^2}{\alpha^2} \right) \alpha^2 R 
   = 
   z^2 \frac{\partial^2 R}{\partial z^2}
   + z \frac{\partial R}{\partial z} + z^2 R = 0.

The *Bessel equation*

.. math::
   z^2 \frac{\partial^2 R}{\partial z^2}
   + z \frac{\partial R}{\partial z} + z^2 R = 0

has two independent solutions:

1. :math:`J_0` is the :index:`Bessel function` of the first kind,

2. :math:`Y_0` is the Bessel function of the second kind.

We can verify this using ``sympy``.
Bessel functions in Julia are provided 
by the package ``SpecialFunctions``.

Starting with the Bessel functions of the first kind:

::

   julia> using SymPy
   julia> z = Sym("z")
   julia> J0 = sympy.besselj(0, z)
   julia> d1J0 = J0.diff(z)
   julia> d2J0 = d1J0.diff(z)

   julia> odeJ0 = z^2*d2J0 + z*d1J0 + z^2*J0

   julia> sympy.simplify(odeJ0)
   0

This verifies symbolically that the Bessel functions
of the first kind satisfy the ordinary differential equation

.. math::

   z^2 \frac{\partial^2 R}{\partial z^2}
   + z \frac{\partial R}{\partial z} + z^2 R = 0.

We contine the Julia session
exploring the Bessel functions of the second kind:

::

   julia> Y0 = sympy.bessely(0, z)
   julia> d1Y0 = Y0.diff(z)
   julia> d2Y0 = d1Y0.diff(z)

   julia> odeY0 = z^2*d2Y0 + z*d1Y0 + z^2*Y0

   julia> sympy.simplify(odeY0)
   0

This verifies symbolically that the Bessel functions
of the second kind satisfy the ordinary differential equation

.. math::

   z^2 \frac{\partial^2 R}{\partial z^2}
   + z \frac{\partial R}{\partial z} + z^2 R = 0.

The Bessel function of the first kind is shown
in :numref:`figBesselJ0`.
We can think of :math:`J_0` as a damped version of a cosine.

.. _figBesselJ0:

.. figure:: ./figBesselJ0.png
   :align: center

   The Bessel function of the first kind.

In :numref:`figBesselY0`, we see
the Bessel function of the second kind.
Observe that :math:`Y_0` is no longer bounded near zero.

.. _figBesselY0:

.. figure:: ./figBesselY0.png
   :align: center

   The Bessel function of the second kind.

Consider the first initial condition:

1. :math:`u(r, 0) = f(r)`, the profile of the drum
   at :math:`t = 0`, :math:`r < 1`,

The solution must be bounded, 
so we use only :math:`J_0`, not :math:`Y_0`.

Consider the boundary condition:

2. :math:`u(1, t) = 0`, for :math:`t \geq 0`,
   the drum is fixed at end, at :math:`r = 1`.

The shape of the solution is :math:`U(r, t) = T(t) R(r)`, with

* :math:`T(t) = a \cos(\alpha t)`,

* :math:`R(r) = J_0(z = \alpha r)`,

thus:

.. math::

   U(r, t) = a \cos(\alpha t) J_0(\alpha r).

The condition :math:`u(1, t) = 0` gives, setting :math:`r = 1`:

.. math::

   \underbrace{a \cos(\alpha t)}_{\not= 0} J_0(\alpha) = 0
   \quad \Rightarrow \quad J_0(\alpha) = 0.

The first ten roots are computed as follows.

::

   julia> using FastGaussQuadrature

   julia> roots = approx_besselroots(0.0, 10)
   10-element Vector{Float64}:
    2.404825557695773
    5.520078110286311
    8.653727912911013
   11.791534439014281
   14.930917708487785
   18.071063967910924
   21.21163662987926
   24.352471530749302
   27.493479132040253
   30.634606468431976

   julia>

Now we derive the general solution.
For every root :math:`\alpha_n$ of $J_0`, we have a solution of the form

.. math::

   U_n(r, t) = a_n \cos(\alpha_n t) J_0(\alpha_n r).

Viewing :math:`U_n` as a basis of functions, we have the general solution

.. math::

   U(r, t) = \sum_{n=1}^\infty a_n \cos(\alpha_n t) J_0(\alpha_n r).

To determine the coefficients :math:`a_n`, use the condition

1. :math:`u(r, 0) = f(r)`, the profile of the drum
   at :math:`t = 0`, :math:`r < 1`.

As before, we are using an :index:`orthogonal basis`.

Applying :math:`u(r, 0) = f(r)` to
:math:`\displaystyle U(r, t) 
= \sum_{n=1}^\infty a_n \cos(\alpha_n t) J_0(\alpha_n r)` yields

.. math::

   U(r, 0) = \sum_{n=1}^\infty a_n  J_0(\alpha_n r) = f(r).

With the inner product

.. math::
   \big\langle f(r), g(r) \big\rangle = \int_0^1 f(r)~\!g(r)~\!r~\!dr.

we have

.. math::

   \big\langle J_0(\alpha_n r), J_0(\alpha_m r) \big\rangle
    = 0, \quad \mbox{for } n \not= m.

By :math:`\displaystyle
\big\langle J_0(\alpha_n r), J_0(\alpha_m r) \big\rangle = 0`,
for :math:`n \not= m`,

.. math::

   \begin{array}{rcl}
   {\displaystyle \big\langle f(r), J_0(\alpha_m r) \big\rangle}
   & = & 
   {\displaystyle \big\langle \sum_{n=1}^\infty a_n  J_0(\alpha_n r) ,
   J_0(\alpha_m r) \big\rangle} \\
   & = & 
   {\displaystyle a_m ~ \big\langle J_0(\alpha_m r) , 
   J_0(\alpha_m r) \big\rangle}.
   \end{array}

The coefficients :math:`a_m` are computed as

.. math::

   a_m = \frac{\big\langle f(r), J_0(\alpha_m r) \big\rangle}
   {\big\langle J_0(\alpha_m r) , J_0(\alpha_m r) \big\rangle}.

Proposal for a Topic of a Project
---------------------------------

*Can you hear the shape of a drum?*
is the title of the paper by Mark Kac,
published in April 1966 by *The American Mathematical Monthly*,
pages 1-23.

This problem is still *open*, according to the paper
**Can't One Really Hear the Shape of a Drum?**
by S. Zuluaga and F. Fonseca, 
in *Acoustical Physics*, Vol. 57, No. 4, pages 465--472, 2011.

* Study the papers and illustrate the importance of this problem
  for industrial applications.

* Summarize the findings of the two papers.

Exercises
---------

1. Use the method of characteristic lines to show that

   .. math::

      u(x, t) = f(x - c t) + g(x + c t)

   satisfies the one dimensional wave equation

   .. math::

      \displaystyle \frac{\partial^2 u}{\partial t^2}
      = c^2 \frac{\partial^2 u}{\partial x^2},
      \quad \mbox{for} \quad x \in [0,1], \quad t \geq 0,

   where :math:`c` is the wave speed.
   The functions :math:`f` and :math:`g` are determined within constants
   that sum to zero.

2. Consider the Bessel functions :math:`J_n`,
   for :math:`n=0,1,2,3,4`.

   Make a plot of :math:`J_n(z)`,
   for :math:`z \in [0,1]`, and for :math:`n=0,1,2,3,4`.

   What do you observe about the roots of those five functions?

3. Verify the orthogonality numerically for :math:`J_0(\alpha_n r)`
   for :math:`n = 1, 2, 3, 4, 5`.

4. For a square drum, for :math:`x \in [0,1]` and :math:`y \in [0,1]`,
   use the proposed solution

   .. math::

      u(x, y, t) = \sum_{m=1}^\infty \sum_{n=1}^\infty
      c_{m,n} \cos\left(\pi t \sqrt{m^2 + n^2}\right)
      \sin(m \pi x) \sin(n \pi y). 

   1. Define the partial differential equation,
      with all initial and boundary conditions for the square drum.

   2. Interpret the solution and argue why this proposal makes sense.

   3. How are the coefficients :math:`c_{m,n}` computed?

Bibliography
------------

1. Charles R. MacCluer, problem 5 of chapter 11 of
   *Industrial Mathematics.  Modeling in Industry, Science, and Government*,
   Prentice Hall, 2000.