New Approach for Calculate Exponential Integral Function

This manuscript presents a new approach to accurately calculating exponential integral function that arises in many applications such as contamination, groundwater flow, hydrological problems and mathematical physics. The calculation is obtained with easily computed components without any restrictive assumptions . A detailed comparison of the execution times is performed. The calculated results by the suggested approach are better and faster accuracy convergence than those calculated by other methods. Error analysis of the calculations is studied using the absolute error and high convergence is achieved. The suggested approach out-performs all previous methods used to calculate this function and this decision is based on the getting results


Introduction
The exponential integral is an important integral because it is used in many physical problems such as groundwater flow models, reactor physics and others [1]. In addition, many integrals and complex functions can be represented by exponential integrals. The important ISSN: 0067-2904 point is how to estimate this function which denoted (well function), W(u) = -Ei(-u). Theis [2] was the beginning of the appearance of exponential integration in groundwater models. Properly calculating the well function is very important in the calculations of the pumping tests which is used to determine the coefficients of confined or leakage aquifers. Formerly, the well function are based on the approximate values in the form of tables for small values or used numerical integrations. Recently, approximations of the well function are valid for the entire range of arguments that have been proposed for both confined and leaky aquifers by Swamee and Ojha see [3][4][5]. Finding an analytical solution for exponential integration is not possible, so numerical methods must be used. Our purpose in this paper is to find a reliable and easy implementation approximate for E1 (u). This approximate is appropriate for the aquifer systems, such approximate will be involved only common mathematical functions and be less computationally intensive to evaluate as well as simple and quick to program, it will also involve errors that are less than would occur in hydrological data sets such that using the formula does not consist any additional errors.
The importance of the calculated exponential integral function is appeared in solving Theis equation which studied the groundwater flow model in confined aquifers in the unsteady state case and it is defined as follows [6][7][8]: where Q is discharge or volume flow rate, Tr is the transmissivities in the r direction.
This article has been arranged as follows: In section 2, we present the exponential integral function. In section 3, we briefly review the most approximate methods. A suggested approach with implementation and discussions for the results will be given in section 4. Finally, the conclusion is given in section 5.

Exponential Integral Function
The exponential integral function is also called Theis function or First order E1(u) which is applying to confined aquifers, it is denoted by W(u) and defined as follows [9][10][11]: Where u is a dimensionless variable and y is a dummy variable. The integral in eq. (1) is a function of the lower integration limit u. This function tends to be infinity, thus it fails to be defined when u approaches zero. The integral in eq. (1) can be expanded in a convergent series that is given in most groundwater textbooks.
As we mentioned earlier, the well function cannot be analytically solved. Therefore, different methods and techniques have been used to approximate this function which varies according to the required of accuracy. For u < 1, we can use the Maclaurin series to the exponential function in equation (1) and then integrate this series term by term to get: where  = 0.5772156649015328606 is the Euler constant. There are some efficient approaches that have been successfully used to calculate Euler constant such as those was suggested by Lebedev (1965) see [14][15][16] and Andrews (1985) see [17][18][19]. But the accurate value of this constant is up to 50 decimal places proposed by Murnaghan and Wrench (1963) see [20][21][22][23].
Here by using MATLAB code, the authors calculate the value of well function for any number of the series terms in equation (2). In the interval (0, 1], it is found that the values of the function are matching with 50 decimal places when calculated for 20 or more terms of the series in equation (2). Table 1 illustrates the matching of the values when calculating the function with 20-terms and 200-terms for different values within the interval (0, 1]. Therefore, we will approximate the following function instead of approximate equation (2)

Approximate Methods
There is various approximation methods have been used and developed to approximate the exponential integral because their nonanalytic nature, each one of them is suitable in some range for the values of u. It will be beneficial to have a simple and handy algorithm for occasional use Barry et al., 2000, see [24]. In the following, most of these methods are reviewed briefly.

Series expansion
The power series in equation (2) is converging with any < 1, that is But for large values of u, the rate of convergence becomes slow. For example, in 1942, Coulson and Duncanson see [25][26] showed to compute up to 75 terms to get accuracy with two significant figures, we need to take u = 20. Therefore, Harris in 1957, see [27] showed that equation (2) is valid if "u < 4", or preferably if "u < 1" (Cody and Thacher, 1968; Stegun and Zucker, 1974, see [28,29]). But, if "u < 0.01" Cooper and Jacob, in 1946, see [30] used the approximation ( ) ≅ − − ln( ) for the practical pump-test problems.

