Computation of the confluent hypergeometric function U(a,b,x) and its derivative for positive arguments

An algorithm and a MATLAB implementation for computing the Kummer function U(a,b,x) and its derivative is given in this paper. The algorithm is efficient and accurate. Numerical tests show that the MATLAB algorithm allows the computation of the function with ∼10−14\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\sim 10^{-14}$\end{document} relative accuracy in the parameter region (a,b,x) ∈ (0,500) × (0,500) × (0,1000) in double-precision floating point arithmetic.


Introduction
Many applications in physics and engineering involve the computation of confluent hypergeometric functions or some of their particular cases (Airy and Bessel functions, Laguerre polynomials, parabolic cylinder functions etc.) See [1,2,8,12] or [11, §13.28] for some examples of applications in physics. As an additional example, in recent work we have obtained expansions for the relativistic Fermi-Dirac integral and its derivatives (of great importance in stellar astrophysics) where confluent hypergeometric functions play a key role [3,5]. Amparo Gil amparo.gil@unican.es 1 In spite of their importance, few algorithms are available in double-precision floating point arithmetic for the computation of any of the standard solutions of the Kummer's equation in the case of real or complex parameters. For the function M(a, b, x), some algorithms can be found in the literature (see [9,10]); however, for the second solution U (a, b, x), no double-precision algorithms with good uniform accuracy seem to be available, as far as we know. For example, the SLATEC library [16] includes the Fortran 77 function dchu.f which implements the confluent Kummer function U(a, b, x) for real parameters a, b and argument x. However, one has to be careful when using dchu.f to compute U(a, b, x) when a or b are large in comparison to x because a significant loss of accuracy occurs in the computations. Another example of a double-precision implementation for the Kummer U(a, b, x) function which is not free of accuracy problems is the function hyperu included in the SCIPY [7] module scipy.special. Having accurate efficient doubleprecision implementations is of interest in, for example, applications where a large number of function evaluations are needed (such as in the condensed matter or astrophysical calculations). On the other hand, extended precision implementations are given in Mathematica, Maple or Arb [6].
In this paper, we describe an efficient and accurate algorithm implemented in MATLAB for computing the Kummer U(a, b, x) function and its derivative. The function U(a, b, x) satisfies the following integral representation valid for a > 0, b ∈ C and x > 0.
In the algorithm, we consider real positive values of the function arguments a, b, x. The algorithm uses series for small values of the parameters, asymptotic expansions for large values of the parameters and a double-precision implementation of the algorithm described in [14] in the rest of cases. Our numerical tests show that the resulting algorithm is both accurate and efficient.

Computation for small values of the parameters
For 0 < a, b < 1 2 and 0 < x < 1, we use the series given in [4] for the Kummer function U(a, b, z) and its derivative where the coefficients w d m can be written in terms of the w m coefficients The coefficients w m satisfy To compute w m for m = 1, 2, ... using (2.4), it is convenient to obtain a stable recursion for u m when b is small. It is possible to verify that The starting u 0 value of the recursion given in (2.5) is computed as sin(πb) and for calculating w 0 , we use For computing G(a, b), we use the series  Table 1 of [4].

Algorithm based on the use of recurrences
One of the key ingredients in the algorithm described in [14] to evaluate U(a, b, x) is the use of three-term recurrence relations satisfied by the function with respect to the parameters a and b.
With respect to the a parameter, we have and with respect to the b parameter, As discussed for example in [13], the Kummer function U(a, b, x) is the minimal solution of the three-term recurrence relation (3.1). Therefore, the forward computation of the recursion is ill-conditioned and backward recursion should be applied. On the contrary, U(a, b, x) is a dominant solution of (3.2) and forward recursion is possible.
Briefly, the algorithm described in [14] uses a sophisticated version of the backward a-recursion to compute U(a, b, x) for 0 < b < 1. This recursion is computed using starting values given by asymptotic expansions in terms of Bessel functions (for small x) or Miller algorithm's when x is not small.
2) is applied. Technical details can be seen in [14]. This algorithm works very well in extended precision; however, in double-precision floating point arithmetic, some loss of accuracy could appear when using the recursions for large values of the parameter a. This is the reason why it is better to use an alternative method of computation for large values of the parameters, such as the expansions described in the next section.

