Two-dimensional integral of a function containing a logarithmic singularity

Question

I have a two-dimensional integral that I'd like to evaluate numerically. I would like to be able to do this as quickly as possible (as I have to calculate many similar integrals). It takes the form $$ \int_{0}^{1}\int_{0}^{1}f(x,y)H_{0}^{(1)}(\sqrt{4(cos(x)-cos(y))^{2}+(sin(x)-sin(y))^{2}})dxdy $$ where $H_{0}^{(1)}(z)$ is the zero order Hankel function of the first kind and $f(x,y)$ is bounded, but is quite complicated itself: $$ f(x,y)=\left(-i\sqrt{4sin^{2}(x)+cos^{2}(x)}+\frac{3cos(x)sin(x)}{\sqrt{4sin^{2}(x)+cos^{2}(x)}}\right)\left(i\sqrt{4sin^{2}(y)+cos^{2}(y)}+\frac{3cos(y)sin(y)}{\sqrt{4sin^{2}(y)+cos^{2}(y)}}\right)e^{-ix}e^{iy} $$ The Hankel function has a logarithmic singularity for small values of its argument. The singularity lies along the line $y=x$ of the integrand. I start from the NIntegrate command as follows:

NIntegrate[f[x,y] HankelH1[0,Sqrt[4(Cos[x]-Cos[y])^2+(Sin[x]-Sin[y])^2]],{x,0,1},{y,0,1}]

How should I modify this to evaluate the integral most efficiently? i.e. what settings associated with NIntegrate should I use? All the examples I've read about online apply to slightly different problems. As it stands the integral takes too long and complains of errors about convergence being too slow. Thanks in advance for any help.

Answer 1

Update: If we tell Mathematica about the singularity along $y=x$ it's much faster and the convergence issues go away with the default method:

NIntegrate[
  f[x, y] HankelH1[0, 
    Sqrt[4 (Cos[x] - Cos[y])^2 + (Sin[x] - Sin[y])^2]],
  {x, 0, 1}, {y, 0, 1}, Exclusions -> (y == x)] // Timing
(* result: {0.375, 2.08433 - 1.89171 I} *)

---

The "LocalAdaptive" method is a bit faster than the "GlobalAdaptive" it chooses by default:

g[a_?NumericQ, b_?NumericQ] := Sqrt[4 Sin[a]^2 + Cos[b]^2]
h[a_?NumericQ, b_?NumericQ] := 3 Cos[a] Sin[b]
f[x_?NumericQ, y_?NumericQ] :=
 (-I g[x, x] + h[x, x]/g[x, x]) (I g[y, y] + h[y, y]/g[y, y]) Exp[I y - I x]

NIntegrate[
 f[x, y] HankelH1[0, 
   Sqrt[4 (Cos[x] - Cos[y])^2 + (Sin[x] - Sin[y])^2]],
   {x, 0, 1}, {y, 0, 1}, Method -> "LocalAdaptive"]//Timing
(* result: {22.75, 2.08433 - 1.89169 I} *)

If you need quick and rough estimates try "AdaptiveQuasiMonteCarlo" which completes in under half a second:

(* result: {0.390625, 2.08553 - 1.88752 I} *)

You may want to take a look at all the options available: http://reference.wolfram.com/language/tutorial/NIntegrateIntegrationStrategies.html