I am having some trouble trying to get Mathematica to do a numerical integral over three dimensions which contains a singularity of dimension 1, and I would like some pointers to solid resources on how to handle this case.
To be more specific, I am interested in integrals of the type $$ \int\frac{\rho(\mathbf r)\:\mathrm d\mathbf r}{\sqrt{(\mathbf r-\mathbf r_1-i\mathbf r_2)^2}}, $$ for which the integrable Coulomb singularity, which reduces to a point when $\mathbf r_2=0$, is now a circle of radius $r_2=|\mathbf r_2|$ centered at $\mathbf r_1$ and with normal vector $\mathbf r_2$. This assumes the standard branch cut on the square root, and a representative charge distribution is $\rho(\mathbf r)=\exp(-r^2)$. (As an aside, this integral also turns out to be completely the wrong way to calculate the analytical continuation of the gaussian charge distribution's electrostatic potential to the complex coordinate space, but I still need to calculate it to show that it doesn't work.)
The naive attempt at doing this is something like
Block[{x = 0, y = 0, z = 1 + I},
NIntegrate[
Exp[-(xx^2 + yy^2 + zz^2)]/Sqrt[(xx - x)^2 + (yy - y)^2 + (zz - z)^2]
, {xx, -∞, ∞}, {yy, -∞, ∞}, {zz, -∞, ∞}
]
]
and it produces slow-convergence warnings:
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.
The integral is not zero or oscillatory, and increasing WorkingPrecision doesn't make the warning go away. I would like to help Mathematica be aware of the singularity, so it can work its magic on it, but I'm having some trouble doing it.
In particular, the documentation details the singularity handling of NIntegrate in NIntegrate Integration Strategies > Singularity Handling, but this seems to congregate almost primarily on (i) single-dimensional integrals, (ii) multi-dimensional integrals where the singularity is an $(N-1)$-dimensional hypersurface of an $N$-dimensional integration domain, or (iii) integrals with singularities at edges or corners of the integration domain.
So, to kick off:
- Is there a way to specify a curve as a singularity, either parametrically or implicitly?
As far as I can tell this is not possible. Certainly, setting
Exclusions -> And[(xx - Re[x])^2 + (yy - Re[y])^2 + (zz - Re[z])^2 ==
Im[x]^2 + Im[y]^2 + Im[z]^2,
Im[x] (xx - Re[x]) + Im[y] (yy - Re[y]) + Im[z] (zz - Re[z]) == 0]
results in an unevaluated integral and the error message
NIntegrate::nexcl: Incorrect Exclusions specification xx^2+yy^2+(-1+zz)^2==1&&-1+zz==0.
The other natural thing to try is a paired exclusion,
Exclusions -> {(xx - Re[x])^2 + (yy - Re[y])^2 + (zz - Re[z])^2 ==
Im[x]^2 + Im[y]^2 + Im[z]^2,
Im[x] (xx - Re[x]) + Im[y] (yy - Re[y]) + Im[z] (zz - Re[z]) == 0}
and this seems to work for some arguments. However, if you try, say, x = 3, y = 0, z = 1 + 5 I, it again slows down and gives a slow-convergence message.
I find it weird that there isn't a built-in way to specify sub-dimensional singularities (so to speak), but then again that does sound like a difficult problem. Luckily, for a ring singularity it is possible to parametrize the integral in a way that turns the singularity into a line segment parallel to a coordinate axis, by taking the integral in spherical coordinates about $\mathbf r_1$ with the poles along $\mathbf r_2$. This gets kind of messy for an arbitrary axis, but if you set $x_1=y_1=0$, with $\mathbf r_2$ along the $z$ axis, then you get a relatively simple expression:
Block[{x1 = 0, y1 = 0, z1 = 1, x2 = 0, y2 = 0, z2 = 1},
NIntegrate[
Exp[-Total[({x1, y1, z1} + {r Sin[θ] Cos[ϕ], r Sin[θ] Cos[ϕ], r Cos[θ]})^2]]
/Sqrt[r^2 - z2^2 - 2 I z2 r Cos[θ]] r^2 Sin[θ]
, {r, 0, z2, ∞}, {θ, 0, π/2, π}, {ϕ, 0, 2 π}
]
]
To start off, this approach enables one to give the coordinates of the singularity as midpoints of the integral domains for r and θ, but that still doesn't help. In fact, it gets worse: you get a slwcon message, a value for the integral that's almost certainly wrong, and the error message
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 16.5826 -2.65363 I and 0.00015976992207685935` for the integral and error estimates.
Similarly, adding the explicit exclusions
Exclusions -> {r == z2, θ == π/2}
which you'd hope would improve things, only slows down the integral a fair bit, keeps the error messages, and keeps the suspect value for the integral.
Similarly, playing around with the integration strategies in the Singularity Handling documentation hasn't yielded anything that behaves any better. So:
How can I make Mathematica aware of my ring singularity, and how can I get it to compute this integral quickly and accurately?
If the built-in integrator simply can't handle situations this complicated, is there a good third-party package that will fill the role?
