MATLAB GUI for computing Bessel functions using continued fractions algorithm ( GUI Matlab para o cálculo de funções de Bessel usando frações continuadas )

Higher order Bessel functions are prevalent in physics and engineering and there exist different methods to evaluate them quickly and efficiently. Two of these methods are Miller’s algorithm and the continued fractions algorithm. Miller’s algorithm uses arbitrary starting values and normalization constants to evaluate Bessel functions. The continued fractions algorithm directly computes each value, keeping the error as small as possible. Both methods respect the stability of the Bessel function recurrence relations. Here we outline both methods and explain why the continued fractions algorithm is more efficient. The goal of this paper is both (1) to introduce the continued fractions algorithm to physics and engineering students and (2) to present a MATLAB GUI (Graphic User Interface) where this method has been used for computing the Semi-integer Bessel Functions and their zeros.


Introduction
Bessel functions arise when using separation of variables to solve some partial differential equations in cylindrical and spherical coordinates [1].They appear naturally when dealing with Laplace's and Helmholtz's equations.The Bessel functions that result as solutions when solving these problems can be applied to various fields, including electricity and magnetism, heat conduction, acoustical vibrations, signal processing, and the radial Schrödinger equation [2].
Students are usually introduced to Bessel functions in their partial differential equations class.Attention is focused on the differential equation to obtain Bessel functions, with, usually, very little to the application of such functions.Due to the many applications of Bessel functions in several scientific fields, a curriculum should be designed to touch on the subject of actually using Bessel functions to solve real-world problems.
While Bessel functions are extremely useful, few algorithms exist to calculate them quickly and efficiently.A common method is Miller's algorithm [3].Another method is the continued fractions algorithm developed by Ratis and Fernández de Córdoba [4].These algorithms are necessary for various problems in physics.For example, when solving the inhomogeneous Bessel equation, Lommel functions arise.Lommel functions of two variables are superpositions of ordinary Bessel functions, and frequently appear when solving problems in diffraction [1].When studying scattering problems at high frequencies [5], one must use high order Hankel functions, which are linear combinations of Bessel functions of the first and second classes.A specific example of using numerical methods to evaluate Hankel functions, and therefore Bessel functions, can be found in an article of Havemann and Baran [6].The continued fractions algorithm can also be use to compute modified Bessel functions and Fresnel integrals [7,8].For problems when a numerically reliable and quick algorithm is needed, the continued fractions algorithm (CFA from now on) is a perfect candidate.In all of these examples, we need to know the values of many functions for a given argument at once.While Bessel functions can be easily evaluated using the built-in functions in Mathematica or MATLAB, these programs present some limitations when several orders and points are needed at the same time.
The goal of this paper is to present the CFA in a pedagogical manner for the use of physics and engineering students and for professors to implement in a classroom.We present both algorithms and give the advantages of the CFA, with the help of a MATLAB GUI that we have developed.This GUI can be downloaded from http://www.intertech.upv.es and it includes all the algorithms used to compute the Bessel functions and their zeros.
The structure of the paper is as follows: In the second section, we explain the concept of a continued fraction and how it is constructed.The third section presents the two previously mentioned algorithms for evaluating Bessel functions of high order.We give detailed descriptions of both Miller's algorithm and the CFA, and compare the two methods.The fourth section is devoted to applications and examples.Finally, some conclusions are outlined.

Continued fractions
In the mid 17th century, a large number of infinite methods were developed for directly computing π, including the method of continued fractions proposed by William Brouncker [9], the president of the Royal Society at the time.However, the theory of continued fractions goes back earlier.In Italy, Pietro Antonio Cataldi had already expressed square roots using this method [10].
Take, for example, √ 2. If we decompose it into the form √ 2 = 1 + x, we see that If we substitute Eq. (1) into itself, and continue this substitution indefinitely, we get a continued fraction Following the notation in Wall [11], we define a linear fractional transformation as where b 0 is an integer, and a n are positive integers, which allows us to construct a continued fraction by inserting z n (x) as the argument to the previous term and repeating indefinitely Consider the sequence given by the successive compositions of the transformations defined in Eq. ( 3) . . .
From this, we see that denoting the n-th convergent of the continued fraction, with z 0 (0) = b 0 denoting the 0-th convergent.
It is useful to use the equivalent notation for Eq. ( 5), following the notation of Wallis [12] where which allows us to rewrite the n-th convergent in Eq. ( 7) For a continued fraction to have convergence, the limit should exist and be finite.An efficient algorithm for calculating continued fractions is the Steed algorithm [13].The method of continued fractions explained in the next section uses the Steed algorithm to calculate a continued fraction.

