# Numerical_integration

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!

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

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


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]


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...

• 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]]]? Apr 18 at 15:21
• @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. Apr 18 at 16:40

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]