Integral involving InverseErf: closed-form expression and asymptotics

Consider the integral

f[n_] := Integrate[s^(n - 1) InverseErf[1 - 2 s], {s, 0, 1},
Assumptions -> n > 0]


a) find a closed form expression of $f$ if possible
b) calculate the exact values of $f(n)$ for $n = 1 .. 6$
c) calculate the asymptotic behaviour of $f(n)$ for $n\to \infty$

J. M. pointed out in a comment that the integral $f$ is equivalent to this one:

f1[n_] :=  2^(1 - n)/Sqrt[\[Pi]]*
Integrate[u Exp[-u^2] Erfc[u]^(n - 1), {u, -\[Infinity], \[Infinity]}]


Remark: $f1$ was my original expression before I found the shorter $f$. But $f1$ has the advantage of containing familiar functions in the integrand rather than the "exotic" InverseErf[].

• for n = 2 wil be: $-\frac{1}{2 \sqrt{2 \pi }}$ Jan 28, 2017 at 16:40
• FWIW, your integral is equivalent to $$\frac{2^{1-n}}{\sqrt{\pi}}\int_{-\infty}^\infty u\exp(-u^2)\operatorname{erfc}(u)^{n-1}\mathrm du$$ Jan 28, 2017 at 16:57
• @J. M. that's where it came from ;-) But maybe it is better to use this form as this Integrand contains better know functions. Jan 28, 2017 at 19:34
• Yes, this one is better behaved than the one with the inverse error function, numerically speaking. Jan 28, 2017 at 19:41
• I'd also recommend asking someone on math.stackexchange.com. There some real integration wizards on there who may certainly point you in the right direction. Jan 28, 2017 at 21:39

We provide a (partial) solution in these sections

• closed-form expressions
• numerical results
• asymptotic behaviour
• discussion

Closed-form expressions

Since the original integral is evaluated by Mathematica only for the first two values we shall look here for equivalent expressions for the integral and see if Mathematica performs better with these.

Equivalent expressions

First, for easier reading, we give the formulas in Latex, then we write down the Mathematica expressions.

Equivalent expressions for $f(n)$ are (the proof of the equivalences is left to the reader. I confine myself to brief hints, the complete proofs can be provided on request).

$$f1(n)={\frac{1}{\sqrt{\pi }}\int_{-\infty }^{\infty } u \exp \left(-u^2\right) \left(\frac{1}{2}\text{erfc}(u)\right)^{n-1} \, du}$$

Noticing that $\text{erfc} = 1 - \text{erf}$ and expanding the power into a binomial sum gives

$$f2(n) = - \frac{2^{1-n}}{\sqrt{\pi }} \sum _{m=1}^{\left\lfloor n/2\right\rfloor } g(2 m-1) \binom{n-1}{2m-1}$$

where

$$g(k)={\int_{-\infty }^{\infty } u \exp \left(-u^2\right) \text{erf}(u)^{k} \, du}$$

Notice that $g(k) = 0$ for even $k$. We have already taken this into account in $f2$ where only odd terms contribute.

The following forms are equivalent to $g$ for odd $k$ but they will be considered also for even $k$ in the following

$$g1(k)={2\int_{0 }^{\infty } u \exp \left(-u^2\right) \text{erf}(u)^{k} \, du}$$

Letting $u^2 \to t$ gives

$$g2(k)={\int_{0 }^{\infty } \exp \left(-t \right) \text{erf}(\sqrt{t})^{k} \, du}$$

The next form is a multiple integral over the $k$-dimensional hypercube obtained by inserting the integral representation of $\text{erf}$ into $g2$ and performing the $t$-Integration:

$$g3(k)=\left(\frac{2}{\sqrt{\pi }}\right)^k \Gamma \left(\frac{k}{2}+1\right) \int _0^1 ... \int _0^1\frac{1}{\left(\sum _{i=1}^k y(i)^2+1\right){}^{\frac{k}{2}+1}}dy(1) ... dy(k)$$

The corresponding Mathematica expressions (with the same names) are:

f1[n_] :=
1/Sqrt[\[Pi]]
Integrate[
u Exp[-u^2] (1/2 Erfc[u])^(n - 1), {u, -\[Infinity], \[Infinity]}]

