Accurate approximations for the complex error function with small imaginary argument

In this paper we present two efficient approximations for the complex error function $w \left( {z} \right)$ with small imaginary argument $\operatorname{Im}{\left[ { z } \right]}<<1$ over the range $0 \le \operatorname{Re}{\left[ { z } \right]} \le 15$ that is commonly considered difficult for highly accurate and rapid computation. These approximations are expressed in terms of the Dawson's integral $F\left( x \right)$ of real argument $x$ that enables their efficient implementation in a rapid algorithm. The error analysis we performed using the random input numbers $x$ and $y$ reveals that in the real and imaginary parts the average accuracy of the first approximation exceeds ${10^{ - 9}}$ and ${10^{ - 14}}$, while the average accuracy of the second approximation exceeds ${10^{ - 13}}$ and ${10^{ - 14}}$, respectively. The first approximation is slightly faster in computation. However, the second approximation provides excellent high-accuracy coverage over the required domain.

None of the integrals above can be taken analytically in closed form. Consequently, the complex error function must be computed numerically.
It should be noted that from identity [14,15] w (−z) = 2e −z 2 − w (z) , it follows that only positive arguments x and y are sufficient in order to cover the entire complex plane. Therefore, we will imply further that both input parameters x and y are always positive.
When the argument z is large enough by absolute value, say |x + iy| > 15, the truncation of the Laplace continued fraction [16,17] (see also [4,5]) can be effectively used for high-accuracy and rapid computation of the complex error function w (z).
However, the highly accurate and simultaneously rapid computation of w (z) at y << 1 and |x + iy| ≤ 15 still remains problematic [23,24]. Different approaches have been implemented to overcome this problem. For example, Zaghloul and Ali developed Algorithm 916 that can cover accurately this domain [15]. In the rapid algorithm developed by Karbach et al. [22] this domain is covered by using the exponential series approximation (equation (14) in [20]), the Taylor expansion series near singularities at z n = nπ/τ m (n is the index integer and τ m is the integration cutoff) and the Laplace continued fraction. In our recent work we proposed to cover this domain by using the master-slave algorithm [18] (see also Matlab source code in [25], where the master-slave approach has been implemented to cover it). In this paper we report two new approximations that also effectively resolve this problem for accurate and rapid computation.
2 Approximations for small y A simplest way to approximate the complex error function w (x, y << 1) is to use the Maclaurin expansion series near the point y = 0 (see for example two approximations (A.4) and (A.5) in Appendix A). However, this method results in a moderate accuracy that does not exceed 10 −6 in the range of a concern |x + iy| ≤ 15 at y << 1.
In order to obtain more efficient approximations, we may attempt to find an appropriate representation of the complex error function w (z). Let us rewrite the equation of the complex error function (4) as Change of the variable u = t + 2y in equation above excludes the parameter y from the integrand: Consequently, this equation can be rearranged as The first integral in the equation above is is the Dawson's integral. As the argument x in F (x) is real, its computation is not difficult and several efficient approximations that can provide rapid and highly accurate computation are reported in literature [26,27,28]. Furthermore, the latest versions of Matlab support the built-in function dawson(x) that computes very rapidly the Dawson's integral of real argument. The advantage of the equation (5) is that the integration range in the second integral at y << 1 is very narrow. Consequently, the second integral in equation (5) makes only a minor contribution to the complex error function w (z). As a result, this minimizes the error in computation.
The representation of the complex error function in form of (5) provides enormous flexibility to approximate it. Here we show two efficient derivations that directly follow from equation (5).
As the upper limit of the second integral of equation (5) is finite, its evaluation is not straightforward (see Appendix B). However, this problem can be overcome by Maclaurin expansion of the exponential function: The first approximation is obtained by substituting only the first term of this expansion into the second integral of equation (5) that yields We will refer to this equation as the basic approximation. This equation is highly accurate while x 3. The second approximation is obtained by substituting the first two terms of the expansion above into second integral of the equation (5): Further, we will refer to this equation as the main approximation. Although the basic approximation (6) is slightly faster in computation, the main approximation (7) provides essentially higher accuracy as it will be shown in the section 4.

