Integral of product of error function difference

Question

In the course of my research I came across the following integral:

$$\int_{-\infty}^{\infty}\,\left(\operatorname{erf}\left(ax-b\right) -\operatorname{erf}\left(\frac{a}{\gamma} x-b- \dfrac{ar}{\gamma}\right) \right) * \left(\operatorname{erf}\left(cx-d\right)-\operatorname{erf}\left(\frac{c}{\gamma}x-d-\frac{cr}{\gamma}\right)\right) \, dx$$
, where $\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^2}dt$ is the error function

Does anyone know if the integration given above can be calculated into closed-form?

Thank you in advance.

Answer 1

Let us take $a \ge 0$, $b \ge 0$, $c \ge 0$, $d \ge 0$, $\gamma \ge 0$ and $r \ge 0$ and let us define:

\begin{eqnarray} {\mathfrak I}^{(a,b)}_{c,d}(\gamma,r):= \int_{-\infty}^{\infty}\,\left(\operatorname{erf}\left(ax-b\right) -\operatorname{erf}\left(\frac{a}{\gamma} x-b- \dfrac{ar}{\gamma}\right) \right) * \left(\operatorname{erf}\left(cx-d\right)-\operatorname{erf}\left(\frac{c}{\gamma}x-d-\frac{cr}{\gamma}\right)\right) \, dx \end{eqnarray} Then by differentiating with respect to the parameter $b$ and using the identity: \begin{equation} \int\limits_{-\infty}^\infty e^{-u^2} \text{erf}(A u + B) du=\sqrt{\pi} \text{erf} \left( \frac{B}{\sqrt{1+A^2}}\right) \end{equation} we get: \begin{eqnarray} &&\frac{\partial }{\partial b} {\mathfrak I}^{(a,b)}_{c,d}(\gamma,r) = \frac{2}{a} \left( \right. \\ && -\text{erf}\left(\frac{\frac{c (a r-b)}{a g}+d}{\sqrt{\frac{c^2}{a^2 g^2}+1}}\right)+\text{erf}\left(\frac{d-\frac{b c}{a}}{\sqrt{\frac{c^2}{a^2}+1}}\right)+g \left(\text{erf}\left(\frac{d-\frac{b c}{a}}{\sqrt{\frac{c^2}{a^2}+1}}\right)-\text{erf}\left(\frac{d-\frac{c (a r+b g)}{a}}{\sqrt{\frac{c^2 g^2}{a^2}+1}}\right)\right) \\ &&\left.\right) \end{eqnarray}

Now all we need to do is to integrate with respect to $b$ from minus infinity to $b$. Since $\int \text{erf}(x) dx = \exp(-x^2)/\sqrt{\pi} + x \text{erf}(x)$ and since in the expression for the partial derivative the parameter $b$ enters the error function in a linear manner only it is clear that it is possible to construct the anti-derivative. Having done this we just take the values at $b$ and at minus infinity . It turns out that the later value is equal to zero and therefore we are left with the value at $b$ only which we simplified by hand. Now the final result reads: \begin{eqnarray} &&{\mathfrak I}^{(a,b)}_{c,d}(\gamma,r) = \frac{2}{a c} \left( \right. \\ && -\frac{(g+1) \sqrt{a^2+c^2} e^{-\frac{(b c-a d)^2}{a^2+c^2}}}{\sqrt{\pi } } + \frac{\sqrt{a^2 g^2+c^2} e^{-\frac{(a c r+a d g-b c)^2}{a^2 g^2+c^2}}+\sqrt{a^2+c^2 g^2} e^{-\frac{(a c r-a d+b c g)^2}{a^2+c^2 g^2}}}{\sqrt{\pi } }+\\ && (g+1) (a d-b c) \text{erf}\left(\frac{b c-a d}{\sqrt{a^2+c^2}}\right)+\\ && (a c r-a d+b c g) \text{erf}\left(\frac{a c r-a d+b c g}{\sqrt{a^2+c^2 g^2}}\right)+(a c r+a d g-b c) \text{erf}\left(\frac{a c r+a d g-b c}{\sqrt{a^2 g^2+c^2}}\right)\\ \left. \right) \end{eqnarray}

For[count = 1, count <= 100, count++,
  {a, b, c, d, g, r} = RandomReal[{0, 1}, 6, WorkingPrecision -> 50];
  I1 = NIntegrate[(Erf[(a x - b)] - 
       Erf[(a/g x - b - a r/g)]) (Erf[c x - d] - 
       Erf[c/g x - d - c r/g]), {x, -Infinity, Infinity}, 
    WorkingPrecision -> 15];
  2/Abs[a] NIntegrate[ (-Erf[(xi c - d a)/Sqrt[a^2 + c^2]] + 
      Erf[(-d a g + c (xi - a r))/Sqrt[a^2 g^2 + c^2]] + 
      Abs[g] (-Erf[(xi c - d a)/Sqrt[a^2 + c^2]] + 
         Erf[(-d a + c (xi g + a r))/Sqrt[
          a^2 + c^2 g^2]])), {xi, -Infinity, b}];
  I2 = 2/(
    Abs[a] c) (-E^(-((b c - a d)^2/(a^2 + c^2))) Sqrt[a^2 + c^2] /
       Sqrt[\[Pi]] (1 + g) + 
      1/ Sqrt[\[Pi]] (E^(-((-b c + a d g + a c r)^2/(c^2 + a^2 g^2)))
           Sqrt[c^2 + a^2 g^2] + 
         E^(-((-a d + b c g + a c r)^2/(a^2 + c^2 g^2))) Sqrt[
          a^2 + c^2 g^2]) + (a d - b c) (1 + g) Erf[(b c - a d)/Sqrt[
        a^2 + c^2]] + (-a d + b c g + a c r) Erf[(-a d + b c g + 
         a c r)/Sqrt[
        a^2 + c^2 g^2]] + (-b c + a d g + a c r) Erf[(-b c + a d g + 
         a c r)/Sqrt[c^2 + a^2 g^2]]);
  If[Abs[I2/I1 - 1] > 10^(-3), 
   Print["results do not match", {a, b, c, d, g, r, {I1, I2}}]; 
   Break[]];
  If[Mod[count, 10] == 0, PrintTemporary[count]];
  ];