6
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Related, on the side of $(N-1)$-dimensional singularities, but inconclusive: How can I handle curve singularity in this NIntegrate integration? $\endgroup$ – Emilio Pisanty Jul 6 '16 at 20:30
  • $\begingroup$ Constructive feedback appreciated =). $\endgroup$ – Emilio Pisanty Jul 6 '16 at 21:28
  • $\begingroup$ I think this discussion is relevant. $\endgroup$ – Anton Antonov Jul 7 '16 at 0:30
  • $\begingroup$ It is good to know your original goal: "I still need to calculate it to show that it doesn't work." I guess you claim that your analytic continuation does not work. It works, however, in at least one case when $\rho(\vec{r})\equiv\delta(\vec r-\vec r_0)$. $\endgroup$ – yarchik Jul 7 '16 at 7:00
  • $\begingroup$ @yarchik For a continuous distribution it doesn't work, because the integrand is not analytic as a function of $\mathbf r_1+i\mathbf r_2$ for $\mathbf r$ at the singularity. The full details are in my thesis, which will be available in a few months. The short story, though, is that if you graph, say, $V(3,0,z)$ as a function of complex $z$, you get a function that's continuous and bounded, which is in violation of Liouville's theorem and the maximum principle. $\endgroup$ – Emilio Pisanty Jul 7 '16 at 9:14
5
$\begingroup$

If you set exclusions more directly, you won't get failure to achieve the precision goal. The NIntegrate::slwcon will remain though, but it's pretty harmless. Here goes the code:

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, -∞, ∞}
  , Exclusions -> (xx - x)^2 + (yy - y)^2 + (zz - z)^2 == 0
  ]
 ]

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.

2.96809041831953 - 1.72434400654356 I

When I set PrecisionGoal->7, the message repeats twice, but no other message like NIntegrate::eincr appears, and the result has small difference from the one with default precision goal. So I suppose it's more or less reliable.

If I set PrecisionGoal->8, I get NIntegrate::eincr, and the error estimate is given as 3.92*^-8, which is just a little short of the requested precision goal. The resulting number is again compatible with previous attempts.


Still, the best results appear when you just switch integration strategy. Namely, with "LocalAdaptive" I was able to get the result without any messages at all, even with precision goal of 12:

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, -∞, ∞}, 
  Method -> "LocalAdaptive", PrecisionGoal -> 12]]

2.9680903889399 - 1.72434400670741 I

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.