Algorithmic implementation
Approximations (6) and (7) cannot be employed directly in the computation flow since at x → 0 the denominators x in approximation (6) and x 3 in approximation (7) blow up the ratios resulting to computational overflow.
However, this problem is very easy to resolve by using, for example, the following supplement: In particular, the computation flow for the basic and the main approximations can be maintained as Thus, according to this scheme if x > 10 −4 the complex error function is computed either by equation (6) or by equation (7). Otherwise, if x ≤ 10 −4 it is computed by equation (8).
The approximation (8) had been used already in our recently published Matlab source code [25] to resolve a similar problem that occurs at |x + iy| → 0. The derivation of the approximation (8) Figure 1a shows the logarithm of relative error for the real part of the basic approximation (6) in the domain 10 −4 ≤ x ≤ 15 and y ≤ 10 −6 . As we can see, at x 3 the approximation (6) is highly accurate and provides accuracy better than 10 −12 . Although the accuracy of the approximation (6) deteriorates as x increases, it, nevertheless, still remains high and better than 10 −9 .    (6) in the domain 10 −4 ≤ x ≤ 15 and y ≤ 10 −6 . As one can see, the imaginary part of the approximation (6) is highly accurate. Specifically, the accuracy is better than 10 −14 almost over all of the domain. However, at small x 10 −3 the accuracy deteriorates as it can be seen by red color along the left edge of y-axis.   One can see that the real part over this domain is highly accurate and better than 10 −13 . Figure 2b depicts the logarithm of relative error log 10 ∆ Im for the imaginary part of the supplementary approximation (8) over the small domain 0 ≤ x ≤ 10 −4 and y ≤ 10 −6 . The accuracy of the imaginary part is also high and better than 10 −12 .
From Figs. 2a and 2b we can conclude that the basic approximation (6) alone may be used practically as the accuracy better than 10 −9 is more than enough for the most applications.  The accuracy of the complex error function w (x, y << 1) can be further improved by using the main approximation (7). Figure 3a depicts the logarithm of relative error for the real part of the approximation (7) in the domain 10 −4 ≤ x ≤ 15 and y ≤ 10 −6 . As it can be seen from this figure, in the real part the accuracy of the approximation (7) is better than 10 −13 . Figure 3b shows the logarithm of relative error for the imaginary part of the main approximation (7) in the domain 10 −4 ≤ x ≤ 15 and y ≤ 10 −6 . Comparing Figs. 3b with 1b one can observe that the graphs are nearly same. This signifies that the behavior of the basic and the main approximations resemble in their imaginary parts. Similar to Fig. 1b, the accuracy in the imaginary part is better than 10 −14 almost over all domain.
To determine the average accuracy for both approximations, we performed calculations with randomly chosen input variables x and y. It has been found that in the real and imaginary parts the basic approximation (6) provides average accuracy exceeding 10 −9 and 10 −14 , while the main approximation (7) provides average accuracy exceeding 10 −13 and 10 −14 , respectively.

Conclusion
We present two efficient approximations for the complex error function that can be used for computation at small y << 1. As these approximation are expressed in terms of the Dawson's integral F (x) of real argument x, they are rapid in computation. Error analysis performed with randomly taken input variables x and y reveals that in the real and imaginary parts the basic approximation (6) provides average accuracy exceeding 10 −9 and 10 −14 , while the main approximation (7) provides average accuracy exceeding 10 −13 and 10 −14 , respectively. Although the basic equation (6) is slightly faster in computation, the main approximation (7) provides excellent high-accuracy coverage.

Appendix A
Consider the following Maclaurin expansions at y = 0:

Substituting (A.1) and (A.2) into identity (1) we obtain
The approximations (A.4) and (A.5) are rapid in computation because the error function of imaginary argument erf (ix) can be expressed in terms of the Dawson's integral of real argument F (x) as follows However, the approximations (A.4) and (A.5) are moderately accurate and should be used when the requirement for accuracy does not exceed 10 −6 . Same technique can be used to obtain approximation for w (x << 1, y << 1) from the expansions near the point z = 0: In particular, substituting these expansions into identity (1) leads to the approximation (8).

Appendix B
The second integral of equation (5) can be expressed as The argument in the first error function erf (ix) is purely imaginary and, therefore, it can be expressed in term of the Dawson's integral of real argument. However, since the argument in the second error function erf (ix − y) is complex, it needs a numerical solution (see for example equations (10) and (11) in [15]). When the high accuracy is not required, the rapid approximation for the error function can obtained, for example, by rearranging the expansion series (A.2) as We can also approximate this integral by taking into account that y << 1. In particular, for small interval we can write du ≈ 2y and take the mid-point in that interval equal to y: Alternatively, according to the trapezoidal rule we get 2y 0 e −u 2 /4 e ixu du ≈ y 1 + e −y 2 e 2ixy .
Once again, these approximations should be applied only when the highaccuracy is not a concern in computation.
The second integral from equation (5) can be approximated more precisely by taking integration by parts and then substituting the first few terms from the following expansion into the right integral of equation (B.1). We will consider the simplest case when only the first term from this expansion is substituted. Following this procedure and using equation (5) we have w (x, y << 1) ≈ e (ix−y) 2 1 + ie x 2 √ π × 2F (x) − 1 − e 2ixy (1 + ix ( √ π erf (y) − 2y)) x .

(B.2)
The consistency between approximations (6) and (B.2) becomes evident from the fact that erf (y << 1) ≈ 2y/ √ π. Although approximation (B.2) is more accurate than the basic approximation (6), it is slightly slower in computation at large size of the input array due to presence of the error function erf (y). Since the approximation (B.2) provides nearly same accuracy as the main approximation (7), both of them can be used when the high-accuracy computation is required.