Application of the Product of Binomials to the Computation of Bessel Functions
Abstract Generated abstract
This paper develops product of binomials approximations for the hypergeometric series 0F1, which forms the basis of several Bessel and related special functions. The method rewrites the logarithmic derivative of the series using a continued fraction, derives explicit numerator and denominator polynomials, and converts the resulting integral representation into an exponential times a finite product with computable exponents. For parameters with nu greater than minus one, remainder estimates are obtained in the complex sector with argument bounded by pi over two. The formulas are then specialized to approximations for Bessel, modified Bessel, Kelvin, and related functions, with numerical examples illustrating the accuracy for selected values.
Full Text
APPLICATION OF THE PRODUCT OF BINOMIALS TO THE COMPUTATION OF BESSEL FUNCTIONS
V. E. KORNILOV
(Presented by the Department of Higher Mathematics of the Tomsk Polytechnic Institute)
In the article, for the series (1), which is a component part of the Bessel functions, approximations in the form of a product of binomials are obtained. These approximations are obtained on the basis of continued fractions and may be applied to the computation of Bessel functions.
An estimate of the remainder term in modulus is obtained in the complex domain \(|\arg z| < \dfrac{\pi}{2}\).
1. Power series
\[ {}_0F_1(\nu+1;z)=\sum_{k=0}^{\infty}\frac{z^k}{(\nu+1)_k k!}, \qquad |z|<N,\quad \nu\ne -1,-2,\ldots; \tag{1} \]
where
\[ (a)_k=\frac{\Gamma(a+k)}{\Gamma(a)}, \qquad (a)_0=1, \]
is transformed as follows:
\[ {}_0F_1(\nu+1;z)=e^{\ln\,{}_0F_1(\nu+1;z)} =\exp\left\{\int_0^z \frac{d[{}_0F_1(\nu+1;t)]}{{}_0F_1(\nu+1;t)}\right\}, \]
i.e.
\[ {}_0F_1(\nu+1;z) = \exp\left[ \int_0^z \frac{{}_0F_1(\nu+2;t)\,dt}{(\nu+1)\,{}_0F_1(\nu+1;t)} \right]. \tag{2} \]
Next we replace the ratio of two series under the integral sign in equality (2) by its known expansion into a continued fraction ([1], p. 136) and represent the power series (1) in the form of the following expansion:
\[ {}_0F_1(\nu+1;z) = \exp\left[ \int_0^z \frac{P_{2n+1}(t)}{Q_{2n+1}(t)}\,dt \right] \exp\left[ -\int_0^z R_{2n+1}(t)\,dt \right], \tag{3} \]
where
\[ \rho_{2n+1}(z)=-\int_0^z R_{2n+1}(t)\,dt, \tag{4} \]
\[ \frac{P_\kappa(t)}{Q_\kappa(t)} = \frac{1}{a_1+\dfrac{1}{a_2+\cdots+\dfrac{1}{a_\kappa}}}, \]
\[ a_1=\nu+1,\qquad a_{2m}=\frac{\nu+2m}{t},\qquad a_{2m+1}=\nu+2m+1,\qquad m=1,2,\ldots \tag{5} \]
\[ R_{2n+1}(t)=\sum_{\kappa=n+1}^{\infty} \frac{a_{2\kappa+1}}{Q_{2\kappa-1}(t)Q_{2\kappa+1}(t)} . \tag{6} \]
For the numerator and denominator of the convergent fraction (5), the following equalities and relation are known ([2], pp. 13, 14):
\[ P_{2n+1}(t)=1+\cdots+\frac{(\nu+2)_{2n}}{t^n} = \frac{(\nu+2)_{2n}}{t^n}\,p_{2n+1}(t), \tag{7} \]
\[ Q_{2n+1}(t)=(n+1)(\nu+n+1)+\cdots+\frac{(\nu+1)_{2n+1}}{t^n} = \frac{(\nu+1)_{2n+1}}{t^n}\,q_{2n+1}(t), \tag{8} \]
\[ \varphi_{2n+1}(z) = \left[ 1+\frac{2z}{(\nu+2n+1)(\nu+2n-1)} \right]\varphi_{2n-1}(z) - \frac{z^2}{(\nu+2n-2)_3(\nu+2n-1)}\, \varphi_{2n-3}(z), \tag{9} \]
where \(p_{2n+1}(z)\) and \(q_{2n+1}(z)=\varphi_{2n+1}(z)\).
The numerator and denominator of the convergent fraction (5) can be expressed explicitly in terms of the parameter \(\nu\) and the variable \(z\), namely:
\[ p_{2n+1}(z) = \sum_{m=0}^{n} \frac{C_{2n-m}^{m}z^m}{(\nu+2)_m(\nu+2n+2-m)_m}, \qquad n=0,1,\ldots; \tag{10} \]
\[ q_{2n+1}(z) = \sum_{m=0}^{n} \frac{C_{2n+1-m}^{m}z^m}{(\nu+1)_m(\nu+2n+2-m)_m}, \qquad n=0,1,\ldots . \tag{11} \]
By direct verification we establish the validity of equality (11) for \(n=0,1\). Suppose that equality (11) is valid for \(q_{2n-3}(z)\) and \(q_{2n-1}(z)\); then, substituting the values of these denominators into the right-hand side of relation (9), we obtain
\[ \begin{aligned} q_{2n+1}(z) &= \left[ 1+\frac{2z}{(\nu+2n+1)(\nu+2n-1)} \right]\times \\[2mm] &\quad \times \left[ \sum_{m=0}^{n-1} \frac{C_{2n-1-m}^{m}z^m}{(\nu+1)_m(\nu+2n-m)_m} \right] - \frac{z^2}{(\nu+2n-2)_3(\nu+2n-1)} \times \\[2mm] &\quad \times \left[ \sum_{m=0}^{n-2} \frac{C_{2n-3-m}^{m}z^m}{(\nu+1)_m(\nu+2n-2-m)_m} \right] \\[2mm] &= 1+\frac{2nz}{(\nu+1)(\nu+2n+1)} + \sum_{m=2}^{n-1} \frac{C_{2n-1-m}^{m-2}z^m}{(\nu+1)_m(\nu+2n-m)_m(\nu+2n+1)} \times \\[2mm] &\quad \times \left\{ \frac{2(2n-2m+1)(n-m)(\nu+2n+1)}{m(m-1)} + \right. \end{aligned} \]
\[ + \frac{(\nu+m)\left[(4n-3m+1)\nu+8n^2-8mn+4n+m^2-2m+1\right]} {(m-1)(\nu+2n)}\right\}+ \]
\[ + \frac{(n+1)z^n}{(\nu+1)_n(\nu+n+2)_n} = \sum_{m=0}^{n} = \frac{C_{2n+1-m}^{m} z^m}{(\nu+1)_m(\nu+2n+2-m)_m}, \]
thereby identity (11) has been proved; equality (10) is proved analogously.
- For \(\nu>-1\) and the independent variable \(z\) in the domain \(|\arg z|\le \dfrac{\pi}{2}\), all \(\alpha_1,\ldots,\alpha_{2m},\ldots\) are positive, and the roots of \(Q_{2n+3}(z)\) are negative and are separated by the roots of \(Q_{2n+1}(z)\) ([2], p. 14); therefore, also in view of equality (8), and ([5], p. 23, (12)),
\[ |q_{2n+1}(z)|<|q_{2n+3}(z)|<\cdots;\quad n=0,1,\ldots \tag{12} \]
Starting from equalities (5), (8), and (12), it is not difficult to obtain an estimate for the functional series (6):
\[ |R_{2n+1}(t)|< \frac{|t|^{2n+1}(\nu+2n+3)_2\left[(\nu+2n+3)_2-|t|^2\right]^{-1}} {(\nu+1)_{2n+1}(\nu+1)_{2n+2}|q_{2n+1}(t)|^2}, \tag{13} \]
where
\[ |t|^2<(\nu+2n+3)_2,\quad |\arg t|\le \frac{\pi}{2}. \]
With respect to the modulus of integral (4), in view of inequality (13), one can establish the upper bound (see [6], (12))
\[ |\rho_{2n+1}(z)|< \frac{|z|^{2n+2}(\nu+2n+3)_2\left[(\nu+2n+3)_2-|z|^2\right]^{-1}} {(\nu+1)_{2n+1}(\nu+1)_{2n+2}|q_{2n+1}(z)|^2}, \]
\[ |\arg z|\le \frac{\pi}{2}. \tag{14} \]
Equality (3), taking into account (7) and (8), can be represented in the form of the following sum:
\[ {}_0F_1(\nu+1;z)= \exp\left[\int_0^z \frac{p_{2n+1}(t)\,dt}{(\nu+1)q_{2n+1}(t)}\right] + \]
\[ +\,r_{2n+1}(z)=F_{2n+1}(z)+r_{2n+1}(z), \tag{15} \]
where
\[ r_{2n+1}(z)=F_{2n+1}(z)\{\exp[\rho_{2n+1}(z)]-1\}. \]
It is not difficult to compute that
\[ |F_{2n+1}(z)|<\left|e^{\frac{z}{\nu+1}}\right|\exp|\rho_1(z)|, \]
therefore
\[ |r_{2n+1}(z)|< \left|e^{\frac{z}{\nu+1}}\right| \exp|\rho_1(z)|\,[\exp|\rho_{2n+1}(z)|-1]. \tag{16} \]
- Theorem. The power series (1) is represented by the following product of binomials:
\[ {}_0F_1(\nu+1;z) = e^{b_0 z}\prod_{m=1}^{n}\left(1+\frac{z}{a_m}\right)^{b_m} +r_{2n+1}(z), \tag{17} \]
where
\[ \left\{ \begin{aligned} \prod_{m=1}^{n}\left(1+\frac{z}{a_m}\right)&=q_{2n+1}(z),\qquad b_0=\frac{1}{(n+1)(\nu+n+1)},\\[4pt] b_m&=\frac{a_m A_{2n+1}(-a_m)} {(\nu+1)^2(\nu+2)\dfrac{d}{dz}\,[q_{2n+1}(-a_m)]},\qquad n=1,2,\ldots; \end{aligned} \right. \tag{18} \]
moreover,
\[ A_{2n+1}(z)=\sum_{m=0}^{n-1} \frac{C_{2n-1-m}^{m}z^m}{(\nu+3)_m(\nu+2n+2-m)}. \tag{19} \]
The polynomial \(q_{2n+1}(z)\) is determined according to equality (11).
An upper bound for the modulus of the remainder term is established by means of inequalities (16) and (14).
Proof. Applying the polynomial
\[ M_{2n+1}(t)=\frac{1}{\nu+1}p_{2n+1}(t) -\frac{1}{(n+1)(\nu+n+1)}q_{2n+1}(t) \tag{20} \]
and equalities (7), (8), and (18), we transform equality (15):
\[ {}_0F_1(\nu+1;z)\exp\left[ \frac{1}{(n+1)(\nu+n+1)}\int_{0}^{z}dt+ \int_{0}^{z}\frac{M_{2n+1}(t)\,dt}{q_{2n+1}(t)} \right]+ \]
\[ +\,r_{2n+1}(z)= \exp\left( b_0\int_{0}^{z}dt+ \sum_{m=1}^{n}\frac{b_m}{a_m} \int_{0}^{z}\frac{dt}{1+\dfrac{t}{a_m}} \right)+r_{2n+1}(z), \]
where
\[ b_0=\frac{1}{(n+1)(\nu+n+1)},\qquad b_m=\frac{M_{2n+1}(-a_m)} {\dfrac{d}{dz}\,[q_{2n+1}(-a_m)]}, \tag{21} \]
and also
\[ {}_0F_1(\nu+1;z)=e^{b_0z}\prod_{m=1}^{n} \left(1+\frac{z}{a_m}\right)^{b_m}+r_{2n+1}(z), \tag{17} \]
besides, on the basis of (15) and (17), the upper bound for \(|r_{2n+1}(z)|\) is found according to inequalities (16) and (14).
In view of the first equality (18), \(q_{2n+1}(-a_m)=0\); therefore, applying equality (20), we obtain
\[ b_m= \frac{ M_{2n+1}(-a_m)+ \left[ \dfrac{1}{(n+1)(\nu+n+1)}-\dfrac{1}{\nu+1} \right]q_{2n+1}(-a_m) }{ \dfrac{d}{dz}\,[q_{2n+1}(-a_m)] } = \]
\[ = \frac{ \dfrac{1}{\nu+1}\,[p_{2n+1}(-a_m)-q_{2n+1}(-a_m)] }{ \dfrac{d}{dz}\,[q_{2n+1}(-a_m)] }. \]
and, applying equalities (10) and (11), we obtain the following equality
\[ \frac{1}{\nu+1}\{p_{2n+1}(z)-q_{2n+1}(z)\}= \]
\[ = \frac{-z}{(\nu+1)^2(\nu+2)} \sum_{m=0}^{n-1} \frac{C_{2n-1-m}^{m} z^m} {(\nu+3)_m(\nu+2n+2-m)_m} = \frac{-z A_{2n+1}(z)} {(\nu+1)^2(\nu+2)}, \]
from which equality (19) follows, as well as the last of equalities (18); thus the theorem is completely proved.
The coefficients \(b_m\) can also be computed from formula (21); then, according to equalities (10), (11), and (20), we compute
\[ M_{2n+1}(z)=\sum_{m=0}^{n-1}\times \]
\[ \times \frac{(n-m)[(2n+1)\nu+2n^2+5n+2-m]\,C_{2n-m}^{m} z^m} {(2n-2m+1)(n+1)(\nu+n+1)(\nu+1)_{m+1}(\nu+2n+2-m)_m}, \tag{22} \]
in particular,
\[ M_1(z)=0,\qquad M_3(z)=\frac{\nu+3}{2(\nu+1)(\nu+2)}, \]
\[ M_5(z)=\frac{2(\nu+4)}{3(\nu+1)(\nu+3)} +\frac{(\nu+4)(5\nu+19)}{3(\nu+1)_5}\,z, \]
\[ M_7(z)=\frac{3(\nu+5)}{4(\nu+1)(\nu+4)} +\frac{(\nu+3)(7\nu+34)}{2(\nu+1)_4(\nu+7)} +\frac{(\nu+5)(7\nu+33)z^2}{2(\nu+1)_7}. \]
- The relation between the Bessel functions and the power series (1) is expressed by the following formulas ([3], pp. 138–149), ([4], p. 150):
\[ J_\nu(x)=\frac{1}{\Gamma(\nu+1)} \left(\frac{x}{2}\right)^\nu \,{}_0F_1\left(\nu+1;-\frac{x^2}{4}\right), \tag{23} \]
\[ I_\nu(x)=\frac{1}{\Gamma(\nu+1)} \left(\frac{x}{2}\right)^\nu \,{}_0F_1\left(\nu+1;\frac{x^2}{4}\right), \tag{24} \]
\[ Y_\nu(x)=\frac{1}{\sin \nu\pi} [\cos \nu\pi\,J_\nu(x)-J_{-\nu}(x)], \qquad \nu\ne 0,\ \pm 1,\ldots; \tag{25} \]
\[ K_\nu(x)=\frac{\pi}{2\sin \nu\pi} [I_{-\nu}(x)-I_\nu(x)], \qquad \nu\ne 0,\ \pm 1,\ldots; \tag{26} \]
\[ \operatorname{bei}_\nu(x)=\operatorname{Im}[J_\nu(i\sqrt{i}\,x)], \tag{27} \]
\[ \operatorname{ber}_\nu(x)=\operatorname{Re}[J_\nu(i\sqrt{i}\,x)]. \tag{28} \]
Applying formula (17) in the case \(n=1,2\), we obtain the following approximate formulas for the power series (1):
\[ {}_0F_1(\nu+1;z) = e^{\frac{z}{2(\nu+2)}} \left[ 1+\frac{2z}{(\nu+1)(\nu+3)} \right]^{\frac{(\nu+3)^2}{4(\nu+2)}}+r_3(z), \tag{29} \]
\[ {}_0F_1(\nu+1;z) = e^{\frac{z}{3(\nu+3)}} \left( 1+ \frac{ 2-\sqrt{\frac{\nu^2+6\nu+17}{\nu^2+6\nu+8}} } {\nu^2+6\nu+5}\,z \right)^{b_1} \times \]
\[ \times \left( 1+ \frac{ 2+\sqrt{\frac{\nu^2+6\nu+17}{\nu^2+6\nu+8}} } {\nu^2+6\nu+5}\,z \right)^{b_2} +r_5(z), \tag{30} \]
where
\[ b_{1,2}=\frac{(\nu+4)\left[(5\nu+19)\sqrt{\nu^2+6\nu+17}\pm4(\nu+2)\sqrt{\nu^2+6\nu+8}\right]}{18(\nu+3)\sqrt{\nu^2+6\nu+17}}. \tag{31} \]
For the functions (27) and (28), we obtain approximations according to formulas (23) and (29):
\[ \operatorname{bei}_{\nu}(x)= \frac{1}{\Gamma(\nu+1)} \left(\frac{x}{2}\right)^{\nu} \left\{ \left[1+\frac{x^4}{4(\nu+1)^2(\nu+3)^2}\right]^{\frac{(\nu+3)^2}{8(\nu+2)}} \times \right. \]
\[ \left. {}\times \sin\left[ \frac{x^2}{8(\nu+2)}+\frac{3\nu\pi}{4} +\frac{(\nu+3)^2}{4(\nu+2)} \operatorname{arctg} \frac{x^2}{2(\nu+1)(\nu+3)} \right] + \right. \]
\[ \left. {}+\operatorname{Im}\left[ r_3\left(i\frac{x^2}{4}\right) \right] \right\}, \tag{32} \]
\[ \operatorname{ber}_{\nu}(x)= \frac{1}{\Gamma(\nu+1)} \left(\frac{x}{2}\right)^{\nu} \left\{ \left[1+\frac{x^4}{4(\nu+1)^2(\nu+3)^2}\right]^{\frac{(\nu+3)^2}{8(\nu+2)}} \times \right. \]
\[ \left. {}\times \cos\left[ \frac{x^2}{8(\nu+2)}+\frac{3\nu\pi}{4} +\frac{(\nu+3)^2}{4(\nu+2)} \operatorname{arctg} \frac{x^2}{2(\nu+1)(\nu+3)} \right]\times \right. \]
\[ \left. {}+\operatorname{Re}\left[ r_3\left(i\frac{x^2}{4}\right) \right] \right\}, \tag{33} \]
where, for \(\nu=0\) and \(x=2\), we obtain
\[ \operatorname{bei}_0(2)\approx 0,97205\;(\Delta=0,00024),\quad \operatorname{ber}_0(2)\approx 0,75332\;(\Delta=0,00058), \]
\(\Delta\) is the absolute error.
We apply formulas (24), (29) to the function (26) and pass to the limit as \(\nu\to0\); as a result we obtain:
\[ K_0(x)\approx e^{\frac{x^2}{16}} \left(1+\frac{x^2}{6}\right)^{\frac{9}{8}} \left[ \psi(1)-\ln\frac{x}{2}+\frac{x^2}{32} \right] - \]
\[ {}-\frac{3}{16}\ln\left(1+\frac{x^2}{6}\right) +\frac{6x^2}{4(6+x^2)} \Bigg], \quad \psi(1)=-0,5772157, \tag{34} \]
where
\[ K_0(1)\approx 0,421049\;(\Delta=-0,000025),\quad K_0(2)\approx 0,1186\;(\Delta=-0,0047). \]
Applying formula (30) to the function (24), we obtain, for \(\nu=0\),
\[ I_0(x)= e^{\frac{x^2}{36}} \left(1+\frac{8-\sqrt{34}}{80}x^2\right)^{\frac{64+38\sqrt{34}}{27\sqrt{34}}} \left(1+\frac{8+\sqrt{34}}{80}x^2\right)^{\frac{64-38\sqrt{34}}{27\sqrt{34}}} + \]
\[ {}+r_5\left(\frac{x^2}{4}\right). \tag{35} \]
We have the table
| \(x\) | \(I_0(x)\) | (35) |
|---|---|---|
| 4 | 11,3019 | 11,3045 |
| 8 | 427,6 | 443,86 |
References
-
A. N. Khovanskii. Application of continued fractions and their generalizations to problems of approximate analysis. Moscow, GITTL, 1956.
-
T. I. Stieltjes. Investigations on continued fractions. Moscow, ONTI, 1936.
-
N. N. Lebedev. Special functions and their applications. Moscow, GITTL, 1953.
-
V. A. Ditkin, A. P. Prudnikov. Integral transforms and operational calculus. Moscow, Fizmatgiz, 1961.
-
V. E. Kornilov. Application of continued fractions to the computation of integrals of binomial differentials. “Izv. TPI,” vol. 131, pp. 21–25, 1965.
-
V. E. Kornilov. Application of the product of binomials to the computation of the Gaussian function. Izv. TPI, vol. 205, 1969.