Quasi-Rational Analytic Approximation for the Modiﬁed Bessel Function I 1 ( x ) with High Accuracy

: A new simple and accurate expression to approximate the modiﬁed Bessel function of the ﬁrst kind I 1 ( x ) is presented in this work. This new approximation is obtained as an improvement of the multi-point quasi-rational approximation technique, MPQA. This method uses the power series of the Bessel function, its asymptotic expansion, and a process of optimization to ﬁt the parameters of a ﬁtting function. The ﬁtting expression is formed by elementary functions combined with rational ones. In the present work, a sum of hyperbolic functions was selected as elementary functions to capture the ﬁrst two terms of the asymptotic expansion of I 1 ( x ) , which represents an important improvement with respect to previous research, where just the leading term of the asymptotic series was captured. The new approximation function presents a remarkable agreement with the analytical solution I 1 ( x ) , decreasing the maximum relative error in more than one order of magnitude with respect to previous similar expressions. Concretely, the relative error was reduced from 10 − 2 to 4 × 10 − 4 , opening the possibility of applying the new improved method to other Bessel functions. It is also remarkable that the new approximation is valid for all positive and negative values of the argument.


Introduction
The Bessel equation and its modified form play an important role for modeling common problems in science and engineering. These equations appear in different contexts spanning areas like electrodynamics [1,2], plasma-physics [3][4][5], heat transfer [6][7][8], chemical engineering [9], fluid mechanics [10,11], and others [12,13]. Bessel functions with integer order appear in systems with cylindrical symmetry, and those with half integer order in spherical symmetries. The first kind solution of the modified Bessel equation is called I ν (x), and it corresponds to a power series which is used to obtain the solution of a variety of problems. However, the numerical calculation for large and intermediate values of x requires a large number of terms of the series to get good accuracy, which has impulsed a range of different techniques for computing the modified Bessel functions properly in a reasonable time frame. In general, these techniques pursue two kinds of goals. A first group focuses on formulating a simple approximation of the modified Bessel function with a low number of fitting parameters. They can be used in correlations originated from the Bessel equation or from more complex equations where the Bessel equation is involved. These approximations allow a simple and direct computation of the solution, without using the modified Bessel function itself. An example of this kind of technique is the model proposed in [14], where the modified Bessel function is approximated by a finite sum of exponential functions, which considers eight fitting parameters established for different intervals of the domain. A more extensive analytic approximation for the specific case of I 1 (x) has also been published in [15]. The method to obtain this approximation was based on the so-called multi-point quasi-rational approximation technique (MPQA) [16]. That approximation is valid for any positive value of x, and its maximum relative error was about one percent, considering just five parameters for the entire domain.
On the other hand, there is a second group of approximation techniques focused exclusively on the accuracy of the fitting function, which implies a larger number of fitting parameters, truncated series, or algorithms based on conditional loops. In this case, the Bessel approximation can be coded in a programming language to be used in computers as a function that can be easily called by a user. An important initial attempt to achieve a high level of accuracy was made by polynomial approximations in [17], and it was coded in different programming languages [18]. Here, eight or ten fitting parameters are used in an expression for I 1 (x), depending on the argument interval. The fitting function resembles the power series and the asymptotic expansion of the modified Bessel function, for low and large values of the argument, respectively. A simple method in this group considers to compute the power series and the asymptotic expansion of the modified Bessel function directly, defining a criterion for limiting the scope of each series and the number of terms to be used. This was proposed in [19], where a code written in FORTRAN language is presented for computing I ν (x). A remarkable fast technique to compute the modified Bessel function with a high level of precision is presented in [20]. This technique is based on an algorithm that evaluates the ratio of the (k + 1)st term to the kth term of the power series of I ν (x) in a loop. The iterative code is very simple and short, and surprisingly faster than some well-established software like MATHEMATICA and MAPLE in several orders of magnitude for large arguments, which was verified in the same work by running time experiments made in a common laptop.
In the present work, a simple and very precise analytic approximation for I 1 (x) has been determined using only six fitting parameters. The methodology corresponds to the MPQA as in [15], but the structure of the new approach differs substantively from the previous one. Now, the new structure is built in such a way that includes more information from the asymptotic expansion of I 1 (x). Specifically, a sum of two hyperbolic functions (combined with rational ones) was introduced. The direct effect of this change on the new approach was the capture of the two first leading terms of the asymptotic expansion, producing a strong increment in accuracy for intermediate values of the x variable. The new approximation is valid for any real value of x, positive or negative.
The theoretical treatment of this study, including the equations to determine the fitting parameters of the analytic approximation of I 1 (x), will be presented in Section 2. The results and discussion will be performed in Section 3. An analysis to determine the convergence of the new approximation to the exact function has also been included at the end of Section 3. Finally, Section 4 will be devoted to the conclusions of this work.

Theoretical Treatment
The general procedure to find an approximate expression for the modified Bessel function consists of: (a) selecting a suitable structure for a fitting function; (b) to establish the constraints on the fitting parameters; (c) to write these parameters as a function of just one parameter, which is taken as an independent variable; and (d) to optimize the fitting parameters minimizing the error.

