AEP 4830 HW3 Root Finding and Special Functions

Plotting of Bessel Functions

The Bessel Functions of the first kind \(J_{\nu}(x)\) and the second kind \(Y_{\nu}(x)\) are two functions of interest in this program. The program incorporates two header files ”nr3.h” and ”bessel.h” and outputs the function values between \(x=0\) to \(x=20\). The first two Bessel functions \(J_{0}\), \(J_{1}\), \(Y_{0}\) and \(Y_{1}\) are built in within ”bessel.h” and we can use the recurrence formula to construct higher order Bessel functions.

\begin{equation} J_{\nu+1}(x)=\frac{2\nu}{x}J_{\nu}(x)-J_{\nu-1}(x)\\ \end{equation} \begin{equation} Y_{\nu+1}(x)=\frac{2\nu}{x}Y_{\nu}(x)-Y_{\nu-1}(x)\\ \end{equation}

The first three \(J_{\nu}\) and \(Y_{\nu}\) are plotted by MATLAB as shown in Fig. 1 and Fig. 2.

The first three Bessel functions of the first kind, \(J_{\nu}(x)\).

The first three Bessel functions of the second kind, \(Y_{\nu}(x)\).

Root Finding

The next part of the program is to find the solution for the following equation.

\begin{equation} J_{0}(x)Y_{0}(x)=J_{2}(x)Y_{2}(x)\\ \end{equation}

First of all, I constructed a new function \(f(x)=J_{2}(x)Y_{2}(x)-J_{0}(x)Y_{0}(x)\) and the question now becomes the root finding for \(f(x)=0\). The function \(f(x)\) is plotted in Fig. 3. Secondly, two root finding algorithms are constructed by two template routines called bisect and falsepo for bisection and false position methods(W. H. Press 2007) respectively. The templates take in an interval bracketed by \(x_{a}\) and \(x_{b}\), the calculation tolerance and the max number of iterations (say, \(10^{4}\) to prevent infinite loops). Given an interval, the program finds the root \(x_{r}\) within the interval to 7 significant digits such that the function value is less than the specified tolerance \((10^{-7})\).

\begin{equation} |f(x_{r})|<10^{-7}\\ \end{equation}

In addition, the numbers of iteration of both algorithms are also shown in Table I.

The function \(f(x)\). We can choose the root finding interval by the plot.

The first five smallest values satisfying the equation \(J_{0}(x)Y_{0}(x)=J_{2}(x)Y_{2}(x)\) and the numbers of iteration of two root finding methods are tabulated.
Interval Root Bisection Iteration False Position Iteration
\((0.1,2.0)\) 0.694391 23 28
\((2.0,4.0)\) 2.924562 23 13
\((4.0,6.0)\) 4.574827 21 9
\((6.0,7.0)\) 6.181686 17 6
\((7.0,8.0)\) 7.773379 15 5


From the numbers of iterations of both methods, the root finding process converges fast. For a reasonable tolerance \(10^{-7}\), each takes less than 30 iterations to complete. However, in general, the false position method is much more efficient because it weighs the sample points differently by taking the linear approximation while bisection method merely keeps sampling the middle point of a given interval.
Interestingly, false position method actually takes more iterations to find the first root than bisection method does. Since the magnitude of the slope around the first root is large, even though \(f(x)\) changes much in each iteration, the \(x_{3}\) just moves a small step (say, \(10^{-}6\)) closer to the real root \(x_{r}\). As a result, the false position method need much more iterations during the final stage to get below the tolerance for an interval with large slope.


  1. S. A. Teukolsky W. H. Press, B. P. Flannery. Numerical Recipes, The Art of Scientific Computing, 3rd edit.. Camb. Univ. Press, 2007.

[Someone else is editing this]

You are editing this file