# Integral involving numerous erf functions

As a part of a bigger problem, I am puzzled with computing:

$$\int_0^\infty e^{-x^2}\cdot \operatorname{erf}(s_1 x)\cdot \operatorname{erf}(s_2 x)\cdot \operatorname{erf}(s_3 x)\cdot \ldots \cdot \operatorname{erf}(s_n x) \; \mathrm{d}x.$$

I have found it done by Briggs, 2003 for one $\operatorname{erf}$ function, and here: http://mathworld.wolfram.com/Erf.html eq.34 for two $\operatorname{erf}$ functions.

Decent approximation would be also appreciated if closed form is absent. Thanks in advance.

• Welcome to MSE. Use mathjax for your math formatting. – Jakobian Aug 18 '18 at 23:01
• Looks like equation $35$ in Math World does it for three error functions. – saulspatz Aug 18 '18 at 23:47
• @saulspatz It looks to me more like $34$, since $35$ has a factor of $x$ in the integral. – Theo Bendit Aug 19 '18 at 1:12
• Thank you, initially I took a LaTex formatted expression, but could not see preview button here, so I went to the plain text. Will check the above format options. – TeamK Aug 19 '18 at 7:42
• Oops, I overlooked that. My bad. – saulspatz Aug 19 '18 at 12:44

Define: $$I\left( {{s}_{1}},...,{{s}_{n}} \right)=\int_{0}^{\infty }{{{e}^{-{{x}^{2}}}}erf\left( {{s}_{1}}x \right)\cdots erf\left( {{s}_{n}}x \right)dx}$$

Show that: \begin{align} & \frac{\partial I}{\partial {{s}_{n}}...\partial {{s}_{1}}}={{\left( \frac{2}{\sqrt{\pi }} \right)}^{n}}\int_{0}^{\infty }{{{x}^{n}}{{e}^{-\left( 1+s_{1}^{2}+s_{2}^{2}+\cdots +s_{n}^{2} \right){{x}^{2}}}}dx} \\ & \quad \quad \quad \ \ ={{\left( \frac{2}{\sqrt{\pi }} \right)}^{n}}\frac{\Gamma \left( \left( n+1 \right)/2\ \right)}{2\sqrt{{{(1+s_{1}^{2}+s_{2}^{2}+\cdots +s_{n}^{2})}^{1+n}}}} \\ \end{align}

Can you find $I\left( {{s}_{1}},...,{{s}_{n}} \right)$?It is not easy.

• Thanks for posting, but I can't follow where this derivative leads. Owen's T function was handy with one erf function and I tried to put that in use for couple of erf functions. – TeamK Aug 19 '18 at 7:53
• you have a point here, but this derivative implies that for $n>2$ there is no closed form because repeated integration can't be expressed in terms of elementary functions – user547564 Aug 19 '18 at 18:58
• I see, then the better approximation remains... – TeamK Aug 19 '18 at 19:25
• indeed, unfortunately – user547564 Aug 19 '18 at 19:28
• Guys, would you recommend me a proper approximation of the erf function, or better idea to solve it numerically if required precision of the solution is given? – TeamK Aug 24 '18 at 8:37

Here we give the result for $$n=3$$. Let $$(s_1,s_2,s_3) \in {\mathbb R}_+^3$$. We note the following result:

$$\begin{eqnarray} &&\int \text{arccos}(\frac{s}{r \sin(\theta)}) \sin(\theta) d\theta = \\ && \frac{s \text{arctan}\left(\frac{\sqrt{r^2} \cos (\theta)}{\sqrt{r^2 \sin ^2(\theta)-s^2}}\right)}{r}-\text{arctan}\left(\frac{\sqrt{s^2} \cos (\theta)}{\sqrt{r^2 \sin ^2(\theta)-s^2}}\right)-\cos (\theta) \text{arccos}\left(\frac{s \csc (\theta)}{r}\right) \quad (i) \end{eqnarray}$$