Fitting Function
The modified Bessel function of order ν and its asymptotic expansion are some well known series [17]; they are given respectively by: If ν is a positive integer, the Gamma function Γ can be evaluated by using the rule: Γ(n + 1) = n!. Specifically, for ν = 1, the following series are obtained: An interesting challenge is given by building a fitting function that captures the behavior of the modified Bessel series for small values of x, and the behavior of the asymptotic expansion for large values of x. The complexity for the conciliation of these limiting behaviors is rooted in the fact that the modified Bessel series is formed by even powers, while its expansion has inverse powers of x and also a prefactor e x / √ x. A simple option to achieve this corresponds to the following approximation function: where λ, p 0 , p 1 , p 2 , p 3 , and q correspond to fitting parameters. The functionĨ 1 (x) captures the prefactor e x / √ x and the first two terms of the asymptotic expansion for big values of x through: (1) the exponential behavior of sinh(x) and cosh(x), and (2) the limits of the different components ofĨ 1 (x) when x → ∞. At the same time, this analytic approximation can be equalized for small values of x with the first tree terms of the modified Bessel function, as it will be described in the next sub-section. The parameter λ constitutes the optimization variable, which means that all the others parameters depend on λ for this purpose. Thus, λ will vary freely to find the optimal fitting functionĨ 1 (x) for which the error will be minimized.
The approximation function shown in Equation (5) also allows for being expanded to capture more terms of the asymptotic series of I 1 (x). Indeed, including successive sums of the hyperbolic functions sinh(x) and cosh(x) multiplied by suitable polynomial fractions, it is possible to capture more terms of Equation (4) for large values of the argument x.
It is important to point out some features ofĨ 1 (x). A remarkable characteristic of the functionĨ 1 (x) is that it is an odd function. The power series shows this symmetry, but this characteristic is broken for the asymptotic expansion. The new functionĨ 1 (x) keeps the symmetry of the power series and also represents two terms of the asymptotic expansion without breaking this symmetry. This is done by obtaining the exponential behavior of expansion by means of x cosh x and sinh x. In the previous work presented in [15], this symmetry is also kept, but only the leading term of the asymptotic expansion is captured. In the present work, the leading term and the next one are used, but additionally more asymptotic terms can be incorporated by increasing the number of even powers in the numerator and denominator.

Fitting Parameters
The fitting parameters will be subjected to different constraints to mimic the behavior of the modified Bessel function I 1 (x) and its corresponding asymptotic expansion. In addition, a condition to guarantee continuity by avoiding singularities will be imposed.

Constraints for Large Values of x
By equating the first two terms of the Equation (4) with the approximation function given in (5) for high values of the variable x, we arrive to: The following expressions are obtained from Equation (6): It is important to emphasize that p 2 is coming from sinh x and p 3 from cosh x. If only one hyperbolic function is used, then only the information coming from the leading term of the asymptotic expansion is considered. To introduce the information from two consecutive terms of the asymptotic expansion, two hyperbolic functions have to be used.

Constraints for Low Values of x
For small values of x, it is possible to use the Taylor series around x = 0 to approximate the expressions contained inĨ 1 (x); this is: Using the Equations (9)- (11), and then matching the approximation functionĨ 1 (x) with the first two terms of the Bessel series, the following polynomial equation is obtained: Truncating the terms of degree 6 or bigger in Equation (12) and equating the polynomial coefficients, it is relatively straightforward to arrive to the following equation system: p 0 120 The equation system for the fitting parameters obtained before can by solved by adding Equations (7) and (8) from the analysis for large values of x. As we said before, the parameter λ will be a free parameter, and then all the rest of parameters will be a function of λ.

Non-Zero Point Constraint
An important requirement forĨ 1 (x) corresponds to avoid artificial zero values in the denominator, which are called defects according to the well known Padé Approximant technique. These points can be removed by restricting q to a positive value. The expression for q corresponds to: where λ is considered to be positive. Figure 1 shows the parameter q as a function of λ. Note that q is positive between the two vertical asymptotes and for values of λ located a little after the asymptote placed in λ = 0.7. In addition, note that all q values that are positive are suitable to use, but values of q (and λ) very close to a vertical asymptote should be determined considering many significant figures due to the big slope in these zones. In the present work, we consider a limit of four significant figures for all the parameters to keep the simplicity of the fitting function.