f2[n_] :=
2^(1 - n)/Sqrt[\[Pi]]
Sum[Binomial[n - 1, k] (-1)^k g[k], {k, 1, n - 1,2}]

g[k_] := Integrate[u Exp[-u^2] Erf[u]^k, {u, -\[Infinity], \[Infinity]}]

g1[k_] := 2 Integrate[u Exp[-u^2] Erf[u]^k, {u, 0, \[Infinity]}]

g2[k_] := Integrate[Exp[-t] Erf[Sqrt[t]]^k, {t, 0, \[Infinity]}]

g3[k_] :=
(2/Sqrt[\[Pi]])^k Gamma[1 + k/2] Integrate[1/(1 + Sum[y[i]^2, {i, 1, k}])^(1 + k/2), Sequence @@ Table[{y[i], 0, 1}, {i, 1, k}]]


Results

We have shown that the calculation of $f$ can be traced back to evaluating the function $g$.

It turns out that $g3$ gives closed-form results for the highest $k$:

Table[{k, g3[k]}, {k, 1, 3}]


$$\left( \begin{array}{cc} 1 & \frac{1}{\sqrt{2}} \\ 2 & \frac{2 \sqrt{2} \tan ^{-1}\left(\frac{1}{\sqrt{2}}\right)}{\pi } \\ 3 & \frac{3 \sqrt{2} \cot ^{-1}\left(2 \sqrt{2}\right)}{\pi } \\ \end{array} \right)$$

Unfortunately, this "highest $k$ is still very modest. Already $g3(4)$ remains unevaluated.

Looking for a rule in the first 3 terms I came up with

gguess[k_] := k Sqrt[2] ArcTan[1/Sqrt[k^2-1]])/\[Pi]}


but it doesn't work for $k=4$ as can be seen by numerical comparison.

Finally, using $f2$ with $g3$ (called $f23$) we obtain the closed-form expressions for $f(n)$