Bessel function evaluation
Bessel functions are the canonical solutions, ω(z), of Bessel's differential equation where α, the order of the Bessel function, is an arbitrary number that can be either real or complex.Bessel functions of integer order have α = n = 0, 1, 2, . . .Semi-integer order Bessel functions, more commonly known as spherical Bessel functions, have α = n + 1/2, n = 0, 1, 2, . . .When z is a purely imaginary argument, we get modified Bessel functions [14].
Bessel functions are split into two different classes: first class and second class.First class Bessel functions are the solutions to Bessel's differential equation that are finite at the origin.Second class Bessel functions are the solutions to Bessel's differential equation that have a singularity at the origin.Usual methods of evaluating Bessel functions use recurrence relations [14].The class of the Bessel function determines whether we use ascending or descending recurrence relations to obtain the next term in the sequence by either increasing or decreasing n.The direction the recurrence relations take serve to maintain numerical stability [15].If we were to go the opposite way of the defined relations, the loss of significant figures would skew the results significantly [15].
Second class Bessel functions use an ascending recurrence relation to maintain stability, meaning we can use the first two known values to calculate all higher terms up to some value n = N .First class Bessel functions use a descending recurrence relation to maintain numerical stability.This means that we cannot start from the first two well known values, but must instead find a clever way around this hindrance.A commonly used method to compute Bessel functions is known as Miller's algorithm [3].

Miller's algorithm
Let us consider the first and second class spherical Bessel functions of semi-integer order, j n (z) and y n (z).
The recurrence relations for these two functions are given by As you can see, Eq. ( 13) is a descending recurrence relation, and Eq. ( 14) is an ascending recurrence relation.As we discussed before, using an ascending relation to evaluate j n would lead to absurd results, but we do not know j n nor j n−1 .Miller's algorithm serves to "guess"the initial values and re-normalize at a later step.
For a desired value n, we use N > n and assume that ĵN+1 = 0 and ĵN = 1 and then use the recurrence relation in Eq. ( 13) to obtain the sequence ĵN−1 , ĵN−2 , . . ., ĵ1 , ĵ0 .If we have chosen N large enough, the terms of this sequence up to n (i.e.ĵn , ĵn−1 , . . ., ĵ1 , ĵ0 ) will be proportional to the corresponding term in the sequence j n , j n−1 , . . ., j 1 , j 0 of desired values.The proportionality factor p can be obtained by comparing the value of ĵ0 with the known value of j 0 .
The terms of the sequence p ĵ0 , p ĵ1 , . . ., p ĵn−1 , p ĵn reproduce the required values j 0 , j 1 , . . ., j n−1 , j n .If the precision obtained is not sufficient, you can repeat the process with a larger value of N .
To obtain the second class Bessel functions, we simply use the initial values, y 0 (z) and y 1 (z), and apply the ascending recurrence relation until desired n.
Below is an example of the application of Miller's algorithm.We have included the same numerical example as illustrated in the work of Abramowitz and Stegun so that the reader may complete the details in the above work [14].

Miller's algorithm example
Suppose we want to evaluate the value of j 15 (x) for x = 24.6 using Miller's algorithm.Let us start, for example, at N = 39 and suppose ĵ40 = 0, ĵ39 = 1. (15) Using the recurrence The value p ĵ0 coincides with the value of j 0 (24.6) with 8 correct significant figures [14].Fortran can be used to evaluate Bessel functions using subroutines of the IMSL library.Numerical analysis for these routines can be found in the work of Ratis and Fernández de Córdoba [4].These routines are based on Miller's algorithm.

