Как преодолеть особенности в численном интегрировании (в Matlab или Mathematica)

Question

Я хочу численно интегрировать следующее:

, где

и a, bβ - это константы, которые для простоты могут быть установлены в 1.

Ни Matlab, использующий dblquad, ни Mathematica, использующий NIntegrate, не могут справиться с сингулярностью, созданной знаменателем.Поскольку это двойной интеграл, я не могу указать, где находится особенность в Mathematica.

Я уверен, что он не бесконечен, поскольку этот интеграл основан на теории возмущений и без

был найден ранее (только не мнойтак что я не знаю, как это делается).

Есть идеи?

Answer 1

(1) Было бы полезно, если вы предоставите явный код, который вы используете. Таким образом, другим (читай: мне) не нужно кодировать его отдельно.

(2) Если интеграл существует, он должен быть равен нулю. Это потому, что вы отрицаете коэффициент n (y) -n (x), когда вы меняете местами x и y, но оставляете все остальное. Тем не менее, симметрия диапазона интеграции означает, что это просто переименование ваших переменных, следовательно, оно должно остаться прежним.

(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

Последовательное разложение в cy, центрированное в cx, указывает на полюс порядка 1. Таким образом, оно представляется сингулярным интегралом.

Даниэль Лихтблау

Answer 2

Интеграл выглядит как интеграл типа принципа Коши (например, имеет сильную особенность). Вот почему вы не можете применять стандартные квадратурные методы.

Вы пробовали PrincipalValue-> True в Mathematica's Integrate?

Answer 3

Может быть, я что-то здесь упускаю, но интеграл 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 [- х, у] = F [х, у]. Поэтому его интеграл по любой области [-u, u] x [-v, v] равен нулю. Кажется, здесь нет необходимости в численном интегрировании. Результат равен нулю.

Answer 4

В дополнение к наблюдению Даниэля об интеграции нечетного подынтегрального выражения в симметричном диапазоне (так что симметрия указывает, что результат должен быть нулевым), вы также можете сделать это, чтобы лучше понять его сходимость (я буду использовать латекс, выписав это с ручка и бумага должны облегчать чтение, написание заняло гораздо больше времени, чем на самом деле, это не так сложно):

Сначала 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 и спросить, что именно они сделали. Возможно, здесь применяется процедура физической регуляризации. Или вы могли бы дать больше информации о физическом происхождении этого (своего рода теория возмущений второго порядка для некоторой бозонной системы?), Если бы это не было не по теме здесь ...