{#, f23[#]} & /@ {1, 2, 3, 4, 5} // Simplify


$$\left( \begin{array}{cc} 1 & 0 \\ 2 & -\frac{1}{2 \sqrt{2 \pi }} \\ 3 & -\frac{1}{2 \sqrt{2 \pi }} \\ 4 & -\frac{3 \cos ^{-1}\left(-1/3\right)}{4 \sqrt{2} \pi ^{3/2}} \\ 5 & -\frac{3 \sqrt{2}}{4 \pi ^{3/2}}\cos ^{-1}\left(\sqrt{\frac{2}{3}}-\frac{1}{6}\right)\\ \end{array} \right)$$

numerically

% // N

(* Out[797]= {{1., 0.}, {2., -0.199471}, {3., -0.199471}, {4., -0.18197}, {5., -0.164468}} *)


Numerical results

Define the numerical analogue to our functions as

fn[n_] := NIntegrate[s^(n - 1) InverseErf[1 - 2 s], {s, 0, 1}]

f1n[n_] :=  1/Sqrt[\[Pi]]
NIntegrate[u Exp[-u^2] (1/2 Erfc[u])^(n - 1), {u, -\[Infinity], \[Infinity]}]


As expected $f1n$ is faster than $fn$;

AbsoluteTiming[fn[10^4]]

(* Out[430]= {1.88997, -0.00027235} *)

AbsoluteTiming[f1n[10^4]]

(* Out[431]= {0.022427, -0.00027235} *)


The first few discrete values are (excluding the trivial case $n = 1$), using

tfn[nn_] := Table[{n, fn[n] // Chop}, {n, 2, nn}]

tfn[10]

(* Out[676]=
{
{2, -0.199471}, {3, -0.199471}, {4, -0.18197}, {5, -0.164468}, {6,-0.149342}, {7, -0.136591}, {8, -0.12583}, {9, -0.116674}, {10, -0.108806}
}


Notice the agreement with the numerical values of the closed-form expressions for $n = 2..4$ which is also a verification of these.

Plotting it

ListPlot[tfn[20],
PlotLabel -> "Numerical values of f(n) for n = 2 .. 20",
AxesLabel -> {"n", "f(n)"}]


Asymptotic behaviour

The asymptotics for large $n$ is most conveniently calculated in the original form of $f$.

The integrand of $f$ is given by

fi[n_, s_] := s^(n - 1) InverseErf[1 - 2 s]


The "exotic" function InverseErf[] does not seem to have an official name (please correct my if I am wrong), hence for brevity we shall call it fre(x) and define

fre[x_] := InverseErf[x]


The function $\text{fre}$ is, of course, just $\text{erf}$ rotated by $\pi /4$ and looks like this

Plot[fre[x], {x, -1, 1}, PlotLabel -> "The function fre(x) = InverseErf[x]",
AxesLabel -> {"x", "fre(x)"}, PlotRange -> All]


The complete integrand for some small values of $n$ is shown here

Plot[{fi[2, s], fi[3, s], fi[5, s]}, {s, 0, 1},
PlotLabel -> "The integrand for some n\nn = 2 (blue), 3 (yellow), 5 (green)",
AxesLabel -> {"s", "fi(n,s)"}, PlotRange -> {-2, .2}]


We see that, as expected, for increasing $n$ the appreciable contribution to the integral concentrates around $s = 1$

Close to $s = 1$ the integrand is

Simplify[Series[s^(n - 1) fre[1 - 2 x], {x, 1, 1}] // Normal, 0 < x < 1]

(* Out[733]=
-s^(-1 + n) Sqrt[-Log[1 - x] -
1/2 Log[2 \[Pi] (Log[1/(2 \[Pi])] - 2 Log[1 - x])]]
*)


First let us look at the plot of $f$ for very large $n$ and try to find analytic expression suggested by the series expansion

Plot[{1, - f1n[n] n/Sqrt[(-Log[1/n] - Log[-Log[1/n]])], -f1n[n] n/Sqrt[
ProductLog[n]], - f1n[n] n/Sqrt[Log[n] ]}, {n, 2, 10^6},
PlotRange -> {0.85, 1.05}, AxesLabel -> {"n", "f/fappr"},
PlotLabel ->
"Behavior of f for large n\nrelative to various approximations\nfappr ~ \
asympt(ProductLog)/n (yellow)\nfappr ~ \!$$\*SqrtBox[\(ProductLog[n]$$]\)/n \
(green)\nfappr ~ \!$$\*SqrtBox[\(Log[n]$$]\)/n (red)"]


We see that, up to a factor close to unity the asymptotics is best described by

$$f(n)\sim \frac{1}{n}\sqrt{W(n)}$$

where $W(n)$ is the LambertW function (Lambert 1758, http://mathworld.wolfram.com/LambertW-Function.html ) which is the inverse of $n = W \exp (W)$.

An analytic derivation of the asymptotics will be given later.

Discussion

Relation to Laplace transform

The function $g2$ is in fact a Laplace transform

$$glp(n) = \mathcal{L}_t\left[\text{erf}\left(\sqrt{t}\right)^n\right](s)$$

at $s\to 1$.

The function $g$ in the form $g3$ is easily generalized to the Laplace form by replacing the $1$ before the sum by $s$, and it can be calculated for the powers 1, 2, and 3. Adding the recent discovery of Mariusz Iwaniuk for $g(4)$ (which, however, contains an unevaluated integral) we have

$$\left( \begin{array}{cc} 1 & \frac{1}{s \sqrt{s+1}} \\ 2 & \frac{4 \tan ^{-1}\left(\frac{1}{\sqrt{s+1}}\right)}{\pi s \sqrt{s+1}} = \frac{4}{\pi s \sqrt{s+1}}\sin ^{-1}\left(\frac{1}{\sqrt{s+2}}\right) \\ 3 & \frac{6 \tan ^{-1}\left(\frac{1}{\sqrt{s+1} \sqrt{s+3}}\right)}{\pi s \sqrt{s+1}} = \frac{6}{\pi s \sqrt{s+1}}\sin ^{-1}\left(\frac{1}{s+2}\right)\\ 4 & \frac{4 \sqrt{\frac{1}{s+1}}}{s}-\frac{96 \int_0^{\sqrt{s+1}} \frac{\cot ^{-1}\left(\sqrt{a^2+2}\right)}{\left(a^2+1\right) \sqrt{a^2+2}} \, da}{\pi ^2 s \sqrt{s+1}}\\ \end{array} \right)$$

In the published tables of Laplace transformations (see e.g. http://authors.library.caltech.edu/43489/1/Volume%201.pdf 4.12 (4)) I have found only the case $n = 1$, hence, as a byproduct of this study, we can add two closed-form expressions and one with an integral left to the tables.

EDIT 06.02.2017

Mariusz Iwaniuk pointed out in a comment that

$$g2(4,s) = \frac{4 \sqrt{\frac{1}{s+1}}}{s}-\frac{96 \int_0^{\sqrt{s+1}} \frac{\cot ^{-1}\left(\sqrt{a^2+2}\right)}{\left(a^2+1\right) \sqrt{a^2+2}} \, da}{\pi ^2 s \sqrt{s+1}}$$

and that the remaining integral is related to Ahmed's integral (https://arxiv.org/pdf/1411.5169.pdf).

The numerical check $g2(4,s=1) = g2n(4)$ is ok.

The integral we need is not Ahmed's integral but a generalization of it. The upper integration limit is not unity but $\sqrt{1+s}$.

EDIT 07.02.2017

In an attempt to "crack" the remaining integral I have defined some integrals appearing naturally in this context, and have derived relations between these (they include the relation of my comment).

The integrals are

$$M(s)=\int_0^{\sqrt{s+1}} \frac{\cot ^{-1}\left(\sqrt{y^2+2}\right)}{\left(y^2+1\right) \sqrt{y^2+2}} \, dy\tag{7.1}$$

$$B(s)= \int_0^1 \frac{\cot ^{-1}\left(\frac{\sqrt{y^2+2}}{\sqrt{s+1}}\right)}{\left(y^2+1\right) \sqrt{y^2+2}} \, dy\tag{7.2}$$

$$C(s) = \frac{1}{2} \int_0^s \frac{\cot ^{-1}\left(\sqrt{z+3}\right)}{\sqrt{z+1} (z+2) \sqrt{z+3}} \, dz\tag{7.3}$$

While trying to understand the recent EDIT of Mariusz Iwaniuk I found that the integral C can be transformed into a simpler form by the substitution $z\to \frac{1}{\sin (w)}-2$

$$C_1(s) = \frac{1}{2} \int_{\sin ^{-1}\left(\frac{1}{s+2}\right)}^{\frac{\pi }{6}} \cot ^{-1}\left(\sqrt{\csc (w)+1}\right) \, dw\tag{7.3.1}$$

The relations are

$$M(s)=\frac{1}{4} \pi \tan ^{-1}\left(\sqrt{s+1}\right)-B(s) \tag {7.4}$$

$$B(s)=\frac{1}{4} \pi \tan ^{-1}\left(\sqrt{s+1}\right)-\frac{\pi ^2}{32} -C(s)\tag{7.5}$$

and, combining these two gives the interesting relation

$$M(s)=\frac{\pi ^2}{32} + C(s)\tag{7.6}$$

It turns out that, for the case of the OP, i.e. $s=1$, the main part of $M(1)$ is given by the explicit fraction, and only about 13% are due to integral $C(1)$.

In fact, numerically for $s=1$ the equation reads (in the above order)

0.3540215 = 0.3084251+0.0455963, and the relation is 0.0455963/0.3540215 = 0.128795.

We can find a good approximate expression for $C$. In fact, since $\cot^{-1}$ is a decreasing function we can replace it in the integrand by its maximum value $\cot^{-1}\left(\sqrt{3}\right)$ to obtain an upper bound (i.e. the function in question is smalller). Now the integral can be performed, and the result for $M$ becomes

$$M(s) \lesssim \frac{\pi ^2}{288}+\frac{1}{6} \pi \tan ^{-1}\left(\sqrt{\frac{s+1}{s+3}}\right)$$

Notice: I have the derivation and some corollaries ready in draft state, but I have to beg for your patience as it takes some time for me to type it in.

EDIT 08.02.17: Generating function and difference-differential equation

Remember that we wish to calculate the Laplace-transformation of the n-th power of the Erf-function of a square root

$$glp(n,s)=\mathcal{L}_t\left[\text{erf}\left(\sqrt{t}\right)^n\right](s)$$

Here's a new idea: instead of looking into single powers we try to "untie the Gordian knot" and cover them all at once.

This is generally done using generating functions. It is convenient here to use the exponential generating function

$$h(z, s)=\int_0^{\infty } \exp (-s t) \exp \left(-z\; \text{erf}\left(\sqrt{t}\right)\right) \, dt\tag{8.1}$$

After some manipulations including partial integration combined with differentiation I came up with this equation for the generating function

$$ddeq=s (s+1) \frac{\partial h(z,s)}{\partial s}+\left(\frac{3 s}{2}+1\right) h(z,s) -\frac{z^2 h(z,s+2)}{\pi }=\frac{1}{2}\tag{8.2}$$

This is a linear inhomogeneous partial difference-differential equation of first order.

This type of equation is new to me. Hence I don't know yet how to approach a solution.

Some first observations

1) the "generating" parameter $z$ appears only in one position as a factor

2) if $z = 0$, the eqation reduces to a simple ODE which can readily be solved with a simple result

DSolve[ddeq /. z -> 0, h[0, s], s]

(* {{h[0, s] -> 1/s + C[1]/(s Sqrt[1 + s])}} *)


The first part is simply $glp(n=0,s)$, as it should be. It is interesting that the next term $glp(n=1,s)$ also appears automatically.

EDIT 18.10.17: Derivation of the difference-differential equation

It's my pleasure to comply with Mariusz Iwaniuk's recent request and present the details of the derivation of equ. (8.2).

The text needs some formatting, sorry.

The derivation consists of a chain of partial integrations and derivatives.

§1) We consider the expression

D[Exp[-s t] Exp[-z Erf[Sqrt[t]]], t] // Expand

-E^(-s t - z Erf[Sqrt[t]]) s - (E^(-t - s t - z Erf[Sqrt[t]]) z)/(
Sqrt[\[Pi]] Sqrt[t])


Integrating over t gives

Exp[-s t] Exp[-z Erf[Sqrt[t]]] /. t -> 0

1

Limit[Exp[-s t] Exp[-z Erf[Sqrt[t]]], t -> \[Infinity], Assumptions -> s > 0]

0


Hence

1 == s Integrate[E^(-s t - z Erf[Sqrt[t]]) , {t, 0, \[Infinity]}] +
z/Sqrt[\[Pi]]
Integrate[E^(-t - s t - z Erf[Sqrt[t]])/ Sqrt[t], {t, 0, \[Infinity]}]

(* 1 *) 1 ==
s h[z, s] +
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]])/ Sqrt[t], {t, 0, \[Infinity]}]


