# Two-dimensional integral of a function containing a logarithmic singularity

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.

• Post all Mathematica code no LaTex? Where is:f[x,y] ?.It's impolite of you to expect an answer without providing the code to reproduce. Commented Jul 7, 2020 at 16:49

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

• The integral of f[x, y] HankelH1[0, Sqrt[4 (Cos[x] - Cos[y])^2 + (Sin[x] - Sin[y])^2]] over the square is an improper double integral. I don't know any numeric method to calculate it. The work of the Exclusions option is not explained in reference.wolfram.com/language/tutorial/… . Commented Jul 7, 2020 at 19:55
• @user64494 sorry, 1) who are you? Are you @Chris? 2) There's nothing wrong with the integral, 3) It is explained - it's literally right there: wolfram.com/xid/0far7me5ywb23ren50sgk-bs9aib Commented Jul 7, 2020 at 19:58
• This is not an improper integral. Series shows, the integrand goes with I Log[x-y]  towards the singularity. This can be integrated regularily Integrate[I Log[x - y], {x, 0, 1}, {y, 0, 1}]  . Commented Jul 8, 2020 at 6:30
• Ha, I'm certainly not user64494. The integral is not improper. Commented Jul 8, 2020 at 6:58