Approximate Closed-Form Formulas for the Zeros of the Bessel Polynomials

Abstract

We find approximate expressions 𝑥(𝑘,𝑛,𝑎) and 𝑦(𝑘,𝑛,𝑎) for the real and imaginary parts of the 𝑘th zero 𝑧_𝑘=𝑥_𝑘+𝑖𝑦_𝑘 of the Bessel polynomial 𝑦_𝑛(𝑥;𝑎). 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 𝑘, 𝑛 and 𝑎 is obtained. It is shown that the resulting complex number 𝑥(𝑘,𝑛,𝑎)+𝑖𝑦(𝑘,𝑛,𝑎) is 𝑂(1/𝑛²)-convergent to 𝑧_𝑘 for fixed 𝑘.

1. Introduction

The polynomial solutions of the differential equation 𝑧²𝑦^′′(𝑧)+(𝑎𝑧+2)𝑦^′(𝑧)−𝑛(𝑛+𝑎−1)𝑦(𝑧)=0,𝑎>0,𝑧∈ℂ, (1.1) were studied systematically in [1] by the first time. They are named (generalized) Bessel polynomials and are given explicitly by 𝑦_𝑛(𝑧;𝑎)=^𝑛_𝑘=0𝑛!(𝑛+𝑎−1)_𝑘()𝑧𝑛−𝑘!𝑘!2^𝑘, (1.2) as it can be shown in [2]. Here, (𝑥)_𝑘 is the Pochhammer symbol and 𝑛=0,1,…. Many properties as well as applications are associated to this equation; the traveling waves in the radial direction which are solutions of the wave equation in spherical coordinates can be written in terms of the polynomial solutions of (1.1). Also, this equation has application in 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] and in [12], explicit expressions for sum rules and for the homogeneous product sum symmetric functions of zeros of these polynomials are given. On the other hand, the electrostatic interpretation of these zeros as the equilibrium configuration in the complex plane with a logarithmic electric potential and a dipole at the origin has been given in [13], and in [14] it is shown that this equilibrium configuration is not stable. Thus, these cases show that it is desirable to acquire new analytical knowledge about the location of the zeros of the Bessel polynomials.

In this paper we give approximate explicit formulas for both the real and imaginary parts of the 𝑘th zero 𝑧_𝑘=𝑥_𝑘+𝑖𝑦_𝑘 of 𝑦_𝑛(𝑧;𝑎) and show that the approximation order of these new formulas to the exact zeros of the Bessel polynomials is 𝑂(1/𝑛²) for fixed 𝑘.

The approach followed in this paper 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 𝑘 whose coefficients depend on 𝑛 and 𝑎.

2. Asymptotic Expressions for the Zeros

Let 𝑧_𝑘=𝑥_𝑘+𝑖𝑦_𝑘, 𝑘=1,2,…,𝑛, be the zeros of the Bessel polynomial 𝑦_𝑛(𝑧;𝑎), ordered according to the imaginary part. Then, from (1.1) it follows that (A procedure for obtaining this kind of nonlinear equations for the zeros of a polynomial satisfying second and higher order differential equations is given in [19].) ^𝑛_𝑘=11𝑧_𝑗−𝑧_𝑘+𝑎𝑧_𝑗+22𝑧²_𝑗=0, (2.1) where 𝑗=1,2,…,𝑛, that is, the real and imaginary parts of the zeros should satisfy the electrostatic equations ^𝑛_𝑘=1𝑥_𝑗−𝑥_𝑘𝑥_𝑗−𝑥_𝑘²+𝑦_𝑗−𝑦_𝑘²+𝑎𝑥³_𝑗+2𝑥²_𝑗+𝑎𝑥_𝑗𝑦²_𝑗−2𝑦²_𝑗2𝑥²_𝑗+𝑦²_𝑗²=0,_𝑛𝑘=1𝑦_𝑗−𝑦_𝑘𝑥_𝑗−𝑥_𝑘²+𝑦_𝑗−𝑦_𝑘²+𝑦_𝑗𝑎𝑥²_𝑗+4𝑥_𝑗+𝑎𝑦²_𝑗2𝑥²_𝑗+𝑦²_𝑗²=0. (2.2) This set of nonlinear equations can be solved by standard methods. We have used a Newton method to solve them up to 𝑛=500 and 𝑎=100.

