wolfram-mathematica - 具有许多奇点的 Mathematica 积分

Question

让 Mathematica 7 或 8 进行积分的最佳方法是什么

NIntegrate[Exp[-x]/Sin[Pi x], {x, 0, 50}]

每个整数都有极点 - 我们需要柯西原理值。这个想法是为了获得从 0 到无穷大的积分的良好近似值。

Integrate 有选项 PrincipleValue -> True。

使用 NIntegrate 我可以给它选项Exclusions -> (Sin[Pi x] == 0)，或者手动给它极点

NIntegrate[Exp[-x]/Sin[Pi x], Evaluate[{x, 0, Sequence@@Range[50], 50}]]

原始命令和上面的两个 NIntegrate 技巧给出了结果 60980 +/- 10。但是他们都吐出错误。在 Mathematica 不想给出错误的情况下快速获得该积分的可靠结果的最佳方法是什么？

Answer 1

Simon，有理由相信你的积分是收敛的吗？

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

看起来问题出在 x==0。对于 k 的整数值，将被积函数 k+eps 拆分为 k+1-eps:

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

如您所见，存在对数奇点。

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}

编辑上面代码的 Out[79] 给出了 eps->0 的级数展开，如果这两个对数项合并，我们得到

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

很明显，-Log[eps]/Pi 来自 x==0 处的极点。因此，如果减去这个，就像主值方法对其他极点执行此操作一样，您最终会得到一个有限值:

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

当然，这个结果很难用数值来验证，但你可能比我对你的问题了解得更多。

编辑 2

此编辑是为了证明计算原始正则化积分的 In[65] 输入是正确的。我们在计算

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

在第三行中，使用了整数 k 的 Sin[Pi*(k+x)] == (-1)^k*Sin[Pi*x]。