In this post we will explore different approaches for the evaluation of
the Goodwin-Staton integral, which is one of the generalizations of the
complementary error function. The Goodwin-Staton integral is defined as
$$\begin{equation}G(z) = \int_0^{\infty} \frac{e^{-t^2}}{t+z} \; dt \quad
|\mathrm{arg}\, z| < \pi \end{equation}$$
From now on, we consider the case where $z$ is real positive number. Note
that this special functions can be expressed in terms of other elementary
special functions; the error function and the exponential integral. Thus,
$$\begin{equation}G(x) = -\frac{e^{-x^2}}{2}(\pi i \mathrm{erf}(ix) + E_i(x^2))
\end{equation}$$
Additionally, the previous representation can be slightly simplified by
considering the following connection formulas: $E_i(x) = -E_1(-x)$ and
$\mathrm{erfi}(x) = -i\mathrm{erf}(ix)$, where $\mathrm{erfi}(x)$ is the
imaginary error function, both valid when $x > 0$. In particular, $E_1(x)$
for real negative argument is not currently supported in Scipy, which will
be used for implemented each of the proposed methods.
from scipy.special import expi, erfi
from scipy import exp, pi
def goodwin_staton_1(x):
mx2 = -x*x
ans = 0.5 * exp(mx2) * (pi * erfi(x) - expi(-mx2))
return ans
When $x$ is small the previous equation can be safely evaluated, but for
$x \gtrsim 20$ the computation will overflow. This problem is avoidable if we
consider some different approaches. A sensible approach consists
of substituting the first term using a well-known relation between the
imaginary error function and Dawson integral $F(x)$,
$$\begin{equation}\frac{\pi \mathrm{erfi}(x)}{2e^{x^2}} = \sqrt{\pi}F(x) \end{equation}$$
On the other hand, the second term involving the exponential integral $E_1(-x^2)$
can be replaced by its asymptotic expansion as follows,
$$\begin{equation} \frac{e^{-x^2} E_1(-x^2)}{2} \sim -\frac{1}{2x^2}
\sum_{k=0}^{\infty} \frac{k!}{x^{2k}},\quad x \to \infty \end{equation}$$
Therefore, for large argument the Goodwin-Staton integral can be computed as
$$\begin{equation} G(x) \sim \sqrt{\pi} F(x) -\frac{1}{2x^2}
\sum_{k=0}^{\infty} \frac{k!}{x^{2k}},\quad x \to \infty \end{equation}$$
from scipy.special import dawsn
from scipy import pi, exp, sqrt
def goodwin_staton_2(x):
# compute Dawson's integral
term1 = sqrt(pi) * dawsn(x)
# compute series
maxiter = 1000
eps = 2.2e-16
acc = False
s = 1.0
sp = s
n = 1
d = 1
x2 = x*x
for k in range(1, maxiter):
n *= k
d *= x2
s += n / d
if abs((s - sp) / sp) < eps:
if acc:
break
else:
acc = True
else:
sp = s
ans = term1 - 0.5 * s / x2
return ans
Another alternative is representing the second term as definite integral,
$$\begin{equation} \frac{e^{-x^2} E_1(-x^2)}{2} = - \frac{1}{2}\int_0^1
\frac{dt}{x^2 + \log(t)}
\end{equation}$$
Therefore, the Goodwin-Station is represented as follows,
$$\begin{equation} G(x) = \sqrt{\pi}F(x) - \frac{1}{2}\int_0^1 \frac{dt}
{x^2 + \log(t)} \end{equation}$$
from scipy.integrate import quad
from scipy.special import dawsn
from scipy import log, inf, sqrt, pi
def integrand(t, x):
return 1 / (x * x + log(t))
def term2_integral(x):
eps = 2.2e-16
return quad(integrand, 0, 1, args=x, epsabs=eps, epsrel=eps, limit=100)
def goodwin_staton_3(x):
# compute Dawson's integral
term1 = sqrt(pi) * dawsn(x)
# compute integral
term2, err = term2_integral(x)
return term1 - 0.5 * term2, err
Finally, we can evaluate the integral using an adaptive quadrature algorithm
implemented in scipy.integrate.quad, based on the Fortran library QUADPACK.
from scipy.integrate import quad
from scipy import exp, inf
def integrand(t, x):
return exp(-t**2) / (t + x)
def goodwin_staton_4(x):
eps = 2.2e-16
return quad(integrand, 0, inf, args=x, epsabs=eps, epsrel=eps, limit=100)
In order to perform a check of correctness, we evaluate the integral definition
using the tanh-sinh or double exponential integration algorithm implemented in
mpmath with 20 digits of precision. For example,
>>> from mpmath import *
>>> mp.dps = 20
>>> z = 10
>>> quad(lambda t: exp(-t*t)/(t + z), [0, inf])
mpf('0.084021593706602168771586')
