# How to calculate this integral?

I am trying to Integrate the following Integral :

$\int_1^{\infty } \dfrac{\left(x^2-1\right)^{13/2} e^{-ax} }{x^{10}} \, dx \,\, \,\,\,\,\,\,\,\,\,\,\,\,(a=\textrm{real>0})$

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

In Mathematica input form:

Integrate[Exp[-a x] (x^2 - 1)^(13/2)/x^10, {x,1,Infinity}]


If I simply enter that into Mathematica, it instantly returns the same expression. How do I go about this? Is there any tricks that can be applied? Thank's.

• I'm also unable to get MMA to calculate the integral symbolically. However, you can calculate it numerically for a given value of a: NIntegrate[(Exp[-a x] (x^2 - 1)^(13/2)/x^10) /. {a -> 0.1}, {x, 1, [Infinity]}] – Cassini Jun 4 '14 at 16:27
• check this one mathematica.stackexchange.com/questions/48766/… you probably have to replace Infinity with a large number to use these schemes. – Sumit Jun 4 '14 at 16:51
• You should define contours appropriately, -1 and 1 are the branch points, thus approaching 1 you should know how to tackle this integral. As a reference examine this answer How to calculate contour integrals with Mathematica? – Artes Jun 4 '14 at 18:13
• is there a way to get the analytical result? – Betatron Jun 5 '14 at 5:47

Here's a way to get an exact symbolic expression for the integral.

Use the notation fi[m,n,a] to denote Integrate[(Exp[-a x](x^2 - 1)^(m/2))/x^n, {x, 1, Infinity}].

Define integration by parts, where the boundary term vanishes.

fi[m_, n_, a_] /; m > 0 && n > 1 :=
(m fi[-2 + m, -2 + n, a] - a fi[m, -1 + n, a])/(-1 + n);


See how we're doing so far.

fi[13, 10, a] // Expand

(*
143/3 fi[3, 0, a] + 9295/384 a^2 fi[5, 0, a] -
3003/128 a fi[5, 1, a] + (3289 a^4 fi[7, 0, a])/2688 -
143/32 a^3 fi[7, 1, a] + (143 a^6 fi[9, 0, a])/10368 -
143/960 a^5 fi[9, 1, a] + (13 a^8 fi[11, 0, a])/362880 - (
13 a^7 fi[11, 1, a])/10080 - (a^9 fi[13, 1, a])/362880
*)


Define how the fi[m_, 0, a_] cases evaluate.

fi[m_, 0, a_] /; m > 0 =
Assuming[Re[m] > 0 && Re[a] > 0, Integrate[Exp[-a x] (x^2 - 1)^(m/2), {x, 1, ∞}]]

(* (2^((1 + m)/2) a^(1/2 (-1 - m)) BesselK[(1 + m)/2, a] Gamma[1 + m/2])/Sqrt[π] *)


See how we are doing so far.

fi[13, 10, a] // Expand

(*
(143 BesselK[2, a])/a^2 + (46475 BesselK[3, a])/(128 a) +
16445/128 BesselK[4, a] + 5005/384 a BesselK[5, a] +
143/384 a^2 BesselK[6, a] - 3003/128 a fi[5, 1, a] -
143/32 a^3 fi[7, 1, a] - 143/960 a^5 fi[9, 1, a] - (
13 a^7 fi[11, 1, a])/10080 - (a^9 fi[13, 1, a])/362880
*)


Define how the fi[m_, 1, a_] cases evaluate, using differentiation w.r.t. a under the integral to get rid of the 1/x factor, then integrating w.r.t. a afterwards.

