How to compute this integral?

I'm trying to calculate this integral :

$$I(z,k,a)= \displaystyle\int_{1}^{\infty} t^2\, \text{ArcTanh} \left(\sqrt{\frac{t^2-1}{t^2}} \dfrac{k}{z}\right)\, e^{-a\,t} \, dt$$

Where :

• $\textrm{ArcTanh(...)}$ : is the inverse hyperbolic tangent function.

• $z=x+i\,y$ : is a complex with the assymption $y<<<x$.

• $k$ and $a$ are reals $>0$.

Mathematica didn't calculate it. Maybe it is too complicated to be done.

In Mathematica input form :

Assuming[{a>0&&k>0},Integrate[t^2 (ArcTanh[Sqrt[(t^2-1)/t^2](k/z)])Exp[-a t],{t,1,Infinity}]]


If I simply enter that into Mathematica, it instantly returns the same expression.

Please, how do I go about this ? Is there any tricks that can be applied ?

Is there a way to get the symbolic result ?

Thank's.

• i like to start by considering a particular case, say a=k=z=1. In this case even that does not yield a closed form result, so I'd call to prospect of obtaining a result for the general case doubtful. May 9 '16 at 19:26
• Is a symbolic solution absolutely necessary, or could you get away with calculating values of this integral numerically after assigning specific values to the parameters? May 9 '16 at 20:40
• @MarcoB Numerical solutions are already calculated, but I would like to calculate the symbolic solution for solving a non linear equation and discuss the analytical solutions. May 9 '16 at 20:53
• @Betatron Then this may not be a Mathematica problem, but a mathematical one, i.e. the existence of a closed-form antiderivative is doubtful. Notice also that, for what it's worth, your Assuming[...] expression evaluates for close to 30 s then returns the expression unevaluated in MMA 10.4.0. May 9 '16 at 22:32
• If you are sure that analytic solution exists please provide some background and preceding steps of your derivation. May 10 '16 at 7:45

2 Answers

Preliminary post

For k/z = 1 there is a closed form result:

I[1,1,a]=((4 + a^2) BesselK[0, a] + a (4 BesselK[1, a] + a BesselK[2, a]))/(2 a^3)


The derivation and the extension to k/z != 1 requires some manual interaction to help Mathematica which we will show in the following.

Solution

Summary

As the integral to be calculated is returned unevaluated by Mathematica we show here to get useful partial results by using a method I like to call "man-machine-interaction". This means we carry out a joint work with paper and pencil on one side and Mathematica on the other side.

The case m = k/z == 1 can be solved symbolically. If m < 1 we show how to get the symbolic soluton to aribtrary precision as a power series in (1-m).

Preliminary remarks

1) The integral to be calculated is

fI := Integrate[
t^2 (ArcTanh[m Sqrt[(t^2 - 1)/t^2]]) Exp[-a t], {t, 1, Infinity},
Assumptions -> {a > 0, 0 < m < 1}]


The numerical integral presents, however, no problem

fN[a_, m_] :=
NIntegrate[t^2 (ArcTanh[m Sqrt[(t^2 - 1)/t^2]]) Exp[-a t], {t, 1, Infinity}]


2) Here we have set m = k/z, because k and z enter only as a quotient.

3) Also we have introduced natural restrictions on the parameters "a" and "m". The restriction for m comes from the requirement that fI be a real quantity. This in turn requires the argument of the ArcTanh to be between 0 and 1 if we keep the upper limit of integration of t at infinity.

4) The factor t^2 in the integrand can be generated by differentiating twice with respect to "a". Hence it is sufficient to study the simpler intergral without the factor t^2.

Hence the integrand to be studied becomes

fi = Exp[-a t] ArcTanh[m Sqrt[(-1 + t^2)/t^2]];


Step by step solution

The first step will be a partial integration. This step is not done automatically by Mathematica, hence we do it manually.

fi can be written as

fi = fid + fi1


where

fid1 = -(1/a) Exp[-a t] ArcTanh[m Sqrt[(-1 + t^2)/t^2]];

fid = D[fid1, t];


and

fi1 = Simplify[ (E^(-a t) m (2/t - (2 (-1 + t^2))/t^3))/(
2 a Sqrt[(-1 + t^2)/t^2] (1 - (m^2 (-1 + t^2))/t^2)), {a > 0, t > 0,
0 < m < 1}]

(*
Out[99]= (E^(-a t) m)/(a Sqrt[-1 + t^2] (t^2 - m^2 (-1 + t^2)))
*)


Indeed

Simplify[fi == fid + fi1, {a > 0, t > 0, 0 < m < 1}]

(* Out[101]= True *)


Because

fid1 /. t -> 1

(* Out[102]= 0 *)

