3
$\begingroup$

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!

$\endgroup$

4 Answers 4

4
$\begingroup$

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]

enter image description here

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

alphanumerical[0.5]

-3.14159 - 2.57874 I

$\endgroup$
1
$\begingroup$

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]
$\endgroup$
1
$\begingroup$

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]

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

$\endgroup$
2
  • $\begingroup$ Integrate[IntegrandEq10[x, x0], {x, x0, \[Infinity]}, Assumptions -> x0 > 3/2] - Pi returns the input. How do you come to -Im[((2*Sqrt[x0])/((-1 + x0)^(1/4)*(3 + x0)^(1/4)))* EllipticF[ ArcSin[((-1 + x0)^(1/4)*Sqrt[Sqrt[-1 + x0] + Sqrt[3 + x0]])/ Sqrt[-3 + x0 + Sqrt[-1 + x0]*Sqrt[3 + x0]]]? $\endgroup$
    – user64494
    Apr 18 at 15:21
  • $\begingroup$ @user64494 substitute x->1/y, do antiderivative wrt y, then transform F[z,k^2] to inverse modulus F[...,1/k^2] and try real or imaginary part of the result for match.I don't know where the checkabort and defer statements come from(Internal?). Copy the Plot command into a new cell and simplify the Im[...] , then the structure of the result is seen better. $\endgroup$
    – Andreas
    Apr 18 at 16:40
1
$\begingroup$

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]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.