Vrftjkry

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

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-1right)+4 log (r)+1right)$$

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.