Numerical_integration

Question

I have this integrand

Clear[IntegrandEq10]; IntegrandEq10[x_,x0_] := (x*Sqrt[((x/x0)^2)*(1 - (1/x0)) - (1 - (1/x))])^(-1)

and I need to solve it numerically. I tried to solve it as follows

alphanumerical[x0_] := 2*NIntegrate[IntegrandEq10[x, x0], {x, x0, \[Infinity]}, 
Method -> {"TrapezoidalRule", "RombergQuadrature" -> False}, 
MaxRecursion -> 100] -1*\[Pi]

but I get some errors for certain values. Is there anything I can do to get ride of errors like

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>

Any suggestions is appreciated. Thank you!

Answer 1

There is an integrable singularity of the integrand at x==x0 as

IntegrandEq10[x_, x0_] := (x*Sqrt[((x/x0)^2)*(1 - (1/x0)) - (1 - (1/x))])^(-1);
Normal[Series[(x*Sqrt[((x/x0)^2)*(1 - (1/x0)) - (1 - (1/x))])^(-1), {x, x0, 
1}]] // Simplify

((x - x0) (x (6 - 5 x0) + 3 x0 (-4 + 3 x0)))/(2 x0^4 (((x - x0) (-3 + 2 x0))/x0^2)^(3/2))

demonstrates. The Method -> {"TrapezoidalRule", "RombergQuadrature" -> False} option is doubtely applicable for improper integrals. The following does the job.

alphanumerical[x0_] := 2*NIntegrate[IntegrandEq10[x, x0], {x, x0, \[Infinity]}, 
Method -> "GlobalAdaptive", MaxRecursion -> 100] - 1*\[Pi];
Plot[alphanumerical[x0], {x0, 1, 4}, PlotRange -> All, PlotStyle -> Thick]

If x0 <1, the one deals with complex numbers, e.g.

alphanumerical[0.5]

-3.14159 - 2.57874 I

Answer 2

Some values of x0 may be giving you integrands that you may not have anticipated. For example:

iEq10[x_, x0_] := (x*Sqrt[((x/x0)^2)*(1 - (1/x0)) - (1 - (1/x))])^(-1)

Note the behavior of the integrand for small x0:

With[{expr = ReIm[FullSimplify[iEq10[x, x0]]]},
 Manipulate[
  Plot[expr, {x, x0, Infinity}],
  {x0, .1, 10}
  ]
 ]

When do things become complex for real x and x0?

rootTerm = Denominator[FullSimplify[iEq10[x, x0]]][[2, 1]]

rootPos = Simplify[Reduce[rootTerm > 0, x0, Reals], Assumptions -> x > x0 > 0]

RegionPlot[rootPos && x > x0, {x, 0, 10}, {x0, 0, 10}]

Your numerical integration:

nint[x0_] := NIntegrate[iEq10[x, x0], {x, x0, Infinity}]

behaves well when x0 is within the region where the integrand is real:

nint[4]

but not when the integrand is not always real:

nint[1.5]

Answer 3

Your integral can be expressed as the negative imaginary part of an EllipticF integral like

LogPlot[{NIntegrate[1/(x*Sqrt[-1 + 1/x + 
(x^2*(-1 + x0))/x0^3]), {x, x0, Infinity}], 
 -Im[(2*Sqrt[x0]*EllipticF[
      ArcSin[((-1 + x0)^(1/4)*Sqrt[Sqrt[-1 + x0] + Sqrt[3 + x0]])/
        Sqrt[-3 + x0 + Sqrt[-1 + x0]*Sqrt[3 + x0]]], 
        (-3 + x0 + Sqrt[-1 + x0]*Sqrt[3 + x0])/(2*
        Sqrt[-1 + x0]*Sqrt[3 + x0])])/
     ((-1 + x0)^(1/4)*(3 + x0)^(1/4))]}, {x0, 3/2, 10}, 
PlotStyle -> {Blue, Dashed}, PlotRange -> All]

It evaluates as a real number for x0 > 3/2. The numerical integration doesn't complain...

Answer 4

You can get anylytical solution via indefinite integration and taking limits.

Clear[IntegrandEq10]; 
IntegrandEq10[x_, x0_] = 
 Map[Together, (x*Sqrt[((x/x0)^2)*(1 - (1/x0)) - (1 - (1/x))])^(-1), 
  Infinity]

(*   1/(x Sqrt[(-x^3 + x^3 x0 + x0^3 - x x0^3)/(x x0^3)])   *)

ii[x_, x0_] = 
 Integrate[IntegrandEq10[x, x0], x, Assumptions -> 3/2 < x0 < x]

