Numerical precision - dawson(x) vs erfcw(real(x+iy))

_{In honor of John's retirement, I shall resurrect this 10 yr old thread. (ha!) I have an option to choose between a Dawson integral form or an error function form of an expression in numerical calculations. My colleague prefers  the former from a first-principled derivation perspective, and I prefer the latter both for ease in making a larger expression dimensionless by removing common factors and for commonality in numerical analysis (i.e. error functions are more "colloquial" than Dawson's integrals).}

_{Just for my own information, I referenced the link below. I am interested in the D+ case.}

_{https://en.wikipedia.org/wiki/Dawson_function}

_{As near as I can now determine, for exact translation between the Dawson integral and error function forms, I must input a complex number to the erf(...) function in Igor Pro and pull off the imaginary component. Hence, the two functions below give the same report for any input.}

// @brief calculate Dawson integral form
Function cal_Dawson(variable kappa)

	variable vrtn
	vrtn = Dawson(kappa)
	return vrtn
end

// @brief calculate erfi form
Function cal_erfi(variable kappa)

	variable vrtn
	vrtn = (sqrt(Pi)/2)*exp(-kappa^2)*imag(erf(cmplx(0,kappa)))
	return vrtn
end

My question is whether the second function can be written directly to "keep it all real"?? By example, just guessing, I tried the function below.

// @brief calculate erfr form
Function cal_erfr(variable kappa)

	variable vrtn
	vrtn = (sqrt(Pi)/2)*exp(kappa^2)*real(erfcw(kappa))
	return vrtn
end

For additional background, I am interested in capturing the equivalence only in the range of positive, real kappa. And I am satisfied with a relative accuracy at the ppm level.

Hello JJ,

I wrote the original Dawson function using code that I wrote in the 1980's.  It used quadrature to compute the integral over the range that was good enough for my PhD work.  Some 20 years later JW replaced the code with a function from the Faddeeva package that boasts higher accuracy.  You can read more about the package here:  http://ab-initio.mit.edu/Faddeeva.

Igor's various error functions are calculated using different algorithms with varying levels of accuracy:

Function   Argument   Algorithm
erf              real                 Incomplete gamma: P(1/2, x^2)
erfc             real                  Incomplete gamma: Q(1/2, x^2), or 1+P(1/2, x^2) for x<0
erf              complex          (2z/sqrt(pi)) * 1F1(1/2; 3/2; -z^2)
erfc            complex          1 - erf(z) via the hypergeometric series
erfcw          complex          Poppe & Wijers, ACM TOMS 16(1), 38-46 (1990)
erfcx            real                  Polynomial approximation (fma-based)

For accuracy reasons I'd stay with Dawson based in the Faddeeva package especially in view of the fact that some of the hypergeometric functions are sometimes difficult to push to high accuracy.

HTH,

AG

Thanks AG. Just curious though. As near as I see, if I input a complex number to the error function, I am tapping into the Faddeeva routine just as in the Dawson function. So, the two functions that I report are equivalent ways in Igor Pro to calculate the same value and should exhibit nearly the same accuracy and speed??

Just for further background, my colleague arrived at the Dawson integral form using derivations by hand. The input was Z = sqrt(X) whereby X is constrained to be positive even though foundations require X to be negative. It seems perhaps that Dawson(Z(positive real)) = Dawson(-i*|Z|(negative complex)) (symmetry). My derivations with symbolic math (Maple in this case) resulted in the error function form, whereby Maple's error function input allowed X to be negative, resulting in a complex Z, resulting in a final agreement. I attach a report from within Maple showing three equivalent ways to obtain the same value in the function we are investigating.

Finally, had you asked me about three weeks ago whether I knew anything about Dawson's integral, I would have said ... Isn't there a TV show with that name? But, everyone in college calculus has been introduced to the error function. Hence my desire to provide future readers of our collaborative publication something that is familiar (error function) even as I now better respect the equivalence (and potential benefits) of using the Dawson's integral.

--
JJW

Hello JJ,

I think you misunderstood the table in my last comment.  If you execute, e.g., 

•print/c erf(cmplx(1,1))

Igor computes the erf() via the HyperGeometric function.

AG

Apologies. I mistook 1F1 as the Faddeeva package approach.

All is now clear.

Thanks!