How to find a more precise value of integral?

Question

Mathematica 11 produces

NIntegrate[ArcCot[x]*Sin[x]/(5/4 + Cos[x]), {x, 0, Infinity}, 
 Method -> "ExtrapolatingOscillatory", AccuracyGoal -> 10, 
 WorkingPrecision -> 20]
0.746733605129940959

This is 4.573.3 from Gradshteyn & Ryzhik. Up to G&R the result is

 N[Pi*Log[3/2/(1 + 1/2*Exp[-1])], 10]
   0.7433557514

The question arises: how to obtain it in Mathematica?

Answer 1

Here's an extrapolatory approach based on Richardson extrapolation of the integrals over each "period." (Richardson[] code from Anton Antonov's answer here):

Clear[Richardson]
Richardson[A_, n_, N_] := 
 Total@Table[(A[n + k]*(n + k)^N*If[OddQ[k + N], -1, 1])/(k! (N - k)!), {k, 0, N}]

(* the integral from 2 n Pi to 2 (n+1) Pi *)
Clear[aa];
aa[n_Integer /; n >= 0] := 
  aa[n] = NIntegrate[ArcCot[x]*Sin[x]/(5/4 + Cos[x]), {x, 2 n Pi, 2 (n + 1) Pi}, 
    WorkingPrecision -> 50];

We get around machine precision with 100 terms plus 6 extrapolation terms:

sf[n_] := Sum[aa[i], {i, 0, n}]; (* partial sum = integral over {0, 2 (n+1) Pi} *)
Richardson[sf, 100, 6]
% - Pi*Log[3/2/(1 + 1/2*Exp[-1])]
(*
  0.743355751361227474780881457479018423835    
  -2.97535713733265290609713*10^-16
*)

Without extrapolation:

sf[100]
% - Pi*Log[3/2/(1 + 1/2*Exp[-1])]
(*
  0.74207790286093203962531860751097243182349263171074
  -0.00127784850029573269127658323333660172377741318898
*)

---

Update

More extrapolatory terms yields an approximation accurate to machine precision with many fewer terms overall:

Richardson[sf, 1, 16]
% - Pi*Log[3/2/(1 + 1/2*Exp[-1])]
(*
  0.74335575136122767678505611095516483469029
  -9.553153907978914419885698*10^-17
*)

Answer 2

You didn't say how accurate you wanted it to be.

You can convert the expression to exponential form, which can be easier to NIntegrate[].

expr = ArcCot[x]*Sin[x]/(5/4 + Cos[x]) // TrigToExp;

NIntegrate[expr, {x, 0, ∞}, 
 Method -> "ExtrapolatingOscillatory", AccuracyGoal -> 80, 
 WorkingPrecision -> 160, MaxRecursion -> 100, MinRecursion -> 50]

This produces as answer 0.74363... with the warnings you anticipated. Further increasing values shows incrementally more accurate results.

Not really great.

So, let's design our intersections manually:

AbsoluteTiming[
 Total[Table[
   NIntegrate[expr, {x, i - Pi, i}, AccuracyGoal -> 10, 
    WorkingPrecision -> 20], {i, Pi, 10000 Pi, Pi}]]]

{57.0974, 0.7433299386508148184}

Takes a while, but it's a lot more accurate, naturally the further you take it, the more accurate it becomes, but I don't have infinite computing time right now.

Answer 3

Although the question was most probably meant to address accuracy problems in NIntegrate it might be of interest to see how to use Mathematica to calculate the exact symbolic result of the integral.

We shall derive the exact expression for the more general integral of Gradshteyn/Ryshik 4.573.3:

fi := Integrate[
  ArcCot[r x] Sin[p x]/(1 + 2 q Cos[p x] + q^2), {x, 0, \[Infinity]}]

and then take the numerical value in the end.

fi reduces to the integral of the OP by letting

rep = {r -> 1, p -> 1, q -> 2};
fop = q^2 fi /. rep;

As Mathematica does not evaluate the integral we make a series expansion with respect to the parameter q around q = 0 and integrate term by term which gives us for the first few terms (for simplicity in list format)

(2 p)/\[Pi] Integrate[
    ArcCot[r x] List @@ 
      Normal[Series[Sin[p x]/(1 + 2 q Cos[p x] + q^2), {q, 0, 3}]], {x, 
     0, \[Infinity]}, Assumptions -> {q^2 < 1, p > 0, r > 0}] // 
  Expand // Simplify

(* Out[2]= {1 - E^(-(p/r)), 1/2 (-1 + E^(-((2 p)/r))) q, 1/3 (1 - E^(-((3 p)/r))) q^2, 1/4 (-1 + E^(-((4 p)/r))) q^3} *)

It is easy to guess the general scheme and see that the these terms are generated by the formula

-1/q Table[(-q)^k (1 - Exp[-k p/r])/k, {k, 1, 4}] // Expand // Simplify

(* Out[3]= {1 - E^(-(p/r)), 1/2 (-1 + E^(-((2 p)/r))) q, 1/3 (1 - E^(-((3 p)/r))) q^2, 
 1/4 (-1 + E^(-((4 p)/r))) q^3} *)

Extending the range of k of the sum to infinity and assuming q^2<1, thus taking care of convergence, we obtain

filt1 = -\[Pi]/(2 p q) Sum[(-q)^k (1 - Exp[-k p/r])/k, {k, 1, \[Infinity]}] //
   Simplify

(* Out[35]= (\[Pi] (Log[1 + q] - Log[1 + E^(-(p/r)) q]))/(2 p q) *)

This formula holds for q^2<1. For q^2 > 1 we write w=1/q in the definition of fi and arrive at the formula

figt1 = w^2 filt1 /. q -> w /. w -> 1/q

(* Out[36]= (\[Pi] (Log[1 + 1/q] - Log[1 + E^(-(p/r))/q]))/(2 p q) *)

Finally, inserting the values of the parameters we find for the integral of the OP

fop = q^2 figt1 /. rep

(* Out[38]= \[Pi] (Log[3/2] - Log[1 + 1/(2 E)]) *)

And its numerical value is

N[%, 20]

(* Out[39]= 0.74335575136122777232 *)