Integral involving numerous erf functions

Question

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.

Answer 1

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.

Answer 2

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 \begin{equation} F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + Li_2\left( \frac{z+a}{a-b}\right) \end{equation} 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]];
  ];