Yet another improper integral which involves a Gaussian and an error function.

Question

Let $ \alpha \ge 0 $ and let $ A \in {\mathbb R} $ and $ B \in {\mathbb R} $ and $ u_* \in {\mathbb R} $. Consider a following integral:

\begin{equation} {\mathfrak I}_{u_*}^{\alpha,A,B}:= \int\limits_{u_*}^\infty u^\alpha \operatorname{erf}(A+B u) \frac{e^{-u^2/2}}{\sqrt{2\pi}} du \end{equation}

We have compute this quantity for in the special case when $\alpha=0$. We have:

\begin{eqnarray} &&{\mathfrak I}_{u_*}^{0,A,B}=\left(\right. \\ && \left. 2 \pi T\left(\frac{A}{\sqrt{B^2+\frac{1}{2}}},\frac{2 A B+2 B^2 u_*+u_*}{\sqrt{2} A}\right)+2 \pi T\left(u_*,\frac{\sqrt{2} (A+B u_*)}{u_*}\right) \right. \\ % && \left. -sign(A \cdot u_*) \cdot \frac{\pi}{2} \right. \\ % && \left. \right)\frac{1}{\pi}+ \frac{1}{2} \text{erf}\left(\frac{A}{\sqrt{2 B^2+1}}\right) \end{eqnarray}

Here $T(,)$ is the Owen's T function.

The result above has been obtained by differentiating the definition with respect to $A$, then doing the integration over $u$ and then integrating the result with respect to $A$ from $A$ to plus infinity and using the result from here. As always, the Mathematica code snippet verifies the result.

In[5169]:= 
ll = {};
For[II = 1, II <= 100, II++,
  {A, B, us} = RandomReal[{-10, 10}, 3];
  x1 = NIntegrate[Erf[A + B u] phi[u], {u, us, Infinity}];
x2 = 1/2 Erf[A /Sqrt[1 + 2 B^2]] +

    1/ [Pi] (-Sign[A us] Pi/2 + 
       2 [Pi] OwenT[A/Sqrt[1/2 + B^2], (2 A B + us + 2 B^2 us)/(
         Sqrt[2] A)] + 2 [Pi] OwenT[us, (Sqrt[2] (A + B us))/us]);
  ll = Join[ll, {x1 - x2}]
  ];
ll // Chop
Out[5171]= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 5.0484210^-10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1.2706610^-8, 

0, 0, 0, -1.1089510^-9, 0, -2.8541410^-10, 0, 0, 0, 0, 0, 0, 0, 0, 

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

0, 2.0869810^-10, 0, -1.8031910^-9, 0, 0, 0, -1.23154*10^-8, 0, 0, 

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

0, 0, 0}

Now, the question is pretty obvious, how do we derive the result in the generic case $\alpha \neq 0 $.

Answer 1

I think I have an answer, although, it isn't pretty...

$${1 \over {\sqrt{2\pi}}}\int_{u_{*}}^{\infty}{u^{\alpha}}{\left(erf(A+Bu)\right)}e^{{-u^2}\over{2}}du$$

$$={\sqrt{2}\over \pi}\int_{u_{*}}^{\infty}{u^{\alpha}}\left(\int_{0}^{A+Bu}e^{-t^2}dt\right)e^{{-u^2}\over{2}}du$$

$$={\sqrt{2}\over \pi}\int_{u_{*}}^{\infty}{u^{\alpha}}\left(\sum_{n=0}^{\infty}{(-1)^n \over n!}\int_{0}^{A+Bu}{t^{2n}}dt\right)e^{{-u^2}\over{2}}du$$

$$={\sqrt{2}\over \pi}\int_{u_{*}}^{\infty}{u^{\alpha}}\left(\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}{(A+Bu})^{2n+1}\right)e^{{-u^2}\over{2}}du$$

So now we have;

$${\sqrt{2}\over \pi}\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}\int_{u_{*}}^{\infty}{u^{\alpha}}\cdot {(A+Bu})^{2n+1}\cdot e^{{-u^2}\over{2}}du$$

We can expand ${(A+Bu})^{2n+1}$ thusly;

$${\sqrt{2}\over \pi}\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}\int_{u_{*}}^{\infty}{u^{\alpha}}\cdot\sum_{k=0}^{2n+1}{2n+1 \choose k} (A)^{2n+1-k}(Bu)^{k}\cdot e^{{-u^2}\over{2}}du$$

