ASSIGNMENT 4 Finding the Roots of Bessel Functions

MACM 316- D100 COMPUTING ASSIGNMENT 4

• Please follow the Guidelines for Computing Assignments found on the Canvas FAQ.
• Submit a one-page report via Crowdmark, grading scheme is also found on the FAQ.

Finding the Roots of Bessel Functions
The so-called Bessel Functions, denoted as Jm(x) where m = 0, 1, 2, ... is an integer, are a class of special
functions with many applications in mathematics and science. A plot of the first three Bessel Functions
can be found in figure 1. In this figure we have marked the first three roots of each, and in general we
denote the k
th positive root of Jm(x) as zm,k. However, since Jm(x) has no closed form representation for
any m, there is no way to find zm,k analytically for general m and k. Your mission, as you shall do in
this assignment, is to find the roots zm,k using both Newton’s method and Matlab’s fzero, as the
basis for comparing the two in terms of accuracy, efficiency, and robustness. Ultimately we will
apply our root-finding efforts to visualize solutions to the wave equation, a physically important PDE.
0 2 4 6 8 10 12
x
-0.5
0
0.5
1
y
First 3 Bessel Functions
J
0
J
1
J
2
Despite the fact that the Jm(x) are challenging to work with analytically, we do know some special
properties that can help us. First and foremost, we can evaluate Jm(x) numerically in Matlab, for any m
and x, using the function call besselj(m,x). Second, we have the following three location properties about
the zeros zm,k that you should in your numerical search:
• A guess for the 1
st zero of J0(x) is z0,1 ≈ 2.5.
• The k
th root of Jm(x) is larger than the (k − 1)th root by a value that is around π
• The 1
st root of Jm(x) is bracketted by the first two roots of Jm−1(x); that is zm−1,1 < zm,1 < zm−1,2.
The other piece of information that we need is how to compute derivatives for the Bessel Functions, in
particular for implementing Newton’s method. There is an incredibly useful recursive relationship for
this:
J
0
m(x) = 1
2
(Jm−1(x) − Jm+1(x))
Because we can evaluate any of the Bessel Functions numerically at any time using besselj, we have an
immediate expression for the derivatives.
DJM 1
MACM 316- D100 COMPUTING ASSIGNMENT 4
Your goal in this computation will be to produce a matrix whose entries are roots of the Bessel functions
[AM,K] =





z0,1 z0,2 · · · z0,K
z1,1 z1,2 · · · z1,K
.
.
.
.
.
.
zM,1 zM,2 · · · zM,K





where Jm(zm,k) = 0 (with k = 1, 2, · · · , K) over the range of the Bessel Functions {J0(x), · · · , JM(x)}.
What do you do with this matrix? We have supplied you with a movie-making script that needs [A16,16] as
its input. If you run the script without modification, it produces a nonsense surface plot and an animated
.gif file. But when you have your matrix of Bessel zeros properly calculated, the script CA4 BesselMovie.m,
then rewards you with a numerical movie that looks a lot better! You do not need to know what this
script is doing, other than it takes your input matrix [A16,16] and produces graphical plots.
To complete this assignment, you should do the following:
• Download the script CA4 bessel.m from Canvas. This code is already set up to find entries of the
matrix [A2,2] using both fzero AND Newton’s method. This code also tracks the number of function
evaluations needed for both methods. However, the initial conditions for Newton’s method have not
been chosen intelligently, so it fails to properly find the roots.
• Your first goal should be to generate 2 codes: one that uses only fzero to compute [A2,2], and the other
uses only Newton’s method. You should use bracketting initial guesses for the fzero code. The three
location properties from the first page will help with finding initial guesses/brackets that
lead to convergence.
• Next, we want to introduce some for loops into our codes to find A16,16
• When you have script working well, produce [A16,16], and modify the CA4 BesselMovie.m by removing the indicated clear command at the start. It should then produce a moving surface, but whose
outer circular edge is stationary with the values z = 0 (blue circle, not moving). If this blue circle is
moving, then you have some wrong entries. Also, there is a plot (figure 3000) of the outer values to
check how close they really are to being zero — these are related to your root-finding tolerance.
• Once you can correctly produce the matrix [A16,16], you can start comparing how the convergence
tolerance and choice of initial guesses/brackets effects performance of each method. You should think
about efficiency in terms of how many function evaluations are required by each method, accuracy
in terms of how close the boundary values are to numerical zero (figure 3000), and robustness. For
this last point, recall that the CA4 Bessel.m was initialized with values that cause failures.
• Your final report must include your discussion of experimental parameters and results. As well, you
need to include a figure of the final frame of the movie generated by CA4 BesselMovie.m, using your
computed [A16,16].
Optional Reading: If you are interested in learning more about the wave equation PDE being solved
in CA4 BesselMovie.m, we will direct you to the Wikipedia page for Vibrations of Circular Membranes:
https://en.wikipedia.org/wiki/Vibrations_of_a_circular_membrane
Recall from lecture that you should aim to communicate critical ideas with clarity. To this end, please
underline the three sentences in your report that you feel convey the biggest impact to the reader.
DJM 2