Differentiating the equation with respect to s gives

0 = D[s h[z, s], s] +
z/Sqrt[\[Pi]]
Integrate[
D[E^(-t - s t - z Erf[Sqrt[t]])/ Sqrt[t], s], {t, 0, \[Infinity]}]


but

D[E^(-t - s t - z Erf[Sqrt[t]])/ Sqrt[t], s]

-E^(-t - s t - z Erf[Sqrt[t]]) Sqrt[t]


inserting this gives

(* 2 *) 0 =
D[s h[z, s], s] -
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}]


§2 2nd partial integration

Consider (notice that we have s+1 now)

D[E^(-t (1 + s) - z Erf[Sqrt[t]]) Sqrt[t], t] // Expand

E^(-(1 + s) t - z Erf[Sqrt[t]])/(2 Sqrt[t]) -
E^(-(1 + s) t - z Erf[Sqrt[t]]) Sqrt[t] -
E^(-(1 + s) t - z Erf[Sqrt[t]]) s Sqrt[t] - (
E^(-t - (1 + s) t - z Erf[Sqrt[t]]) z)/Sqrt[\[Pi]]


Integrating this expression gives zero since

Limit[E^(-t (1 + s) - z Erf[Sqrt[t]]) Sqrt[t], t -> 0]

0

