# Integral of product of error function difference

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?

• I doubt it. Have you tried Wolfram|Alpha? By the way, you can get proper formatting for functions like $\operatorname{erf}$ using the syntax \operatorname{erf}; you can get displayed equations (that have nicer integrals and fractions) by enclosing them in double instead of single dollar signs; and you can get proper parentheses that adapt to their content by preceding them with \left and \right. More generally, here's a basic tutorial and quick reference on how to typeset math on this site: math.meta.stackexchange.com/questions/5020. – joriki Jul 5 '18 at 8:36
• @joriki Thank you for your advice. I tried this on Wolframa Alpha and also on Mathematica, but it failed – Yunho Choi Jul 5 '18 at 12:14

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: $$$$\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)$$$$ 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]];
];