fi[m_, 1, a_] /; m > 0 =
Assuming[Re[m] > 0 && Re[a] > 0,
Integrate[-Exp[-a x] (x^2 - 1)^(m/2), {x, 1, ∞}] //
Integrate[#, a] & //
Simplify //
(# - (Limit[#, a -> ∞]//FullSimplify)&)]

(*
-(1/2) π Csc[(m π)/2] +
1/4 a^-m Sqrt[π]
Gamma[1 + m/2] (a^(1 + m) Sqrt[π]
HypergeometricPFQRegularized[{1/2}, {(3 + m)/2, 3/2}, a^2/4]
- 2^(1 + m) Gamma[-(m/2)] HypergeometricPFQRegularized[{-(m/2)}, {1/2 - m/2, 1 - m/2}, a^2/4]) Sec[(m π)/2]
*)


The last step fixes the constant of integration to ensure that the result goes to 0 (as it should) as a goes to infinity.

However, if m is an odd integer then the Sec[(m π)/2] factor blows up, and the factor containing the difference between the two HypergeometricPFQRegularized terms goes to zero, but these factors have a well-behaved product. So we need to manually intervene to make the expression manifestly finite when m is an odd integer.

It is sufficient to do a series expansion about odd integer values of m to extract the relevant coefficients. Here is how I did it.

factor1 = Assuming[m0 ∈ Integers && m0 >= 0,
SeriesCoefficient[
a^(1 + m) Sqrt[π]
HypergeometricPFQRegularized[{1/2}, {(3 + m)/2, 3/2}, a^2/4] -
2^(1 + m) Gamma[-(m/2)] HypergeometricPFQRegularized[{-(m/2)}, {1/2 - m/2,
1 - m/2}, a^2/4], {m, 2 m0 + 1, 1}] // FullSimplify];

factor2 = Assuming[m0 ∈ Integers && m0 >= 0,
SeriesCoefficient[Sec[(m π)/2], {m, 2 m0 + 1, -1}]];


Gather the pieces together, and rewrite the definition for the fi[m_, 1, a_] cases.

fi[m_, 1, a_] /; m > 0 =
-(1/2) π Csc[(m π)/2] + 1/4 a^-m Sqrt[π] Gamma[1 + m/2] (factor1 factor2 /. m0 -> (m - 1)/2)

(*
(-(1/2))*Pi*Csc[(m*Pi)/2] - (Gamma[1 + m/2]*Sec[(1/2)*(-1 + m)*Pi]*
(a^(1 + m)*(2*
HypergeometricPFQ[{1/2}, {3/2, 3/2 + m/2},
a^2/4]*(-Log + 2*Log[a] +
PolyGamma[0, -(1/2) + (1 - m)/2]) +
Sqrt[Pi]*Gamma[2 + (1/2)*(-1 + m)]*
Derivative[{0}, {1, 0}, 0][
HypergeometricPFQRegularized][{1/2}, {2 + (1/2)*(-1 + m),
3/2}, a^2/4]) +
4^(1 + (1/2)*(-1 + m))*Gamma[-(1/2) + (1 - m)/2]*
Gamma[2 + (1/2)*(-1 + m)]*
(Derivative[{0}, {0, 1}, 0][
HypergeometricPFQRegularized][{-(1/2) + (1 - m)/
2}, {(1 - m)/2, 1/2 + (1 - m)/2},
a^2/4] +
Derivative[{0}, {1, 0}, 0][
HypergeometricPFQRegularized][{-(1/2) + (1 - m)/2},
{(1 - m)/2, 1/2 + (1 - m)/2}, a^2/4] +
Derivative[{1}, {0, 0}, 0][HypergeometricPFQRegularized][
{-(1/2) + (1 - m)/2}, {(1 - m)/2, 1/2 + (1 - m)/2},
a^2/4])))/(a^m*(4*Sqrt[Pi]*Gamma[2 + (1/2)*(-1 + m)]))
*)


Numerically verify that this gives you the right answer.

With[{m = 13, n = 10, a = RandomReal[{10^-6, 3}]},
{NIntegrate[(Exp[-a x] (x^2 - 1)^(m/2))/x^n, {x, 1, ∞}], fi[m, n, a]}]


The required integral is then

fi[13, 10, a] // Expand

(* Lots of HypergeometricPFQ and HypergeometricPFQRegularized *)


The symbolic result looks incredibly messy to me, so maybe there are some further simplifications that could compress it down somewhat ...

• Cool. +1. Your polynomial part agrees with my polynomial part. I guess the rest of it simplifies to the MeijerG part of my answer. – Michael E2 Jun 6 '14 at 20:15
• Indeed, my messy BesselK + HypergeometricPFQ + HypergeometricPFQRegularized expression is numerically the same as your nice MeijerG term. So, combining my power series piece with your MeijerG piece solves the problem nicely. I also noticed that you can directly obtain my symbolic power series from your numerical power series by evaluating Pi Rationalize[Expand[<your power series>/N[Pi]],10^-10], but that needs a bit of hindsight to pull off with confidence! – Stephen Luttrell Jun 6 '14 at 22:00
• @ Stephen Luttrell, I am very grateful for your help. Thank you again for everything you’ve done. – Betatron Jun 8 '14 at 20:43

Here is a way to get an approximate symbolic expression for the integral. Some of the coefficients are approximate because at some points in the process we need the value of a definite integral at a = 1, and unfortunately those integrals have to be computed numerically.

First, a lemma: Let $f(a) = \int_1^\infty g(x, a) \; dx$. Then $f(a) = f(1) + \int_1^a \int_1^\infty {\partial g \over \partial a}(x, a) \; dx \; da$ under suitable conditions for $g$.

The gist is that we can differentiate with respect to a first, then integrate with respect to x, and finally integrate with respect to a. This can be helpful since ${\partial g / \partial a}$ has a different form than $g$; it might turn into something that can be integrated. In the OP's case, if we have to differentiate ten times (which gets eliminates x from the denominator), we get to an integrate Mathematica will solve.

i1 = Simplify[
Integrate[D[Exp[-a x] (x^2 - 1)^(13/2)/x^10, {a, 10}], {x, 1, Infinity}],
a > 0]
(* 135135 BesselK[7, a] / a^7 *)


Next we'll calculate the constants (the $f(1)$ values) for each antiderivative. These are the integrals that do not evaluate symbolically.

ClearAll[ni];
ni[n_Integer, a0_?NumericQ, opts___] :=
NIntegrate[D[Exp[-a x] (x^2 - 1)^(13/2)/x^10, {a, n}] /. a -> a0, {x, 1, Infinity}, opts]

nii = Table[ni[k, 1, WorkingPrecision -> 50], {k, 0, 9}];


Next we'll integrate ten times. There are numerical issues, so we converted the numbers above to infinite precision representations.

i2 = Block[{k = 0, n = SetPrecision[Reverse@nii, Infinity]},
Nest[
With[{int = Integrate[#, a]}, n[[++k]] + int - (int /. a -> 1)] &,
i1, 10]];


This is gives the answer but there are many exact MeijerG expressions from the repeated integrations. For practical use, it would be faster to convert them once to approximate numbers. Here is a way to get 50-digit reals:

integral50 = Block[{nh, $MaxExtraPrecision = 100}, SetAttributes[nh, NHoldAll]; N[i2 /. g_MeijerG /; Not[FreeQ[g, a]] :> nh[g], 50] /. nh -> Identity ]  The even-power coefficients are in the range$10^{-48}$to$10^{-45}$, which suggests they might be zero. If we recompute nii and then the integral with 200 digits of precision, the coefficients drop to the range$10^{-198}$to$10^{-196}$which strongly suggests that they are zero. nii = Table[ni[k, 1, WorkingPrecision -> 200, MaxRecursion -> 20], {k, 0, 9}]; i2 = Block[{k = 0, n = SetPrecision[Reverse@nii, Infinity]}, Nest[With[{int = Integrate[#, a]}, n[[++k]] + int - (int /. a -> 1)] &, i1, 10] ]; integral200 = Block[{nh,$MaxExtraPrecision = 500},
SetAttributes[nh, NHoldAll];
N[i2 /. g_MeijerG /; Not[FreeQ[g, a]] :> nh[g], 200] /. nh -> Identity
];

integral = Rationalize[integral200] // Chop
(*
36.85235444816464... a   -  7.01949608536469... a^3 +
0.23398320284548... a^5 -  0.00202582859606... a^7 +
4.32869358133514...*10^-6 a^9 -
135135/256 MeijerG[{{11/2}, {}}, {{-2, 0}, {1/2}}, a/2, 1/2]
*)


(Note: I replaced the coefficient of MeijerG by means of Rationalize, since it is in fact exact. The coefficients of the polynomial part are approximate. The relatively small coefficient of a^9 did not change with the change in precision.)

Plot:

Plot[integral, {a, 1, 10}, PlotRange -> All, AxesOrigin -> {0, 0}] • I am very grateful for your help. Thank you again for everything you’ve done. – Betatron Jun 8 '14 at 20:44
• @betatron You're welcome. – Michael E2 Jun 8 '14 at 22:00

You can have an analytical result, though not exact. I have several times addressed the analogous questions. You may look here and here where the approach is described within another problem depending upon a single parameter. I will consider your particular case fixing the domain of variation of the parameter a between 0.5 and 2. It is just for shortness. One may choose whatever he wants. Let us make a table like the following:

J = Table[{a, NIntegrate[Exp[-a x] (x^2 - 1)^(13/2)/x^10, {x, 1, Infinity}]}, {a, 0.5, 2,0.1}];


Now we have a list of points with the structure {a, J(a)} and we can fit it with any reasonable function within this interval:

    ff = FindFit[int, k1*Exp[-k2*a] + k3, {k1, k2, k3}, a]
Show[{
ListPlot[int, PlotRange -> All,
AxesLabel -> {Style["a", Italic, 16], Style["J", Italic, 16]}],
Plot[k1*Exp[-k2*a] + k3 /. ff, {a, 0.5, 2}, PlotStyle -> Red,
PlotRange -> All]
}]


This yields the parameters of the fitting function:

(*   {k1 -> 3504.27, k2 -> 7.62848, k3 -> 0.650277}  *)


and the plot to check the fitting quality visually. It should look like the following: where the blue points show the values of the integral, and the red line shows the fit.

If you need to get a better quality you may think about the modification of the function. For example, like this:

    Manipulate[
ff = FindFit[int, (k1*Exp[-k2*a] + k3)/a^\[Alpha], {k1, k2, k3}, a];
im2 = Show[{
ListPlot[int, PlotStyle -> {Blue, PointSize[0.01]},
PlotRange -> All,
AxesLabel -> {Style["a", Italic, 16], Style["J", Italic, 16]}],
Plot[(k1*Exp[-k2*a] + k3)/a^\[Alpha] /. ff, {a, 0.5, 2},
PlotStyle -> Red, PlotRange -> All]
}], {\[Alpha], 0.1, 3}]


You should see the following: on the screen. The fit is somewhat better. Of course, the larger is the fitting interval, the more challenging is the task of guessing a good fitting function.

As I said in my comment, you can calculate the integral for a specific value of a. Here's a plot of the integral as a function of a:

 Plot[NIntegrate[(Exp[-a x] (x^2 - 1)^(13/2)/x^10), {x, 1, Infinity}], {a, 1, 10}, PlotRange -> All, AxesOrigin -> {0, 0}, AxesLabel -> {"a", "Integral"}] 