What's the best way of getting Mathematica 7 or 8 to do the integral
NIntegrate[Exp[-x]/Sin[Pi x], {x, 0, 50}]
There are poles at every integer - and we want the Cauchy principle value. The idea is to get a good approximation for the integral from 0 to infinity.
With Integrate there is the option PrincipleValue -> True.
With NIntegrate I can give it the option Exclusions -> (Sin[Pi x] == 0), or manually give it the poles by
NIntegrate[Exp[-x]/Sin[Pi x], Evaluate[{x, 0, Sequence@@Range[50], 50}]]
The original command and the above two NIntegrate tricks give the result 60980 +/- 10. But they all spit out errors. What is the best way of getting a quick reliable result for this integral without Mathematica wanting to give errors?
Simon, is there reason to believe your integral is convergent ?
In[52]:= f[k_Integer, eps_Real] :=
NIntegrate[Exp[-x]/Sin[Pi x], {x, k + eps, k + 1 - eps}]
In[53]:= Sum[f[k, 1.0*10^-4], {k, 0, 50}]
Out[53]= 2.72613
In[54]:= Sum[f[k, 1.0*10^-5], {k, 0, 50}]
Out[54]= 3.45906
In[55]:= Sum[f[k, 1.0*10^-6], {k, 0, 50}]
Out[55]= 4.19199
It looks like the problem is at x==0. Splitting integrand k+eps to k+1-eps for integer values of k:
In[65]:= int =
Sum[(-1)^k Exp[-k ], {k, 0, Infinity}] Integrate[
Exp[-x]/Sin[Pi x], {x, eps, 1 - eps}, Assumptions -> 0 < eps < 1/2]
Out[65]= (1/((1 +
E) (I + \[Pi])))E (2 E^(-1 + eps - I eps \[Pi])
Hypergeometric2F1[1, (I + \[Pi])/(2 \[Pi]), 3/2 + I/(2 \[Pi]),
E^(-2 I eps \[Pi])] +
2 E^(I eps (I + \[Pi]))
Hypergeometric2F1[1, (I + \[Pi])/(2 \[Pi]), 3/2 + I/(2 \[Pi]),
E^(2 I eps \[Pi])])
In[73]:= N[int /. eps -> 10^-6, 20]
Out[73]= 4.1919897038160855098 + 0.*10^-20 I
In[74]:= N[int /. eps -> 10^-4, 20]
Out[74]= 2.7261330651934049862 + 0.*10^-20 I
In[75]:= N[int /. eps -> 10^-5, 20]
Out[75]= 3.4590554287709991277 + 0.*10^-20 I
As you see there is a logarithmic singularity.
In[79]:= ser =
Assuming[0 < eps < 1/32, FullSimplify[Series[int, {eps, 0, 1}]]]
Out[79]= SeriesData[eps, 0, {(I*(-1 + E)*Pi -
2*(1 + E)*HarmonicNumber[-(-I + Pi)/(2*Pi)] +
Log[1/(4*eps^2*Pi^2)] - 2*E*Log[2*eps*Pi])/(2*(1 + E)*Pi),
(-1 + E)/((1 + E)*Pi)}, 0, 2, 1]
In[80]:= Normal[
ser] /. {{eps -> 1.*^-6}, {eps -> 0.00001}, {eps -> 0.0001}}
Out[80]= {4.191989703816426 - 7.603403526913691*^-17*I,
3.459055428805136 -
7.603403526913691*^-17*I,
2.726133068607085 - 7.603403526913691*^-17*I}
EDIT Out[79] of the code above gives the series expansion for eps->0, and if these two logarithmic terms get combined, we get
In[7]:= ser = SeriesData[eps, 0,
{(I*(-1 + E)*Pi - 2*(1 + E)*HarmonicNumber[-(-I + Pi)/(2*Pi)] +
Log[1/(4*eps^2*Pi^2)] - 2*E*Log[2*eps*Pi])/(2*(1 + E)*
Pi),
(-1 + E)/((1 + E)*Pi)}, 0, 2, 1];
In[8]:= Collect[Normal[PowerExpand //@ (ser + O[eps])],
Log[eps], FullSimplify]
Out[8]= -(Log[eps]/\[Pi]) + (
I (-1 + E) \[Pi] -
2 (1 + E) (HarmonicNumber[-((-I + \[Pi])/(2 \[Pi]))] +
Log[2 \[Pi]]))/(2 (1 + E) \[Pi])
Clearly the -Log[eps]/Pi came from the pole at x==0. So if one subtracts this, just like principle value method does this for other poles you end up with a finitely value:
In[9]:= % /. Log[eps] -> 0
Out[9]= (I (-1 + E) \[Pi] -
2 (1 + E) (HarmonicNumber[-((-I + \[Pi])/(2 \[Pi]))] +
Log[2 \[Pi]]))/(2 (1 + E) \[Pi])
In[10]:= N[%, 20]
Out[10]= -0.20562403655659928968 + 0.*10^-21 I
Of course, this result is difficult to verify numerically, but you might know more that I do about your problem.
EDIT 2
This edit is to justify In[65] input that computes the original regularized integral. We are computing
Sum[ Integrate[ Exp[-x]/Sin[Pi*x], {x, k+eps, k+1-eps}], {k, 0, Infinity}] ==
Sum[ Integrate[ Exp[-x-k]/Sin[Pi*(k+x)], {x, eps, 1-eps}], {k, 0, Infinity}] ==
Sum[ (-1)^k*Exp[-k]*Integrate[ Exp[-x]/Sin[Pi*x], {x, eps, 1-eps}],
{k, 0, Infinity}] ==
Sum[ (-1)^k*Exp[-k], {k, 0, Infinity}] *
Integrate[ Exp[-x]/Sin[Pi*x], {x, eps, 1-eps}]
In the third line Sin[Pi*(k+x)] == (-1)^k*Sin[Pi*x] for integer k was used.