Approximate closed-form formulas for the zeros of the Bessel Polynomials

We find approximate expressions x(k,n) and y(k,n) for the real and imaginary parts of the kth zero z_k=x_k+i y_k of the Bessel polynomial y_n(x). To obtain these closed-form formulas we use the fact that the points of well-defined curves in the complex plane are limit points of the zeros of the normalized Bessel polynomials. Thus, these zeros are first computed numerically through an implementation of the electrostatic interpretation formulas and then, a fit to the real and imaginary parts as functions of k and n is obtained. It is shown that the resulting complex number x(k,n)+i y(k,n) is O(1/n^2)-convergent to z_k for fixed k


Introduction
The polynomial solutions of the differential equation x 2 y ′′ (x) + 2(x + 1)y ′ (x) − n(n + 1)y(x) = 0 (1) were studied systematically in [1] by the first time. They are named Bessel polynomials and are given explicitly by where n = 0, 1, · · · . Many properties as well as applications are associated to these polynomials: the solution of the wave equation in spherical coordinates, network and filter design, isotropic turbulence fields, and more (see the monograph [2] or [3]- [14] and references therein for some other results). Among these, several results about the important problem concerning the location of its zeros have been obtained [8]- [11]. Just for instance, the use of the reverse Bessel polynomials θ n (x) = x n y n (1/x) in filter design is known since the time when these polynomials began to be studied [2,14]. Here, the poles of the transfer function are essentially the zeros of θ n (x). Thus, it is desirable to acquire new analytical knowledge about the location of the zeros of the Bessel polynomials.
In this note we give approximate explicit formulas for both the real and imaginary parts of the kth zero z k = x k + iy k of y n (x) and show that the approximation order of these new formulas to the exact zeros of the Bessel polynomials is O(1/n 2 ) for fixed k.
The approach followed in this note is simple and based on three items. The first is the electrostatic interpretation of the zeros of polynomials satisfying second order differential equations [15]- [17], the second is a simple curve fitting of numerical data and the third is the known fact that the points of well-defined curves in the complex plane are limit points of the zeros of the normalized Bessel polynomials [8]- [11]. The formulas yielded by the electrostatic interpretation of the zeros of Bessel polynomials are used to find them numerically as it has been done previously with these and other sets of points [7]- [19]. Several sets of zeros are computed in this way and the sets of real and imaginary values are fitted by polynomials depending on the index k whose coefficients depend on n. Finally, it is found that the approximate expression for the kth zero of y n (x) is O(1/n 2 )-convergent to a limit point of the zeros of the Bessel polynomials. Since the exact zero z k is O(1/n 2 )convergent to its limit point [8], we conclude that the approximation order of our approximate expression to the exact kth zero is also O(1/n 2 ).

Asymptotic expressions for the zeros
Let z k = x k + iy k , k = 1, 2, · · · , n, be the zeros of the Bessel polynomial y n (x). Then, from (1) follows that where j = 1, 2, · · · , n, i.e., the real and imaginary parts of the zeros should satisfy the electrostatic equations This set of nonlinear equations can be solved by standard methods. We have used a Newton method to solve them up to n = 500. As it is shown in Fig. 1, the piecewise linear interpolation of the real and imaginary parts of the zeros of the normalized Bessel polynomials y n (x/n) can be fitted by polynomials of the second and third degree in the index k. Thus, we propose the following expressions x(k, n) = a 2 (n)k 2 + a 1 (n)k + a 0 (n), to fit our data. To find the relationship between the coefficients of these polynomials and n, we take into account the numerical behavior of the data at the middle and end points. According to this,x(k, n) andỹ(k, n) can be determined byx and These conditions lead to the following coefficients a 2 (n) = 6 n 2 (n + 2) , a 1 (n) = − 6(n + 1) n 2 (n + 2) , a 0 (n) = 0, b 3 (n) = − 8(n − 2) n 2 (n 3 − 3n 2 + 2) , b 2 (n) = 12(n − 2)(n + 1) (n − 1)n 2 (n 2 − 2n − 2) , b 1 (n) = − 2 (n 3 + 6n 2 − 12n − 4) (n − 1)n 2 (n 2 − 2n − 2) , b 0 (n) = − n 3 − 5n 2 + 6 n (n 3 − 3n 2 + 2) .
The substitution of (7) in (4) yields approximate closed-form expressions k = 1, 2, · · · , n, that converge to the zeros z k of the Bessel polynomial y n (x), as we will show in the following.

Convergence
Following [9], we define and denote by Γ the curve defined by Γ = {z ∈ C : |W (z)| = 1 and |argz| ≥ π 2 }, which contains the limit pointsω k of the zeros of the normalized Bessel polynomial y n (x/n). Then, it has been proved in [8] that the zero ω k of y n (x/n) approaches to order O(1/n) the limit valueω k , i.e, as n → ∞. Thus, if we show that |ω k −ω k | = O(1/n), we will have proved that and therefore, taking into account that ω k = nz k , the explicit expression (8) approaches to order O(1/n 2 ) the zero z k of the Bessel polynomial y n (x).
To this purpose, we simply substituteω k = nz k in (9) to obtain, after a lengthily calculation, that the expansion of W (ω k ) in terms of 1/n is and this implies that for fixed k. Thus,ω k approaches to order O(1/n) the Γ curve and (11) follows. From here we have that as n → ∞. Numerical calculations confirm and extend this result. Figure 2 shows the behavior of the maxima of |z k −z k | over k as they depend on n. The numbers computed by (3) are taken as the exact zeros z k . The displayed data shows numerical uniform convergence on the values of k, not only for fixed k. A fit of these data gives 1/n a with a = 1.7.

Some few tests
Just to give examples of the application of the approximate expression (8), we consider the following cases.

The real zero
An closed-form formula for the unique real zero α(n) of the Bessel polynomial y n (x) can be obtained by the substitution of k = (n + 1)/2 in the real part of (8),x(k, n). This gives α(n) = − 3(n + 1) 2 2n 2 (n + 2) as our new result. In [2,11] asymptotic expressions for α(n) are given. Particularly, the following formula can be obtained from a more general expression given in [11]. Deleting the O-terms and expanding both results in powers of 1/n we find that indicating good agreement between the two approaches.

Power sums
Here we carry out the corresponding multiplications and use some cases of Faulhaber's formula. Then we compare our results with the exact ones. 1. Sum of the zeros. The simple sum ofz k [cf. (8)], gives The exact result is s 1 (n) = −1, as can be seen from (2). 2. Sum of the squares of the zeros. In this case we obtaiñ The exact sum s 2 (n) = 1 2n − 1 , can be found elsewhere [12]. Expanding both expressions in powers of 1/n we find thats 2 (n) = 11 21n + O(1/n 2 ), s 2 (n) = 1 2n + O(1/n 2 ).

Final comment
Note that the approximate formula for z k given above is not unique. There exist other functions to fit the zeros obtained through the electrostatic equations (3), and there are other conditions to impose at the extreme and middle points of the fitting interval. For instance, the imaginary partỹ(k, n), can be fitted by a polynomial of degree 5, but this does not improve the rate of convergence and, on the other hand, the calculations become more lengthy.