Limit[fid1, t -> \[Infinity], Assumptions -> a > 0]

(* Out[103]= 0 *)


the integral becomes

fI2 := Integrate[(E^(-a t) m)/(
a Sqrt[t^2 - 1] (t^2 + m^2 (1 - t^2))), {t, 1, \[Infinity]},
Assumptions -> {a > 0, 0 < m < 1}]


This integral is also returned unevaluated by Mathematica. Hence we look first at a special case.

The case m = 1

If m = 1 the integrand simplifies and the integral can be done expicitly

(E^(-a t) m)/(a Sqrt[t^2 - 1] (t^2 + m^2 (1 - t^2))) /. m -> 1

(* Out[125]= E^(-a t)/(a Sqrt[-1 + t^2]) *)

fI3 = Integrate[%, {t, 1, \[Infinity]}, Assumptions -> a > 0]

(* Out[126]= BesselK[0, a]/a *)


In order to find the original integral we need to differentiate this expression twice

fI4 = D[fI3, {a, 2}] // Simplify

(*
Out[127]= ((4 + a^2) BesselK[0, a] + a (4 BesselK[1, a] + a BesselK[2, a]))/(2 a^3)
*)


The result in the original designation

$$I(1,1,a)=\frac{\left(a^2+4\right) K_0(a)+a (4 K_1(a)+a K_2(a))}{2 a^3}$$

The numerical agreement is excellent:

Plot[fN[a, 1]/fI4, {a, 0.1, 2}, PlotRange -> {{0, 2}, {0.99, 1.01}},
PlotLabel ->
"Integral for m = 1\nComparison of symbolical and numerical solution",
AxesLabel -> {"a", "fN/fI4"}]


Series expansion

Now we can proceed further developing the integrand of fi1 into a power series about m = 1.

Remark: Although it leads to cumbersome expressions I took four terms in order to possibly find a rule in the result (but, alas, I didn't)

Series[fi1, {m, 1, 4}] // Normal

(*
Out[122]= E^(-a t)/(a Sqrt[-1 + t^2]) + (E^(-a t) (-1 + m) (-1 + 2 t^2))/(
a Sqrt[-1 + t^2]) + (E^(-a t) (-1 + m)^2 (1 - 5 t^2 + 4 t^4))/(
a Sqrt[-1 + t^2]) + (
E^(-a t) (-1 + m)^4 (-1 + t^2)^(3/2) (1 - 12 t^2 + 16 t^4))/a + (
E^(-a t) (-1 + m)^3 (-1 + 9 t^2 - 16 t^4 + 8 t^6))/(a Sqrt[-1 + t^2])
*)

fis = List @@ %;

Timing[Integrate[fis, {t, 1, \[Infinity]}, Assumptions -> {a > 0}]]

(*
{1233.4843069, {BesselK[0, a]/a, ((-1 + m) BesselK[2, a])/a, (
3 (-1 + m)^2 (4 a BesselK[0, a] + (8 + a^2) BesselK[1, a]))/a^4, (
15 (-1 + m)^4 (a (112 + a^2) BesselK[2, a] +
4 (168 + 5 a^2) BesselK[3, a]))/
a^6, ((-1 +
m)^3 (24 a (20 + a^2) BesselK[0, a] + (960 + 168 a^2 + a^4) BesselK[1,
a]))/a^6}}
*)


Finally, we need the second derivative with respect to a:

D[%[[2]], {a, 2}] // Simplify

(*
{((4 + a^2) BesselK[0, a] + a (4 BesselK[1, a] + a BesselK[2, a]))/(
2 a^3), ((-1 + m) (a^2 BesselK[0, a] + 4 a BesselK[1, a] +
8 BesselK[2, a] + 2 a^2 BesselK[2, a] + 4 a BesselK[3, a] +
a^2 BesselK[4, a]))/(4 a^3), (
3 (-1 + m)^2 (16 a (20 + a^2) BesselK[0,
a] + (640 + 144 a^2 + 3 a^4) BesselK[1, a] +
a (8 + a^2) (16 BesselK[2, a] + a BesselK[3, a])))/(4 a^6), (1/(
4 a^8))15 (-1 + m)^4 (a^3 (112 + a^2) BesselK[0, a] +
32 a^2 (91 + a^2) BesselK[1, a] + 29568 a BesselK[2, a] +
592 a^3 BesselK[2, a] + 2 a^5 BesselK[2, a] +
112896 BesselK[3, a] + 5184 a^2 BesselK[3, a] +
52 a^4 BesselK[3, a] + 16128 a BesselK[4, a] +
432 a^3 BesselK[4, a] + a^5 BesselK[4, a] +
672 a^2 BesselK[5, a] + 20 a^4 BesselK[5, a]), (1/(
4 a^8))(-1 + m)^3 (8 a (10080 + 600 a^2 + 7 a^4) BesselK[0, a] +
3 (53760 + 11840 a^2 + 368 a^4 + a^6) BesselK[1, a] +
a (8 (2880 + 456 a^2 + 7 a^4) BesselK[2, a] +
a (960 + 168 a^2 + a^4) BesselK[3, a]))}
*)


Observation

Related simpler example where Mathematica needs some help

This integral is returned unevaluated

Integrate[Exp[-a Cosh[u]], {u, 0, \[Infinity]}, Assumptions -> a > 0]

(*
Out[124]= Integrate[E^(-a Cosh[u]), {u, 0, \[Infinity]}, Assumptions -> a > 0] *)


But with the substitution Cosh[u] -> t we obtain an eqivalent expression of the integral which now is deone by Mathematica:

Integrate[Exp[-a t]/Sqrt[t^2 - 1], {t, 1, \[Infinity]}, Assumptions -> a > 0]

(* Out[142]= BesselK[0, a] *)

• "[...] I'll show it later." -- Looking forward to see it! :) May 15 '16 at 17:58
• @Dr. Wolfgang Hintze, If I simply enter into Mathematica I[1,1,a] it instantly returns the same expression. Please, How did you get the closed form I[1,1,a]=((4 + a^2) BesselK[0, a] + a (4 BesselK[1, a] + a BesselK[2, a]))/(2 a^3) ? Thanks May 15 '16 at 20:30

