A new application methodology of the Fourier transform for rational approximation of the complex error function

This paper presents a new approach in application of the Fourier transform to the complex error function resulting in an efficient rational approximation. Specifically, the computational test shows that with only $17$ summation terms the obtained rational approximation of the complex error function provides the average accuracy ${10^{ - 15}}$ over the most domain of practical importance $0 \le x \le 40,000$ and ${10^{ - 4}} \le y \le {10^2}$ required for the HITRAN-based spectroscopic applications. Since the rational approximation does not contain trigonometric or exponential functions dependent upon the input parameters $x$ and $y$, it is rapid in computation. Such an example demonstrates that the considered methodology of the Fourier transform may be advantageous in practical applications.


Introduction
The forward and inverse Fourier transforms can be defined as [1,2] and respectfully. Approximation theory based on the Fourier trigonometric series for functions or signals remains a topical subject in mathematical analysis and many new efficient methodologies have been reported in the recent scientific literature (see for example [3,4,5]). In our recent publication [6] we have shown that a sampling with the Gaussian function of the kind he −(t/c) 2 / (c √ π) leads to the trigonometric approximations for the forward and inverse Fourier transforms where h is the step between two adjacent sampling points and c is the fitting parameter, e −2πiνnh = cos (2πiνnh) − i sin (2πiνnh) and e 2πitnh = cos (2πitnh)+i sin (2πitnh). The parameters h, c and N in the equations (2a) and (2b) may be the same in the forward and inverse Fourier transforms only when we imply the most favorable conditions h << 1, c << 1 and N >> 1.
In practical tasks, however, these conditions may be compromised in order to reduce the number of the summation terms. As a result, these parameters may not be necessarily equal to each other in the forward and inverse Fourier transforms. Consequently, it is convenient to rewrite two equations above in form where h f , c f , M and h i , c i , N are the steps, the fitting parameters and the integers corresponding to the forward and inverse Fourier transforms, respectively.
The presence of the damping functions e −(πc f ν) 2 and e −(πc i t) 2 in the equations above excludes periodicity of the approximated functions f (t) and F (ν). Consequently, a solitary wavelet (or non-periodic pulse) can be effectively approximated in the Fourier transform. However, when we take c f = c i = 0, the right side of these equations become periodic with corresponding periods 1/h f , 1/h i and represent the Fourier-type expansion series as follows and It should be noted that if the integral (1a) is not analytically integrable, then the function f (t) can be approximated by substituting equation (3a) into (3b). This substitution yields In this work we show a new application methodology of the Fourier transform to the complex error function. Due to representation of the complex error function as a rational approximation, it is rapid in computation. Furthermore, with only 17 summation terms the obtained rational approximation of the complex error function provides accuracy 10 −15 over the most domain of practical importance 0 ≤ x ≤ 40, 000 ∩ 10 −4 ≤ y ≤ 10 2 required for applications utilizing the HITRAN molecular spectroscopic database [7].

Function overview
The complex error function, also known as the Faddeeva function or the Kramp function, can be defined as [8,9,10,11,12] where z = x + iy is the complex argument. The complex error function is a solution of the differential equation [12] w with initial condition w (0) = 1. The complex error function is closely related to a family of the special functions. Among them the most important one is the complex probability function [11,12,13] This principal value integral implies that the complex probability function has no discontinuity at y = 0 and x = t. In particular, where daw (x) is the Dawson's integral that will be briefly introduced later. There is a direct relationship between complex error function and complex probability function [11,12] The real part of the complex probability function, denoted as K (x, y), is known as the Voigt function. Mathematically, the Voigt function represents a convolution integral of the Gaussian and Lorentzian distributions [11,12,14,15,16] K where the principal value integral also implies that it has no discontinuity at y = 0 and x = t. Specifically, from equation (4) it follows that At non-negative argument y the real part of the complex error function is also the Voigt function in accordance with identity (5). The Voigt function is widely used in many spectroscopic applications as it describes the line broadening effects [17,18,19,20,21]. Therefore, the application of the complex error function is very significant in quantitative spectroscopy.
Other closely related functions are the error function of complex argument [12] w the plasma dispersion function [22] the Dawson's integral [23,24,25,26,27] and the normal distribution function [29] Φ It is not difficult to show that the complex error function can be represented in an alternative form (see equation (3) in [30] and [31], see also Appendix A in [27] for derivation) This representation of the complex error function will be used for derivation of a rational approximation.