Results and Discussion
The criterion for determining λ is by comparing the maximum relative error of the fitting function for different λ values. In this sense, two kinds of error definitions were established: (1) an error in a point x given the λ value ε p (x, λ), and (2) a global error ε(λ) evaluated on an interval [a, b] of the x variable. This last error is given by the maximum error ε p (x, λ) in the interval [a, b]. The expressions for the error (1) and (2) are: Figure 2 shows the error ε in the interval x ∈ [0, 500] as a function of λ. In Figure 2a, four local minimum points of the ε(λ) function can be localized. Figure 2b zooms in the zone where the absolute minimum value of the error is located. This absolute minimum value called ε min corresponds to 0.0003938, and it occurs for λ min = 0.4800. Note that, when λ = λ min , the value of the parameter q corresponds to 1.297, which satisfied the non-singular point constraint. Then, the optimal values of each parameter of the fitting function can be determined, and they are the ones shown in Table 1.  Using the values from Table 1, the analytic fitting function is expressed as: Figure 3a shows the relative error of the fitting functionĨ 1 (x) (Equation (19)). The maximum error, which also corresponds to ε min , is 0.0004 or 0.04%. This maximum error is reached for x ≈ 14, and it is more than one order of magnitude lower than previous research, where the error reaches a value of 1 percent [15]. Additionally, in Figure 3b, the analytic approximation functionĨ 1 (x) is shown in comparison with the actual modified Bessel function, across the zone where the error reaches its maximum value of 0.0004. Note the strong similarity between both functions. It is interesting to point out that the relative error of the approximation decreases with x, which shows the convergence between the approximation and the real function for large values of x. This is a consequence of the convergence between the asymptotic behaviors of the two functions. On the other side, there is also a convergence in the values near zero. The points with the higher errors will be always for intermediate x values.
It is important to note that the approximation here obtained is an odd function of x, and therefore is good for positive as well as negative values of x. To clarify more this matter, it is useful to compare the present approximation with that in Ref. [17], where two approximations have been determined, one valid for x in the interval −3.75 ≤ x ≤ 3.75 and the second one for positive x > 3.75. The total number of parameters to be determined there are 16, six for the first interval and ten for the second one. Now, in the present work, only one approximation has been determined valid for every value of x, positive or negative, using only six parameters. In our analysis, relative errors are analyzed instead of absolute ones because, in most of applications, the function I 1 (x) appears as a factor, and, for these cases, the relative error is the simplest option to be evaluated.
The approximations by finite sum of exponential functions are also important in [14]. The characteristics of these approximations are that each one is limited to an interval. Using four intervals, they obtain fitting functions for all positive values of x. However, the maximum relative error corresponds to 0.002 at x = 0, which is larger by a factor of five than the one here presented with only one function, and also valid for any positive and negative value of x.
As it has been mentioned before, there are two different goals when an approximation of the modified Bessel function is built: simplicity or extreme accuracy. Our goal was to produce a simple expression suitable for using in correlations or rapid computations with a calculator, but, at the same time, providing an excellent level of accuracy. Methods based on loops or finite series achieve in general more accuracy, but they need to be coded, and they also have a computational cost intrinsic, which is negligible for our new approximation function. The relative error in a loop routine can be of order 10 −15 , but the running time in a common two core laptop can fluctuate from a fraction of second to several minutes to compute I 1 (x) for large values of x, depending on the routine method [20].
To study the convergence of the present approximationĨ 1 (x) with respect to the exact function, the following analysis has been performed. First, the asymptotic behavior of I 1 (x) is analyzed. After that, we look for the differences with the asymptotic expansion of I 1 (x). Finally, the relative error ε p (x) for large values of x is determined, showing the actual convergence of ε p (x) for large values of x.
4·0.4800 3 1.297 0.2289 − 0.08585 It is clear from this result that the relative error for large values of x will be 2.8 · 10 −5 , lower than the maximum relative error here found of 0.0004. In the present approximation, all the parameters are given with four decimals since that is enough to keep the largest relative error as low as 0.0004. The calculation we have done shows that, if the parameters are given with more digits, then the relative error for large x can be decreased, but the maximum relative error will be the same. In case that, for any reason, a researcher needs higher approximation for large values of x, the parameters given in Table 1 should be improved to obtain an error of O 1 x or even of O 1 x 2 .

Conclusions
A simple and accurate analytic approximation for the first kind modified Bessel function I 1 (x) has been determined in the present work. This approximation is much better than in a previous version [15], since the relative error is more than one order of magnitude lower than the previous one, decreasing from 1% to 0.04%, and also lower (by a factor of five) than other method of the same category based on exponential functions. This substan-tial improvement is rooted in the new structure of the fitting function, which captures two terms of the asymptotic series of I 1 (x), instead of just one as in the previous approximation version. The structural change in the new fitting function consists of the incorporation of the sum of two hyperbolic functions weighted through rational functions, which can be extended to capture more terms of the asymptotic expansion of I 1 (x). An important point is that the new improvement can be applied to all Bessel functions, opening the door to new studies and perspectives. The approximation here presented is valid not only for any positive values of x, but also for negative ones. It can also be differentiated and integrated as the real I 1 (x), since only one fitting function is good for all the values of the argument. Furthermore, the relative error of the new approximation is about 4 · 10 −4 , which will be right for most of the applications. In addition, the calculation can be performed with a simple pocket calculator. Finally, one last analysis has also been performed showing the convergence ofĨ 1 (x) to I 1 (x) for large values of x. It is clear that the relative error decreases for large values of x, and the maximum relative error is always for intermediate values of the variable x.