具有许多奇点的 Mathematica 积分

发布于 2024-10-30 18:10:21 字数 582 浏览 1 评论 0原文

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

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

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

对于Integrate，有一个选项PrincipleValue ->确实如此。

使用NIntegrate，我可以给它选项Exclusions -> (Sin[Pi x] == 0)，或者通过原始命令和上述两个 NIntegrate 技巧手动给出极点，

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

得到结果 60980 +/- 10< /代码>。但他们都吐出了错误。在 Mathematica 不想给出错误的情况下，获得快速可靠的积分结果的最佳方法是什么？

原文

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?

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

绮烟 2024-11-06 18:10:21

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

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

第三行 Sin[Pi*(k+x)] == (-1)^k*Sin[Pi*x] 对于整数 k 使用。

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.

回复收藏 0 原文

下雨或天晴 2024-11-06 18:10:21

西蒙，我没有花太多时间研究你的积分，但你应该尝试看看稳态相位近似。您拥有的是平滑函数 (exp) 和高度振荡函数 (sine)。现在涉及的工作是将 1/sin(x) 转换成 exp(if(x)) 的形式。

另外，您也可以使用余割（在极点无效）：

In[1]:=Series[Csc[x], {x, 0, 5}]
(formatted) Out[1]=1/x + x/6 + 7/360 x^3 + 31/15120 x^5 +O[x]^6

请注意，对于所有 m>-1，您有以下结果：

In[2]:=Integrate[x^m Exp[-x], {x, 0, Infinity}, Assumptions -> m > -1]
Out[2]=Gamma[1+m]

但是，将级数与余割系数相加（来自维基百科），不包括 1/x Exp[-x] 情况，它不收敛于 [0,Infinity]。

c[m_] := (-1)^(m + 1) 2 (2^(2 m - 1) - 1) BernoulliB[2 m]/Factorial[2 m];
Sum[c[m] Gamma[1 + 2 m - 1], {m, 1, Infinity}]

也不收敛...

所以，我不确定您是否可以计算出无穷大积分的近似值，但如果您对达到某个大 N 的解决方案感到满意，我希望这些有所帮助。

Simon, I haven't spent much time with your integral, but you should try looking at stationary phase approximation. What you have is a smooth function (exp), and a highly oscillatory function (sine). The work involved is now in brow-beating the 1/sin(x) into the form exp(if(x))

Alternatively, you could use the series expansion of the cosecant (not valid at poles):

In[1]:=Series[Csc[x], {x, 0, 5}]
(formatted) Out[1]=1/x + x/6 + 7/360 x^3 + 31/15120 x^5 +O[x]^6

Note that for all m>-1, you have the following:

In[2]:=Integrate[x^m Exp[-x], {x, 0, Infinity}, Assumptions -> m > -1]
Out[2]=Gamma[1+m]

However, summing the series with the coefficients of cosecant (from wikipedia), not including 1/x Exp[-x] case, which doesn't converge on [0,Infinity].

c[m_] := (-1)^(m + 1) 2 (2^(2 m - 1) - 1) BernoulliB[2 m]/Factorial[2 m];
Sum[c[m] Gamma[1 + 2 m - 1], {m, 1, Infinity}]

does not converge either...

So, I'm not sure that you can work out an approximation for the integral to infinity, but I if you're satisfied with a solution upto some large N, I hope these help.

回复收藏 0 原文

慈悲佛祖 2024-11-06 18:10:21

我必须同意 Sasha，积分似乎不是收敛。但是，如果排除 x == 0 并将积分分解为

Integrate[Exp[-x]/Sin[Pi x], {x, n + 1/2, n + 3/2}, PrincipalValue -> True]

n >= 0 && 的片段。 Element[n, Integers]，那么看来你可能会得到一个交替的序列

I Sum[ (-1/E)^n, {n, 1, Infinity}] == - I / (1 + E )

，现在我只取出到n == 4，但看起来很合理。然而，对于上面的积分假设 ->元素[n, 整数] && n >= 0 Mathematica 给出的

If[ 2 n >= 1, - I / E, Integrate[ ... ] ]

结果并不符合具体情况。作为补充说明，如果极点位于积分区域的边界，即您的极限是 {x, n, n + 1}，您只能得到 DirectedInfinity 。快速查看该图意味着您在限制 {x, n, n + 1} 的情况下只有严格的正或负被积函数，因此无限值可能是由于缺乏补偿{x, n + 1/2, n + 3/2} 为您提供。使用 {x, n, n + 2} 检查，但它只输出未计算的积分。

I have to agree with Sasha, the integral does not appear to be convergent. However, if you exclude x == 0 and break the integral into pieces

Integrate[Exp[-x]/Sin[Pi x], {x, n + 1/2, n + 3/2}, PrincipalValue -> True]

where n >= 0 && Element[n, Integers], then it seems you may get an alternating series

I Sum[ (-1/E)^n, {n, 1, Infinity}] == - I / (1 + E )

Now, I only took it out to n == 4, but it looks reasonable. However, for the integral above with Assumptions -> Element[n, Integers] && n >= 0 Mathematica gives

If[ 2 n >= 1, - I / E, Integrate[ ... ] ]

which just doesn't conform to the individual cases. As an additional note, if the pole lies at the boundary of the integration region, i.e. your limits are {x, n, n + 1}, you only get DirectedInfinitys. A quick look at the plot implies that you with the limits {x, n, n + 1} you only have a strictly positive or negative integrand, so the infinite value may be due to the lack of compensation which {x, n + 1/2, n + 3/2} gives you. Checking with {x, n, n + 2}, however it only spits out the unevaluated integral.

回复收藏 0 原文

~没有更多了~