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 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.

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.

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