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

Authors: Rafael G. Campos; Marisol L. Calderón
Journal: International Journal of Mathematics and Mathematical Sciences
 (2012)

---

## Abstract

We find approximate expressionsx̃(k,n,a) and ỹ(k,n,a) for the real and imaginary parts of the kth zero zk=xk+iyk of the Bessel polynomial yn(x;a). 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, n and a is obtained. It is shown that the resulting complex number x̃(k,n,a)+iỹ(k,n,a) is O(1/n2)-convergent to zk for fixed k.


---

## Body

## 1. Introduction

The polynomial solutions of the differential equation(1.1)z2y′′(z)+(az+2)y′(z)-n(n+a-1)y(z)=0,a>0,z∈ℂ,
were studied systematically in [1] by the first time. They are named (generalized) Bessel polynomials and are given explicitly by
(1.2)yn(z;a)=∑k=0n‍n!(n+a-1)k(n-k)!k!(z2)k,
as it can be shown in [2]. Here, (x)k is the Pochhammer symbol and n=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 thekth zero zk=xk+iyk of yn(z;a) and show that the approximation order of these new formulas to the exact zeros of the Bessel polynomials is O(1/n2) for fixed k.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 k whose coefficients depend on n and a.


## 2. Asymptotic Expressions for the Zeros

Letzk=xk+iyk, k=1,2,…,n, be the zeros of the Bessel polynomial yn(z;a), 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].) (2.1)∑k=1n‍1zj-zk+(azj+2)2zj2=0,
where j=1,2,…,n, that is, the real and imaginary parts of the zeros should satisfy the electrostatic equations
(2.2)∑k=1n‍xj-xk(xj-xk)2+(yj-yk)2+axj3+2xj2+axjyj2-2yj22(xj2+yj2)2=0,∑k=1n‍yj-yk(xj-xk)2+(yj-yk)2+yj(axj2+4xj+ayj2)2(xj2+yj2)2=0.
This set of nonlinear equations can be solved by standard methods. We have used a Newton method to solve them up to n=500 and a=100.Letωk=μk+iνk be the kth zero of the normalized Bessel polynomials yn(2z/(2n+a-2);a), that is, μk=(2n+a-2)xk/2 and νk=(2n+a-2)xk/2. As it is shown in Figure 1, the piecewise linear interpolation of the real and imaginary parts of ωk can be fitted by polynomials of the second and third degree in the index k.Real and imaginary parts of the zeros of the normalized Bessel polynomialsyn(2z/(2n+a-2);a) for a=2,40,100. They are plotted as functions of k for n=100,200,300,400,500, in gray-level intensity, from lower to higher, according to the value of n.
(a)
(b)
(c)
(d)
(e)
(f)Thus, we propose the following expressions(2.3)μ~(k,n,a)=a2(n,a)k2+a1(n,a)k+a0(n,a),ν~(k,n,a)=b3(n,a)k3+b2(n,a)k2+b1(n,a)k+b0(n),
for the approximate zero ω~k=μ~(k,n,a)+iν~(k,n,a) to fit our data. To find the dependence of the coefficients of these polynomials on n and a, 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μ~(k,n,a) at k=1 and k=n/2. In Figure 2 we show the dependence of μ~(1,n,a) and μ~(n/2,n,a) on n for some values of a.Dependence ofμ~(1,n,a) and μ~(n/2,n,a) on n for a=10,20,30,40,100. Plots are shown in gray-level intensity, from lower to higher, according to the value of a.
(a)
(b)A fit of these data to the models-A/(n+B) and -A/(n+B)-3/2 yields
(2.4)μ~(1,n,a)=-54a+860100n+50a+715,μ~(n2,n,a)=-75n-2a+40050n+11a+220.
These conditions and the symmetry of μ~(k,n,a) with respect to the middle point lead to the following coefficients:
(2.5)a2(n,a)=p2(n,a)r(n,a),a1(n,a)=p1(n,a)r(n,a),a0(n,a)=p0(n,a)r(n,a),
where
(2.6)p2(n,a)=4(7500n2+50625n+850an-694a2-2770a+96800),p1(n,a)=-4(n+1)(7500n2+50625n+850an-694a2-2770a+96800),p0(n,a)=-100(130n2-993n-7656)-20(135n2+627n-1580)a-2(297n+794)a2,r(n,a)=5n(n-2)(50n+11a+220)(20n+10a+143).
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,n,a) and ∂ν~(k,n,a)/∂k|k=n/2 on n for some values of a is shown in Figure 3.Dependence ofν~(1,n,a) and ∂ν~(k,n,a)/∂k|k=n/2 on n for a=10,20,30,40,100. Plots are shown in gray-level intensity, from lower to higher, according to the value of a.
(a)
(b)Again, a fit of these data to the modelsA/(n+B)-1 and A/n yields
(2.7)ν~(1,n,a)=-25n-a-5025(n-1),∂ν~(k,n,a)∂k|k=n/2=96n(a+25).
In addition, we have that
(2.8)ν~(n2,n,2)=0,ν~(n,n,2)=-ν~(1,n,2),
therefore, the coefficients are given by
(2.9)b3(n,a)=q3(n,a)s(n,a),b2(n,a)=q2(n,a)s(n,a),b1(n,a)=q1(n,a)s(n,a),b0(n,a)=q0(n,a)s(n,a),
where
(2.10)q3(n,a)=-200(23n3-92n2+90n+4)+200(n3-5n2+8n-6)a-8(n2-2n+2)a2,q2(n,a)=100(69n4-255n3+186n2+92n+8)-100(3n4-12n3+9n2+4n-12)a+4(3n3-3n2+4)a2,q1(n,a)=50(21n4+54n3-522n2+748n-176)-50(3n4-9n3-6n2+26n-24)a+2(3n3-6n+8)a2,q0(n,a)=-25(25n4-217n3+618n2-660n+184)-25(n4-2n3-7n2+16n-12)a+(n3+n2-4n+4)a2s(n,a)=25(n-1)2(n-2)2(a+25).
Thus, the substitution of (2.10), (2.9), (2.6), and (2.5), respectively, in (2.3) yields the approximate closed-form expressions
(2.11)z~k=x~(k,n,a)+iy~(k,n,a)=2μ~(k,n,a)2n+a-2+i2ν~(k,n,a)2n+a-2,
where k=1,2,…,n, for the zeros of the unnormalized Bessel polynomial yn(z;a). The expressions given in (2.11) converge to the zeros zk of these polynomials, as we will show in the following.


