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 ansWhen $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 ansAnother 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, errFinally, 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')