For k/z = 1, and integrating by parts:

ClearAll["Global*"]
Iv[t_] := Exp[-a*t]*t^2;
u[t_] := ArcTanh[Sqrt[(t^2 - 1)/t^2]];
v = Integrate[Iv[t], t];
Du = Simplify@D[u[t], t];
Int == u[t]*v - Integrate[Du*v, t]

HoldForm[Integrate[Exp[-a*t]*t^2*ArcTanh[Sqrt[(t^2 - 1)/t^2]], {t, 1, Infinity}] ==
Limit[-((E^(-a t) (2 + 2 a t + a^2 t^2) ArcTanh[Sqrt[(-1 + t^2)/
t^2]])/a^3), t -> Infinity] -
Limit[-((
E^(-a t) (2 + 2 a t + a^2 t^2) ArcTanh[Sqrt[(-1 + t^2)/t^2]])/
a^3), t -> 1] +
Integrate[(E^(-a t) (2 + 2 a t + a^2 t^2))/(
Sqrt[1 - 1/t^2] t), {t, 1, Infinity}]/a^3]


$\int_1^{\infty } \exp (-a t) t^2 \tanh ^{-1}\left(\sqrt{\frac{t^2-1}{t^2}}\right) \, dt=\lim_{t\to \infty } \, -\frac{e^{-a t} \left(2+2 a t+a^2 t^2\right) \tanh ^{-1}\left(\sqrt{\frac{-1+t^2}{t^2}}\right)}{a^3}-\lim_{t\to 1} \, -\frac{e^{-a t} \left(2+2 a t+a^2 t^2\right) \tanh ^{-1}\left(\sqrt{\frac{-1+t^2}{t^2}}\right)}{a^3}+\frac{\int_1^{\infty } \frac{e^{-a t} \left(2+2 a t+a^2 t^2\right)}{\sqrt{1-\frac{1}{t^2}} t} \, dt}{a^3}$

Limit[u[t]*v, t -> Infinity, Assumptions -> a > 0] -
Limit[u[t]*v, t -> 1, Assumptions -> a > 0]


$0$

ExpandAll@Simplify[Du*v, t > 0]

in1 = Assuming[{a > 0, t > 0}, Integrate[-((2 E^(-a t))/(a^3 Sqrt[-1 + t^2])), {t, 1, Infinity}]];
in2 = Assuming[{a > 0, t > 0}, Integrate[-((2 E^(-a t) t)/(a^2 Sqrt[-1 + t^2])), {t, 1,
Infinity}]];
in3 = Assuming[{a > 0, t > 0}, Integrate[(-E^(-a t) t^2)/(a Sqrt[-1 + t^2]), {t, 1, Infinity}]];
HoldForm[Integrate[Exp[-a*t]*t^2*ArcTanh[Sqrt[(t^2 - 1)/t^2]], {t, 1,
Infinity}]] == -(in1 + in2 + in3)


$\int_1^{\infty } \exp (-a t) t^2 \tanh ^{-1}\left(\sqrt{\frac{t^2-1}{t^2}}\right) \, dt=\frac{2 K_0(a)}{a^3}+\frac{2 K_1(a)}{a^2}+\frac{a K_0(a)+K_1(a)}{a^2}$

where: $K_0(a)$ and $K_1(a)$ the modified Bessel function of the second kind.