(*   (1/(Sqrt[-1 + 1/x + (x^2 (-1 + x0))/x0^3] x0))2 Sqrt[2] x Sqrt[(
 x^2 (-1 + x0) + x (-1 + x0) x0 - x0^2)/(
 x^2 (-3 + 2 x0 + x0^2))] Sqrt[-3 + 2 x0 + x0^2] Sqrt[(-x + x0)/(
 x (-3 + x0 + Sqrt[-3 + 2 x0 + x0^2]))]
  EllipticF[
  ArcSin[Sqrt[(-1 + x0 - (2 x0)/x + Sqrt[-3 + 2 x0 + x0^2])/
   Sqrt[-3 + 2 x0 + x0^2]]/Sqrt[2]], (
  2 Sqrt[-3 + 2 x0 + x0^2])/(-3 + x0 + Sqrt[-3 + 2 x0 + x0^2])]   *)

lim2 = Limit[ii[x, x0], x -> Infinity, Assumptions -> x0 > 3/2]

(*   2 I Sqrt[2] Sqrt[x0/(-3 + x0 + Sqrt[-3 + 2 x0 + x0^2])]
  EllipticF[
  ArcSin[Sqrt[-1 + x0 + Sqrt[-3 + 2 x0 + x0^2]]/(
   Sqrt[2] (-3 + 2 x0 + x0^2)^(1/4))], (
  2 Sqrt[-3 + 2 x0 + x0^2])/(-3 + x0 + Sqrt[-3 + 2 x0 + x0^2])]   *)

ComplexExpand has to be used to help Limit. lim2 has a imaginary part, which must be subtracted at the end.

cereii = ComplexExpand[
   Re[ii[x, x0] // FullSimplify[#, Assumptions -> 3/2 < x0 < x] &], 
   TargetFunctions -> {Re, Im}] // 
  FullSimplify[#, Assumptions -> 3/2 < x0 < x] &

(*   -((2 Sqrt[x0]
   Im[EllipticF[
    ArcTan[Sqrt[(x + 2 x0 - x x0 - x Sqrt[-3 + 2 x0 + x0^2])/(
     2 x - 2 x0)]], (-3 + x0 + Sqrt[(-1 + x0) (3 + x0)])/(
    2 Sqrt[(-1 + x0) (3 + x0)])]])/((-1 + x0) (3 + x0))^(1/4))   *)

ceimii = ComplexExpand[
   Im[ii[x, x0] // FullSimplify[#, Assumptions -> 3/2 < x0 < x] &], 
   TargetFunctions -> {Re, Im}] // 
  FullSimplify[#, Assumptions -> 3/2 < x0 < x] &

(*  (2 Sqrt[x0]
  Re[EllipticF[
   ArcTan[Sqrt[(x + 2 x0 - x x0 - x Sqrt[-3 + 2 x0 + x0^2])/(
    2 x - 2 x0)]], (-3 + x0 + Sqrt[(-1 + x0) (3 + x0)])/(
   2 Sqrt[(-1 + x0) (3 + x0)])]])/((-1 + x0) (3 + x0))^(1/4)    *)

lim1re = Limit[cereii, x -> x0, Direction -> -1, 
  Assumptions -> x0 > 3/2]

(*   0   *)

lim1im = Limit[ceimii, x -> x0, Direction -> -1, 
  Assumptions -> x0 > 3/2]

(*   (2 Sqrt[x0]
  EllipticK[(-3 + x0 + Sqrt[-3 + 2 x0 + x0^2])/(
  2 Sqrt[-3 + 2 x0 + x0^2])])/((-1 + x0) (3 + x0))^(1/4)   *)

int[x0_] = 
 lim2 - I lim1im // FullSimplify[#, Assumptions -> 3/2 < x0] & // 
  PowerExpand[#, Assumptions -> 3/2 < x0] &

(*   (2 I Sqrt[x0]
   EllipticF[
   ArcSin[((-1 + x0)/(3 + x0))^(1/4) Sqrt[(
     3 + x0 + Sqrt[(-1 + x0) (3 + x0)])/(-3 + x0 + 
      Sqrt[(-1 + x0) (3 + x0)])]], (-3 + x0 + 
    Sqrt[(-1 + x0) (3 + x0)])/(
   2 Sqrt[(-1 + x0) (3 + x0)])])/((-1 + x0) (3 + x0))^(1/4) - (
 2 I Sqrt[x0]
   EllipticK[(-3 + x0 + Sqrt[(-1 + x0) (3 + x0)])/(
   2 Sqrt[(-1 + x0) (3 + x0)])])/((-1 + x0) (3 + x0))^(1/4)   *)

Plot[int[x0], {x0, 3/2, 20}, PlotRange -> {0, 10}, 
 AxesOrigin -> {0, 0}, GridLines -> {{3/2}, Automatic}, 
 PlotStyle -> Thick]