# How to find a more precise value of integral?

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?

• I guess G&R is Gradshteyn & Ryzhik – Michael E2 Aug 31 '16 at 10:46
• @Michael E2:Thank you for your interest to the question. You are right. – user64494 Aug 31 '16 at 10:50
• , AccuracyGoal -> 80, WorkingPrecision -> 160 yields 0.74394... Other than that, standard techniques don't do much. – Feyre Aug 31 '16 at 11:46
• @Feyre: Thank you for your interest to the question. I am troubled by no warning and error estimate. – user64494 Aug 31 '16 at 11:56
• Your can get closer approximation if you split range {0,Infinity} into say {0,300 Pi} and {300 Pi, Infinity}. For the second interval "ExtrapolatingOscillatory" and high working presicion is ok. The first range is more complicated and even using high working precision and increased MaxRecursion hardly helps to get correct 5-th digit. – user18792 Aug 31 '16 at 12:43

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
*)

• This is effectively a manual reimplementation of Longman's ("ExtrapolatingOscillatory") method, with the only possible difference of Michael using Richardson extrapolation, and the internal routines using the Wynn $\varepsilon$ algorithm (via SequenceLimit[]). – J. M. will be back soon Oct 5 '16 at 12:23

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.

• @ Feyre: Many thanks from me to you for your work. I think the accuracy obtained by you is sufficient in most cases. The purpose of the question is the following. The integrand is not very complicated and oscillating. Because of this a quite accurate answer done without tricks is expected, not with the relative error 0.005. I believe this can be reached so I forbear to accept your answer as yet. I voted up your answer. – user64494 Aug 31 '16 at 12:35
• @user64494 Don't worry, I understand you wanting a more formal answer, which has now come. – Feyre Aug 31 '16 at 12:43

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 *)

• @ Dr. Wolfgang Hintze : Many thanks from me to you for your work well done. – user64494 Sep 4 '16 at 16:53