¿Cómo superar las singularidades en la integración numérica (en Matlab o Mathematica)

Question

Quiero integrar numéricamente los siguientes:

donde

y a, b y β son constantes que, por simplicidad, pueden establecerse en 1.

Ni Matlab usando dblquad, ni Mathematica usando NIntegrate pueden tratar con la singularidad creada por el denominador. Como es una integral doble, no puedo especificar dónde está la singularidad en Mathematica.

Estoy seguro de que no es infinita ya que esta integral se basa en la teoría de las perturbaciones y sin la

se ha encontrado antes (no por mí, así que no sé cómo es hecho).

¿Alguna idea?

Answer 1

(1) Sería útil si proporciona el código explícito que utiliza. De esa forma, otros (léase: yo) no necesitan codificarlo por separado.

(2) Si existe la integral, tiene que ser cero. Esto se debe a que anulas el factor n (y) -n (x) cuando intercambias x e y, pero mantienes el resto igual. Sin embargo, la simetría del rango de integración significa que equivale simplemente a cambiar el nombre de las variables, por lo tanto, debe permanecer igual.

(3) Aquí hay un código que muestra que será cero, al menos si ponemos a cero la parte singular y una pequeña banda a su alrededor.

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] 

Agregamos 1 al integrando solo para evitar problemas de precisión con resultados cercanos a cero.

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) 

obtengo el resultado a continuación.

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 

Este es un buen indicio de que el resultado no adulterado será cero.

(4) Sustituyendo cx por cos (x) y cy por cos (y), y eliminando factores extraños para propósitos de evaluación de convergencia, da la expresión a continuación.

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

Una expansión de serie en cy, centrada en cx, indica un polo de orden 1. Por lo tanto, parece ser una integral singular.

Daniel Lichtblau

Answer 2

La integral se ve como una integral de Tipo de valor principal de Cauchy (es decir, tiene una singularidad fuerte). Es por eso que no puedes aplicar técnicas estándar de cuadratura.

¿Has probado PrincipalValue-> True in Mathematica's Integrate?

Answer 3

Además de la observación de Daniel acerca de la integración de un integrando extraña en un rango simétrico (por lo que la simetría indica el resultado debe ser cero), también se puede hacer esto para entender su convergencia mejor (voy usar látex, escribir esto con lápiz y papel debería facilitar la lectura; tardó mucho más tiempo escribir que hacerlo, no es tan complicado):

Primero, epsilon(x)-\epsilon(y)\propto\cos(y)-\cos(x)=2\sin(\xi_+)\sin(\xi_-) donde he definido \xi_\pm=(x\pm y)/2 (entonces ' he girado los ejes en pi/4). La región de integración es \xi_+ entre \pi/\sqrt{2} y -\pi/\sqrt{2} y \xi_- entre \pm(\pi/\sqrt{2}-\xi_-). Entonces el integrando toma la forma \frac{1}{\sin^2(\xi_-)\sin^2(\xi_+)} veces términos sin divergencias. Entonces, evidentemente, hay polos de segundo orden, y esto no es convergente como se presentó.

Quizás podría enviar un correo electrónico a las personas que obtuvieron una respuesta con el término cos y preguntar qué es exactamente lo que hicieron. Tal vez hay un procedimiento de regularización física que se emplea. ¿O podría haber dado más información sobre el origen físico de esto (algún tipo de teoría de perturbación de segundo orden para algún tipo de sistema bosónico?), Si no hubiera estado fuera de tema aquí ...

Answer 4

Puede ser que me esté faltando algo aquí, pero el 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]]) y eps [x] = 2 (ab * Cos [x]) es de hecho una función simétrica en x e y: f [x, -y] = f [-x, y] = f [x, y]. Por lo tanto, su integral sobre cualquier dominio [-u, u] x [-v, v] es cero. No parece ser necesaria una integración numérica aquí. El resultado es solo cero.