# Symbolic integration of potential over a disc : branch cut problem?

Context

I am trying to explore the geometry of a crystal made of irregular bubbles.

See animation here.

very vaguely in the spirit of this post (it is in fact motivated by cosmology and galaxy formation).

So I give myself an interaction potential (which is both attractive and repulsive at large and small distances resp.)

pot[r_] = 1/r^2 + r^2


looking like this

Plot[pot[r], {r, 0.1, 5}]


and I integrate it over a Disk

int= Integrate[ pot[Sqrt[(x - x0)^2 + (y - y0)^2]], {x0, y0} ∈
Disk[{0, 0}, 1]]

(* π (x^2+y^2+1/2)  *)


which incidentally looks suspicious, because it is lacking a repulsion near the disc.

But if I take a specific value for {x,y}

rxy = Thread[{x, y} -> {2, 3}]


and carry out the integration numerically

NIntegrate[
pot[Sqrt[(x - x0)^2 + (y - y0)^2]] /. rxy, {x0, y0} ∈
Disk[{0, 0}, 1], PrecisionGoal -> 6]

(* 42.663 *)


I get a different answer from

  int /. rxy

(* 42.4115 *)


Indeed if I do the replacement First

Integrate[pot[Sqrt[(x - x0)^2 + (y - y0)^2]] /. rxy, {x0, y0} ∈ Disk[{0, 0}, 1]]

(* π (27/2+log(13/12)) *)

N[%]

(* 42.663 *)


So mathematica seems to be doing the general integration wrong.

Questions

Is this a bug? Any workaround?

Check

Indeed I can check by integrating numerically radially away from the edge of the disk that the potential generated by the disc is repulsive at close distance:

dat = ParallelTable[
NIntegrate[
pot[Sqrt[(x - x0)^2 + (y - y0)^2]] /. {x -> r Cos[t],
y -> r Sin[t]} /. t -> Pi/4, {x0, y0}∈
Disk[{0, 0}, 1], PrecisionGoal -> 8],
{r, 1.01, 2, 0.025}];
dat // ListLinePlot


• Indeed, Integrate appears to have a problem with the repulsive part; Integrate[ 1/((x - x0)^2 + (y - y0)^2), {x0, y0} \[Element] Disk[{0, 0}, 1]] returns 0 which is obviously wrong. I'd say, this is a bug. Please inform Wolfram Support. – Henrik Schumacher Nov 28 '18 at 8:16

It's worth noting that the integrals will evaluate separately!

totalPot[x_]=Integrate[(x-x0)^2+(y0)^2,{x0,y0}∈ Disk[{0,0},1]]+
Integrate[1/((x-x0)^2+(y0)^2),{x0,y0}∈ Disk[{0,0},1],Assumptions->{x>1}];
N[totalPot[Sqrt[2^2 + 3^2]]]
(* 42.663 *)


The exact form of the potential being

$$\frac{1}{2} \pi \left(2 r^2-2 \log \left(r^2-1\right)+4 \log (r)+1\right)$$

Where I made sure to use the manifest rotational symmetry to put y=0, and also added an assumption that x is greater than 1 to avoid any issues with divergences in the 1/r^2 case.

Since it's of physical interest, to put units back in, if I take the potential to be an energy density $$k_1 r^2+k_2/r^2$$ and the disk is of radius $$R$$, I find:

$$E(r)=k_1 \frac{\pi}{2}(R^4+2 R^2 r^2)-k_2 \pi \log(1-\frac{R^2}{r^2})$$

As noted by Henrik in the comments, this looks like a bug & should be reported to wolfram support.

• thanks! Would you know how to do the integral over an elliptic disc? int= Integrate[ pot[Sqrt[(x - x0)^2 + (y - y0)^2]], {x0, y0} ∈ Disk[{0, 0}, {1,2}]] – chris Nov 28 '18 at 8:31
• @chris hmm I don't. In 3D for the 1/r potential and an ellipsoid I know that the result is pretty complicated with no nice answer. But maybe it's easier here. – David Nov 28 '18 at 10:33