Limit[E^(-t (1 + s) - z Erf[Sqrt[t]]) Sqrt[t], t -> \[Infinity],
Assumptions -> s > 0]

0

0 == Integrate[E^(-(1 + s) t - z Erf[Sqrt[t]])/(
2 Sqrt[t]), {t, 0, \[Infinity]}] - (1 + s) Integrate[
E^(-(1 + s) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}] -
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 2) t - z Erf[Sqrt[t]]) , {t, 0, \[Infinity]}]


or, noticing that the last integral is h(z,s+2), i.e. h(z,s) with s shifted by two units,

(* 3 *) Integrate[E^(-(1 + s) t - z Erf[Sqrt[t]])/(
2 Sqrt[t]), {t,
0, \[Infinity]}] == (1 + s) Integrate[
E^(-(1 + s) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}] +
z/Sqrt[\[Pi]] h[z, s + 2]


Now the left hand side of (* 3 ) can be replaced using ( 1 *)

(* 1 rep. *) 1 ==
s h[z, s] +
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]])/ Sqrt[t], {t, 0, \[Infinity]}]


or

Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]])/(
2 Sqrt[t]), {t, 0, \[Infinity]}] == Sqrt[\[Pi]]/(2 z) (1 - s h[z, s])


giving

Sqrt[\[Pi]]/(
2 z) (1 - s h[z, s]) == (1 + s) Integrate[
E^(-(1 + s) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}] +
z/Sqrt[\[Pi]] h[z, s + 2]