Rational Chebyshev approximation
This method is suggested by Cody and Thacher in 1968, see [31] and the idea of the method is based on Chebyshev approximation polynomial to find a new formula for evaluation 1 ( )as follows: and the suitable result when n = 6 is got in [32].

Asymptotic expansion
Integration by parts can be used to equation (1) to get the following asymptotic expansion of the exponential integral (1): Coulson and Duncanson in 1942 [33] used equation (5) to approximate E1(u) when u >15 , while Harris in 1957, see [34] used equation (5) when u >50. Bleistein and Handelsman in 1975 see [35] got the best approximation of equation (5) if n (greatest integer) ≤ u.

Continued Fraction
Gautschi and Cahill in 1964 see [36] suggest another useful representation of E1(u) which was the continued fraction as: Van der Laan and Temme in 1984 see [37] discussed the properties of this representation in detail. The generalized of En(u) when u > 1 has been approximated as the even form of the continued fraction by Stegun and Zucker in 1974. Stegun and Zucker in 1974 see [38] showed that the continued fraction converges more rapidly than the asymptotic expansion for large u.

Polynomial, Rational and Chebyshev polynomial approximations
Polynomials are one of the most common methods of approximating functions in a given interval. Nevertheless, the optimal rational function, which is the ratio of two polynomials, is able to achieve higher accuracy than the optimal polynomial with the same number of coefficients. The best rational approximation is difficult to be obtained using easy methods; therefore, many techniques have been proposed based on minimizing some measure of errors (Hamming, 1973, see [39]). The widely used rational approximation of functions is the socalled minimax or Chebyshev approximation, which minimizes the maximum error. Many authors used this technique to approximate the exponential integral; Gautschi and Cahill in 1964 see [40] approximated E1(u) for u[0, 1] by a polynomial of degree five. Rational functions with the same degree in numerator and denominator have been developed to approximate the expression "ue u El(u)" in the Chebyshev sense by Hastings in 1955, see [41]. Cody and Thacher in 1968 see [42] introduced rational Chebyshev approximations in three intervals, later in 1969 see [43]; they used the method to approximate the exponential integral in four intervals.

Other Approximate Methods
There are some other methods used to approximate the exponential integral such stehfest algorithm of the inverse Laplace transform, finite difference and numerical quadrature formulas such as trapezoidal, Simpson, and Gauss-Laguerre. Moreover, Prodanoff et al., in 2006 see [44] and Nadarajah in 2007 see [21] used smoothing techniques to approximate the well function which is employed whenever the integrand gradient is high around a point, or even when it is singular. Some approximations of Theis' solution are reviewed and compared with the numerical evaluation of the exponential integral by means of a short FORTRAN code suggested by Giao, in 2003, see [45]. While Swamee with Ojha, in 1990, see [41] and Vatankhah in 2014 see [24] proposed full-range approximations valid for all values but got the result with less accuracy.
The results obtained by previous methods were even less accurate than required, so we suggest the following approach.

Suggested Approach
Here, we suggest the following algorithm to approximate equation (3) based on optimal rational Chebyshev approximation and to compute the best approximation.
iv.Repeat steps (ii) and (iii) at the k th time using the set of points { (1) }.

Discussion
From the previous studies, we note that the results obtained by using rational Chebyshev approximation which is suggested by Cody and Thacher, is the best compared with those obtained by using other methods which are studied in section three. However, the results of our proposed method are more accurate and have fewer calculations than others as illustrated in Figure 1 and Table 3.

Conclusion
In this article, we present a comprehensive review and study of the numerical calculations for exponential integral function (also called the Theis well function) that is frequently appearing in well-hydraulic literature. We suggest an efficient approach to calculate this function. We see that all previous studies until now do not give the best approximations to exponential integral that are valid for the full range of u. The suggested method in this work achieves this aim. For small u, the series representation is still the most widely used technique among all these methods which is used to calculate E1(u). For large u, an effort and special techniques are required. The practical results show the suggested approach can be easily evaluated this value in software algorithms with acceptable accuracy.