Rational approximation
In our recent publications we have shown a new technique to obtain a rational approximation for the integrals of kind [32,33] ∞ 0 e −t 2 f (t) dt.
We apply this approach together with the Fourier transform methodology discussed above in the Introduction. We can use either of equation (3a) or (3b). For example, we may choose the equation (3b) corresponding to the inverse Fourier transform. Consider the function f (t) = e −t 2 /4 . Let us find first its forward Fourier transform by substituting f (t) = e −t 2 /4 into equation (1a). These leads to Now substituting 2 √ πe −(2πν) 2 into equation (3b) yields the following approximation for the exponential function Taking into account that the approximation (7) can be simplified as given by The right side limitation t ≤ 1/ (2h i ) along the positive t-axis in equation (8) can be readily excluded by multiplying both its sides to exp (−σt) if a constant σ is positive and sufficiently large. This can be explained by considering Fig. 1 that shows two functions computed according to right side of equation (8) at σ = 0.1 (blue curve) and σ = 0.2 (red curve). For example, at σ = 0.1 we can observe two additional peaks at 1/h i and 2/h i (blue curve). However, as σ increases the additional peaks are suppressed stronger to zero due to multiplication to the damping exponential function exp (−σt). As a result, at σ = 0.2 only a single additional peak at 1/h i remains visible (red curve). By σ 1 all additional peaks completely vanish and, therefore, do not contribute to error in integration. Consequently, if the constant σ is large enough, say approximately equal or greater than 1, we can write the approximation that remains always valid without any limitation along the positive t-axis.
Assuming y ≥ 0 we, therefore, can write now Since e −t 2 /4 e −yt = e σ 2 e −(t−2σ) 2 /4 e −(y+σ)t from approximation (9) we obtain Once again, due to presence of the rapidly damping exponential multiplier e −(y+σ)t this approximation is valid without any limitation along the positive t-axis. As the peak of the function e −(t−2σ) 2 /4 is shifted towards right with respect to the origin, we may regard to the value σ as the shift constant. Finally, substituting approximation (10) into integral (6) yields where and C n = 2πh i n.
As the expansion coefficients A n , B n and C n are independent of the argument z, the obtained equation (11) is a rational approximation.
In algorithmic implementation it is more convenient to use ψ-function defined as

