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

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.

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.

• I tried this out, and in my excitement thought that it worked, but it doesn't seem to. Even with this method, I simply get an integral result of "zero", which makes me think things were so small during the process, that they got rounded off to zero? here's the wolfram link: wolframalpha.com/input/?i=Integrate%5Be%5E%28x%5E2%29*e%5E%28-y%5E2%29%2C+%7Bx%2C+-16.85%2C+-16.08%7D%2C+%7By%2C+-16.85%2C+x%7D%5D – user89 Oct 15 '15 at 23:53
• fred, here's the link fixed up; I also tried this out in Matlab, and Matlab gives a small number (on the order of 1e-7), so perhaps Matlab doesn't suffer from Wolfram Alpha's issues? – user89 Oct 16 '15 at 0:40
• I think it was a Wolfram issue. Matlab works fine. – user89 Oct 16 '15 at 2:34
• I am a little nervous about this. I had assumed that $a$ and $b$ were both positive. Can you explain what you need this for? You probably need to have some sanity checks before you take Matlab's answer as the gospel truth. I gave the integral to Mathematica and it said it was 0.0293917959... – fred Oct 16 '15 at 3:33
• That's what Matlab gave it as too! Wolfram Alpha can't handle it though. I am not sure if I can really explain what I need this for, because it may take more space than the comments provide ("I am calculating the mean time it takes for something to happen, given a force"). $a$ and $b$ are not necessarily positive, but they are necessarily negative -- they depend on "force". So, basically I have a mean time given force function, and for various values of the force, I can plot out the mean time. The derivative of this function then tells me how the mean time changes as the force changes. – user89 Oct 16 '15 at 3:41