$$={\sqrt{2}\over \pi}\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}\sum_{k=0}^{2n+1}{2n+1 \choose k} (A)^{2n+1-k}(B)^{k}\int_{u_{*}}^{\infty}(u)^{k+\alpha}\cdot e^{{-u^2}\over{2}}du$$

Applying the substitution $j={u^2 \over 2}$ gives us;

$$={1\over \pi}\sum_{n=0}^{\infty}{(-1)^n \over (2n+1)\cdot n!}\sum_{k=0}^{2n+1}{2n+1 \choose k} (A)^{2n+1-k}(B)^{k}\cdot {(\sqrt{2}})^{k+\alpha}\int_{u_{*}^{2}\over2}^{\infty}(j)^{{k+\alpha-1}\over 2}\cdot e^{-j}dj$$

The remaining integral is simply a special case of the upper incomplete gamma function $\small \Gamma(s,x)= \int_{x}^{\infty}{t^{(s-1)}e^{-t}dt}$;

So, finally, we have the following;

$$\mathfrak{I}_{u_{*}}^{{\alpha,A,B}}={A\over \pi}\sum_{n=0}^{\infty}\sum_{k=0}^{2n+1}{(-1)^{n}\cdot {(\sqrt{2}})^{k+\alpha} \cdot (2n)! \over n! \cdot k! \cdot (2n+1-k)!} (A)^{2n-k}(B)^{k}\cdot \Gamma\left({k+\alpha+1\over2},{u_{*}^{2}\over2}\right)$$

I'm not sure if you specifically wanted a final result that made use of Owen's $T$ function but maybe this disaster of a derivation will prove somewhat useful to you...

Answer 2

If $\alpha \in {\mathbb N}$ then we can obtain the answer readily by inserting an exponential with a parameter into the integrand and differentiating by that parameter. This gives:

\begin{equation} {\mathfrak J}^{n,A,B}_{u_*} = \left. \frac{d^n}{d \theta^n} e^{\frac{\theta^2}{2}} {\mathfrak J}^{0,A+B \theta,B}_{u_*-\theta} \right|_{\theta=0} \end{equation}

Then we can simplify the above further using the chain rule.

In[5216]:= ll = {};
For[II = 1, II <= 100, II++,
  {A, B, us} = RandomReal[{-10, 10}, 3];
  n = RandomInteger[{0, 10}]; th =.;
  x1 = NIntegrate[u^n Erf[A + B u] phi[u], {u, us, Infinity}];
  J0[A_, B_, us_] := 
   1/2 Erf[A /Sqrt[1 + 2 B^2]] +  
    1/ \[Pi] (-Sign[A us] Pi/2 + 
       2 \[Pi] OwenT[A/Sqrt[1/2 + B^2], (2 A B + us + 2 B^2 us)/(
         Sqrt[2] A)] + 2 \[Pi] OwenT[us, (Sqrt[2] (A + B us))/us]);
  ex = D[Exp[th^2/2] (J0[A + B th, B, us - th]), {th, n}] /. th :> 0;
  x2 = ex /. Derivative[n_][Sign][a_] :> 0 /; a != 0;
  ll = Join[ll, {x1 - x2}];
  ];
ll // Chop
Out[5218]= {-1.8604210^-8, 0, -3.8605810^-9, 0, 
 4.9047310^-10, 0, 0, 0, 0, 0, 0, 0, 0, 
 1.5492510^-9, 0, 0, 0, 0, 0, 0, -6.6543710^-9, -1.4187810^-10, 0, 

0, 0, 0, -8.847510^-10, 1.7997810^-10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
 3.1364310^-8, 0, -2.9410^-8, 0, -2.4390410^-9, 0, 0, 0, 0, 
 4.5184710^-8, 0, 0, 0, 0, 1.0278410^-9, 0, 0, 
 3.8722910^-9, 0, 0, -1.5508410^-8, 
 1.0679110^-10, 0, 0, 0, 0, 0, 0, 0, 0, -5.5492810^-9, 0, 0, 0, 

-2.6521910^-9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

-1.7591610^-10, 0, 0, 0, 1.0023510^-9, 0, 0, 0, 0, 0, 0, 0}