Computational procedure and error analysis
Due to a remarkable identity of the complex error function [28,34] it is sufficient to consider only I and II quadrants in order to cover the entire complex plane. This can be seen explicitly by representation of the identity (13) in form Thus, if the parameter y is negative we can simply take it by absolute value and then compute the complex error function according to right side of this equation. Therefore, further we will always assume that y ≥ 0. When the argument z is large enough by absolute value, say |x + iy| > ∼ 15, we can truncate the Laplace continued fraction [9,35] Approximation based on the Laplace continued fraction is rapid in computation. However, its accuracy deteriorates as the argument z decreases by absolute value.
There are different approximations for computation of the narrow-band domain 0 ≤ x ≤ 15 and 0 ≤ y < 10 −6 [32,36,37]. We can apply, for example, an approximation proposed in our recent work [32] w (x, y << 1) ≈ 1 − y y min e −x 2 + y y min K (x, y min ) + iL (x, y min ) , y min << 1, where L (x, y min ) = Im [w (x, y min )] and y min can be taken equal to 10 −5 . It has been shown that this approximation can provide accuracy better than 10 −9 over the narrow-band domain 0 ≤ x ≤ 15 and 0 ≤ y < 10 −6 . The domain |x + iy| ≤ 15 ∩ y ≥ 10 −6 is the most difficult for computation. Nevertheless, with only 17 summation terms (at N = 16) the proposed rational approximation (12) covers this domain providing high-accuracy and rapid computation. In computational procedure we have to choose properly the margin value ν m for the exponential function e −(2πν) 2 that appears from the forward Fourier transform F (ν) = 2 √ πe −(2πν) 2 . As it has been justified by Melone et al. [38], the margin value for integration involving the exponential function e −t 2 can be taken as t = t m = 6. We can use this result in order to determine the required value by solving the following equation with respect to the variable ν as follows There are two solutions for this equation ν 1,2 = ±6/ (2π). Consequently, the margin value for the exponential function e −(2πν) 2 can be taken as ν m = 6/ (2π). As a parameter h i is the step between two adjacent sampling points along positive ν-axis (see [6] for details), its value can be calculated as h i = ν m /N . Taking N = 16 we can find that h i = ν m /16 ≈ 5.968310365946075 × 10 −2 .
In order to quantify the accuracy of the rational approximation (12) we may define the relative errors where w ref (x, y) is the reference, for the real and imaginary parts, respectively. The highly accurate reference values can be generated by using, for example, Algorithm 680 [39,40], recently published Algorithm 916 [34] or C++ code from the RooFit package, CERN's library [41]. Figure 2 shows log 10 ∆ Re for the real part of the complex error function computed over the domain 0 ≤ x ≤ 15 and 10 −6 ≤ y ≤ 15 at N = 16, σ = 1.5 and h i = 5.968310365946075 × 10 −2 . As we can see from this figure, the rational approximation (12) provides accuracy 10 −15 (blue color) over the most of this domain. Although accuracy deteriorates with decreasing y, it remains better than 10 −9 (red color) in the range 10 −4 ≤ y ≤ 10 −6 . This indicates that at the same N = 16 the accuracy of the rational approximation (12) is by several orders of the magnitude higher than the accuracy of the Weideman's rational approximation (see equation 38(I) in [13]).  One can see that in the imaginary part the accuracy is also highly accurate 10 −15 (blue color) over the most domain. There is only a small area 0 ≤ x < 1 and 10 −6 ≤ y ≤ 10 −4 near the origin where the accuracy deteriorates as the parameters x and y both tend to zero. Nevertheless, the accuracy in this area still remains high and better than 10 −9 (red color). The computational test reveals that with only 17 summation terms (at N = 16) the rational approximation (12) alone can cover the entire HITRAN domain 0 ≤ x ≤ 40, 000 ∩ 10 −4 ≤ y ≤ 10 2 providing average accuracy 10 −15 for an input array consisting of 3×10 7 elements. Algorithmic implementation of the rational approximation (12) results to the same computational speed as that of described in our recent work where we proposed a sampling by incomplete cosine expansion of the sinc function to approximate the complex error function [33].
A Matlab subroutine code that covers the HITRAN domain with highaccuracy is presented in Appendix A.

Conclusion
We present a new efficient rational approximation to the complex error function by application of the Fourier transform that provides computationally rapid and highly accurate results. The computational test we performed with only 17 summation terms shows that the accuracy of the rational approximation of the complex error function is 10 −15 over the most domain of practical importance. In particular, the proposed rational approximation of the complex error function alone can cover with high accuracy the entire domain 0 ≤ x ≤ 40, 000 ∩ 10 −4 ≤ y ≤ 10 2 required for the HITRAN-based spectroscopic applications. Cn = 2*pi*hi*n; z = z + 1i*sig; % redefine input z (see formula (12) for the psi-function) zz = z.^2; VF = 1i*(2*hi*exp(sig^2))./z; % define first term for n = 1:num VF = VF + (An(n) -1i*z*Bn(n))./(Cn(n)^2 -zz); end end