Let 𝜔_𝑘=𝜇_𝑘+𝑖𝜈_𝑘 be the 𝑘th zero of the normalized Bessel polynomials 𝑦_𝑛(2𝑧/(2𝑛+𝑎−2);𝑎), that is, 𝜇_𝑘=(2𝑛+𝑎−2)𝑥_𝑘/2 and 𝜈_𝑘=(2𝑛+𝑎−2)𝑥_𝑘/2. As it is shown in Figure 1, the piecewise linear interpolation of the real and imaginary parts of 𝜔_𝑘 can be fitted by polynomials of the second and third degree in the index 𝑘.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

Real and imaginary parts of the zeros of the normalized Bessel polynomials 𝑦𝑛(2𝑧/(2𝑛+𝑎−2);𝑎) for 𝑎=2,40,100. They are plotted as functions of 𝑘 for 𝑛=100,200,300,400,500, in gray-level intensity, from lower to higher, according to the value of 𝑛.

Thus, we propose the following expressions 𝜇(𝑘,𝑛,𝑎)=𝑎₂(𝑛,𝑎)𝑘²+𝑎₁(𝑛,𝑎)𝑘+𝑎₀(),𝜈()𝑛,𝑎𝑘,𝑛,𝑎=𝑏₃(𝑛,𝑎)𝑘³+𝑏₂(𝑛,𝑎)𝑘²+𝑏₁(𝑛,𝑎)𝑘+𝑏₀(𝑛), (2.3) for the approximate zero 𝜔_𝑘=𝜇(𝑘,𝑛,𝑎)+𝑖𝜈(𝑘,𝑛,𝑎) to fit our data. To find the dependence of the coefficients of these polynomials on 𝑛 and 𝑎, we take into account the numerical behavior of the data at the middle and end points.

We begin by finding the coefficients of the second-order polynomial giving the real part by fitting the values of 𝜇(𝑘,𝑛,𝑎) at 𝑘=1 and 𝑘=𝑛/2. In Figure 2 we show the dependence of 𝜇(1,𝑛,𝑎) and 𝜇(𝑛/2,𝑛,𝑎) on 𝑛 for some values of 𝑎.

(a)

(b)

Figure 2 

Dependence of 𝜇(1,𝑛,𝑎) and 𝜇(𝑛/2,𝑛,𝑎) on 𝑛 for 𝑎=10,20,30,40,100. Plots are shown in gray-level intensity, from lower to higher, according to the value of 𝑎.

A fit of these data to the models −𝐴/(𝑛+𝐵) and −𝐴/(𝑛+𝐵)−3/2 yields 𝜇(1,𝑛,𝑎)=−54𝑎+860,𝜇𝑛100𝑛+50𝑎+7152,𝑛,𝑎=−75𝑛−2𝑎+400.50𝑛+11𝑎+220 (2.4) These conditions and the symmetry of 𝜇(𝑘,𝑛,𝑎) with respect to the middle point lead to the following coefficients: 𝑎₂()=𝑝𝑛,𝑎²(𝑛,𝑎)𝑟(𝑛,𝑎),𝑎₁()=𝑝𝑛,𝑎¹(𝑛,𝑎)𝑟(𝑛,𝑎),𝑎₀()=𝑝𝑛,𝑎⁰(𝑛,𝑎)𝑟(),𝑛,𝑎 (2.5) where 𝑝₂(𝑛,𝑎)=47500𝑛²+50625𝑛+850𝑎𝑛−694𝑎²,𝑝−2770𝑎+96800₁(𝑛,𝑎)=−4(𝑛+1)7500𝑛²+50625𝑛+850𝑎𝑛−694𝑎²,𝑝−2770𝑎+96800₀(𝑛,𝑎)=−100130𝑛²−993𝑛−7656−20135𝑛²+627𝑛−1580𝑎−2(297𝑛+794)𝑎²,𝑟(𝑛,𝑎)=5𝑛(𝑛−2)(50𝑛+11𝑎+220)(20𝑛+10𝑎+143). (2.6) Now, to find the coefficients of the third-order polynomial giving the imaginary part we follow a similar procedure. The dependence of the real part of 𝜈(1,𝑛,𝑎) and 𝜕𝜈(𝑘,𝑛,𝑎)/𝜕𝑘|_𝑘=𝑛/2 on 𝑛 for some values of 𝑎 is shown in Figure 3.

(a)

(b)

Figure 3 

Dependence of 𝜈(1,𝑛,𝑎) and 𝜕𝜈(𝑘,𝑛,𝑎)/𝜕𝑘|𝑘=𝑛/2 on 𝑛 for 𝑎=10,20,30,40,100. Plots are shown in gray-level intensity, from lower to higher, according to the value of 𝑎.

