How can I avoid roundoff error when calculating the difference $\textrm{erfc}(a) - \textrm{erfc}(b)$?

Question

In this excellent answer, it is recommended that one make use of the $\textrm{erfcx}$ function to avoid roundoff error in calculating dealing with $x < 25$ (approximately). So, one scales their answer by $e^{x^2}$, and then scales back, I suppose (?) to recover the number they are really looking for?

What if we have the case where we need to calculate $\textrm{erfc}(a) - \textrm{erfc}(b)$? I could write:

$$\textrm{erfcx}(a) - \textrm{erfcx}(b) = e^{a^2}\textrm{erfc}(a) - e^{b^2}\textrm{erfc}(b) $$

but that's not particularly helpful if I want to numerically integrate the function ($x, a, b \in \mathbb{R}$):

$$ \int_{a}^{b} e^{x^2}(\textrm{erf}(x) - \textrm{erf(a)})\;\textrm{d}x = \int_{a}^{b} e^{x^2}(\textrm{erfc}(a) - \textrm{erfc(x)})\;\textrm{d}x$$

If $a$ and $b$ are close enough/large enough in magnitude (not large in terms of sign) then we get catastrophic cancellation. For example, try $a = -16.85$ and $b = -16.08$. Rescaling inside the integrand would seem to be unhelpful, rather than helpful!

What do?

If you have Python, here's the toy problem script on Github Gist. It implements many of the recommendations provided in this answer to avoid overflow errors, but still runs into underflow problems.

Answer 1

I am not sure exactly what you are trying to do. I am assuming you want to evaluate the integral using some sort of quadrature. This will work for that.

Note that we have

$$ \begin{align*} \int_{x=a}^b \exp(x^2) (\mathrm{erf}(x) - \mathrm{erf}(a))\,\mathrm{d} x &= \frac{2}{\sqrt{\pi}} \int_{x=a}^b \exp(x^2) \int_{\tau = a}^x \exp(-\tau^2) \, \mathrm{d} \tau\, \mathrm{d} x\\ &= \frac{2}{\sqrt{\pi}} \int_{x=a}^b \int_{\tau= a}^x \exp(x^2-\tau^2) \mathrm{d} \tau\, \mathrm{d}x \end{align*} $$

And in the triangle of integration, the integrand is strictly positive so you shouldn't have any cancellation.