Asymptotic expansions for large values of the parameters
We consider the recent expansion given in [15] where with τ = t 0 /μ, β = b/z and μ = (b − a)/z.

Expansion for U (a, b + 1, z)
For obtaining the asymptotic expansion of U (a, b + 1, z) we consider the relation and use the integral (see (4.3) of [15]) We write this in the form where, with μ = λ/z, The function φ(t) is the same as used in [15, §3], and the saddle point t 0 and the coefficients are the same as used in that section, up to a factor (−1) n . We obtain the expansion where and (4.10) The saddle point t 0 is given by with expansion (4.12) The first few f n coefficients are f 0 = 1 and (4.13)

Numerical testing and algorithm
The set of MATLAB functions implementing the methods used are: 1. Uabxsmall(a,b,x). Implementation of the method given in Section 2.
2. Uabxrec(a,b,x). Implementation of the method given in Section 3.
3. Uabxexpan(a,b,x). Implementation of the method given in Section 4.
For testing the computation of U (a, b, x) and U (a, b, x) using the different methods, we consider the recurrence relations Also, for testing the method in Section 2, we consider the following alternative relation When using this relation, values a ∼ b should not be included to avoid the loss of significant digits by cancellation. We have randomly generated 10 6 values in the parameter region (a, b, x) ∈ (0, 0.5) × (0, 0.5) × (0, 1) for testing (5.2) using Uabxsmall (a,b,x). The results obtained show that for more of the ∼ 98% of the points, the obtained accuracy was better than 10 −14 . Also, for the same parameter region, this method is more than a 20% faster than the method implemented in Uabxrec(a,b,x).
In Figs. 1 and 2, we show the accuracy obtained when using the expansion (4.1) and (4.8) implemented in the MATLAB function Uabxexpan(a,b,x) to approximate the values of U(a, b, x) and U (a, b, x). We take n = 5 terms in both  Taking into account the results obtained in the different tests of the methods, the resulting computational scheme is given in Algorithm 1. This algorithm is implemented in the MATLAB function Uabxcomp. 1 The front factors in the expansions (4.1) and (4.8) are used in the full algorithm as an estimation of the function value to avoid possible overflow/underflow problems.
For testing the accuracy of the full algorithm for computing U (a, b, x), we have considered 10 8 random values in the parameter region (a, b, x) ∈ (0, 500) ×  (0, 500) × (0, 1000). Of all the randomly generated points, 2.64 × 10 7 were points not in the underflow or overflow regions. The results show that at approximately 54% of these points, the accuracy is better than 10 −14 , while at around 43% of the points the accuracy is between 10 −14 and 10 −13 . The maximum error obtained was ∼ 10 −11 , at a point for which the function value was close to the underflow limit. See also Fig. 3. Similar results were obtained for the derivative. U(a, b, x) and its derivative.

Algorithm 1 Computation of the Kummer function
As an additional test, we have compared our function Uabxcomp against the intrinsic MATLAB function to evaluate the Kummer U(a.b, x) function which is called kummerU. This function, included in the symbolic math toolbox, provides also floating-point results for numeric arguments. When checking the accuracy obtained in comparison to our function, the tests reveal that for some parameters, the computations obtained with kummerU are not accurate. For example, in Fig. 4, we show that the calculation of U(a, b, x) using kummerU in MATLAB R2022b gives a negative number for certain values of the parameters (U(a, b, x) is non-negative for x > 0, a ≥ 0, b ∈ ). With respect to computational efficiency, the computations with kummerU take much more time than those obtained with our algorithm, but this is expected since kummerU is part of a symbolic package. As an example, the computation of 500 function values with random arguments in the parameter region (a, b, x) ∈ (0, 500) × (0, 500) × (0, 1000) took 9.4 10 −2 s using Uabxcomp and 199 s using kummerU. The results have been obtained using MATLAB R2022b in a PC with 8 GB RAM and Intel Core i5-4310U processor.  U(a, b, x) using the intrinsic function in MATLAB R2022b gives a negative number