Again, a fit of these data to the models 𝐴/(𝑛+𝐵)−1 and 𝐴/𝑛 yields 𝜈(1,𝑛,𝑎)=−25𝑛−𝑎−50(),𝜕𝜈()25𝑛−1𝑘,𝑛,𝑎𝜕𝑘_𝑘=𝑛/2=96𝑛().𝑎+25 (2.7) In addition, we have that 𝜈𝑛2,𝑛,2=0,𝜈(𝑛,𝑛,2)=−𝜈(1,𝑛,2), (2.8) therefore, the coefficients are given by 𝑏₃()=𝑞𝑛,𝑎³(𝑛,𝑎)𝑠(𝑛,𝑎),𝑏₂()=𝑞𝑛,𝑎²(𝑛,𝑎)𝑠(),𝑏𝑛,𝑎₁()=𝑞𝑛,𝑎¹(𝑛,𝑎)𝑠(𝑛,𝑎),𝑏₀()=𝑞𝑛,𝑎⁰(𝑛,𝑎)𝑠(),𝑛,𝑎 (2.9) where 𝑞₃(𝑛,𝑎)=−20023𝑛³−92𝑛²+90𝑛+4+200𝑛³−5𝑛²𝑎𝑛+8𝑛−6−8²−2𝑛+2𝑎²,𝑞₂(𝑛,𝑎)=10069𝑛⁴−255𝑛³+186𝑛²+92𝑛+8−1003𝑛⁴−12𝑛³+9𝑛²+4𝑛−12𝑎+43𝑛³−3𝑛²+4𝑎²,𝑞₁(𝑛,𝑎)=5021𝑛⁴+54𝑛³−522𝑛²+748𝑛−176−503𝑛⁴−9𝑛³−6𝑛²+26𝑛−24𝑎+23𝑛³−6𝑛+8𝑎²,𝑞₀(𝑛,𝑎)=−2525𝑛⁴−217𝑛³+618𝑛²𝑛−660𝑛+184−25⁴−2𝑛³−7𝑛²+16𝑛−12𝑎+𝑛³+𝑛²−4𝑛+4𝑎²𝑠(𝑛,𝑎)=25(𝑛−1)²(𝑛−2)²(𝑎+25). (2.10) Thus, the substitution of (2.10), (2.9), (2.6), and (2.5), respectively, in (2.3) yields the approximate closed-form expressions 𝑧_𝑘=𝑥()𝑦()=2𝜇()𝑘,𝑛,𝑎+𝑖𝑘,𝑛,𝑎𝑘,𝑛,𝑎2𝜈()2𝑛+𝑎−2+𝑖𝑘,𝑛,𝑎,2𝑛+𝑎−2 (2.11) where 𝑘=1,2,…,𝑛, for the zeros of the unnormalized Bessel polynomial 𝑦_𝑛(𝑧;𝑎). The expressions given in (2.11) converge to the zeros 𝑧_𝑘 of these polynomials, as we will show in the following.

3. Convergence

Following [9], we define 𝑊(𝑧)=𝑒^1+1/𝑧2𝑧1+1+1/𝑧², (3.1) and denote by Γ the curve defined by |𝑊(𝑧≥𝜋Γ=𝑧∈ℂ∶)|=1,arg𝑧2, (3.2) which contains the limit points 𝜔_𝑘 of the zeros of the normalized Bessel polynomial 𝑦_𝑛(2𝑧/(2𝑛+𝑎−2);𝑎). Then, it has been proved in [8] that the zero 𝜔_𝑘 of 𝑦_𝑛(2𝑧/(2𝑛+𝑎−2);𝑎) approaches to order 𝑂(1/𝑛) the limit value 𝜔_𝑘, that is, 𝜔_𝑘−𝜔_𝑘1=𝑂𝑛, (3.3) as 𝑛→∞.

Thus, if we show that |𝜔_𝑘−𝜔_𝑘|=𝑂(1/𝑛), we will have proved that 𝜔_𝑘−𝜔_𝑘1=𝑂𝑛, (3.4) and therefore, taking into account that 𝜔_𝑘=(2𝑛+𝑎−2)𝑧_𝑘/2, the explicit expression (2.11) approaches to order 𝑂(1/𝑛²) the zero 𝑧_𝑘 of the Bessel polynomial 𝑦_𝑛(𝑧;𝑎).