We have: $$\begin{eqnarray} &&I(s_1,s_2,s_3)=\\ &&(\frac{2}{\sqrt{\pi}})^n \int\limits_0^{s_1} \int\limits_0^{s_2} \int\limits_0^{s_3} \frac{\Gamma((n+1)/2)} {2 \sqrt{(1+t_1^2+t_2^2+t_3^2)^{1+n}} } dt_1 dt_2 d t_3=\\ &&(\frac{2}{\sqrt{\pi}})^n \int\limits_0^{\sqrt{s_1^2+s_2^2+s_3^2}} \int\limits_{0}^{\pi/2} \int\limits_{0}^{\pi/2} \frac{\Gamma((n+1)/2)} {2 \sqrt{(1+r^2)^{1+n}} } r^2 \sin(\theta) 1_{0 < r \sin(\theta) \cos(\phi) < s_1} 1_{0 < r \sin(\theta) \sin(\phi) < s_2} 1_{0< r \cos(\theta) < s_3} d\phi d\theta dr=\\ &&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!(\frac{2}{\sqrt{\pi}})^n \int\limits_0^{\sqrt{s_1^2+s_2^2+s_3^2}} \int\limits_{0}^{\pi/2} \frac{\Gamma((n+1)/2)} {2 \sqrt{(1+r^2)^{1+n}} } r^2 \sin(\theta) \cdot \text{max}\left(\text{arccos}(\frac{s_1}{r \sin(\theta)}) - \text{arcsin}(\frac{s_2}{r \sin(\theta)}) ,0\right) 1_{0< r \cos(\theta) < s_3} d\theta dr=\\ && \frac{1}{\sqrt{\pi}}\left( \arctan(s_1)+\arctan(s_2)+\arctan(s_3) \right)+\\ &&\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\frac{2 \imath }{\pi^{3/2}} \cdot \sum\limits_{1 \le i < j \le 3}\frac{s_i^2+s_j^2}{s_i s_j} \int\limits_0^{\sqrt{\frac{\sqrt{s_1^2+s_2^2+s_3^2}-\sqrt{s_i^2+s_j^2}} {\sqrt{s_1^2+s_2^2+s_3^2}+\sqrt{s_i^2+s_j^2}}}} \log\left( \frac{1-t^2+\imath (1+t^2) \sqrt{s_i^2+s_j^2}}{1-t^2-\imath (1+t^2) \sqrt{s_i^2+s_j^2}}\right) \cdot \frac{(1-t^2)(1+6 t^2+t^4)}{((1+t^2)^2+4 t^2 \frac{s_j^2}{s_i^2})((1+t^2)^2+4 t^2 \frac{s_i^2}{s_j^2})} dt \end{eqnarray}$$

In the top line we used the result for the partial derivative given in the first answer above, in the second line we went to spherical coordinates, in the third line we integrated over $$\phi$$ and in the final line we used $$(i)$$. Note that the integral above can be expressed through di-logarithms by decompositing the rational function in the integrand into partial fractions. We will complete that mundane task asap.

Update: Define $$S:=\sqrt{s_1^2+s_2^2+s_3^2}$$ and $$k_{i,j}=3 \cdot 1_{{i,j}={1,2}} + 2 \cdot 1_{{i,j}={1,3}} + 1 \cdot 1_{{i,j}={2,3}}$$. Then also define:

$$\begin{eqnarray} {\mathfrak F}^{(A,B)}_{a,b} &:=& \int\limits_A^B \frac{\log(z+a)}{z+b} dz\\ &=& F[B,a,b] - F[A,a,b] + 1_{t^* \in (0,1)} \left( -F[A+(t^*+\epsilon)(B-A),a,b] + F[A+(t^*-\epsilon)(B-A),a,b] \right) \end{eqnarray}$$ where $$\begin{eqnarray} t^*:=-\frac{Im[(A+b)(b^*-a^*)]}{Im[(B-A)(b^*-a^*)]} \end{eqnarray}$$ and $$$$F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right)$$$$ for $$a$$,$$b$$,$$A$$,$$B$$ being complex.

Then the final result is: $$\begin{eqnarray} &&I(s_1,s_2,s_3)=\\ &&\frac{1}{\sqrt{\pi}} \sum\limits_{i=1}^3 \arctan(s_i) -\frac{2}{\pi^{3/2}}\sum\limits_{1 \le i < j \le 3} \left(\pi - \arctan(\sqrt{s_i^2+s_j^2}) \right)\left( \arctan(\frac{s_j s_{k_{i,j}}}{s_i S}) + \arctan(\frac{s_i s_{k_{i,j}}}{s_j S})\right)+\\ &&-\frac{1}{2\pi^{3/2}} \sum\limits_{1 \le i < j \le 3} \sum\limits_{\xi=1}^4 \sum\limits_{\eta=1}^8 (-1)^{\left\lfloor \frac{\eta -1}{4}\right\rfloor +\eta +\left\lfloor \frac{\xi -1}{2}\right\rfloor } {\mathfrak F}^{(0,\frac{S-\sqrt{s_i^2+s_j^2}}{s_{k_{i,j}}})}_{(-1)^{\xi } \exp \left(i (-1)^{\left\lfloor \frac{\xi -1}{2}\right\rfloor } \tan ^{-1}\left(\sqrt{s_i^2+s_j^2}\right)\right),\frac{i (-1)^{\eta +1} \left((-1)^{\left\lfloor \frac{\eta -1}{4}\right\rfloor } s_{\frac{1}{2} (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor +1} (j-i)+\frac{i+j}{2}}+\sqrt{s_i^2+s_j^2}\right)}{s_{\frac{1}{2} (-1)^{\left\lfloor \frac{\eta -1}{2}\right\rfloor } (j-i)+\frac{i+j}{2}}}} \end{eqnarray}$$

