NIntegrate can't solve this complicated multiple integration with absolute value in Exp and Cos functions

Question

I'm trying to use NIntegrate to calculate this complicated multiple integration, which has absolute value in the Exp and Cos functions. $$ \int_{-\infty}^{\infty} dx_1 \int_{-\infty}^{\infty} dx_2 \int_{-\infty}^{\infty} dy_1 \int_{-\infty}^{\infty} dy_2 \quad \frac{1}{r_1 r_2} \exp(-|x_1-x_2|-2|y_1-y_2|) \cos(3 |x_1-x_2|) $$ where $r_1=\sqrt{x_1^2+y_1^2}$ and $r_2=\sqrt{x_2^2+y_2^2}$

I've tried,

eqn = (Exp[-RealAbs[x1 - x2] - 2 RealAbs[y1 - y2]] Cos[3 RealAbs[x1 - x2]])/(Sqrt[x1^2 + y1^2] Sqrt[x2^2 + y2^2])
NIntegrate[eqn, {x1, -Infinity, Infinity}, {x2, -Infinity, 
  Infinity}, {y1, -Infinity, Infinity}, {y2, -Infinity, Infinity} ]

but NIntegrate can't converge and returns with result of 43.04 which is apparently wrong for the large accompanying error estimate 106.

NIntegrate: 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::ncvb: NIntegrate failed to converge to prescribed accuracy after 36 recursive bisections in x2 near {x1,x2,y1,y2} = {0.844124,0.5,0.844124,0.844124}. NIntegrate obtained 43.04088365833965` and 106.01214916269602` for the integral and error estimates.

Other methods, such as QuasiMonteCarlo, also can't work and returns with wrong value.

NIntegrate[eqn, {x1, -Infinity, Infinity}, {x2, -Infinity, 
  Infinity}, {y1, -Infinity, Infinity}, {y2, -Infinity, Infinity}, 
 Method -> "QuasiMonteCarlo"]

NIntegrate:The integral failed to converge after 50000 integrand evaluations. NIntegrate obtained 0.5832445898563868` and 4.830661305937795` for the integral and error estimates.

I've been also thinking about possible ways to simplify this integration analytically, which seems impossible due to the absolute symbol and r terms in the denominator.

I would appreciate a lot if any advice or thought is given on either numerical or analytical side.

Answer 1

If I am not mistaken, the improper multivariate integral under consideration diverges. Here are my arguments.

We begin from the change of variables {x1 == (t + r)/2, x2 == (r - t)/2, y1 == (s + p)/2, y2 == (p - s)/2}

eqn = (Exp[-RealAbs[x1 - x2] - 2  RealAbs[y1 - y2]]  Cos[
3  RealAbs[x1 - x2]])/(Sqrt[x1^2 + y1^2]  Sqrt[x2^2 + y2^2]);
IntegrateChangeVariables[Inactive[Integrate][eqn, {x1, -Infinity, Infinity}, 
{x2, -Infinity, Infinity}, {y1, -Infinity, Infinity}, 
{y2, -Infinity,   Infinity}], {t, r, s, p}, {x1 == (t + r)/2, 
x2 == (r - t)/2,   y1 == (s + p)/2, y2 == (p - s)/2}]

Inactive[ Integrate][(E^(-2 RealAbs[s] - RealAbs[t]) Cos[3 RealAbs[ t]])/(\[Sqrt]((r^2 + (-p + s)^2 - 2 r t + t^2) (r^2 + (p + s)^2 + 2 r t + t^2))), {t, -\[Infinity], \[Infinity]}, {r, -\[Infinity], \ \[Infinity]}, {s, -\[Infinity], \[Infinity]}, {p, -\[Infinity], \ \[Infinity]}]

As we see, the numerator of the integrand (E^(-2 RealAbs[s] - RealAbs[t])* Cos[3 RealAbs[t]]) does not depend on p and r.

Making use of that circumstance, let us integrate over p and r, switching to the polar coordinates.

IntegrateChangeVariables[Inactive [Integrate][
1/ Sqrt[(r^2 + (-p + s)^2 - 2*r*t + t^2)*(r^2 + (p + s)^2 + 2*r*t +  t^2)],  
  {r, -Infinity, Infinity}, {p , -Infinity, Infinity}],
 {\[Rho], \[Theta]}, "Cartesian" -> "Polar", Assumptions -> {s, t} \[Element] Reals]

Inactive[ Integrate][\[Rho]/(\[Sqrt](s^4 + 2 s^2 t^2 + t^4 + \[Rho]^4 + 2 (s^2 - t^2) \[Rho]^2 Cos[2 \[Theta]] - 4 s t \[Rho]^2 Sin[2 \[Theta]])), {\[Rho], 0, \[Infinity]}, {\[Theta], -\[Pi], \[Pi]}]

Let us consider the integral over \[Rho]

Integrate[\[Rho]/(\[Sqrt](s^4 + 2  s^2  t^2 + t^4 + \[Rho]^4 + 
  2  (s^2 - t^2)  \[Rho]^2  Cos[2  \[Theta]] - 
  4  s  t  \[Rho]^2  Sin[2  \[Theta]])), {\[Rho], 0, \[Infinity]},
  Assumptions -> {\[Theta], s, t} \[Element] Reals]

Integrate::idiv: Integral of [Rho]/Sqrt[s^4+2 s^2 t^2+t^4+[Rho]^4+2 (s^2-t^2) [Rho]^2 Cos[2 [Theta]]-4 s t [Rho]^2 Sin[2 [Theta]]] does not converge on {0,[Infinity]}.