Complex Error Function in Mathematica

Question

The complex error function w (z) is defined as e^(-x^2) erfc(-ix) . The problem with using w (z), as defined above, is that erfc tends to break for large x (complemented by an exponent going to 0, so everything remains small), so Mathematica goes back to arbitrary precision calculations that make life VERY Slow. This function is used to implement the voigt profile — the line shape commonly used in spectroscopy and other related fields. Right now, I am returning to calculating the line at once and using interpolation to speed up the process, but this does not allow me to easily change the parameters of the line (or match them).

scipy has a nice and fast w (z) implementation like scipy.special.wofz , and I was wondering if there is an equivalent in Mathematica.

Answer 1

The complex error function can be written in terms of the Hermite “polynomial” H_{-1}(x) :

 In[1]:= FullSimplify[2 HermiteH[-1,I x] == Sqrt[Pi] Exp[-x^2] Erfc[I x]] Out[1]= True 

And the score does not suffer so many overflows and overflows

 In[68]:= 2 HermiteH[-1, I x] /. x -> 100000. Out[68]= 6.12323*10^-22 - 0.00001 I In[69]:= Sqrt[Pi] E^-x^2 Erfc[I x] /. x -> 100000. During evaluation of In[69]:= General::unfl: Underflow occurred in computation. >> During evaluation of In[69]:= General::ovfl: Overflow occurred in computation. >> Out[69]= Indeterminate 

However, some quick tests show that the rate of estimation of the Hermite function should be slower than the speed of the exponent and the error function ...

Answer 2

The expansion in a series at infinity shows that the real and imaginary parts have very different scales. I would suggest to calculate them separately and not add. Below, I use the first few members of the series extension to get the imaginary part.

 In[186]:= w[x_?NumericQ] := {N[Exp[-SetPrecision[x, 25]^2], 20], N[(3 /(4 Sqrt[\[Pi]] x^5) + 1/(2 Sqrt[\[Pi]] x^3) + 1/( Sqrt[\[Pi]] x))]} In[187]:= w[11] Out[187]= {2.8207700884601354011*10^-53, 0.05150453151309212} In[188]:= w[1000] Out[188]= {3.296831478088558579*10^-434295, 0.0005641898656429712} 

Not sure how much you want this very small real part. If you can drop it, that will keep numbers within a reasonable range. In some ranges (or if higher accuracy of the machine is required), you can use more terms from the extension to this imaginary part.

Daniel Lichtblow Wolfram Research

Answer 3

The real and imaginary parts of the complex error function on a real line can be explicitly and efficiently calculated in Mathematica using the Dawson integral :

 In[9]:= Sqrt[Pi] Exp[-x^2] Erfc[I x] == E^-x^2 Sqrt[\[Pi]] - 2 I DawsonF[x] // FullSimplify Out[9]= True 

This is about 4 times faster than using HermiteH[-1,z] .

 In[10]:= w1[x_] := E^-x^2 Sqrt[\[Pi]] - 2 I DawsonF[x] w2[x_] := 2 HermiteH[-1, I x] In[15]:= AbsoluteTiming[w1 /@ Range[-5.0, 5.0, 0.001];] Out[15]= {2.3272327, Null} In[16]:= AbsoluteTiming[w2 /@ Range[-5.0, 5.0, 0.001];] Out[16]= {10.2400239, Null} 

Answer 4

A C program for the complex error function (the so-called Faddeeva function), which can be run from Mathematica, is also available in RooFit . Read an article by Karbach et al. arXiv: 1407.0748 for more information.

Answer 5

Just wrap the libcerf C library .