如何克服数值积分中的奇点(在Matlab或Mathematica中) [英] How to overcome singularities in numerical integration (in Matlab or Mathematica)

问题描述

我想对以下内容进行数字积分:

I want to numerically integrate the following:

其中

和a，b和β是常量，为简单起见，可以全部设置为1.

and a, b and β are constants which for simplicity, can all be set to 1.

使用dblquad的Matlab或使用NIntegrate的Mathematica都无法处理分母创建的奇异性.由于它是一个双整数，因此我无法指定Mathematica中的奇点.

Neither Matlab using dblquad, nor Mathematica using NIntegrate can deal with the singularity created by the denominator. Since it's a double integral, I can't specify where the singularity is in Mathematica.

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

I'm sure that it is not infinite since this integral is based in perturbation theory and without the

之前被发现过(只是不是我一个人，所以我不知道它是如何完成的.)

has been found before (just not by me so I don't know how it's done).

有什么想法吗?

推荐答案

(1)如果提供使用的显式代码，将很有帮助.这样一来，其他人(阅读:我)就无需对其进行单独编码.

(1) It would be helpful if you provide the explicit code you use. That way others (read: me) need not code it up separately.

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

(2) If the integral exists, it has to be zero. This is because you negate the n(y)-n(x) factor when you swap x and y but keep the rest the same. Yet the integration range symmetry means that amounts to just renaming your variables, hence it must stay the same.

(3)这是一些代码，表明它将为零，至少如果我们将奇异部分和其周围的一小段零置零.

(3) Here is some code that shows it will be zero, at least if we zero out the singular part and a small band around it.

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加到被积数上只是为了避免结果接近零的准确性问题.

We add 1 to the integrand just to avoid accuracy issues with results that are near zero.

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

这很好地表明了纯净结果将为零.

This is a good indication that the unadulterated result will be zero.

(4)将cx替换为cos(x)，将cy替换为cos(y)，并为了进行收敛性评估而删除无关的因素，得出以下表达式.

(4) Substituting cx for cos(x) and cy for cos(y), and removing extraneous factors for purposes of convergence assessment, gives the expression below.

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

在cy上以cx为中心的级数展开表示一个1阶极点.因此，它的确是一个奇异积分.

A series expansion in cy, centered at cx, indicates a pole of order 1. So it does appear to be a singular integral.

丹尼尔·里奇布劳

这篇关于如何克服数值积分中的奇点(在Matlab或Mathematica中)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！