Complex Error Function in Mathematica

6
votes

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 the erfc tends to explode out for larger x (complemented by the exponential going to 0 so everything stays small), so that Mathematica reverts to arbitrary precision calculations that make life VERY slow. The function is used in implementing the voigt profile - a line shape commonly used in spectroscopy and other related areas. Right now I'm reverting to calculating the lineshape once and using an interpolation to speed things up, however this doesn't let me alter the parameters of the lineshape (or fit to them) easily.

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

5
wolfram-mathematica
Shouldn't it be the other way around in that the error function blows up and the exponential decays for x large (and real)? - Heike
I tried playing with SetSystemOptions["CatchMachineUnderflow"->False], which, however, results in getting 0. +0. I for large arguments. Then I tried defining the function as Exp[Log[-z^2+Log@Erfc[-I*z]]], but this turns out to not be any faster than with automatic switching. So, it seems hard to speed this up, except as @Daniel Lichtblau suggests - acl
@Heike you are right, a bit of a slip on my part. - crasic

5 Answers

5
votes

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 evaluation does not suffer as many underflows 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

That said, some quick tests show that the evaluation speed of the Hermite function to be slower than that of the product of the exponential and error function...

4
votes

A series expansion at infinity shows that the real and imaginary parts are of very different scales. I'd suggest computing them separately and not adding them. Below I use the first few terms of the series expansion 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 badly you want that very small real part. If you can drop it that will keep the numbers in a reasonable range. In some ranges (or if higher than machine precision is desired) you may want to use more terms from the expansion on that imaginary part.

Daniel Lichtblau Wolfram Research

4
votes

The real and imaginary parts of the complex error function on the real line can be explicitly and efficiently computed in Mathematica using 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}
1
votes

A program in language C for the complex error function (aka the Faddeeva function) that can be run from Mathematica is also available in RooFit. Read the article by Karbach et al. arXiv:1407.0748 for more information.

0
votes

Just wrap the C library libcerf.