# Finding closed or numeric value for $\int_0^{\infty } \frac{x \csc (a x)}{x^2+b^2} \, dx$

According to Gradshteyn 3.747-3 $$\int_0^{\infty } \frac{x \csc (a x)}{x^2+b^2} \, dx=\frac{\pi}{2\sinh(ab)}$$ $$b>0$$

so I'm trying to get some numeric results from Mathematica 12.

a = 6; b = 12;

NIntegrate[(x/((x^2 + b^2) Sin[a x])), {x, 0, Infinity}, MaxRecursion -> 12]


However this does not match the values given by

Pi/(2*Sinh[a*b])


So it is not clear where I'm going wrong.

• How did you evaluate the NIntegrate for undefined parameters a,b ? Please show some code! Aug 3, 2020 at 15:44
• @Ulrich Neumann I just did a = 6; b = 12; but I also just put them in manually too to be sure. Aug 3, 2020 at 15:46
• There are singularities in the integral when a*x = Pi *k for some integer k, and when the intergral goes through the point x = b Aug 3, 2020 at 15:46
• @flinty is there any code to deal with these? Aug 3, 2020 at 15:49
• ^ 3-747.3 seems to be discussed here on page 9. researchgate.net/publication/… Aug 3, 2020 at 17:04

I want to focus on the numerics here. There is a way to do this integral taking infinite number of singularities into account.

One can split the integral into the domains containing only one singularity at $$a x_n=\pi n$$, i.e. $$x\in[x_n-\frac{\pi}{2a},x_n+\frac{\pi}{2a}]$$ and an integral in the interval $$x\in[0,\frac{\pi}{2a}]$$. On the last step we use NSum with some options. The options are very important. Consider first my failed attempt:

Clear[f]
f[n_?NumericQ,a_,b_]:=NIntegrate[(x/((x^2+b^2) Sin[a x])),{x,(π n)/a-π/(2a),(π n)/a,(π n)/a+π/(2a)},WorkingPrecision->100,Method->PrincipalValue]
f0[a_,b_]:=NIntegrate[(x/((x^2+b^2) Sin[a x])),{x,0,π/(2a)},WorkingPrecision->100]

f0[6,12]+NSum[f[n,6,12],{n,1,∞},Method>"AlternatingSigns",WorkingPrecision->100]
Out[1]= 6.39989549924364176901258523623081516506764870738550937643852237103123602088582134002479513849115*10^-6


The numerical value seems to be highly accurate and well converged, but is different from the analytic one

e[a_,b_]:=N[Pi/(2*Sinh[a*b]),100]
e[6,12]
Out[2]= 1.690235331526788818439805170791473807196429480676031631266783609275102725281157127348346908376558298*10^-31


The discrepancy can be corrected by changing the default NSumTerms value. Let us introduce the relative error as a function of this parameter

relErr[a_,b_,k_]:=Abs[f0[a,b]+NSum[f[n,a,b],{n,1,∞},
Method->"AlternatingSigns",
WorkingPrecision->100,
NSumTerms->k]
-e[a,b]]/Abs[e[a,b]]//N


and try a few values:

relErr[6,12,21]
Out[3]= 3.78639*10^25
relErr[6,12,22]
Out[4]= 2.38766*10^-68
relErr[6,12,25]
Out[5]= 2.3903*10^-68


Thus, it is important to explicitly consider the first 22 terms!

• It should be noticed that f0[6, 12] + NSum[f[n, 6, 12], {n, 1, \[Infinity]}, WorkingPrecision -> 100] performs 0.*10^-7. Aug 3, 2020 at 20:42
• You don't take into account the singularity at the origin in your f0[a_,b_]. Aug 3, 2020 at 20:45
• @user64494 Thank you for checking other options. 1. According to the documentation one needs to "use a method for alternating series to get a very precise sum approximation". 2. This is a removable singularity (the limit $x\rightarrow 0$) is finite. Aug 3, 2020 at 20:51
• -1. In you answer you use e[a_,b_]:=N[Pi/(2*Sinh[a*b]),100] . Imagine you don't know it. Such sort answer does not make a good impression. Aug 4, 2020 at 4:04
• I'd like to add that the NSumTerms->k option means you consider a truncated sum and that sum is a sum for the principal value of the integral over $(0,\frac {\pi k} a)$. This is not the required result. Aug 4, 2020 at 4:30

Using the substitution  u==a x the identity is transformed to Integrate[u/((u^2 + ab^2) Sin[u]), {u, 0, Infinity}] == Pi/(2 Sinh[ab])

The integral is singular at u==k Pi, k=1,2,... and the integration range is splitted into subintervals (similar to @yarchik answer) containing only one singularity (thanks to @flinty and @ ChipHurst comments):

int[ab_?NumericQ, n_?NumericQ] :=NIntegrate[u/((u^2 + ab^2) Sin[u]) , {u, 0, Pi/2}] +
Sum[ NIntegrate[u/((u^2 + ab^2) Sin[u]) , {u, k Pi - Pi/2, k Pi + Pi/2},Method -> PrincipalValue, Exclusions -> {u == k Pi}], {k, 1, n} ]


This finite sum matchs quite well:

Plot[{int[ab, 10], Pi/(2 Sinh[ab])}, {ab, 0, 1},PlotStyle -> {{Thickness[0.01],Lighter[Blue]}, Red},
PlotLabel ->"{int[ab,10],\!$$\*FractionBox[\(Pi$$, $$2\\\ Sinh[ab]$$]\)}",AxesLabel -> {"ab ", None}]


• This is a modification of the @yarchik's approach. Aug 4, 2020 at 7:48
• +1. int[ab, Infinity] works well to me. I think you give a correct numerical answer. Aug 4, 2020 at 8:20
• @user64494 Thanks, the "modification" shows a straightforward approach. Ticky idea int[ab, Infinity] Aug 4, 2020 at 8:22
• I changed my mind. In fact, it does not work: int[72, Infinity] returns the input. I was missed by the plot of  Pi/(2 Sinh[ab]. Aug 4, 2020 at 8:25
• @yarchik Your solution is great! I tried to avoid such high values WorkingPrecision->100 in my answer. You can't see convergence in my comment, but a trend of decaying values. Aug 5, 2020 at 11:38

Following the @flinty's comment, we obtain

Residue[x/(x^2 + b^2)/Sin[a*x], {x, π/a*n},  Assumptions -> n ∈ Integers]
(*((-1)^-n n π)/(a^2 b^2 + n^2 π^2)*)

Sum[%, {n, -∞, ∞}]
(*0*)


and

2*π*I*Residue[x/(x^2 + b^2)/Sin[a*x], {x, I*b}]
(*π Csch[a b]*)


Now, making use of the Jordan's lemma, we conclude that $$PV\int_{-\infty}^\infty\frac x {(x^2+b^2)\sin(ax)}\,dx= \pi \text{csch}(a b).$$ It remains to apply the parity of the integrand.