To this purpose, we simply substitute the expression for 𝜔_𝑘=𝜇(𝑘,𝑛,𝑎)+𝑖𝜈(𝑘,𝑛,𝑎) given by (2.3) in (3.1) to obtain, after a lengthy calculation, that the expansion of 𝑊(𝜔_𝑘) in terms of 1/𝑛 is 𝑊𝜔_𝑘=1+2𝑎²−100𝑎(3𝑘−2)−50(42𝑘−67)()+25𝑎+25ℎ(𝑘,𝑎)25(𝑎+25)cos²(𝑡(𝑘,𝑎))−sin²(𝑡(1𝑘,𝑎))𝑛1+𝑂𝑛^3/2, (3.5) where ℎ(𝑘,𝑎)=(25+𝑎)²(130+27𝑎+300𝑘)²+2𝑎²−100𝑎(3𝑘−2)−50(42𝑘−67)²,𝑡()=()𝑘,𝑎25+𝑎)(130+27𝑎+300𝑘41675−100𝑎+𝑎².−1050𝑘+150𝑎𝑘 (3.6) This implies that 𝑊𝜔_𝑘1=1+𝑂𝑛, (3.7) for fixed 𝑘 and 𝑎. Thus, 𝜔_𝑘 approaches to order 𝑂(1/𝑛) the Γ curve and (3.4) follows. From here we have that 𝑧_𝑘−𝑧_𝑘1=𝑂𝑛², (3.8) as 𝑛→∞. Numerical calculations confirm and extend this result. Figure 4 shows the behavior of the maxima of |𝑧_𝑘−𝑧_𝑘| over 𝑘 as they depend on 𝑛 for the particular case of 𝑎=2. The numbers computed by (2.2) are taken as the exact zeros 𝑧_𝑘. A fit of these data gives 1/𝑛^𝑎 with 𝑎=1.7.

Figure 4 

Plot of the values of max𝑛𝑘=1|𝑧𝑘−𝑧𝑘| against 𝑛 for the case of 𝑎=2.

4. Some Few Tests

Just to give examples of the application of the approximate expression (2.11), we consider the following cases.

4.1. The Real Zero

A closed-form formula for the unique real zero 𝛼_𝑛(𝑎) of the Bessel polynomial 𝑦_𝑛(𝑧;𝑎) can be obtained by the substitution of 𝑘=(𝑛+1)/2 in the real part of (2.11), 𝑥(𝑘,𝑛). This gives 𝛼_𝑛(𝑎)=−4𝑛3+1()+7550−71𝑎−707𝑎²+5900𝑎+3500017500𝑛+𝑂𝑛², (4.1) as our new result. In [2, 11] very accurate expressions for 𝛼_𝑛(𝑎) are given. Particularly, the following formula: 2𝛼_𝑛(𝑎)1=−1.325486838𝑛−1.00628995𝑎+1.349836480+𝑂,2𝑛+𝑎−2 (4.2) is given in [11]. Expanding 2/𝛼_𝑛(𝑎) in powers of 𝑛 we find that 2𝛼_𝑛(𝑎)≃−1.33333𝑛−0.946667𝑎+0.666667, (4.3) indicating good relative agreement between the two results.

4.2. 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 of 𝑧_𝑘 (cf. (2.11)) gives a complicated expression for the real part. However, expanding both the real and imaginary parts of this sum gives 𝑠₁̃(𝑛)=^𝑛_𝑘=1𝑧_𝑘=−1+53𝑎−1007130+𝑖321𝑎+25𝑛1+𝑂𝑛². (4.4) The exact result is 𝑠₁(𝑛)=−1, as can be seen from (1.2).

(2)  Sum of the Squares of the Zeros. In this case we take the particular case of 𝑎=1. For this value we obtain 𝑠₂̃(𝑛)=^𝑛_𝑘=1𝑧²_𝑘=346915915𝑛+𝑂𝑛². (4.5) The exact sum 𝑠₂(𝑛)=1=12𝑛−112𝑛+𝑂𝑛² (4.6) can be found elsewhere [12].

The use of the approximate formula (2.11) for obtaining sums of higher powers of these zeros is not expected to give satisfactory results, since the powers and the sum magnify the total error.

5. Final Comment

The approximate formula for 𝑧_𝑘 given above is far from being unique. There exist many other functions to fit the zeros obtained through the electrostatic equations (2.2), and there are other conditions to impose at the extreme and middle points of the fitting interval. For instance, the imaginary part 𝑦(𝑘,𝑛) 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 complicated.

Acknowledgment

The authors thank Consejo Nacional de Ciencia y Tecnología for the financial support given to this project.