4
$\begingroup$

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$

Gradshteyn

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.

$\endgroup$
11
  • $\begingroup$ How did you evaluate the NIntegrate for undefined parameters a,b ? Please show some code! $\endgroup$ Aug 3, 2020 at 15:44
  • $\begingroup$ @Ulrich Neumann I just did a = 6; b = 12; but I also just put them in manually too to be sure. $\endgroup$
    – onepound
    Aug 3, 2020 at 15:46
  • $\begingroup$ 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 $\endgroup$
    – flinty
    Aug 3, 2020 at 15:46
  • $\begingroup$ @flinty is there any code to deal with these? $\endgroup$
    – onepound
    Aug 3, 2020 at 15:49
  • 1
    $\begingroup$ ^ 3-747.3 seems to be discussed here on page 9. researchgate.net/publication/… $\endgroup$
    – flinty
    Aug 3, 2020 at 17:04

3 Answers 3

5
$\begingroup$

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!

$\endgroup$
6
  • $\begingroup$ It should be noticed that f0[6, 12] + NSum[f[n, 6, 12], {n, 1, \[Infinity]}, WorkingPrecision -> 100] performs 0.*10^-7. $\endgroup$
    – user64494
    Aug 3, 2020 at 20:42
  • $\begingroup$ You don't take into account the singularity at the origin in your f0[a_,b_]. $\endgroup$
    – user64494
    Aug 3, 2020 at 20:45
  • $\begingroup$ @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. $\endgroup$
    – yarchik
    Aug 3, 2020 at 20:51
  • $\begingroup$ -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. $\endgroup$
    – user64494
    Aug 4, 2020 at 4:04
  • $\begingroup$ 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. $\endgroup$
    – user64494
    Aug 4, 2020 at 4:30
4
$\begingroup$

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

enter image description here

$\endgroup$
10
  • $\begingroup$ This is a modification of the @yarchik's approach. $\endgroup$
    – user64494
    Aug 4, 2020 at 7:48
  • $\begingroup$ +1. int[ab, Infinity] works well to me. I think you give a correct numerical answer. $\endgroup$
    – user64494
    Aug 4, 2020 at 8:20
  • $\begingroup$ @user64494 Thanks, the "modification" shows a straightforward approach. Ticky idea int[ab, Infinity] $\endgroup$ Aug 4, 2020 at 8:22
  • $\begingroup$ 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]`. $\endgroup$
    – user64494
    Aug 4, 2020 at 8:25
  • 1
    $\begingroup$ @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. $\endgroup$ Aug 5, 2020 at 11:38
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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