5
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • $\begingroup$ 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]}] $\endgroup$ Commented Jun 4, 2014 at 16:27
  • $\begingroup$ check this one mathematica.stackexchange.com/questions/48766/… you probably have to replace Infinity with a large number to use these schemes. $\endgroup$ Commented Jun 4, 2014 at 16:51
  • $\begingroup$ 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? $\endgroup$ Commented Jun 4, 2014 at 18:13
  • $\begingroup$ is there a way to get the analytical result? $\endgroup$ Commented Jun 5, 2014 at 5:47

4 Answers 4

Reset to default
8
$\begingroup$

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[4] + 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 ...

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Cool. +1. Your polynomial part agrees with my polynomial part. I guess the rest of it simplifies to the MeijerG part of my answer. $\endgroup$ Commented Jun 6, 2014 at 20:15
  • 1
    $\begingroup$ 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! $\endgroup$ Commented Jun 6, 2014 at 22:00
  • $\begingroup$ @ Stephen Luttrell, I am very grateful for your help. Thank you again for everything you’ve done. $\endgroup$ Commented Jun 8, 2014 at 20:43
7
$\begingroup$

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}]

Mathematica graphics

Improve this answer
$\endgroup$
2
  • $\begingroup$ I am very grateful for your help. Thank you again for everything you’ve done. $\endgroup$ Commented Jun 8, 2014 at 20:44
  • $\begingroup$ @betatron You're welcome. $\endgroup$ Commented Jun 8, 2014 at 22:00
5
$\begingroup$

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:

enter image description here 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:

enter image description here

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.

Improve this answer
$\endgroup$
2
$\begingroup$

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"}]

enter image description here

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.