Now on the r.h.s. of (* 3 *) we use (*2 *)

(* 2 rappel *) 0 =
D[s h[z, s], s] -
z/Sqrt[\[Pi]]
Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}]


or

Integrate[E^(-(s + 1) t - z Erf[Sqrt[t]]) Sqrt[t], {t, 0, \[Infinity]}] ==
Sqrt[\[Pi]]/z D[s h[z, s], s]


to get from (* 3 *)

(* 3a *)  Sqrt[\[Pi]]/(
2 z) (1 - s h[z, s]) == (1 + s) Sqrt[\[Pi]]/z D[s h[z, s], s] +
z/Sqrt[\[Pi]] h[z, s + 2]


simplifying and sorting

1/2 (1 - s h[z, s]) == (1 + s) D[s h[z, s], s] + z^2/\[Pi] h[z, s + 2]

1/2 (1 - s h[z, s]) == (1 + s) (h[z, s] + s D[h[z, s], s]) +
z^2/\[Pi] h[z, s + 2]

1/2 (1 - s h[z, s]) == (1 + s) h[z, s] + (1 + s) s D[h[z, s], s] +
z^2/\[Pi] h[z, s + 2]

1/2 - s /2 h[z, s] - (1 + s) h[z, s] -
z^2/\[Pi] h[z, s + 2] == (1 + s) s D[h[z, s], s]


gives finally the equation

(* 4 *) ddeq = (1 + s) s D[h[z, s], s] + (1 + (3 s )/2) h[z, s] -
z^2/\[Pi] h[z, s + 2] == 1/2;


End of derivation of (8.2).

– ciao
Feb 3, 2017 at 8:56
• @ciao Thank you for your kind comment. It encourages me to continue and share my fun. Feb 3, 2017 at 14:25
• $\mathcal{L}_t\left[\text{erf}\left(\sqrt{t}\right)^4\right](s)=\frac{4 \sqrt{\frac{1}{1+s}}}{s}-\frac{96 \int_0^{\sqrt{s+1}} \frac{\cot ^{-1}\left(\sqrt{2+a^2}\right)}{\left(1+a^2\right) \sqrt{2+a^2}} \, da}{\pi ^2 s \sqrt{1+s}}$ and solve Ahmed's integral mathworld.wolfram.com/AhmedsIntegral.html Feb 5, 2017 at 14:39
• @Mariusz Iwaniuk Very interesting result. Numerically it is correct. How did you derive it? The difference to Ahmed's integral is that here we need the indefinite integral, or at least one with a variable upper Integration limit, or at least with the upper limit equal to $\sqrt 2$, corresponding to $s=1$. Feb 5, 2017 at 19:22
• @Mariusz Iwaniuk After some manipulations I end up with this integral to be calculated: $$\int_0^1 \frac{\tan ^{-1}\left(\frac{\sqrt{s+1}}{\sqrt{y^2+2}}\right)}{\left(y^2+1\right) \sqrt{y^2+2}} \, dy$$ This is a generalization of Ahmed's integral differing from it by the factor $\sqrt {1+s}$ under the $\tan ^-1$ Feb 5, 2017 at 21:23

Let $F(t)=\text{erf}\left(\sqrt{t}\right)^4$.Then $F'(t)=\frac{4 e^{-t} \text{erf}\left(\sqrt{t}\right)^3}{\sqrt{\pi } \sqrt{t}}$ and $F(0)=0$. $$\mathcal{L}_t\left[F'(t)\right](s)=s \left(\mathcal{L}_t[F(t)](s)\right)-F(0)$$ $$\mathcal{L}_t\left[\frac{4 e^{-t} \text{erf}\left(\sqrt{t}\right)^3}{\sqrt{\pi } \sqrt{t}}\right](s)=s \left(\mathcal{L}_t\left[\text{erf}\left(\sqrt{t}\right)^4\right](s)\right)-0$$