(*Definitions*)

Clear[F]; Clear[FF];
F[z_, a_, b_] :=
Log[a + z] Log[(b + z)/(-a + b)] + PolyLog[2, (a + z)/(a - b)];
FF[A_, B_, a_, b_] :=
Module[{result, ts, zs, zsp, zsm, eps = 10^(-15)},
(*This is Integrate[Log[z+a]/(z+b),{z,A,B}] where all a,b,A,
and B are complex. *)
result = F[B, a, b] - F[A, a, b];

ts = - (Im[(A + b) (Conjugate[b] - Conjugate[a])]/
Im[(B - A) (Conjugate[b] - Conjugate[a])]);
If[0 <= ts <= 1,
zsp = A + (ts + eps) (B - A);
zsm = A + (ts - eps) (B - A);
result += -F[zsp, a, b] + F[zsm, a, b];
];

result
];

n = 3; Clear[II];
II[s_] :=
NIntegrate[
Exp[-x^2] Product[Erf[s[[j]] x], {j, 1, Length[s]}], {x, 0,
Infinity}];
For[count = 1, count <= 200, count++,
s = Table[0, {3}];
{s[[1]], s[[2]], s[[3]]} =
Sort[RandomReal[{0, 5}, 3, WorkingPrecision -> 50]];
(*s[[3]]=RandomReal[{Sqrt[s[[1]]^2+s[[2]]^2],10},
WorkingPrecision\[Rule]50];*)
mArcSin[x_] := If[x > 1, Pi/2, ArcSin[x]];
mArcCos[x_] := If[x > 1, 0, ArcCos[x]];
x1 = II[s];
(*(2/Sqrt[Pi])^nNIntegrate[Gamma[(n+1)/2]/(2Sqrt[(1+t1^2+t2^2+
t3^2)^(1+n)]),{t1,0,s[[1]]},{t2,0,s[[2]]},{t3,0,s[[3]]}];*)
(*(2/Sqrt[Pi])^nNIntegrate[Gamma[(n+1)/2]/(2Sqrt[(1+r^2)^(1+n)])
r^2If[0<r Sin[th] Cos[phi]<s[[1]]&&0<r Sin[th] Sin[
phi]\[LessEqual]s[[2]]&&0<r Cos[th]<s[[3]],Sin[th],0],{r,0,Sqrt[s[[
1]]^2+s[[2]]^2+s[[3]]^2]},{th,0,Pi/2},{phi,0,Pi/2},
WorkingPrecision\[Rule]15]*)
k[i_, j_] :=
Which[{i, j} == {1, 2}, 3, {i, j} == {1, 3}, 2, {i, j} == {2, 3},
1];
S = Sqrt[s[[1]]^2 + s[[2]]^2 + s[[3]]^2];
x2 = (ArcTan[s[[1]]] + ArcTan[s[[2]]] + ArcTan[s[[3]]])/(\[Pi]^(
1/2))  + -2/\[Pi]^(3/2)
Sum[ (- ArcTan[Sqrt[s[[i]]^2 + s[[j]]^2]] +
Pi) (ArcTan[( s[[k[i, j]]] s[[j]] )/( S s[[i]]  )] +
ArcTan[( s[[k[i, j]]] s[[i]] )/( S s[[j]])]) , {i, 1, 3}, {j,
i + 1, 3}] +
-(1 /(2 \[Pi]^(3/2)))
Sum[ (Sum[(-1)^Floor[(xi - 1)/2] (-1)^(
eta + Floor[(eta - 1)/4])
FF[0, (-Sqrt[s[[i]]^2 + s[[j]]^2] + S)/
s[[k[i, j]]], (-1)^
xi E^(((-1)^Floor[(xi - 1)/2])
I ArcTan[Sqrt[s[[i]]^2 + s[[j]]^2]] ),
I (-1)^(eta +
1) ((-1)^
Floor[(eta - 1)/
4] s[[(-1)^(1 + Floor[(eta - 1)/2]) (j - i)/
2 + (j + i)/2]] + Sqrt[(s[[i]]^2 + s[[j]]^2)])/
s[[(-1)^Floor[(eta - 1)/2] (j - i)/2 + (j + i)/2]]], {xi,
1, 4}, {eta, 1, 8}]) , {i, 1, 3}, {j, i + 1, 3}];
If[Abs[x2/x1 - 1] > 10^(-3),
Print["results do not match ..", count, s, {x1, x2}]; Break[]];

If[Mod[count, 10] == 0, PrintTemporary[count]];
];