如何克服数值积分中的奇点（在 Matlab 或 Mathematica 中）答案

【问题标题】：How to overcome singularities in numerical integration (in Matlab or Mathematica)如何克服数值积分中的奇点（在 Matlab 或 Mathematica 中）
【发布时间】：2011-10-28 01:47:51
【问题描述】：

我想对以下内容进行数值积分：

在哪里

和a、b和β是常量，为简单起见，都可以设置为1。

使用dblquad 的Matlab 和使用NIntegrate 的Mathematica 都不能处理分母产生的奇异性。由于它是一个双积分，我无法指定 Mathematica 中的奇点在哪里。

我确信它不是无限的，因为这个积分是基于微扰理论并且没有

以前发现过（只是不是我发现的，所以我不知道是怎么做到的）。

有什么想法吗？

【问题讨论】：

在 Mathematica 中，Exclusions -> {Cos[x] == Cos[y]} 怎么样，或者您可以手动拆分范围...
或将 vars 更改为 xi_\pm=(x\pm y)/2 （这相当于旋转空间，所以你也必须注意你的限制），这样两极在每个 var 都独立于另一个 var（因为 cos(x)-cos(y) 是 xi_\pm 的可分离 fn，至少我认为是这样）。顺便说一句，这不是一个编程问题（不是我投票关闭它或任何东西）。

标签： matlab wolfram-mathematica numerical-integration

【解决方案1】：

(1) 如果您提供您使用的显式代码，将会很有帮助。这样其他人（阅读：我）不需要单独编码。

(2) 如果积分存在，它必须为零。这是因为您在交换 x 和 y 时否定了 n(y)-n(x) 因子，但其余部分保持不变。然而，积分范围对称意味着仅重命名变量，因此它必须保持不变。

(3) 这里有一些代码显示它将为零，至少如果我们将奇异部分和围绕它的小带归零。

a = 1;
b = 1;
beta = 1;
eps[x_] := 2*(a-b*Cos[x])
n[x_] := 1/(1+Exp[beta*eps[x]])
delta = .001;
pw[x_,y_] := Piecewise[{{1,Abs[Abs[x]-Abs[y]]>delta}}, 0]

我们在被积函数上加 1 只是为了避免结果接近于零的准确性问题。

NIntegrate[1+Cos[(x+y)/2]^2*(n[x]-n[y])/(eps[x]-eps[y])^2*pw[Cos[x],Cos[y]],
  {x,-Pi,Pi}, {y,-Pi,Pi}] / (4*Pi^2)

我得到下面的结果。

NIntegrate::slwcon: 
   Numerical integration converging too slowly; suspect one of the following:
    singularity, value of the integration is 0, highly oscillatory integrand,
    or WorkingPrecision too small.

NIntegrate::eincr: 
   The global error of the strategy GlobalAdaptive has increased more than 
    2000 times. The global error is expected to decrease monotonically after a
     number of integrand evaluations. Suspect one of the following: the
     working precision is insufficient for the specified precision goal; the
     integrand is highly oscillatory or it is not a (piecewise) smooth
     function; or the true value of the integral is 0. Increasing the value of
     the GlobalAdaptive option MaxErrorIncreases might lead to a convergent
     numerical integration. NIntegrate obtained 39.4791 and 0.459541
     for the integral and error estimates.

Out[24]= 1.00002

这很好地表明纯正的结果将为零。

(4) 用 cx 代替 cos(x) 和 cy 代替 cos(y)，并去除无关因素以进行收敛性评估，得到以下表达式。

((1 + E^(2*(1 - cx)))^(-1) - (1 + E^(2*(1 - cy)))^(-1))/
 (2*(1 - cx) - 2*(1 - cy))^2

以 cx 为中心的 cy 中的级数展开式表示一个 1 阶极点。因此它看起来确实是一个奇异积分。

丹尼尔·利希特布劳

【讨论】：

【解决方案2】：

积分看起来像一个 Cauchy Principal Value 类型的积分（即它具有很强的奇异性）。这就是为什么你不能应用标准的正交技术。

您是否在 Mathematica 的 Integrate 中尝试过 PrincipalValue->True？

【讨论】：

感谢@siemann，但 PrincipleValue 只能用于一维积分。

【解决方案3】：

除了 Daniel 关于在对称范围内积分奇数被积函数的观察（因此对称性表明结果应该为零），您还可以这样做以更好地理解它的收敛性（我将使用乳胶，用笔和纸应该更容易阅读；写比写要花更长的时间，没那么复杂）：

首先，epsilon(x)-\epsilon(y)\propto\cos(y)-\cos(x)=2\sin(\xi_+)\sin(\xi_-) 我在其中定义了\xi_\pm=(x\pm y)/2（所以我将轴旋转了 pi/4）。那么，整合区域是\xi_+ 在\pi/\sqrt{2} 和-\pi/\sqrt{2} 之间以及\xi_- 在\pm(\pi/\sqrt{2}-\xi_-) 之间。然后被积函数采用\frac{1}{\sin^2(\xi_-)\sin^2(\xi_+)} 乘以没有分歧的形式。所以，很明显，存在二阶极点，这并不像呈现的那样收敛。

也许您可以通过电子邮件向获得 cos 项答案的人发送电子邮件，并询问他们究竟做了什么。也许正在使用物理正则化程序。或者你可以提供更多关于这个物理起源的信息（某种玻色子系统的某种二阶微扰理论？），如果这不是题外话...

【讨论】：

【解决方案4】：

可能是我在这里遗漏了一些东西，但是被积分 f[x,y]=Cos^2[(x+y)/2]*(n[x]-n[y])/(eps[x]-eps[y]) 其中 n[x]=1 /(1+Exp[Beta*eps[x]]) 和 eps[x]=2(ab*Cos[x]) 确实是 x 和 y 中的对称函数： f[x,-y]= f[- x,y]=f[x,y]。因此它在任何域 [-u,u]x[-v,v] 上的积分为零。这里似乎不需要数值积分。结果只是零。

【讨论】：