$$\mathcal{L}_t\left[\text{erf}\left(\sqrt{t}\right)^4\right](s)=\frac{4 \left(\mathcal{L}_t\left[\frac{\text{erf}\left(\sqrt{t}\right)^3}{\sqrt{t}}\right](1+s)\right)}{\sqrt{\pi } s}$$ Solve Laplace Transform with integral:

$$\int_0^{\infty } \frac{\text{erf}\left(\sqrt{t}\right)^3 \exp (-(s+1) t)}{\sqrt{t}} \, dt$$ Integrating with parts and substituting Sqrt[t]=x

$${\frac {\sqrt {\pi}}{\sqrt {s+1}}}-\int_{0}^{\infty }\!6\,{\frac { {\rm erf} \left(\sqrt {s+1}x\right) \left( {\rm erf} \left(x\right) \right) ^{2}{{\rm e}^{-{x}^{2}}}}{\sqrt {s+1}}}\,{\rm d}x$$

Consider the parametric integral : $$I(a,b) = \int_{0}^{\infty}\!\!\!\text{erf}(ax)\,\text{erf}(bx)\,\text{erf}(bx)\,e^{-x^2}dx \tag{1}$$ and notice that:

$$\frac{\partial }{\partial a}\frac{\partial I(a,b)}{\partial b}=\int_0^{\infty } \frac{8 e^{-x^2-a^2 x^2-b^2 x^2} x^2 \text{erf}(b x)}{\pi } \, dx=\frac{4 b}{\pi ^{3/2} \left(a^2+b^2+1\right) \left(a^2+2 b^2+1\right)}+\frac{4 \tan ^{-1}\left(\frac{b}{\sqrt{a^2+b^2+1}}\right)}{\pi ^{3/2} \left(a^2+b^2+1\right)^{3/2}}$$ integrating over $b\in \{0,1\}$ and $a\in \{0,\sqrt{s+1}\}$

 Integrate[(4 b)/((1 + a^2 + b^2) (1 + a^2 + 2 b^2) \[Pi]^(3/2)) + (
4 ArcTan[b/Sqrt[1 + a^2 + b^2]])/((1 + a^2 + b^2)^(3/2) \[Pi]^(
/2)), {b, 0, 1}, Assumptions -> a > 0]
(*(4 ArcTan[1/Sqrt[2 + a^2]])/((1 + a^2) Sqrt[2 + a^2] \[Pi]^(3/2))*)

Integrate[(4 ArcTan[1/Sqrt[2 + a^2]])/((1 + a^2) Sqrt[2 + a^2] \[Pi]^(3/2)), {a,
0, Sqrt[s + 1]}]
(*----*)


$$\int_0^{\sqrt{s+1}} \frac{4 \tan ^{-1}\left(\frac{1}{\sqrt{2+a^2}}\right)}{\left(1+a^2\right) \sqrt{2+a^2} \pi ^{3/2}} \, da$$ Connect everything together:

$$\mathcal{L}_t\left[\text{erf}\left(\sqrt{t}\right)^4\right](s)=\frac{4 \sqrt{\frac{1}{s+1}}}{s}-96\frac{\int_0^{\sqrt{s+1}} \frac{\tan ^{-1}\left(\frac{1}{\sqrt{2+a^2}}\right)}{\left(1+a^2\right) \sqrt{2+a^2}} \, da}{\pi ^2 s \sqrt{s+1}}$$

.EDIT 15.02.2017.

Finding Integral:

$\int_0^{\sqrt{s+1}} \frac{\tan ^{-1}\left(\frac{1}{\sqrt{2+a^2}}\right)}{\left(1+a^2\right) \sqrt{2+a^2}} \, da=\frac{1}{4} \pi \tan ^{-1}\left(\sqrt{s+1}\right)-\int_0^1 \frac{\tan ^{-1}\left(\frac{\sqrt{s+1}}{\sqrt{a^2+2}}\right)}{\sqrt{a^2+2} \left(a^2+1\right)} \, da$

