Accurate approximations for the complex error function with small imaginary argument

# Accurate approximations for the complex error function with small imaginary argument

S. M. Abrarov111Dept. Earth and Space Science and Engineering, York University, Toronto, Canada, M3J 1P3.  and B. M. Quine222Dept. Physics and Astronomy, York University, Toronto, Canada, M3J 1P3.
November 14, 2014
###### Abstract

In this paper we present two efficient approximations for the complex error function with small imaginary argument over the range that is commonly considered difficult for highly accurate and rapid computation. These approximations are expressed in terms of the Dawson’s integral of real argument that enables their efficient implementation in a rapid algorithm. The error analysis we performed using the random input numbers and reveals that in the real and imaginary parts the average accuracy of the first approximation exceeds and , while the average accuracy of the second approximation exceeds and , respectively. The first approximation is slightly faster in computation. However, the second approximation provides excellent high-accuracy coverage over the required domain.
Keywords: Complex error function, Voigt function, Faddeeva function, Kramp function, Dawson’s integral

## 1 Introduction

The complex error function, also known as the Faddeeva function or the Kramp function, can be defined as [1, 2, 3, 4, 5, 6]

 w(z) =e−z2[1−erf(−iz)] (1) =e−z2⎛⎜⎝1+2i√πz∫0et2dt⎞⎟⎠,

where is the complex argument. Using the Fourier transforms, the real and imaginary parts of the complex error function (1) can be represented as [3, 7, 8]

 K(x,y)=1√π∞∫0exp(−t2/4)exp(−yt)cos(xt)dt,y>0 (2)

and

 L(x,y)=1√π∞∫0exp(−t2/4)exp(−yt)sin(xt)dt,y>0, (3)

respectively. The real part of the complex error function is known as the Voigt function , widely used in many fields of Physics, Astronomy and Chemistry [10, 11, 12, 13]. Combining the real (2) and imaginary (3) parts together as yields

 w(x,y)=1√π∞∫0exp(−t2/4)exp(−yt)exp(ixt)dt,y>0. (4)

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−z2−w(z),

it follows that only positive arguments and are sufficient in order to cover the entire complex plane. Therefore, we will imply further that both input parameters and are always positive.

When the argument is large enough by absolute value, say , the truncation of the Laplace continued fraction [16, 17] (see also [4, 5])

 w(z)=μ0z−1/2z−1z−3/2z−2z−5/2z−⋯,μ0≡i/√π,

can be effectively used for high-accuracy and rapid computation of the complex error function .

There are several approximations such as the Chiarella and Reichel approximation (equation (15) in ), the Weideman’s rational approximation (equation (38-I) in ), the exponential series approximation (equation (14) in  and its modification (3) in ) and the rational approximation (equation (14) in ) that provide highly accurate and rapid calculation for the complex error function within domain and .