Continued fractions algorithm
The method of continued fractions introduced in section 2 can be used to directly evaluate the first class Bessel functions without any normalization.By applying the ascending recurrence relation to the second class Bessel functions, we generate the set {y n (z); n = 0, 1, 2, . . ., N }, for a desired value of the order n = N .We then use the last two values, y N (z) and y N −1 (z), a continued fraction, and the Bessel function Wronskian to solve for j N (z) and j N −1 (z).We then apply the descending recurrence relation to evaluate the first class Bessel functions, {j n (z); n = 0, 1, 2, . . ., N }.This method achieves the correct values without the need for a normalization factor.Relying on normalization relations to maintain stability hinders the performance speed of Miller's algorithm.By disposing of this necessity, the CFA runs faster, while still preserving stability by using the appropriate recurrence relations.
In Fig. 1 we present the flowcharts for Miller's algorithm and the CFA to illustrate each method.
We continue to use the usual notation of Abramowitz and Stegun and formally introduce spherical Bessel functions of the first and second class [14] as solutions to the differential equation, We generate all spherical Bessel functions of the second class, {y n (z), n = 0, 1, 2, . . ., N }, starting with the known values of by using the ascending recurrence relation in Eq. ( 14) until the fixed value N .
To calculate the highest order of the first class spherical Bessel functions, {j n (z), n = 0, 1, 2, . . ., N }, we use the calculated values of y N (z) and y N −1 (z), and the value of the spherical Bessel function Wronskian [14] We can rewrite Eq. () as To evaluate the ratio j n /j n−1 , we rearrange the recurrence relation for j n (z) given in Eq. ( 13) to read allowing us to construct the continued fraction Eq. ( 27) is easily evaluated using the Steed algorithm for a fixed N and z.Using the resulting value, we see that Equations ( 25) and (28) allow us to generate all spherical Bessel functions of the first class, {j n (z), n = 0, 1, 2, . . ., N }, by considering the calculated values of j N (z) and j N −1 (z) and using the descending recurrence relation, The CFA maintains the stability of each function by ensuring the use of the proper recurrence relations.Unlike Miller's algorithm, the CFA directly calculates the first class Bessel functions, rather than using arbitrary values and normalizing.Miller's algorithm relies on the normalization process, which requires more values than needed in order to converge to a reasonable proportionality factor.This necessity introduces not only another source of error, but also longer calculation times.A detailed numerical analysis can be found in the work of Ratis and Fernández de Córdoba [4].

MATLAB GUI development
In Fig. 2 we have shown the MATLAB GUI developed for this article.This GUI is divided in four parts.In the left-most section there are several tools to control the functions to plot, number of points, order and precision (Steed tolerance) of the CFA.The user can plot the Bessel function of order n or the complete set of functions from orders 0 to n.It is also possible to layer the graphics using the hold on option, and change the color of the new plots using the first menu of the second column.The y n (z) functions are plotted using a continuous line, while the j n (z) functions are shown by a dashed line.
The MATLAB's Bessel functions section compares the computation times using CFA and MATLAB's libraries, and checks the realtive error between CFA and MATLAB codes.It is easy to check that, for n = 50 with 500 points in the range (100, 200) of z, CFA code is much faster than MATLAB's.See Fig. 3 for this example.In a modern computer, the difference between the computation times of the two methods can be up to two orders of magnitude.In this figure it is possible to see that the relative error distribution in the j n (z) s is higher than that of y n (z) s.
The last group of the second column and the rightmost column are devoted to computing the roots of the Bessel functions.The algorithm that we have implemented is a combination of a brute force strategy and a bisection method.First, we compute in the desired interval (zmin, zmax) the Bessel functions with (zmax − zmin) * 10 points.Second, using this information, we compute the zeros in the subintervals where the functions change their sign, using a simple bisection method.There also exists another implementation of the code that computes the desired number of roots starting at zmin.The firsts ten zeros found from each function are displayed in the rightmost part of the figure.Beneath the graph, it is possible to save all the data computed.

Example
A nice exercise is to compute the zeros of y n (z) and check how they distribute in space.The same procedure can be used to check the zeros of j n (z), but the dotted plotting line makes higher order functions difficult to follow.One procedure to compare zeros is the following:  It is easy to see in Fig. 5 that the zeros are more dispersed as we increase the order.However, the student can increase the value of zmin to check that the larger z is, the closer the difference between zeros is to π.This is shown in Fig. 6, but one can also prove this statement using asymptotic Bessel function expansions [14].
Of course, the code and GUI can be easily modified in order to show many other interesting properties of Bessel functions.

Conclusions
We have explained how both Miller's algorithm and the continued fractions algorithm can be used to compute Bessel functions of high order in a manner conducive to the understanding of the average physics or engineering student.However, we focused on the benefits of using the method of continued fractions for such computations.
The continued fraction algorithm maintains the stability of each function by ensuring the use of the proper recurrence relations.Unlike Miller's algorithm, the continued fraction algorithm directly calculates the first class Bessel functions, rather than using arbitrary values and normalizing.Miller's algorithm relies on the normalization process, which requires more values than needed in order to converge to a reasonable proportionality factor.This necessity introduces not only another source of error, but also longer calculation times.A detailed numerical analysis can be found in Ratis and Fernández de Córdoba [4].A MATLAB implementation of this algorithm, together with a GUI, can be downloaded from http://www.intertech.upv.es.

Figura 3 - 2 . 3 .
Figura 3 -Relative error between CFA and MATLAB.The points with the largest error are close to the roots of the Bessel function involved.

Figura 5 -
Figura 5 -Variation of the mean difference between zeros in terms of the order.

Figura 6 -
Figura 6 -Difference between zeros of the y n (z) when z is large.As z increases, the difference approaches π.