$\int_0^1 \frac{\tan ^{-1}\left(\frac{\sqrt{s+1}}{\sqrt{a^2+2}}\right)}{\sqrt{a^2+2} \left(a^2+1\right)} \, da = \frac{1}{6} \left(3 \sqrt{s+1} X(s)+\pi \tan ^{-1}\left(\frac{\sqrt{s+1}}{\sqrt{3}}\right)\right)$

where X[s]is:

$X(s)=\int_{\sin ^{-1}\left(\frac{\sqrt{3} \sqrt{2}}{3}\right)}^{\frac{\pi }{2}} \frac{i \left(\tan ^{-1}\left(\frac{2 \cos (x)}{\sin ^2(x)}\right) \cos (x)\right)}{2 \left(i \sqrt{2}+\sqrt{s+1} \sin (x)\right)} \, dx+\int_{\sin ^{-1}\left(\frac{\sqrt{3} \sqrt{2}}{3}\right)}^{\frac{\pi }{2}} \frac{i \left(\tan ^{-1}\left(\frac{2 \cos (x)}{\sin ^2(x)}\right) \cos (x)\right)}{2 \left(i \sqrt{2}-\sqrt{s+1} \sin (x)\right)} \, dx$

Check numerics:

NIntegrate[ArcTan[Sqrt[s + 1]/Sqrt[x^2 + 2]]/(Sqrt[x^2 + 2]*(x^2 + 1)) /.
s -> 1, {x, 0, 1}]
(* 0.396282 *)

X[s_?NumericQ] :=
NIntegrate[(I/2)*(ArcTan[2 *Cos[x]/Sin[x]^2] Cos[x])/(
I Sqrt[2] + Sqrt[s + 1]* Sin[x]), {x, ArcSin[Sqrt[3]*Sqrt[2]/3],
Pi/2}] +
NIntegrate[(I/2)*(ArcTan[2 *Cos[x]/Sin[x]^2] Cos[x])/(
I Sqrt[2] - Sqrt[s + 1]* Sin[x]), {x, ArcSin[Sqrt[3]*Sqrt[2]/3],
Pi/2}]
1/6 (\[Pi] ArcTan[Sqrt[1 + s]/Sqrt[3]] + 3 Sqrt[1 + s] X[s]) /.
s -> 1 // N
(* 0.396282 *)
Ok works.


Finding a integral :

$\int_{\sin ^{-1}\left(\frac{\sqrt{3} \sqrt{2}}{3}\right)}^{\frac{\pi }{2}} \frac{i \left(\tan ^{-1}\left(\frac{2 \cos (x)}{\sin ^2(x)}\right) \cos (x)\right)}{2 \left(i \sqrt{2}-\sqrt{s+1} \sin (x)\right)} \, dx=?$

MMA has problem to caluculate this integral with limits,but without limits is: where a is Sqrt[s+1]

 Integrate[(I/2)*(ArcTan[2 *Cos[x]/Sin[x]^2] Cos[x])/(
I Sqrt[2] -a*Sin[x]), x]
(*the result is a very long*)

• Call the integral B(s). Differentiating with respect to s I found this formula. Don't know if it is useful $$B(s) = \int_0^1 \frac{\arctan\left(\sqrt{\frac{s+1}{y^2+2}}\right)}{\left(y^2+1\right) \sqrt{y^2+2}} \, dy=-\int_0^s \frac{\cot^{-1}\left(\sqrt{z+3}\right)}{2 \sqrt{z+1} (z+2) \sqrt{z+3}} \, dz+\frac{ \pi}{4} \arctan\left(\sqrt{s+1}\right)-\frac{\pi ^2}{32}$$ At $s = 0$ we have $$B(s=0) = \int_0^1 \frac{\cot ^{-1}\left(\sqrt{y^2+2}\right)}{\left(y^2+1\right) \sqrt{y^2+2}} \, dy=\frac{\pi ^2}{32}$$ Feb 6, 2017 at 22:29
• It seems that the formula in your update is incorrect. It looks as if it was meant as a consequence of the relations arctan(x) + arccot(x) = pi/2 and arccot(x) = arctan(1/x) but the argument in the arctan on the right hand side of your formula is wrong. Feb 7, 2017 at 14:36
• @Dr.WolfgangHintze. Yes you are right.Thanks :) Feb 7, 2017 at 14:49