Computing the Goodwin-Staton integral

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')