## 3. Convergence

Following [9], we define
(3.1)W(z)=e1+1/z2z(1+1+1/z2),
and denote by Γ the curve defined by
(3.2)Γ={z∈ℂ:|W(z)|=1,|argz|≥π2},
which contains the limit points ω^k of the zeros of the normalized Bessel polynomial yn(2z/(2n+a-2);a). Then, it has been proved in [8] that the zero ωk of yn(2z/(2n+a-2);a) approaches to order O(1/n) the limit value ω^k, that is,
(3.3)|ωk-ω^k|=O(1n),
as n→∞.Thus, if we show that|ω~k-ω^k|=O(1/n), we will have proved that
(3.4)|ωk-ω~k|=O(1n),
and therefore, taking into account that ωk=(2n+a-2)zk/2, the explicit expression (2.11) approaches to order O(1/n2) the zero zk of the Bessel polynomial yn(z;a).To this purpose, we simply substitute the expression forω~k=μ~(k,n,a)+iν~(k,n,a) given by (2.3) in (3.1) to obtain, after a lengthy calculation, that the expansion of W(ω~k) in terms of 1/n is
(3.5)W(ω~k)=1+[2a2-100a(3k-2)-50(42k-67)25(a+25)+h(k,a)25(a+25)(cos2(t(k,a))-sin2(t(k,a)))]1n+O(1n3/2),
where
(3.6)h(k,a)=(25+a)2(130+27a+300k)2+(2a2-100a(3k-2)-50(42k-67))2,t(k,a)=(25+a)(130+27a+300k)4(1675-100a+a2-1050k+150ak).
This implies that
(3.7)|W(ω~k)|=1+O(1n),
for fixed k and a. Thus, ω~k approaches to order O(1/n) the Γ curve and (3.4) follows. From here we have that
(3.8)|zk-z~k|=O(1n2),
as n→∞. Numerical calculations confirm and extend this result. Figure 4 shows the behavior of the maxima of |zk-z~k| over k as they depend on n for the particular case of a=2. The numbers computed by (2.2) are taken as the exact zeros zk. A fit of these data gives 1/na with a=1.7.Figure 4
Plot of the values ofmaxk=1n|zk-z~k| against n for the case of a=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αn(a) of the Bessel polynomial yn(z;a) can be obtained by the substitution of k=(n+1)/2 in the real part of (2.11), x~(k,n). This gives
(4.1)α~n(a)=-4n3+175(50-71a)+(-707a2+5900a+35000)7500n+O(1n2),
as our new result. In [2, 11] very accurate expressions for αn(a) are given. Particularly, the following formula:
(4.2)2αn(a)=-1.325486838n-1.00628995a+1.349836480+O(12n+a-2),
is given in [11]. Expanding 2/α~n(a) in powers of n we find that
(4.3)2α~n(a)≃-1.33333n-0.946667a+0.666667,
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 z~k (cf. (2.11)) gives a complicated expression for the real part. However, expanding both the real and imaginary parts of this sum gives
(4.4)s~1(n)=∑k=1n‍z~k=-1+(53a100-7130+i32a+25)1n+O(1n2).
The exact result is s1(n)=-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 a=1. For this value we obtain
(4.5)s~2(n)=∑k=1n‍z~k2=34695915n+O(1n2).
The exact sum
(4.6)s2(n)=12n-1=12n+O(1n2)
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.


## 4.1. The Real Zero

A closed-form formula for the unique real zeroαn(a) of the Bessel polynomial yn(z;a) can be obtained by the substitution of k=(n+1)/2 in the real part of (2.11), x~(k,n). This gives
(4.1)α~n(a)=-4n3+175(50-71a)+(-707a2+5900a+35000)7500n+O(1n2),
as our new result. In [2, 11] very accurate expressions for αn(a) are given. Particularly, the following formula:
(4.2)2αn(a)=-1.325486838n-1.00628995a+1.349836480+O(12n+a-2),
is given in [11]. Expanding 2/α~n(a) in powers of n we find that
(4.3)2α~n(a)≃-1.33333n-0.946667a+0.666667,
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 z~k (cf. (2.11)) gives a complicated expression for the real part. However, expanding both the real and imaginary parts of this sum gives
(4.4)s~1(n)=∑k=1n‍z~k=-1+(53a100-7130+i32a+25)1n+O(1n2).
The exact result is s1(n)=-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 a=1. For this value we obtain
(4.5)s~2(n)=∑k=1n‍z~k2=34695915n+O(1n2).
The exact sum
(4.6)s2(n)=12n-1=12n+O(1n2)
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 forzk 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 y~(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 complicated.



---

*Source: 873078-2012-10-23.xml*