Jak pokonać osobliwości w integracji numerycznej (w Matlab czy Mathematica)

Question

Chcę numerycznie integrować następujące:

gdzie

i a, b i β są stałe, które dla uproszczenia można ustawić na 1.

Ani Matlab używający dblquad, ani Mathematica używający NIntegrate nie poradzą sobie ze szczególną osobowością stworzoną przez mianownik. Ponieważ jest to podwójna integracja, nie mogę określić, gdzie znajduje się osobliwość w Mathematica.

Jestem pewien, że to nie jest nieskończony ponieważ integralny z siedzibą w rachunku zaburzeń i bez

stwierdzono wcześniej (nie tylko przeze mnie, więc nie wiem jak to jest Gotowe).

Wszelkie pomysły?

Answer 1

(1) Byłoby pomocne, gdyby podać jawny kod, którego używasz. W ten sposób inni (czytaj: ja) nie muszą kodować go osobno.

(2) Jeśli całka istnieje, musi wynosić zero. Dzieje się tak, ponieważ negujesz czynnik n (y) -n (x), gdy zamieniasz x i y, ale zachowujesz resztę tak samo. Jednak symetria zakresu integracji oznacza, że ​​kwoty tylko zmieniają nazwy zmiennych, a więc muszą pozostać takie same.

(3) Oto kod, który pokazuje, że będzie zero, przynajmniej jeśli wyzerujemy pojedynczą część i mały pas wokół niej.

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] 

Dodajemy 1 do integrandu, aby uniknąć problemów z dokładnością przy wynikach bliskich zeru.

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) 

Otrzymuję wynik poniżej.

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 

Jest to dobry sygnał, że nieskażony wynik będzie wynosił zero.

(4) Zastępując cx dla cos (x) i cy dla cos (y) i usuwając czynniki zewnętrzne dla celów oceny zbieżności, podaje poniższe wyrażenie.

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

Rozszerzenie serii w cy, wyśrodkowany na cx, wskazuje biegun rzędu 1. Tak więc wydaje się być całką pojedynczą.

Daniel Lichtblau

Answer 2

Całka wygląda jak całka typu Cauchy Principal Value (tzn. Ma silną osobliwość). Dlatego nie można stosować standardowych technik kwadraturowych.

Czy próbowałeś PrincipalValue-> True w integracji z Mathematica?

Answer 3

Oprócz obserwacji Daniela o integracji dziwny całki w przedziale symetrycznym (tak, że symetria oznacza, że ​​wynik powinien wynosić zero), można też to zrobić, aby zrozumieć jego zbieżność lepiej (będę użyj lateksu, napisanie tego za pomocą pióra i papieru powinno ułatwić czytanie, pisanie było dużo dłuższe niż pisanie, nie jest to aż tak skomplikowane):

Po pierwsze, epsilon(x)-\epsilon(y)\propto\cos(y)-\cos(x)=2\sin(\xi_+)\sin(\xi_-) gdzie zdefiniowałem \xi_\pm=(x\pm y)/2 (więc ja ' ve obróciły osie o pi/4). Region integracji to \xi_+ między \pi/\sqrt{2} i -\pi/\sqrt{2} i \xi_- między \pm(\pi/\sqrt{2}-\xi_-). Następnie integrand przyjmuje formę razy, bez żadnych rozbieżności. Tak więc, widocznie są bieguny drugiego rzędu, a to nie jest zbieżne, jak przedstawiono.

Być może mógłbyś wysłać wiadomość e-mail do osób, które uzyskały odpowiedź z terminem cos i zapytać, co dokładnie zrobili. Być może jest stosowana procedura fizycznej regulacji. Czy mógłbyś podać więcej informacji na temat fizycznego pochodzenia tego (jakiś rodzaj teorii perturbacji drugiego rzędu dla jakiegoś bosonośnego systemu?), Gdyby to nie było nie na temat tutaj ...

Answer 4

Może być brakujący coś tutaj, ale integran f [x, y] = Cos^2 [(x + y)/2] * (n [x] -n [y])/(eps [x] -eps [y]) z n [x] = 1/(1 + Exp [Beta * eps [x]]) i eps [x] = 2 (ab * Cos [x]) jest rzeczywiście funkcją symetryczną w x i y: f [x, -y] = f [-x, y] = f [x, y]. Dlatego jego integralna część w dowolnej domenie [-u, u] x [-v, v] wynosi zero. Wydaje się, że tutaj nie ma potrzeby integracji numerycznej. Wynik wynosi zero.