However, the highly accurate and simultaneously rapid computation of at and 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 . In the rapid algorithm developed by Karbach et al.  this domain is covered by using the exponential series approximation (equation (14) in ), the Taylor expansion series near singularities at ( is the index integer and 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  (see also Matlab source code in , 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 is to use the Maclaurin expansion series near the point (see for example two approximations (Appendix A) and (Appendix A) in Appendix A). However, this method results in a moderate accuracy that does not exceed in the range of a concern at .

In order to obtain more efficient approximations, we may attempt to find an appropriate representation of the complex error function . Let us rewrite the equation of the complex error function (4) as

 w(x,y)=ey2√π∞∫0e−(t+2y)2/4eixtdt.

Change of the variable in equation above excludes the parameter from the integrand:

 w(x,y) =ey2√π∞∫2ye−u2/4eix(u−2y)du =ey2−2ixy√π∞∫2ye−u2/4eixudu.

Consequently, this equation can be rearranged as

 w(x,y)=ey2−2ixy√π⎛⎜⎝∞∫0e−u2/4eixudu−2y∫0e−u2/4eixudu⎞⎟⎠. (5)

The first integral in the equation above is

 ∞∫0e−u2/4eixudu=e−x2√π+2iF(x),

where

 F(x)=12∞∫0e−u2/4sin(xu)du

is the Dawson’s integral. As the argument in 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 is very narrow. Consequently, the second integral in equation (5) makes only a minor contribution to the complex error function . 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:

 e−u2/4=1−u24+u432−u6384+O(u8).

The first approximation is obtained by substituting only the first term of this expansion into the second integral of equation (5) that yields

 w(x,y<<1)≈e(ix−y)2[1+iex2√π(2F(x)−1−e2ixyx)]. (6)

We will refer to this equation as the basic approximation. This equation is highly accurate while .

The second approximation is obtained by substituting the first two terms of the expansion above into second integral of the equation (5):

 w(x,y<<1) ≈e(ix−y)2[1+iex2√π (7) ×(2F(x)−x2+0.5+e2ixy(x2(y2−1)+ixy−0.5)x3)].

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 .

## 3 Algorithmic implementation

Approximations (6) and (7) cannot be employed directly in the computation flow since at the denominators in approximation (6) and 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:

 w(x<<1,y<<1)≈(1−(x+iy)2)(1+2i(x+iy)√π). (8)

In particular, the computation flow for the basic and the main approximations can be maintained as

 w(x,y<<1)≈{Eq.(6)orEq.(7),10−4

Thus, according to this scheme if the complex error function is computed either by equation (6) or by equation (7). Otherwise, if it is computed by equation (8).

The approximation (8) had been used already in our recently published Matlab source code  to resolve a similar problem that occurs at . The derivation of the approximation (8) is shown in Appendix A.

## 4 Error analysis

Define the relative errors for the real and imaginary parts of the complex error function in forms

 ΔRe=∣∣ ∣∣Re[w(z)]−Re[wref.(z)]Re[wref.(z)]∣∣ ∣∣

and

 ΔIm=∣∣ ∣∣Im[w(z)]−Im[wref.(z)]Im[wref.(z)]∣∣ ∣∣,

respectively, where is the reference. The highly accurate reference values can be obtained according to equation (1) by using the latest versions of Wolfram Mathematica that supports error function of complex argument.

Figure 1a shows the logarithm of relative error for the real part of the basic approximation (6) in the domain and . As we can see, at the approximation (6) is highly accurate and provides accuracy better than . Although the accuracy of the approximation (6) deteriorates as increases, it, nevertheless, still remains high and better than .

Figure 1b illustrates the logarithm of relative error for the imaginary part of the basic approximation (6) in the domain and . As one can see, the imaginary part of the approximation (6) is highly accurate. Specifically, the accuracy is better than almost over all of the domain. However, at small the accuracy deteriorates as it can be seen by red color along the left edge of -axis.

Figure 2a illustrates the logarithm of relative error for the real part of the supplementary approximation (8) over the small domain and . One can see that the real part over this domain is highly accurate and better than .

Figure 2b depicts the logarithm of relative error for the imaginary part of the supplementary approximation (8) over the small domain and . The accuracy of the imaginary part is also high and better than .

From Figs. 2a and 2b we can conclude that the basic approximation (6) alone may be used practically as the accuracy better than is more than enough for the most applications.

The accuracy of the complex error function 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 and . As it can be seen from this figure, in the real part the accuracy of the approximation (7) is better than .

Figure 3b shows the logarithm of relative error for the imaginary part of the main approximation (7) in the domain and . 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 almost over all domain.

To determine the average accuracy for both approximations, we performed calculations with randomly chosen input variables and . It has been found that in the real and imaginary parts the basic approximation (6) provides average accuracy exceeding and , while the main approximation (7) provides average accuracy exceeding and , respectively.

## 5 Conclusion

We present two efficient approximations for the complex error function that can be used for computation at small . As these approximation are expressed in terms of the Dawson’s integral of real argument , they are rapid in computation. Error analysis performed with randomly taken input variables and reveals that in the real and imaginary parts the basic approximation (6) provides average accuracy exceeding and , while the main approximation (7) provides average accuracy exceeding and , 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 :

 e−(x+iy)2=e−x2−2ixye−x2+O(y2),
 erf(−i(x+iy))=−ierfi(x)+2y√πex2+O(y2)

and

 w(x,y) ≡e−(x+iy)2[1−erf(−i(x+iy))] =e−x2(1+ierfi(x))+2y(xe−x2erfi(x)−ixe−x2−1√π)+O(y2),

Substituting (Appendix A) and (Appendix A) into identity (1) we obtain

 w(x,y<<1)≈(1−2ixy)(e−x2(1+erf(ix))−2y√π),

while equation (Appendix A) immediately results in

The approximations (Appendix A) and (Appendix A) are rapid in computation because the error function of imaginary argument can be expressed in terms of the Dawson’s integral of real argument as follows

 erf(ix)≡2i√πex2F(x).

However, the approximations (Appendix A) and (Appendix A) are moderately accurate and should be used when the requirement for accuracy does not exceed .

Same technique can be used to obtain approximation for from the expansions near the point :

 exp(−z2)=1−z2+O(z4),

and

 erf(−iz)=−2iz√π+O(z3),

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

 2y∫0e−u2/4eixudu=√πe−x2[erf(ix)−erf(ix−y)].

The argument in the first error function 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 is complex, it needs a numerical solution (see for example equations (10) and (11) in ). When the high accuracy is not required, the rapid approximation for the error function can obtained, for example, by rearranging the expansion series (Appendix A) as

 erf(ix−y) ≈ierfi(x)−2y√πex2 ≡erf(ix)−2y√πex2,y<<1.

We can also approximate this integral by taking into account that . In particular, for small interval we can write and take the mid-point in that interval equal to :

 2y∫0e−u2/4eixudu≈2ye−y2/4eixy.

Alternatively, according to the trapezoidal rule we get

 2y∫0e−u2/4eixudu≈y(1+e−y2e2ixy).

Once again, these approximations should be applied only when the high-accuracy is not a concern in computation.

The second integral from equation (5) can be approximated more precisely by taking integration by parts

 2y∫0e−u2/4eixudu=√πerf(u2)eixu∣∣∣2y0−√πix2y∫0erf(u2)eixtdu

and then substituting the first few terms from the following expansion

 erf(u2)=u√π−u312√π+u5160√π−u72688√π+O(u9)

into the right integral of equation (Appendix B). 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+iex2√π ×(2F(x)−1−e2ixy(1+ix(√πerf(y)−2y))x)].

The consistency between approximations (6) and (Appendix B) becomes evident from the fact that . Although approximation (Appendix B) 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 . Since the approximation (Appendix B) provides nearly same accuracy as the main approximation (7), both of them can be used when the high-accuracy computation is required.

## Acknowledgments

This work is supported by National Research Council of Canada, Thoth Technology Inc. and York University. The authors wish to thank Prof. Ian McDade and Dr. Brian Solheim for discussions and constructive suggestions.

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters   