Come superare singolarità di integrazione numerica (in Matlab o Mathematica)

Question

Voglio integrare numericamente quanto segue:

dove

e a, b e β sono costanti che per semplicità possono essere impostate su 1.

Né Matlab che utilizza dblquad, né Mathematica che utilizza NIntegrate può gestire la singolarità creata dal denominatore. Dato che è un integrale doppio, non posso specificare dove si trova la singolarità in Mathematica.

Sono sicuro che non è infinita in quanto questo integrale si basa in teoria delle perturbazioni e senza la

è stato trovato prima (non solo da me, quindi non so come sia fatto).

Qualche idea?

Answer 1

(1) Sarebbe utile se fornisci il codice esplicito che usi. In questo modo gli altri (leggi: io) non devono codificarlo separatamente.

(2) Se l'integrale esiste, deve essere zero. Questo perché si nega il fattore n (y) -n (x) quando si scambiano x e y, ma si mantiene il resto uguale. Tuttavia, la simmetria dell'intervallo di integrazione significa semplicemente rinominare le variabili, quindi deve rimanere invariata.

(3) Ecco un codice che mostra che sarà zero, almeno se azzeriamo la parte singolare e una piccola banda attorno ad esso.

Aggiungere 1 all'integrand solo per evitare problemi di precisione con risultati vicini allo 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) 

Ottengo il risultato qui sotto.

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 

Questa è una buona indicazione che il risultato non adulterato sarà zero.

(4) Sostituendo cx per cos (x) e cy per cos (y) e rimuovendo fattori estranei ai fini della valutazione della convergenza, viene fornita l'espressione seguente.

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

Un'espansione serie nel cy, centrata a cx, indica un polo di ordine 1. Quindi sembra essere un integrale singolare.

Daniel Lichtblau

Answer 2

L'integrale ha l'aspetto di un integrale di tipo Valore principale Cauchy (vale a dire ha una forte singolarità). Ecco perché non è possibile applicare le tecniche di quadratura standard.

Hai provato PrincipalValue-> True in Mathematica's Integrate?

Answer 3

Oltre all'osservazione di Daniel di integrare un integrando strana su una gamma simmetrica (in modo che la simmetria indica il risultato dovrebbe essere zero), si può anche fare questo per capire la sua convergenza meglio (io usare lattice, scrivendo questo fuori con carta e penna dovrebbe rendere più facile la lettura, c'è voluto molto più tempo di scrivere che a fare, non è così complicato):

In primo luogo, epsilon(x)-\epsilon(y)\propto\cos(y)-\cos(x)=2\sin(\xi_+)\sin(\xi_-) dove ho definito \xi_\pm=(x\pm y)/2 (così ho' abbiamo ruotato gli assi di pi/4). La regione di integrazione è quindi \xi_+ tra \pi/\sqrt{2} e -\pi/\sqrt{2} e \xi_- tra \pm(\pi/\sqrt{2}-\xi_-). Quindi l'integrando assume il formato \frac{1}{\sin^2(\xi_-)\sin^2(\xi_+)} termini senza divergenze. Quindi, evidentemente, ci sono poli di secondo ordine, e questo non è convergente come presentato.

Forse potreste inviare per e-mail le persone che hanno ottenuto una risposta con il termine cos e chiedere che cosa esattamente hanno fatto. Forse c'è una procedura di regolarizzazione fisica che viene impiegata. Oppure avresti potuto fornire maggiori informazioni sull'origine fisica di questo (una sorta di teoria perturbativa del secondo ordine per una sorta di sistema bosonico?), Se non fosse stato off-topic qui ...

Answer 4

Forse mi manca qualcosa qui, ma l'integrando f [x, y] = Cos^2 [(x + y)/2] * (n [x] -n [y])/(eps [x] -eps [y]) con n [x] = 1/(1 + Exp [Beta * eps [x]]) ed eps [x] = 2 (ab * Cos [x]) è in effetti una funzione simmetrica in x e y: f [x, -y] = f [-x, y] = f [x, y]. Quindi il suo integrale su qualsiasi dominio [-u, u] x [-v, v] è zero. Nessuna integrazione numerica sembra essere necessaria qui. Il risultato è solo zero.